
# The object data is a data set where each row represents a day. For each column in the row, it is represented the flow that was evaluated in that instant of time.
# The object day is a vector that contains the flow of a day evaluated in several instants of time (each 5 minutes).
# The commented instructions show the use of the functions by defect

###################
##### FIGURE 6 ####
###################

#rainbow plots
library(rainbow)

#fts(x, y, start = 1, frequency = 1, xname, yname)
f = fts(x=1:dim(data)[2], y=t(data))
plot(f, xlab = 'Time', ylab = 'Flow', cex.lab = 1.3)

##########################
##### FIGURE 7 (left) ####
##########################

# download the necessary functions to construct the tolerance bands.
#http://www.utdallas.edu/~pkc022000/functional_tolerance_band/TB_functions.R

# construct tolerance bands. 

# TB.calc(data, alpha = 0.05, prob = 0.90, simul = TRUE, n.boot)
result = TB.calc(as.matrix(data), alpha = .05, prob = 0.97, n.boot = 500, simul = T)
    

# plot the data with the tolerance bands
plot(x = 1:dim(data)[2], y = as.vector(data[1,]), xaxt = "n", xlab = "Time", ylab = "Flow", ylim = c(0, 300), main = "Tolerance bands")

for(j in 1:dim(data)[1]){
    lines(x = which(data[j,] != 'NA'), y = data[j,][which(data[j,] != 'NA')], type = 'o', col = "gray80")
}

lines(result[[3]][,1], type = 'l', lty = 2, col = 'red', lwd = 2)
lines(result[[3]][,2], type = 'l', lty = 2, col = 'red', lwd = 2)
lines(result[[4]][,1], type = 'l', lty = 1, col = 'blue', lwd = 2)
lines(result[[4]][,2], type = 'l', lty = 1, col = 'blue', lwd = 2)
lines(result[[5]][,1], type = 'l', lty = 3, col = 'black', lwd = 2)

legend("topright", legend = c("estimated mean", "pointwise band", "simultaneous band"), lty = c(3, 2, 1), lwd = c(2, 2, 2), col = c("black", "red", "blue"), cex = .75)

#################################################
##### FIGURE 7 (right) AND FIGURE 9 (right) ####
#################################################

# smooth the data

library(fda)

#create.bspline.basis(rangeval=NULL, nbasis=NULL, norder=4,
#      breaks=NULL, dropind=NULL, quadvals=NULL, values=NULL,
#      basisvalues=NULL, names="bspl")

bs = create.bspline.basis(rangeval = c(0,length(day)), nbasis = 15)

#Data2fd(argvals=NULL, y=NULL, basisobj=NULL, nderiv=NULL,
#        lambda=3e-8/diff(as.numeric(range(argvals))),
#        fdnames=NULL, covariates=NULL, method="chol",
#        dfscale=1)

aprxbspline = Data2fd(argvals = c(which(!is.na(day))), y = c(na.omit(day)), basisobj = bs) # to create a functional data object from data

#eval.fd(evalarg, fdobj, Lfdobj=0, returnMatrix=FALSE)

vec = c()
for( i in 1:length(day) ){
  vec = c(vec, eval.fd(evalarg = i , fdobj = aprxbspline)) # In the vector vec we have the smoothed data of a day
}

###################
##### FIGURE 8 ####
###################

basis = seq(4,22) 
variance = rep(0,19)

for(k in 1:dim(data)[1]){
	
	j = 0

	for(nbase in basis){

      n = 0
      
      mibspline = create.bspline.basis(rangeval = c(0,length(data[k,])), nbasis = nbase, norder = 4)  
      aprxbspline = Data2fd(argvals = c(which(!is.na(data[k,]))), y = c(na.omit(data[k,])), basisobj = mibspline)
      
      aux = 0
      for(i in 1:length(data[k,])){
        
        evalua = eval.fd(evalarg = i , fdobj = aprxbspline)
        
        if(!is.na(data[k,][i])){
          aux = aux + (data[k,][i] - evalua)^2
          n = n + 1
        }
      }  

      aux = aux/(n - basis)
      variance[j] = variance[j] + aux
      j = j+1  
    }
}
 
plot(variance/dim(data)[1], type = "b", xlab = "Number of basis", ylab = "Pooled variace", xaxt = "n",  cex.lab = 1.3, col = "blue")
axis(1, at = seq(1,19), labels = basis, las=2, cex.axis = 1)


###################
##### FIGURE 10 ####
###################

# outlier detection by FB

library(fda)

#fbplot(fit, x = NULL, method = "MBD", depth = NULL, plot = TRUE,
#	 prob = 0.5, color = 6, outliercol = 2, barcol = 4,
#	 fullout=FALSE, factor=1.5,xlim=c(1,nrow(fit)),
#	 ylim=c(min(fit)-.5*diff(range(fit)),max(fit)+.5*diff(range(fit))),...)

out = fbplot(t(data))
out$outpoint # it returns the numbers of the rows that are detected as outliers

#repeat if necessary (if outliers were found) with cleaned data