############
# SETTINGS #
############
 
rm(list=ls())
 
# Directory di lavoro
workdir <- "/Users/dariosolari/Documents/Jobs/Strategico"
 
# Library
require(tseries)
require(R2HTML)
require(hwriter)
require(forecast)
#require(ast)
 
###############
# IMPORT DATA #
###############
 
setwd(workdir)
 
model <- url("http://rospo.dyndns.org/strategico/projects/web-ltp/1200/600261/V1/model.RData")
 
###################
# OTHER FUNCTIONS #
###################
 
.incSampleTime <- function(now, period.freq = 2, increment = 1) {
	if (now[2] + increment - 1 <= period.freq - 1)
	now[2] = now[2] + increment
	else now = c(now[1] + (now[2] - 1 + increment)%/%period.freq, ((now[2] + increment - 1)%%period.freq) + 1)
	now
	}
 
############################################
# STRATEGICO CODE - .plot.best & .plot.all #
############################################
 
### .PLOT.BEST
 
.plot.best = function(best, plot.trend =TRUE, color.ic,
  color.forecast, fies.name, title,filename="modelBest.png", width = width, height = height) {
 
  y = best$ts.product
  period.freq = frequency(y)
  end_serie = end(y)
  period.start = start(y)
  start_pred = .incSampleTime(period.freq = period.freq, now = end_serie)
 
  pred = best$prediction
  ic.lwr = best$IC$lwr
  ic.upr = best$IC$upr
 
  if (!is.ts(pred)) pred = ts(pred,start=start_pred,frequency=period.freq)
  if (!is.ts(ic.lwr)) ic.lwr = ts(ic.lwr,start=start_pred,frequency=period.freq)
  if (!is.ts(ic.upr)) ic.upr = ts(ic.upr,start=start_pred,frequency=period.freq)
 
  p.best = append(as.vector(window(y,end=end_serie)),pred[1])
  p.best = ts(p.best,start=period.start,frequency=period.freq)
 
  y.best = append(as.vector(window(y,end=end_serie)),pred)
  y.best = ts(y.best,start=period.start,frequency=period.freq)
 
  inf = min(ic.lwr,y,na.rm = TRUE)
  sup = max(ic.upr,y,na.rm = TRUE)
 
  # ORIGINAL CODE
  #bitmap(units="px",filename, width = width, height = height)
  # MODIFIED CODE
  png(units="px",filename, width = width, height = height)
 
  plot(window(p.best, end = start_pred), ylim = c((inf - (inf/4)),
                                           (sup + (sup/2))), xlim = c(period.start[1], end(pred)[1]),
       main = title,ylab="y")
  lines(y.best, col = color.forecast[1], pch = "*", lwd = 2)
  lines(ic.lwr, col = color.ic, lwd = 1)
  lines(ic.upr, col = color.ic, lwd = 1)
  points(ic.upr, col = color.ic, cex = 1, pch = "*")
  points(ic.lwr, col = color.ic, cex = 1, pch = "*")
  lines(y, pch = "*", lwd = 2)
  legend(x = period.start[1], y = (sup + (sup/2)), legend = c("Prediction",
                                                     "Confidence band 95%"), col = c(color.forecast[1], color.ic),
         lty = 1, lwd = 2, horiz = FALSE, x.intersp = 1)
 
  if (plot.trend) {
    trend.best = try(smooth.spline(y.best),TRUE)
    if(!is(trend.best,"try-error"))  lines(trend.best, col = color.forecast[1], lwd = 1)
  }
  dev.off()
}
 
# .PLOT.ALL
 
.plot.all = function(model, color.forecast, plot.trend = TRUE, fies.name, title,filename="modelAll.png",width = width, height = height) {
  y = model[[model$BestModel]]$ts.product
  period.freq = frequency(y)
  end_serie = end(y)
  period.start = start(y)
  start_pred = .incSampleTime(period.freq = period.freq, now = end_serie)
 
  pred=list(pred.lm = model[["LinearModel"]]$prediction,
    pred.arima = model[["Arima"]]$prediction,
	pred.es = model[["ExponentialSmooth"]]$prediction,
    pred.trend = model[["Trend"]]$prediction,
    pred.mean = model[["Mean"]]$prediction)
  names(pred)=c("LinearModel","Arima","ExponentialSmooth","Trend","Mean")
 
  pred = lapply(pred, function(pr){ if (!is.null(pr))  if (!is.ts(pr)) pr = ts(pr, start = start_pred, frequency = period.freq); pr})
 
  p = lapply( pred, function(pr) {pr=append(as.vector(window(y, end = end_serie)), pr)
                                  pr= ts(pr, start = period.start, frequency = period.freq)})
 
  yy=list()
  for(i in which(sapply(p,function(yyy)!is.null(yyy) ))){
    if(!is.null(p[i])){
      yy[[i]]=append(as.vector(window(y, end = end_serie)), p[[i]])
      yy[[i]] = ts(yy[[i]], start = period.start, frequency = period.freq)
    }
  }
 
  inf = min(unlist(pred), y,na.rm = TRUE)
  sup = max(unlist(pred), y,na.rm = TRUE)
 
  #bmp(file=fies.name)
 
  # ORIGINAL CODE
  #bitmap(units="px",filename, width = width, height = height)
  # MODIFIED CODE
  png(units="px",filename, width = width, height = height)
  plot(y, ylim = c((inf - (inf/4)), (sup + (sup/2))), xlim = c(period.start[1], end(p[[model$BestModel]])[1]),
       lwd = 2, main = title,ylab="y")
 
   # ORIGINAL CODE
   #for(i in names(pred)[which(sapply(pred,function(pp)!is.null(pp) ))]){
   # MODIFIED CODE
   for(i in seq_along(names(pred))[which(sapply(pred,function(pp)!is.null(pp) ))]){
    lines(window(p[[i]], start = end_serie) , col = color.forecast[i], pch = "*", cex = 2, lwd = 2)
  }
 
legend(x = period.start[1], y = (sup + (sup/2)),legend = c("Linear Model","Arima" , "Exp. Smooth", "Trend" ,"Mean" )[which(sapply(pred,function(pp)!is.null(pp) ))],col = color.forecast[which(sapply(pred,function(pp)!is.null(pp) ))], lty = 1, lwd = 2, horiz = FALSE, y.intersp=1)  
 
  if (plot.trend) {
    for(i in which(sapply(pred,function(yyy)!is.null(yyy) ))){
      trend = try(smooth.spline(yy[[i]]),TRUE)
      if(!is(trend,"try-error")) lines(trend, pch = "*", col = color.forecast[i], lwd = 1)
    }
  }
dev.off()
  }
 
############
# OUTPUT 1 #
############
 
best <- model[[model$BestModel]]
 
.plot.best(best,filename="BEST_model.png",height=600,width=1500,title="BEST Model",color.forecast=c("darkblue","darkred","darkcyan","red","darkgreen","purple"),plot.trend=FALSE,color.ic="cyan")
 
############
# OUTPUT 2 #
############
 
.plot.all(model,filename="ALL_models.png",height=600,width=1500,title="ALL Models",color.forecast=c("darkblue","darkred","darkcyan","red","darkgreen","purple"),plot.trend=FALSE)