### Clean up the workspace
rm(list = ls())

## additional libraries
library(devtools)
library(splines)
library(quantreg)

## Load the onlineforecast package - this is the development version
load_all(".") # need to be under the package directory

# Load the data
D <- Dsolar
names(D)

# Short Data Anlysis
tstart <- as.POSIXct("2009-01-01", tz = "GMT")
tend <- as.POSIXct("2011-01-01", tz = "GMT")
D <- subset(D, in_range(tstart,D$t,tend))


plotmat <- matrix(1:3, nrow = 3, byrow=TRUE)
layout(plotmat, widths=c(1,1,1), heights=c(1,1,1))
par(mar=c(0, 3.5, 0.5, 1), lwd = 2, cex.lab = 1.5, cex.axis = 1.5)

plot(D$t, D$Ps, type = "l", bty = "l", xlab = "", ylab = "", axes = F)
axis(2, seq(0,ceiling(max(D$Ps)+0.2), by = 0.5), las = 1)

plot(D$t, lagvec(D$I$k1,1), type = "l", col = "gold", bty = "l", xlab = "", ylab = "", axes = F)
lines(D$t, D$Iobs, col = rgb(0, 0, 255, max = 255, alpha = 125))
axis(2, seq(0,ceiling(max(D$I$k1)+0.2), by = 100), las = 1)

plot(D$t, lagvec(D$Ta$k1,1), type = "l", col = rgb(240, 37, 37, max = 255, alpha = 150), bty = "l", xlab = "", ylab = "", axes = F, ylim = c(min(D$Ta$k1) - 4,max(D$Ta$k1)))
axis(2, seq(floor(min(D$Ta$k1)),ceiling(max(D$Ta$k1)+0.2), by = 10), las = 1)
axis.POSIXct(side = 1, x = D$t, xaxt = "s",
        at = seq(D$t[1], D$t[length(D$t)], by = "1 month"),
        format = "%Y-%m-%d", mgp = c(4,-1,-2))


# Make a training set of 3 months, and a then test set with 1 month
tstart <- "2009-06-01"
ttest <- "2009-03-31"
tend <- "2009-07-01"

# Cut only the necessary period
D$scoreperiod <- in_range("2009-06-10", D$t)
Dtrain <- subset(D, in_range(tstart,D$t,tend))



## Create the model, where N1 is the initial number of data points - for the cold start
model <- qmodel$new(N1 = 50)
## The output is the Ps from the Dsolar data
model$output = "Ps"
## The design matrix, here we are using 5 bsplines fixed at 5 and 18 hour of the day - ie when sun is up and down. 
model$add_inputs(I = "bspline(tday, df=5, Boundary.knots=c(5,18)) %**% I")

## Add the regularization parameter, lambda, and the bounds
model$add_regprm("rls_prm(lambda=0.9)")
model$add_prmbounds(lambda = c(0.99, 0.999, 0.9999))

### Setup of the booking keeping matrix - 
##############################################################
####### THIS SHOULD BE DONE IN quantile_fit.R ################
##############################################################
## First 5 columns are the input Matrix
model$IX <- 0:4
## The number of predictors, we could just use the length of the input matrix or the IX
model$K <- 5
## The output column is the 6th one
model$Iy <- 5

## Just quick fun to see the transformation
## maybe have splines that sum to one always..
model$datatr <- model$transform_data(Dtrain)
plot(Dtrain$I$k1[1:100], type = "l")
lines(model$datatr$I.bs1$k1, col = "red")
lines(model$datatr$I.bs2$k1, col = "red")
lines(model$datatr$I.bs3$k1, col = "red")
lines(model$datatr$I.bs4$k1, col = "red")
lines(model$datatr$I.bs5$k1, col = "red")



## Select which quantiles we want to optimise and fit
model$tau <- c(0.05,0.2, 0.5, 0.8,0.95)
## Select which prediction horizon
model$kseq <- c(1,4,8,12,24)
model$kseq <- 1:24

#opt_model <- quantile_optim(model = model, data = Dtrain)

PAR <- c("lambda" = 0.999)
quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[1])
PAR <- c("lambda" = 0.999)
quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[2])
PAR <- c("lambda" = 0.999)
quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[3])
PAR <- c("lambda" = 0.999)
quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[4])
PAR <- c("lambda" = 0.999)
quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[5])


Pred_model <- quantile_predict(model = model, datatr = model$datatr)


plot(Dtrain$Ps, type = "l")
lines(Pred_model$q0.05[,"k24"], type = "l", col = "green")
lines(Pred_model$q0.5[,"k24"], type = "l", col = "blue")
#lines(Pred_model$q0.5[,"k1"], type = "l", col = "green")
lines(Pred_model$q0.95[,"k24"], type = "l", col = "red")



str(Pred_model)

k <- 24
colSeq <- grey(seq(0.9,0.1,len=3))
#colSeq <- colorJet(nPolygons)

require( RColorBrewer )
colSeq <- brewer.pal(11 , "Spectral" )

