source("conv.r")

fraser=scan("fraser-river.dat",nmax=36)
data=list(s=c(1:length(fraser)),y=log(fraser))

# examples of changing kernel spacing
par(mfrow=c(3,1))
w3 <- as.matrix(seq(0,length(fraser)+2,by=2.0))
data.conv.gp1 <- conv.gp.1d(data, w=w3, sdkern=4, iter=1000)
plot(data$s,data$y,pch="*",xlab="s",ylab="z")
lines(data.conv.gp1$s,apply(data.conv.gp1$a$z[501:(dim(data.conv.gp1$a$z))[1],],2,mean),lwd=2)

w3 <- as.matrix(seq(0,length(fraser)+3,by=3.0))
data.conv.gp2 <- conv.gp.1d(data, w=w3, sdkern=4, iter=1500)
plot(data$s,data$y,pch="*",xlab="s",ylab="z")
lines(data.conv.gp2$s,apply(data.conv.gp2$a$z[501:(dim(data.conv.gp2$a$z))[1],],2,mean),lwd=2)

w3 <- as.matrix(seq(0,length(fraser)+4,by=4.0))
data.conv.gp3 <- conv.gp.1d(data, w=w3, sdkern=4, iter=1500)
plot(data$s,data$y,pch="*",xlab="s",ylab="z")
lines(data.conv.gp3$s,apply(data.conv.gp3$a$z[501:(dim(data.conv.gp3$a$z))[1],],2,mean),lwd=2)


w3 <- as.matrix(seq(0,length(fraser)+4,by=5))
data.conv.gp4 <- conv.gp.1d(data, w=w3, sdkern=4, iter=1500)
plot(data$s,data$y,pch="*",xlab="s",ylab="z")
lines(data.conv.gp4$s,apply(data.conv.gp4$a$z[501:(dim(data.conv.gp4$a$z))[1],],2,mean),lwd=2)

w3 <- as.matrix(seq(0,length(fraser)+6,by=6))
data.conv.gp5 <- conv.gp.1d(data, w=w3, sdkern=4, iter=1500)
plot(data$s,data$y,pch="*",xlab="s",ylab="z")
lines(data.conv.gp5$s,apply(data.conv.gp5$a$z[501:(dim(data.conv.gp5$a$z))[1],],2,mean),lwd=2)

w3 <- as.matrix(seq(0,length(fraser)+6,by=7))
data.conv.gp6 <- conv.gp.1d(data, w=w3, sdkern=4, iter=1500)
plot(data$s,data$y,pch="*",xlab="s",ylab="z")
lines(data.conv.gp6$s,apply(data.conv.gp6$a$z[501:(dim(data.conv.gp6$a$z))[1],],2,mean),lwd=2)


# multiscale fit to Fraser River data

m1 <- c(11,41) 
w1 <- list()  # (must be a list)
for(i in 1:2) w1[[i]] <- seq(0,40,length=m1[i])
sdkern <- c(4,1) # kernels of size 4 and 1

# estimate the 1-d ML multiscale convolution model 
data.conv.gp.multi <- conv.gp.multi.1d(data, w=w1, sdkern=sdkern, iter=1000)

#plot multiscale fit
par(mfrow=c(3,1))
plot(data$s,data$y,pch="*",xlab="s",ylab="z")
lines(data.conv.gp.multi$s,apply(data.conv.gp.multi$a$z[501:(dim(data.conv.gp.multi$a$z))[1],],2,mean),lwd=2)

#plot separate components of multiscale fit
partial.fit=c(rep(mean(data.conv.gp.multi$a$x[[1]][501:1000,1]),length(data$s)))
for(i in 1:length(w1[[1]])) {
  partial.fit=partial.fit+mean(data.conv.gp.multi$a$x[[1]][501:1000,i+1])*data.conv.gp.multi$df$K[[1]][,i+1]
}
plot(data$s,data$y,pch="*",xlab="s",ylab="z")
lines(data$s,partial.fit)
partial.fit=c(rep(0,length(data$s)))
for(i in 1:length(w1[[2]])) {
  partial.fit=partial.fit+mean(data.conv.gp.multi$a$x[[2]][501:1000,i])*data.conv.gp.multi$df$K[[2]][,i]
}
lines(data$s,partial.fit+mean(data$y),type="l")

#compare to standard convolution with size 2 kernels
w3 <- as.matrix(seq(0,length(fraser)+2,by=2))
data.conv.gp7 <- conv.gp.1d(data, w=w3, sdkern=2, iter=1000)
plot(data$s,data$y,pch="*",xlab="s",ylab="z")
lines(data.conv.gp7$s,apply(data.conv.gp7$a$z[501:(dim(data.conv.gp7$a$z))[1],],2,mean),lwd=2)


# multiscale synthetic example

s2=seq(0.5,40,by=.5)
s2true=seq(0.5,40,by=.25)
y2true=sin(2*pi*s2true/20)+0.25*cos(2*pi*s2true/4)
y2=sin(2*pi*s2/20)+0.25*cos(2*pi*s2/4)+rnorm(80,sd=0.1)
data2=list(s=s2,y=y2)

