library("waveslim")
### read in oxygen isotope data from Antarctic ice core
### using one of the following two R commands:
temp <- scan("/Users/dbp/R/Lund/oxygen.dat")  ### will need to adjust pathname
temp <- scan("http://faculty.washington.edu/dbp/DATA/oxygen.dat")
### extract times from odd indexed elements of temp
ox.times.0 <- temp[2*(1:352)-1]
### extract actual data from even indexed elements of temp
ox <- temp[2*(1:352)]
N <- length(ox)
par(mfcol=c(2,1))
plot(ox.times.0, ox)
plot(ox.times.0, ox, type="l")
### look at documentation for dwt function
?dwt
### level 1 Haar
ox.dwt.haar <- dwt(ox, "haar", 1)
ox.dwt.haar
names(ox.dwt.haar)
class(ox.dwt.haar)
### check that energy of series is preserved in DWT coefficients
sum(ox^2)
sum(ox.dwt.haar$d1^2) + sum(ox.dwt.haar$s1^2)
N1 <- length(ox.dwt.haar$d1)
ox.times.haar.1 <- (ox.times.0[2*(1:N1)] + ox.times.0[2*(1:N1)-1])/2
plot(ox.times.haar.1,ox.dwt.haar$d1,type="h")
### level 1 LA(8)
ox.dwt.la8 <- dwt(ox, "la8", 1)
names(ox.dwt.la8)
### again,let's check the energy
sum(ox^2)
sum(ox.dwt.la8$d1^2) + sum(ox.dwt.la8$s1^2)
### define function that gives times for LA(8) wavelet coefficients
### (see overhead III-64)
LA8.wc.times <- function(t0, delta.t, level, N)
  {
    ttj <- 2^level
    time.indices <- 0:(N/ttj-1)
    t0 + ((ttj*(time.indices + 1) - 1 - (7*(ttj-1)+1)/2) %% N) * delta.t
  }
###
oxy.delta <- ox.times.0[2] - ox.times.0[1]
ox.times.la8.1 <- LA8.wc.times(ox.times.0[1], oxy.delta, 1, N)
ox.times.la8.1[1:5]
plot(ox.times.la8.1,ox.dwt.la8$d1,type="h")
### function that gives number of boundary DWT coefficients
### (see overhead IV-60)
Lpj <- function(L,j) ceiling((L-2)*(1- 1/2^j))
Lpj(8,1)
plot(ox.times.la8.1,ox.dwt.la8$d1,type="h")
plot(ox.times.la8.1[4:N1],ox.dwt.la8$d1[-(1:3)],type="h")
lines(ox.times.la8.1[1:3],ox.dwt.la8$d1[1:3],type="h",col="red")
### compute level 2 coeffs from level 1 scaling coeffs
ox.dwt.la8.2.only <- dwt(ox.dwt.la8$s1, "la8", 1)
ox.dwt.la8.2.only$d1[1:5]
### compute level 2 dwt for ox data
ox.dwt.la8.2 <- dwt(ox, "la8", 2)
ox.dwt.la8.2$d2[1:5]
### extract DWT coefficient vectors
V2 <- ox.dwt.la8.2$s2
W2 <- ox.dwt.la8.2$d2
W1 <- ox.dwt.la8.2$d1
### plot level 2 wavelet coefficients
N2 <- N1/2
ox.times.la8.2 <- LA8.wc.times(ox.times.0[1], oxy.delta, 2, N)
par(mfrow=c(2,1))
plot(ox.times.la8.2, ox.dwt.la8.2$d2,type="h",xlim=c(1800,2000))
Lpj(8,2)
plot(ox.times.la8.2[6:N2],ox.dwt.la8.2$d2[-(1:5)],type="h",xlim=c(1800,2000))
lines(ox.times.la8.2[1:5],ox.dwt.la8.2$d2[1:5],type="h",col="red")
### level J0 = 2 MODWT
ox.modwt.la8.2 <- modwt(ox, "la8", 2)
class(ox.modwt.la8.2)
names(ox.modwt.la8.2)
### extract MODWT coefficient vectors
tV2 <- ox.modwt.la8.2$s2
tW2 <- ox.modwt.la8.2$d2
tW1 <- ox.modwt.la8.2$d1
sum(ox^2)
sum(tW1^2) + sum(tW2^2) + sum(tV2^2)
sum(ox^2) - (sum(tW1^2) + sum(tW2^2) + sum(tV2^2))
### need to align coefficients (see page 114 of WMTSA);
### first, need to define function that returns width of
### equivalent filter of level j for wavelet filter
### whose width at level 1 is L
Lj <- function(L,j) (2^j - 1) * (L - 1) + 1
###
abs.nuHj.LA8 <- function(j) { Lj(8,j)/2 }
abs.nuHj.LA8(1)
abs.nuHj.LA8(2)
abs.nuGj.LA8 <- function(j) { (Lj(8,j)-1)*6/14 }
### circularly shift coefficients to align with correct times
tV2s <- cshift(tV2, -abs.nuGj.LA8(2))
tW2s <- cshift(tW2, -abs.nuHj.LA8(2))
tW1s <- cshift(tW1, -abs.nuHj.LA8(1))
### define a function to plot boundary markers 
add.boundary.markers <- function (times, j, abs.nu.func=abs.nuHj.LA8, L=8) {
  n.boundary.coeffs <- Lj(L,j) - 1 # see overhead IV-61
  ne <- abs.nu.func(j)             # number of coeffs at end
  nb <- n.boundary.coeffs - ne     # number of coeffs at beginning
  abline(v=mean(times[nb:(nb+1)]), col="red")
  abline(v=mean(ox.times.0[(N-ne):(N-ne+1)]), col="red")
}
### plot MODWT wavlets coefficients
par(mfrow=c(4,1))
plot(ox.times.0, tV2s, typ="l")
add.boundary.markers(ox.times.0, 2, abs=abs.nuGj.LA8)
plot(ox.times.0, tW2s, typ="l")
add.boundary.markers(ox.times.0, 2)
plot(ox.times.0, tW1s, typ="l")
add.boundary.markers(ox.times.0, 1)
plot(ox.times.0, ox,typ="l")
###
par(mfrow=c(1,1))
stack.plot(data.frame(cbind(tV2s,tW2s,tW1s,ox)))
### note that we can extract DWT coeffs from MODWT coeffs
### (see overhead IV-12)
W2[1:5]
2*tW2[4*((0:4)+1)]

