# load the package:

library("TrendLSW")

#---------Section 2.2 Code----------

# create the trend as a function, and the spectrum as a list of functions:
n1 <- 512

trend1 <- function(z) {
  3 * (32 * z^3 - 48 * z^2 + 22 * z - 3)
}

spec1 <- vector(mode = "list", length = log2(n1))

spec1[[2]] <- function(z) {
  2 + 20 * z - 20 * z^2
}
spec1[[4]] <- function(z) {
  3
}

# simulate an example time series:

set.seed(123)
x1 <- TLSWsim(trend = trend1, spec = spec1)

# Figure 1: plot time series and true trend:

plot(x1, ylab = expression(X[t]), type = "l", xlab = "Time")

index1 <- seq(from = 0, to = 1, length = n1)
lines(trend1(index1), col = 2, lwd = 3)

#---------Section 3.3 Code-----------

# use the TLSW function on x1:

x1.TLSW <- TLSW(x1)

# Figure 2: Plot the spectrum estimate of x1:

plot(x1.TLSW, plot.type = "spec")

# calculate the lacf from the TLSW object:

x1.lacf <- TLSWlacf(x1.TLSW)

# Figure 3: Plot the lacf at a few time points:

plot(x1.lacf$lacr[50, ],
  ylim = c(-0.4, 1), type = "l",
  xlab = "Lag", ylab = "ACF", lwd = 2, col = 4
)

lines(x1.lacf$lacr[200, ], col = 2, lwd = 2, lty = 2)

lines(x1.lacf$lacr[350, ], col = 1, lwd = 2, lty = 4)

legend(20, 1, c("t = 50", "t = 200", "t = 350"),
  lty = c(1, 2, 4),
  lwd = c(2, 2, 2), col = c(4, 2, 1), bty = "n", cex = 1.2
)


#---------Section 4.3 Code-----------

# plot the trend estimate of x1:

plot(x1.TLSW, plot.type = "trend")


#---------Section 5 Example-----------

# create a new example with trend given by numeric vector and spectrum given by numeric matrix:

n2 <- 1024
index2 <- seq(from = 0, to = 1, length = n2)

trend2 <- 5 * sin(pi * 6 * index2) + c(seq(from = 0, to = 10, length = 300), seq(from = 10, to = -4, length = 724))

spec2 <- matrix(0, nrow = log2(n2), ncol = n2)

spec2[1, ] <- seq(from = 2, to = 10, len = n2)

spec2[3, ] <- c(
  rep(1, 200), seq(from = 1, to = 6, length = 200),
  seq(from = 6, to = 1, length = 200), rep(1, 424)
)

spec2[5, ] <- 2 + 4 * sin(4 * pi * index2)^2

# simulate an example of this process:

set.seed(1234)
x2 <- TLSWsim(trend = trend2, spec = spec2, filter.number = 4)

# Figure 5: plot the data, true trend and true spectrum. Internal package function is used to convert the
# numeric matrix to a wd object amenable for plotting.

par(mfrow = c(1, 2))

plot.ts(x2, ylab = expression(X[t]), col = "grey60")
lines(trend2, col = 2, lwd = 2)

plot(TrendLSW:::mat.to.spec(spec2),
  ylab = "Scale", sub = "", xlab = "Time",
  main = ""
)

# estimate the trend and spectrum of x2. Seed set for bootstrap replicability.

set.seed(10)
x2.TLSW <- TLSW(x2,
  T.filter.number = 6, T.family = "DaubLeAsymm",
  T.est.type = "nonlinear", T.CI = TRUE, T.reps = 500,
  S.filter.number = 4, S.family = "DaubExPhase",
  S.do.diff = TRUE, S.smooth.type = "median",
  S.binwidth = 128
)

# Figures 6 and 7: plot the trend and spectrum estimates:

par(mfrow = c(1, 1))

plot(x2.TLSW)


#---------Section 5.1 Code-----------

# Figures 8 and 9: code for plotting the TLSW estimates of x2:

plot(x2.TLSW,
  trend.plot.args = list(
    col = "black", type = "p", pch = 16, T.col = "blue", T.lwd = 2,
    poly.col = "grey", CI.col = "purple", CI.lwd = 2, CI.lty = 1
  ),
  spec.plot.args = list(scaling = "by.level", ylab = "Level")
)



#--------------------------Section 6.1 Code--------------