for(i in k:(length(Dtrain$t)-1)){

    par(mfrow = c(2,1))
    par(mar = c(2,4,2,2)) # bottom left top right
    
    with(Pred_model, {

        plot(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], bty = "l", lwd = 2, col = "black", pch=19, cex=0.5,  axes = FALSE, xaxt = "n",
             type="n", ylim = range(Dtrain$Ps[(i-10):(i+k)], q0.5[i,],  q0.95[i,], q0.05[i,], na.rm = T), main = "TimeAdaptive with Weights", xlab = "Time", ylab = "Ps")
        
        axis(2)
        
        
        
        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
                c(Pred_model$q0.05[i,],rev(Pred_model$q0.2[i,])), col= colSeq[3], border = NA)
        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
                c(Pred_model$q0.8[i,],rev(Pred_model$q0.95[i,])), col= colSeq[3], border = NA)

        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
               c(Pred_model$q0.2[i,],rev(Pred_model$q0.5[i,])), col= colSeq[4], border = NA)
        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
               c(Pred_model$q0.5[i,],rev(Pred_model$q0.8[i,])), col= colSeq[4], border = NA)
        
        
        lines(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], col = "black", type = "b", lwd = 2)
        lines(Dtrain$t[(i+1):(i+k)], q0.5[i,], type = "b", col = "grey", lwd = 2)

        axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)], 
                     at = seq(from = Dtrain$t[(i-10):(i+k)][1],
                              to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])],
                              by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a",
                     las = 1, cex.axis = 1, srt = 45)
        
        #lines(Pred_model$q0.05[i,], type = "b", col = "blue")
        #lines(Pred_model$q0.95[i,], type = "b", col = "blue")

        plot(Dtrain$t[(i-10):(i+k)], c(rep(NA,11), as.numeric(Dtrain$I[i, 1:k])), type = "b", col = "steelblue", axes = FALSE, xlab = "Time", ylab = "Temp", lwd = 2, ylim = range(Dtrain$I[i, 1:k],  Dtrain$Iobs[(i+1):(i+k)]))
        #lines(Dtrain$t[i], Dtrain$Iobs[i+1,1], type = "b", col = "red", axes = FALSE)
        lines(Dtrain$t[(i+1):(i+k)], Dtrain$Iobs[(i+1):(i+k)], type = "b", col = "red", lwd = 2)
        axis(2)
        abline(h=0,v=Dtrain$t[i],lty=2, col = "lightgrey", lwd = 2)
        
        legend("topleft", legend=c("I Obs", "I Pred"),
               col=c("red", "steelblue"), lty=1, cex=1, lwd = 2)
    
        
i 
        axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)], 
                     at = seq(from = Dtrain$t[(i-10):(i+k)][1],
                              to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])],
                              by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a",
                     las = 1, cex.axis = 1, srt = 45)
    })


}













library(animation)
k <- 24
colSeq <- grey(seq(0.9,0.1,len=3))
#colSeq <- colorJet(nPolygons)

require( RColorBrewer )
colSeq <- brewer.pal(11 , "Spectral" )

saveHTML({
for(i in k:(60)){

    par(mfrow = c(2,1))
    par(mar = c(2,4,2,2)) # bottom left top right
    
    with(Pred_model, {

        plot(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], bty = "l", lwd = 2, col = "black", pch=19, cex=0.5,  axes = FALSE, xaxt = "n",
             type="n", ylim = range(Dtrain$Ps[(i-10):(i+k)], q0.5[i,],  q0.95[i,], q0.05[i,], na.rm = T), main = "TimeAdaptive with Weights", xlab = "Time", ylab = "Ps")
        
        axis(2)
        
        
        
        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
                c(Pred_model$q0.05[i,],rev(Pred_model$q0.2[i,])), col= colSeq[3], border = NA)
        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
                c(Pred_model$q0.8[i,],rev(Pred_model$q0.95[i,])), col= colSeq[3], border = NA)

        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
               c(Pred_model$q0.2[i,],rev(Pred_model$q0.5[i,])), col= colSeq[4], border = NA)
        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
               c(Pred_model$q0.5[i,],rev(Pred_model$q0.8[i,])), col= colSeq[4], border = NA)
        
        
        lines(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], col = "black", type = "b", lwd = 2)
        lines(Dtrain$t[(i+1):(i+k)], q0.5[i,], type = "b", col = "grey", lwd = 2)

        axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)], 
                     at = seq(from = Dtrain$t[(i-10):(i+k)][1],
                              to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])],
                              by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a",
                     las = 1, cex.axis = 1, srt = 45)
        
        #lines(Pred_model$q0.05[i,], type = "b", col = "blue")
        #lines(Pred_model$q0.95[i,], type = "b", col = "blue")

        plot(Dtrain$t[(i-10):(i+k)], c(rep(NA,11), as.numeric(Dtrain$I[i, 1:k])), type = "b", col = "steelblue", axes = FALSE, xlab = "Time", ylab = "Temp", lwd = 2, ylim = range(Dtrain$I[i, 1:k],  Dtrain$Iobs[(i+1):(i+k)]))
        #lines(Dtrain$t[i], Dtrain$Iobs[i+1,1], type = "b", col = "red", axes = FALSE)
        lines(Dtrain$t[(i+1):(i+k)], Dtrain$Iobs[(i+1):(i+k)], type = "b", col = "red", lwd = 2)
        axis(2)
        abline(h=0,v=Dtrain$t[i],lty=2, col = "lightgrey", lwd = 2)
        
        legend("topleft", legend=c("Ta Obs", "Ta Pred"),
               col=c("red", "steelblue"), lty=1, cex=1, lwd = 2)
    
        
i 
        axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)], 
                     at = seq(from = Dtrain$t[(i-10):(i+k)][1],
                              to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])],
                              by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a",
                     las = 1, cex.axis = 1, srt = 45)
    })


}},  htmlfile = "misc-R/quantile/model.html", ani.height = 800, ani.width = 1600)