library(fda)
library(ggplot2)

# Canadian Weather
data(CanadianWeather)
temp <- CanadianWeather$dailyAv[,,1]
precip <- CanadianWeather$dailyAv[,,2]

daytime <- (1:365) - 0.5

df <- data.frame(daytime = daytime, temp = c(temp), precip = c(precip),
                 place = rep(CanadianWeather$place, each = 365), region = rep(CanadianWeather$region, each = 365))

dfp <- subset(df, place == "St. Johns")


# Define cubic B-splines with full basis

rangeval <- range(dfp$daytime)
breaks <- dfp$daytime
norder <- 4

bbasis <- create.bspline.basis (rangeval, norder = norder, breaks = breaks)

# We'll concentrate on B-splines and second-derivative penalties.
curv.Lfd <- int2Lfd(2)
# Set smoothing parameter
lambda <- 10^(-1)

# Now we can define the fdPar object 
curv.fdPar <- fdPar(bbasis, curv.Lfd, lambda)

precipSmooth <- smooth.basis(dfp$daytime, dfp$precip, curv.fdPar)

### Plot smooth data
plotfit.fd(dfp$precip, dfp$daytime, precipSmooth$fd)

# ggplot
# fdobjSmootheval <- eval.fd(fdobj = precipSmooth$fd, evalarg = dfp$daytime)
# dfp$precipsm <- fdobjSmootheval
# 
# x11()
# p0 <- ggplot (data = dfp, aes (x = daytime, y = precip, color = place))
# p0 <- p0 + geom_point ()
# p0 <- p0 + geom_line (aes(y = precipsm), colour = 4)
# p0 <- p0 + labs (x = "Days", y = "Precipitation") 
# p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns")
# print (p0)
# 
# dev.copy2pdf (file = "precip_estBspl.pdf", width = 5, height = 4)


# The Fourier basis

fbasis <- create.fourier.basis(rangeval, nbasis = 365)
# Set smoothing parameter
lambda <- 10^(-1)
curv.fdPar <- fdPar(fbasis, curv.Lfd, lambda)
precipSmoothF <- smooth.basis(dfp$daytime, dfp$precip, curv.fdPar)

# Plot smooth data
plotfit.fd(dfp$precip, dfp$daytime, precipSmoothF$fd)

# ggplot
# fdobjSmootheval <- eval.fd(fdobj = precipSmoothF$fd, evalarg = dfp$daytime)
# dfp$precipsmF <- fdobjSmootheval
# 
# x11()
# p0 <- ggplot (data = dfp, aes (x = daytime, y = precip, color = place))
# p0 <- p0 + geom_point ()
# p0 <- p0 + geom_line (aes(y = precipsmF), colour = 4)
# p0 <- p0 + labs (x = "Days", y = "Precipitation") 
# p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns")
# print (p0)
# 

# Degrees of Freedom

degfre <- precipSmoothF$df

# Harmonic acceleration

fbasis <- create.fourier.basis(rangeval, nbasis = 365)

# We'll concentrate on Fourier basis and harmonic penalties.
omega <- 2*pi/365
accelvec <- c(0, omega^2, 0, 1)  # vector of coefficients, m-th is by Dx^(m-1)
harmLfd <- vec2Lfd (bwtvec = accelvec, rangeval = rangeval)

lambda <- 10^(-1)
curv.fdPar <- fdPar(fbasis, harmLfd, lambda)

FSmooth <- smooth.basis(dfp$daytime, dfp$precip, curv.fdPar)

# Plot smooth data
plotfit.fd(dfp$precip, dfp$daytime, FSmooth$fd)


# ggplot
# fdobjSmootheval <- eval.fd(fdobj = FSmooth$fd, evalarg = dfp$daytime)
# dfacc <- data.frame(x = dfp$daytime, y = HarmAccel)
# 
# x11()
# p0 <- ggplot (data = dfp, aes (x = daytime, y = precip, color = place))
# p0 <- p0 + geom_point ()
# p0 <- p0 + geom_line (aes(y = fdobjSmootheval), colour = 4)
# p0 <- p0 + labs (x = "Days", y = "Precipitation")
# p0 <- p0 + theme_bw () + guides (colour = FALSE) + ggtitle("St. Johns")
# print (p0)


# Cross - validation

Rmat <- eval.penalty(fbasis, harmLfd)

phi <- t(FSmooth$y2cMap)
phiinv <- solve(t(phi)%*%phi + lambda*Rmat)
# Smoothing matrix S
S <- phi%*%phiinv%*%t(phi)
sii <- diag(S)

# Fitted values
yhat <- S%*%dfp$precip
  
OCV <- mean(((dfp$precip - yhat)/(1-sii))^2)
GCV <- FSmooth$gcv

# ------------------
# Problems
# ------------------

# 1. Melanoma Data
?melanoma

melanoma <- data.frame(year = melanoma[, 2], incidence = melanoma[, 3])

x <- melanoma$year
y <- melanoma$incidence

# 2. Canadian Weather Data
?CanadianWeather

temp <- CanadianWeather$dailyAv[,,1]
daytime <- (1:365) - 0.5
rangeval <- range(daytime)

places <- c("Halifax", "Montreal", "Edmonton", "Ottawa")
temp.4 <- temp[, places]