havok: Procedures for the Analysis of Nonlinear Dynamical Systems

Documented in plot.havok

#' Plotting functions for the havok package of class “havok"
#'
#' @description Generic plotting function for object of class ("havok")
#' @param x A "havok" object.
#' @param what See details.
#' @param ... Other calls to plot.
#' @details Arguments for \code{what} parameter: \itemize{
#' \item{\code{"Interactive"} - }{An interactive plotting function.}
#' \item{\code{"Reconstruction"} - }{Scaled reconstruction of the major component of a time-series.}
#' \item{\code{"Forcing"} - }{Forcing vector derived from HAVOK.}
#' \item{\code{"Both"} - }{A combination of 'reconstruction' and 'forcing'.}
#' \item{\code{"U-modes"} - }{U modes of the reconstructed time series.}
#' \item{\code{"Vembedded"} - }{A 2D reconstruction of the attractor based on Eigen vectors.}
#' \item{\code{"Nonlinear"} - }{A 2D reconstruction of the attractor with nonlinear regions colored red.}
#' \item{\code{"SSmod"} - }{Ouput time series of a HAVOK state space model.}
#' \item{\code{"SSembedded"} - }{A 2D reconstruction of the attractor based a HAVOK state space model.}}
#' @examples
#'
#' \dontrun{
#'#Generate Data
#'##Set Lorenz Parameters
#'parameters <- c(s = 10, r = 28, b = 8/3)
#'n <- 3
#'state <- c(X = -8, Y = 8, Z =27) ##Inital Values
#'
#'#Intergrate
#'dt <- 0.001
#'tspan <- seq(dt, 200, dt)
#'N <- length(tspan)
#'
#'Lorenz <- function(t, state, parameters) {
#'  with(as.list(c(state, parameters)), {
#'    dX <- s * (Y - X)
#'    dY <- X * (r - Z) - Y
#'    dZ <- X * Y - b * Z
#'    list(c(dX, dY, dZ))
#'  })
#'}
#'
#'out <- ode(y = state, times = tspan, func = Lorenz, parms = parameters, rtol = 1e-12, atol = 1e-12)
#'xdat <- out[, "X"]
#'t <- out[, "time"]
#'hav <- havok(xdat = xdat, dt = dt, stackmax = 100, lambda = 0,
#'             rmax = 15, polyOrder = 1, useSine = FALSE)
#'
#'plot(hav)
#'}
###################################

