R/plot_functions.R

In zejiang-unsw/WASP_1.0.0: WAvelet System Prediction

#' Plot function: Plot origianl time series and decomposed frequency components
#'
#' @param y           Original time series (Y).
#' @param y.mra       Decomposed frequency components (d1,d2,..,aJ).
#' @param limits.x    x limit for plot.
#' @param limits.y    y limit for plot.
#'
#' @return A plot with origianl time series and decomposed frequency components.
#' @export
#'
#' @examples
#' ###synthetic example
#' #frequency, sampled from a given range
#' fd <- c(3,5,10,15,25,30,55,70,95)
#' data.SW3 <- data.gen.SW(nobs=512,fp=c(15,25,30),fd=fd)
#'
#' x <- data.SW3$x
#' xx <- padding(x,pad="zero")
#' ###wavelet transfrom
#' # wavelet family, extension mode and package
#' wf <- "d4" # wavelet family D8 or db4
#' boundary <- "periodic"
#' pad <-  "zero"
#' if(wf!="haar") v <- as.integer(as.numeric(substr(wf,2,3))/2) else v <- 1
#'
#' #Maximum decomposition level J
#' n <- length(x)
#' J <- ceiling(log(n/(2*v-1))/log(2)) #(Kaiser, 1994)
#'
#' ###decomposition
#' x.mra <- waveslim::mra(xx,wf=wf, J=J, method="dwt", boundary="periodic")
#' x.mra.m <- matrix(unlist(x.mra), ncol=J+1)
#'
#' print(sum(abs(x-rowSums(x.mra.m[1:n,]))))               #additive check
#' var(x);sum(apply(x.mra.m[1:n,],2,var))                  #variance check
#'
#' limits.x <- c(0,n); limits.y <- c(-3,3)
#' mra.plot(x, x.mra.m, limits.x, limits.y)

mra.plot <- function(y, y.mra, limits.x, limits.y,...){

  J <- ncol(y.mra)-1
  par(mfcol=c(J+2,1), mar=c(0,3,0.2,1),# margin of the plot
      oma = c(2, 1, 1, 1), # move plot to the right and up
      mgp=c(1, 1, 0), # move axis labels closer to axis
      bg = "transparent", pty="m", # maximal plotting region
      ps=8) # text size

  plot(y, type = "l", axes=FALSE, xlab=NA, xaxs="i",
       xlim=limits.x, ylim=limits.y,...);box()

  for(i in 1:J){

    plot(y.mra[,i], type="l", axes=FALSE, xlab=NA, ylab=paste0("d",i), xaxs="i",
         xlim=limits.x, ylim=limits.y);box()
  }

  plot(y.mra[,J+1], type="l", axes=FALSE, xlab=NA, ylab=paste0("a",J), xaxs="i",
       xlim=limits.x, ylim=limits.y);box()
  axis(side=1, at=seq(limits.x[1],limits.x[2],by = 50), labels=seq(limits.x[1],limits.x[2],by = 50))

  return(recordPlot())

}


#' Plot function: Variance structure before and after variance transformation
#'
#' @param dwt.data  Output data from variance transformation function
#'
#' @return A plot with variance structure before and after variance transformation.
#' @export
#'
#' @examples
#' data("data.HL")
#' data("data.SW1")
#' #variance transfrom
#' dwt.SW1<- dwt.vt(data.SW1, wf, J, method, pad, boundary)
#'
#' #plot
#' fig1 <- fig.dwt.vt(dwt.SW1)
#' fig1
#'
#' #variance transfrom
#' dwt.HL<- dwt.vt(data.HL, wf, J, method, pad, boundary)
#'
#' #plot
#' fig2 <- fig.vt(dwt.HL)
#' fig2

fig.dwt.vt <- function(dwt.data){

  x <- dwt.data$x
  dp <- dwt.data$dp;dp.n <- dwt.data$dp.n
  n=nrow(dp);ndim <- ncol(dp)

  wf <- dwt.data$wavelet; method <- dwt.data$method; boundary <- dwt.data$boundary; pad=dwt.data$pad
  if(wf!="haar") v <- as.integer(as.numeric(substr(wf,2,3))/2) else v <- 1
  J <- ceiling(log(n/(2*v-1))/log(2))

  #variance structure of response in spectrum
  mra.x <- waveslim::mra(padding(x,pad), wf = wf, J = J, method = method, boundary = boundary)
  idwt.x <- lapply(mra.x, function(z) z[1:n])

  x.Dj=c(unlist(lapply(idwt.x, var))/var(x))
  x.Dj <- data.frame(matrix(rep(x.Dj,ndim), ncol = ndim))
  colnames(x.Dj) <- paste0("X",1:ndim)
  x.Dj$Group <- 1

  #variance structure of predictors in spectrum
  idwt.dp <- lapply(1:ndim, function(i) waveslim::mra(padding(dp[,i],pad), wf = wf, J = J, method = method, boundary = boundary))
  idwt.dp.n <- lapply(1:ndim, function(i) waveslim::mra(padding(dp.n[,i],pad), wf = wf, J = J, method = method, boundary = boundary))

  dpred.dp <- lapply(1:length(idwt.dp), function(i) lapply(idwt.dp[[i]], function(z) z[1:n]))
  dpred.dp.n <- lapply(1:length(idwt.dp.n), function(i) lapply(idwt.dp.n[[i]], function(z) z[1:n]))

  dp.Dj = data.frame(sapply(1:length(dpred.dp), function(i) c(unlist(lapply(dpred.dp[[i]], var)))/var(dp[,i])))
  dp.Dj$Group <- 2
  dp.Dj.n = data.frame(sapply(1:length(dpred.dp.n), function(i) c(unlist(lapply(dpred.dp.n[[i]], var)))/var(dp.n[,i])))
  dp.Dj.n$Group <- 3

  df <- cbind(Level=rep(1:length(idwt.x),3),rbind(x.Dj, dp.Dj, dp.Dj.n))
  df.n <- gather(df, Predictor, Value, 2:(ndim+1))
  df.n$Predictor <- factor(df.n$Predictor, levels=paste0("X",1:ndim))

  #barplot+lineplot
  fig <- ggplot(df.n[df.n$Group!=1,], aes(factor(Level), Value)) +
    geom_bar(aes(fill = factor(Group)), position = "dodge", stat="identity") +
    geom_line(data=df.n[df.n$Group==1,],aes(x=factor(Level), y=Value),group=1) +
    facet_grid(Predictor~., scales = "free", space = "free") +

    scale_y_continuous(breaks=seq(0,1,0.5), limits=c(0,1))

    scale_fill_manual(values=c("red","blue"))+

    xlab("Decomposition level") +
    ylab("Variance (percent)") +

    theme_bw() +
    theme(text = element_text(size = 8),
          plot.margin = unit(c(1,1,1,1), "cm"),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),

          strip.text.y = element_text(angle = 90),
          legend.position="none",
          legend.title = element_blank()

    )

  return(fig)

}