plot(data2$s,data2$y,pch="*",xlab="s",ylab="z")

#standard convolutions
par(mfrow=c(3,1))
w4 <- as.matrix(seq(-2,43,by=5.0))
data2.conv.gp1 <- conv.gp.1d(data2, pdata=s2true, w=w4, sdkern=5, iter=1000)
plot(data2$s,data2$y,pch="*",xlab="s",ylab="z")
lines(s2true,y2true,col=8,lwd=2)
lines(s2true,apply(data2.conv.gp1$a$z.pred[501:(dim(data2.conv.gp1$a$z.pred))[1],],2,mean),lwd=2)

w4 <- as.matrix(seq(0,42,by=2.0))
data2.conv.gp2 <- conv.gp.1d(data2, pdata=s2true, w=w4, sdkern=2, iter=1000)
plot(data2$s,data2$y,pch="*",xlab="s",ylab="z")
lines(s2true,y2true,col=8,lwd=2)
lines(s2true,apply(data2.conv.gp2$a$z.pred[501:(dim(data2.conv.gp2$a$z.pred))[1],],2,mean),lwd=2)

w4 <- as.matrix(seq(0,42,by=1.0))
data2.conv.gp3 <- conv.gp.1d(data2, pdata=s2true, w=w4, sdkern=1, iter=1000)
plot(data2$s,data2$y,pch="*",xlab="s",ylab="z")
lines(s2true,y2true,col=8,lwd=2)
lines(s2true,apply(data2.conv.gp3$a$z.pred[501:(dim(data2.conv.gp3$a$z.pred))[1],],2,mean),lwd=2)

#standard convolution (again, now with kernels plotted)
w4 <- as.matrix(seq(0,42,by=1.0))
plot(data2$s,data2$y,pch="*",xlab="s",ylab="z")
lines(s2true,y2true,col=8,lwd=2)
lines(s2true,apply(data2.conv.gp3$a$z.pred[501:(dim(data2.conv.gp3$a$z.pred))[1],],2,mean),lwd=2)
for(i in 1:length(w4)) {
  lines(s2true,mean(data2.conv.gp3$a$x[501:1000,i+1])*data2.conv.gp3$df$K.pred[,i+1],col=8,lty=3)
  lines(c(w4[i],w4[i]),c(0,mean(data2.conv.gp3$a$x[501:1000,i+1])*dnorm(0,sd=1)),col=8,lty=2)
}

#multiscale convolutions with kernels of size 5 and 1
m5 <- c(9,41) 
w5 <- list()
for(i in 1:2) w5[[i]] <- seq(0,40,length=m5[i])
sdkern <- c(5,1)
data2.conv.gp.multi1 <- conv.gp.multi.1d(data2, pdata=s2true, w=w5, sdkern=sdkern, iter=1000)
plot(data2$s,data2$y,pch="*",xlab="s",ylab="z")
lines(s2true,y2true,col=8,lwd=2)
lines(s2true,apply(data2.conv.gp.multi1$a$z.pred[501:(dim(data2.conv.gp.multi1$a$z.pred))[1],],2,mean),lwd=2)
for(j in 1:2) {
  for(i in 1:length(w5[[j]])) {
    lines(s2true,mean(data2.conv.gp.multi1$a$x[[j]][501:1000,i+(2-j)])*data2.conv.gp.multi1$df$K.pred[[j]][,i+(2-j)],col=8,lty=j+1)
    lines(c(w5[[j]][i],w5[[j]][i]),c(0,mean(data2.conv.gp.multi1$a$x[[j]][501:1000,i+(2-j)])*dnorm(0,sd=sdkern[j])),col=8,lty=j+1)
    # (2-j) is because there is an intercept on the first level only
  }
}

#multiscale convolution with kernels of size 3 and 1
m5 <- c(13,41) 
w5 <- list()
for(i in 1:2) w5[[i]] <- seq(0,40,length=m5[i])
sdkern <- c(3,1)
data2.conv.gp.multi2 <- conv.gp.multi.1d(data2, pdata=s2true, w=w5, sdkern=sdkern, iter=1000)
plot(data2$s,data2$y,pch="*",xlab="s",ylab="z")
lines(s2true,y2true,col=8,lwd=2)
lines(s2true,apply(data2.conv.gp.multi1$a$z.pred[501:(dim(data2.conv.gp.multi1$a$z.pred))[1],],2,mean),lwd=2)


plot.1d(data2.conv.gp.multi1, main="1-d Multiple Resolution")
plot.1d(data2.conv.gp.multi2, main="1-d Multiple Resolution")


for(j in 1:2) {
  for(i in 1:length(w5[[j]])) {
    lines(s2true,mean(data2.conv.gp.multi2$a$x[[j]][501:1000,i+(2-j)])*data2.conv.gp.multi2$df$K.pred[[j]][,i+(2-j)],col=8,lty=j+1)
    lines(c(w5[[j]][i],w5[[j]][i]),c(0,mean(data2.conv.gp.multi2$a$x[[j]][501:1000,i+(2-j)])*dnorm(0,sd=sdkern[j])),col=8,lty=j+1)
  }
}