## S3 method for class "havok"
#' @export
plot.havok <- function(x, what = "interactive", ...) {

  if (tolower(what) == "interactive"){

    cat("--- Please select a plot type by number ---\n
        Plot Types:
        1 - V1 Reconstruction  5 - V Embedded Attractor
        2 - Forcing            6 - Strongly Nonlinear Region
        3 - Both               7 - SS Model Output
        4 - U-modes            8 - SS Embedded Attractor")

    repeat {
      what <- readline("Please select a number (press esc to exit): ")

      if (!what %in% 1:8){
        stop("Please pick a number between 1 and 8")
      }


      if (what == 1){
        graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
        graphics::plot(x$havokSS$t, scale(x$Vr_aligned[,1]), type = "l", xlab = "Time", ylab = "V1 vs. V1_hat", ...)
        graphics::lines(x$havokSS$t, scale(x$havokSS$y[1,]), col = "red")
        graphics::legend("topleft", legend = c("V1", "Predicted V1"), lwd = c(1,1), col = c("black", "red"))
        graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
      }

      if (what == 2){
        graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
        graphics::plot(x$havokSS$t, x$Vr_aligned[,x$r], type = "l", xlab = "Time", ylab = "Forcing", ...)
        graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
      }

      if (what == 3){
        graphics::par(mfrow=c(2,1), mai = c(0.5, 1.0, 0.1, 0.1))

        graphics::plot(x$havokSS$t, scale(x$Vr_aligned[,1]), type = "l", xlab = "Time", ylab = "V1 vs. V1_hat", ...)
        graphics::lines(x$havokSS$t, scale(x$havokSS$y[1,]), col = "red")
        graphics::legend("topleft", legend = c("V1", "Predicted V1"), lwd = c(1,1), col = c("black", "red"))

        graphics::par(mai = c(1, 1.0, 0.1, 0.1))

        graphics::plot(x$havokSS$t, x$Vr_aligned[,x$r], type = "l", xlab = "Time", ylab = "Forcing", ...)

        graphics::par(mfrow=c(1,1), mai = c(1.02, 0.82, 0.82, 0.42))
      }

      if (what == 4){
        graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
        plotBuff <- (.1 * (max(x$U[, 1:x$r]) - min(x$U[, 1:x$r])))

        graphics::plot(x$U[,1], ylab = "U", xlab = "Time",
                       type = "l",
                       xaxt = "n",
                       col = grDevices::rainbow(x$r)[1],
                       ylim = c(min(x$U[,1:x$r]) - (.1*plotBuff), max(x$U[,1:x$r]) + (.1*plotBuff)))

        graphics::axis(1, at = seq(0, length(x$U[,1]), length.out = 10),
                       labels = round(seq(0, x$params["stackmax",], length.out = 10)*x$params["dt",], 2))

        for (i in 1:x$r){
          graphics::lines(x$U[,i], col = grDevices::rainbow(x$r, alpha = 1/sqrt(i))[i])
        }

        if (x$r <= 6){
          graphics::legend("topleft",
                           fill = grDevices::rainbow(x$r, alpha = 1/sqrt(1:x$r)),
                           legend = paste("r = ", 1:x$r, sep = ""))
        }

        if (x$r > 6){
          graphics::legend("topleft",
                           fill = grDevices::rainbow(x$r, alpha = 1/sqrt(1:x$r))[c(1:6,x$r)],
                           legend = c(paste("r = ", 1:5, sep = ""),
                                      "...",
                                      paste("r = ", x$r, sep = ""))
          )
        }
        graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))

      }


      if (what == 5){
        graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
        graphics::plot(x$Vr[,1], x$Vr[,2], type = "l", xlab = "V1", ylab = "V2", ...)
        graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
      }

      if (what == 6){
        havForce <- active_forcing(x)
        graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
        graphics::plot(x$Vr[,1], x$Vr[,2], type = "l", xlab = "V1", ylab = "V2", ...)
        graphics::segments(utils::head(x$Vr[,1], -1), utils::head(x$Vr[,2], -1), x$Vr[,1][-1], x$Vr[,2][-1], ifelse(havForce$active==1,"red","black"))
        graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
      }

      if (what == 7){
        graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
        graphics::plot(x$havokSS$t, x$havokSS$y[1,], type = "l", xlab = "Time", ylab = "y", ...)
        graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
      }

      if (what == 8){
        graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
        graphics::plot(x$havokSS$y[1,], x$havokSS$y[2,], type = "l", xlab = "y1", ylab = "y2", ...)
        graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
      }

    }
  }

  if (what == 1 | tolower(what) == "reconstruction"){
    graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
    graphics::plot(x$havokSS$t, scale(x$Vr_aligned[,1]), type = "l", xlab = "Time", ylab = "V1 vs. V1_hat", ...)
    graphics::lines(x$havokSS$t, scale(x$havokSS$y[1,]), col = "red")
    graphics::legend("topleft", legend = c("V1", "Predicted V1"), lwd = c(1,1), col = c("black", "red"))
    graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
  }

  if (what == 2 | tolower(what) == "forcing"){
    graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
    graphics::plot(x$havokSS$t, x$Vr_aligned[,x$r], type = "l", xlab = "Time", ylab = "Forcing", ...)
    graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
  }

  if (what == 3 | tolower(what) == "both"){
    graphics::par(mfrow=c(2,1), mai = c(0.5, 1.0, 0.1, 0.1))

    graphics::plot(x$havokSS$t, scale(x$Vr_aligned[,1]), type = "l", xlab = "Time", ylab = "V1 vs. V1_hat", ...)
    graphics::lines(x$havokSS$t, scale(x$havokSS$y[1,]), col = "red")
    graphics::legend("topleft", legend = c("V1", "Predicted V1"), lwd = c(1,1), col = c("black", "red"))

    graphics::par(mai = c(1, 1.0, 0.1, 0.1))

    graphics::plot(x$havokSS$t, x$Vr_aligned[,x$r], type = "l", xlab = "Time", ylab = "Forcing", ...)

    graphics::par(mfrow=c(1,1), mai = c(1.02, 0.82, 0.82, 0.42))
  }

  if (what == 4 | tolower(what) == "u-modes"){
    graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
    plotBuff <- (.1 * (max(x$U[, 1:x$r]) - min(x$U[, 1:x$r])))

    graphics::plot(x$U[,1], ylab = "U", xlab = "Time",
                   type = "l",
                   xaxt = "n",
                   col = grDevices::rainbow(x$r)[1],
                   ylim = c(min(x$U[,1:x$r]) - (.1*plotBuff), max(x$U[,1:x$r]) + (.1*plotBuff)))

    graphics::axis(1, at = seq(0, length(x$U[,1]), length.out = 10),
                   labels = round(seq(0, x$params["stackmax",], length.out = 10)*x$params["dt",], 2))

    for (i in 1:x$r){
      graphics::lines(x$U[,i], col = grDevices::rainbow(x$r, alpha = 1/sqrt(i))[i])
    }

    if (x$r <= 6){
      graphics::legend("topleft",
                       fill = grDevices::rainbow(x$r, alpha = 1/sqrt(1:x$r)),
                       legend = paste("r = ", 1:x$r, sep = ""))
    }

    if (x$r > 6){
      graphics::legend("topleft",
                       fill = grDevices::rainbow(x$r, alpha = 1/sqrt(1:x$r))[c(1:6,x$r)],
                       legend = c(paste("r = ", 1:5, sep = ""),
                                  "...",
                                  paste("r = ", x$r, sep = ""))
      )
    }
    graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))

  }


  if (what == 5 | tolower(what) == "vembedded"){
    graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
    graphics::plot(x$Vr[,1], x$Vr[,2], type = "l", xlab = "V1", ylab = "V2", ...)
    graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
  }

  if (what == 6 | tolower(what) == "nonlinear"){
    havForce <- active_forcing(x)
    graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
    graphics::plot(x$Vr[,1], x$Vr[,2], type = "l", xlab = "V1", ylab = "V2", ...)
    graphics::segments(utils::head(x$Vr[,1], -1), utils::head(x$Vr[,2], -1), x$Vr[,1][-1], x$Vr[,2][-1], ifelse(havForce$active==1,"red","black"))
    graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
  }

  if (what == 7 | tolower(what) == "ssmod"){
    graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
    graphics::plot(x$havokSS$t, x$havokSS$y[1,], type = "l", xlab = "Time", ylab = "y", ...)
    graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
  }

  if (what == 8 | tolower(what) == "ssembedded"){
    graphics::par(mai = c(0.9, 0.8, 0.1, 0.1))
    graphics::plot(x$havokSS$y[1,], x$havokSS$y[2,], type = "l", xlab = "y1", ylab = "y2", ...)
    graphics::par(mai = c(1.02, 0.82, 0.82, 0.42))
  }



}



# Copyright 2020 Robert Glenn Moulder Jr. & Elena Martynova
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.