# Figure 10: load the celegansbio data, use TLSW() for estimation, and plot the trend estimate, along with turning points:

data("celegansbio")
out <- TLSW(celegansbio,
  T.filter.number = 10, T.family = "DaubExPhase",
  T.transform = "nondec"
)
plot(out, plot.type = "trend", trend.plot.args = list(
  ylab = "Luminescence (AU)",
  T.col = "blue", T.lwd = 3, main = ""
), cex.lab = 1.3, col = "black")
delta <- diff(out$trend.est$T)
turns <- which(delta[-1] * delta[-length(delta)] < 0) + 1
points(turns, out$trend.est$T[turns], col = "red", pch = 18, cex = 2)

# Figure 11: comparing different wavelet fits:

out.EP1 <- TLSW(celegansbio, T.filter.number = 1, T.family = "DaubExPhase",
                    T.transform = "nondec")
out.EP4 <- TLSW(celegansbio, T.filter.number = 4, T.family = "DaubExPhase",
                    T.transform = "nondec")
out.LA10 <- TLSW(celegansbio, T.filter.number = 10, T.family = "DaubLeAsymm",
                     T.transform = "nondec")


plot(out, plot.type = "trend", trend.plot.args = list(
  ylab = "Luminescence (AU)",
  T.col = "blue", T.lwd = 2, main = ""),
  cex.lab = 1.3, col = "grey60")

lines(out.EP1$trend.est$T, col = "black", lwd = 2, lty = 2)
lines(out.EP4$trend.est$T, col = "red", lty = 3, lwd = 2)
lines(out.LA10$trend.est$T, col = "green", lty = 4, lwd = 2)

legend(10, 85000, c("EP10", "EP1", "EP4", "LA10"),
       lty = c(1, 2, 3, 4),lwd = c(2, 2, 2, 2),
       col = c("blue", "black", "red", "green"), bty = "n", cex = 1.2)

n <- length(celegansbio)

trends <- list(EP1 = out.EP1$trend.est$T,EP4 = out.EP4$trend.est$T,
               EP10 = out$trend.est$T, LA10 = out.LA10$trend.est$T)

trends.rmse <- lapply(trends,
                      function(y){sqrt(sum((y-celegansbio)^2)/n)})

#------------------- Section 6.2 Code---------------

# snapshot of experiment 3, user 2
# z.acc contains the time series, and z.labels is a data frame containing the
# activities and their start and end times.

data("z.acc")
data("z.labels")

# Figure 12: plot z.acc along with activity start and end times:

plot(z.acc,
  ylab = expression("Z axis acceleration (m/s"^2 * ")"),
  type = "l", xlab = "Time"
)

abline(v = c(z.labels$start[c(1, 3, 5)], z.labels$end[c(1, 3, 5)]), col = "blue", lwd = 1.5)
abline(v = c(z.labels$start[c(2, 4, 6)], z.labels$end[c(2, 4, 6)]), col = "red", lwd = 1.5, lty = 2)

# use the TLSW function on z.acc:

z.acc.TLSW <- TLSW(z.acc, S.smooth.type = "epan", S.binwidth = 150)

# Figure 13: plot the estimated spectrum of z.acc:

plot(z.acc.TLSW, plot.type = "spec", spec.plot.args = list(scaling = "by.level"))

abline(v = c(z.labels$start[c(1, 3, 5)], z.labels$end[c(1, 3, 5)]), col = "blue", lwd = 1.5)
abline(v = c(z.labels$start[c(2, 4, 6)], z.labels$end[c(2, 4, 6)]), col = "red", lwd = 1.5, lty = 2)

# Figure 14: plot the spectrum only at scale 5 and add regions for activity.
# The wavethresh function accessD is used (note scaling convention is different):

plot(1:6000, wavethresh::accessD(z.acc.TLSW$spec.est$S, level = 8)[1:6000],
  ylab = expression(S[5](z)),
  xlab = "Time", type = "l", lwd = 1.5
)

for (i in c(1, 3, 5)) {
  rect(
    xleft = z.labels$start[i], xright = z.labels$end[i], ybottom = -1, ytop = 1,
    border = NA, col = adjustcolor("blue", alpha = 0.3)
  )
}

for (i in c(2, 4, 6)) {
  rect(
    xleft = z.labels$start[i], xright = z.labels$end[i], ybottom = -1, 1,
    border = NA, col = adjustcolor("red", alpha = 0.3)
  )
}




## R Session Information

sessionInfo()