### DWT-based MRA
ox.dwt.mra.la8.2   <- mra(ox, "la8", 2, "dwt")
names(ox.dwt.mra.la8.2)
S2 <- ox.dwt.mra.la8.2$S2
D2 <- ox.dwt.mra.la8.2$D2
D1 <- ox.dwt.mra.la8.2$D1
ox.reconstructed <- D1 + D2 + S2
par(mfrow=c(3,1))
plot(ox.times.0, ox, typ="l")
plot(ox.times.0, ox.reconstructed, typ="l")
plot(ox.times.0, ox - ox.reconstructed, typ="l")
sum((ox - (D1 + D2 + S2))^2)
###
par(mfrow=c(4,1))
plot(ox.times.0, S2, typ="l")
plot(ox.times.0, D2, typ="l")
plot(ox.times.0, D1, typ="l")
plot(ox.times.0, ox, typ="l")
###
par(mfrow=c(1,1))
stack.plot(data.frame(cbind(S2,D2,D1,ox)))

### MODWT-based MRA
ox.modwt.mra.la8.2 <- mra(ox, "la8", 2, "modwt")
tS2 <- ox.modwt.mra.la8.2$S2
tD2 <- ox.modwt.mra.la8.2$D2
tD1 <- ox.modwt.mra.la8.2$D1
ox.reconstructed <- tD1 + tD2 + tS2
par(mfrow=c(3,1))
plot(ox.times.0, ox, typ="l")
plot(ox.times.0, ox.reconstructed, typ="l")
plot(ox.times.0, ox - ox.reconstructed, typ="l")
sum((ox - (tD1 + tD2 + tS2))^2)
###
par(mfrow=c(4,1))
plot(ox.times.0, tS2, typ="l")
plot(ox.times.0, tD2, typ="l")
plot(ox.times.0, tD1, typ="l")
plot(ox.times.0, ox, typ="l")
### compare DWT and MODWT MRAs
par(mfrow=c(3,1))
plot(ox.times.0, tS2, typ="l", col="red")
lines(ox.times.0, S2, typ="l", col="blue")
plot(ox.times.0, tD2, typ="l", col="red")
lines(ox.times.0, D2, typ="l", col="blue")
plot(ox.times.0, tD1, typ="l", col="red")
lines(ox.times.0, D1, typ="l", col="blue")

### a peak at some improvements ...
old.par <- par()
library(wavelets)
plot(dwt(ox, "la8", 2))
plot(modwt(ox, "la8", 2))
detach("package:wavelets")
par(old.par)