# data - time series observation
x <- c(40.3,40.7,38.5,37.9,38.6,41.1,45.2,45.7,46.7,46.5,
        45.2,45.1,45.8,46.3,47.5,48.5,49.1,51.7,50.6,48.0,
        44.7,41.2,40.0,40.3)
t <- 1:length(x)

# auxilliary functions to create design matrices F, V and the projection matrix M_F
makeF <- function(times, freqs) {

        n <- length(times)
        f <- length(freqs)
        c1 <- matrix(1, n)
        F <- matrix(c1, n, 1)

        for (i in 1:f) {
                F <- cbind(F, cos(2 * pi * freqs[i] * times))
                F <- cbind(F, sin(2 * pi * freqs[i] * times))
        }

        return(F)
}

makeV <- function(times, freqs) {

        V <- makeF(times, freqs)
        V <- V[, -1]

        return(V)
}

makeM_F <- function(F) {

        n <- nrow(F)
        c1 <- rep(1, n)
        I <- diag(c1)
        M_F <- I - F %*% solve((t(F) %*% F)) %*% t(F)

        return(M_F)
}

# model parameters
n <- 24 
k <- 3 
l <- 4

# model - design matrices F, V
F <- makeF(t, c(1/24))
V <- makeV(t, c(3/24, 4/24))

library(nlme)

REMLE_nlme <- function(X, F, V, method) {
    g <- rep(1,length(X))
    d <- data.frame(g,F,V,X)
    colnames(d) <- c("g","f1","f2","f3","v1","v2","v3","v4","X")
    
    result <- NULL
    
    if(method == "ML") {
        result <- lme(fixed=X~f2+f3,random=list(g=pdDiag(~-1+v1+v2+v3+v4)),data=d,method="ML")
    } else {
        result <- lme(fixed=X~f2+f3,random=list(g=pdDiag(~-1+v1+v2+v3+v4)),data=d,method="REML")
    }
    
    return(as.vector(c(result$sigma^2, diag(getVarCov(result)))))
}

# auxilliary function to calculate the norm of a vector
norm_vec <- function(x) sqrt(sum(x^2))

MLEnlme <- REMLE_nlme(x, F, V, "ML")
print(MLEnlme)
norm_vec(MLEnlme)

# loading all fdslrm functions as an R script from GiHub
devtools::source_url("https://github.com/fdslrm/fdslrmAllinOne/blob/master/fdslrmAllinOne.R?raw=TRUE")
initialFDSLRM()

# fit particular FDSLRM employing ML estimation of variances (by 'lme' function)
MLEfitnlme <- fitDiagFDSLRM(x, t, freq_mean = c(1/24), freq_random = c(3/24, 4/24), 
                            fit_function = "lme", var_estim_method = "ML")

print(as.vector(c(MLEfitnlme$error_variance, diag(MLEfitnlme$rand_eff_variance))))
norm_vec(as.vector(c(MLEfitnlme$error_variance, diag(MLEfitnlme$rand_eff_variance))))

# load function 'MMEinR'
devtools::source_url("https://github.com/fdslrm/MMEinR/blob/master/MMEinR.R?raw=TRUE")

MLEmixed <- mixed(x, F, V, c(1,1,1,1), c(1,1,1,1,1), 1)
print(MLEmixed$s2)
print(norm_vec(MLEmixed$s2))

# fit particular FDSLRM employing ML estimation of variances (by 'MMEinR (mixed)' function)
MLEfitmixed <- fitDiagFDSLRM(x, t, freq_mean = c(1/24), freq_random = c(3/24, 4/24), 
                             fit_function = "mixed", var_estim_method = "ML")

print(as.vector(c(MLEfitmixed$error_variance, MLEfitmixed$rand_eff_variance)))
norm_vec(as.vector(c(MLEfitmixed$error_variance, MLEfitmixed$rand_eff_variance)))

REMLEnlme <- REMLE_nlme(x, F, V, "REML")
print(REMLEnlme)
norm_vec(REMLEnlme)

# fit particular FDSLRM employing REML estimation of variances (by 'lme' function)
REMLEfitnlme <- fitDiagFDSLRM(x, t, freq_mean = c(1/24), freq_random = c(3/24, 4/24), 
                              fit_function = "lme", var_estim_method = "REML")

print(as.vector(c(REMLEfitnlme$error_variance, diag(REMLEfitnlme$rand_eff_variance))))
norm_vec(as.vector(c(REMLEfitnlme$error_variance, diag(REMLEfitnlme$rand_eff_variance))))

REMLEmixed <- mixed(x, F, V, c(1,1,1,1), c(1,1,1,1,1), 2)
print(REMLEmixed$s2)
print(norm_vec(REMLEmixed$s2))

# fit particular FDSLRM employing REML estimation of variances (by 'MMEinR (mixed)' function)
REMLEfitmixed <- fitDiagFDSLRM(x, t, freq_mean = c(1/24), freq_random = c(3/24, 4/24), 
                               fit_function = "mixed", var_estim_method = "REML")

print(as.vector(c(REMLEfitmixed$error_variance, REMLEfitmixed$rand_eff_variance)))
norm_vec(as.vector(c(REMLEfitmixed$error_variance, REMLEfitmixed$rand_eff_variance)))

library(sommer)

REMLE_sommer <- function(X, F, V) {
    f1 <- F[,1]
    f2 <- F[,2]
    f3 <- F[,3]
    v1 <- V[,1]
    v2 <- V[,2]
    v3 <- V[,3]
    v4 <- V[,4]
    
    DT <- data.frame(X, f1, f2, f3, v1, v2, v3, v4)

    suppressWarnings(result_new <- mmer(fixed = X~f2+f3, random = ~ vs(ds(v1),1)+vs(ds(v2),1)+vs(ds(v3),1)+vs(ds(v4),1), 
                                        data = DT, verbose = FALSE))
    
    output <- as.vector(unlist(result_new$sigma))
    
    return(output)

}

REMLEsommer <- REMLE_sommer(x, F, V)
print(REMLEsommer)
print(norm_vec(REMLEsommer))

# fit particular FDSLRM employing REML estimation of variances (by 'mmer' function)
REMLEfitsommer <- suppressWarnings(fitDiagFDSLRM(x, t, freq_mean = c(1/24), freq_random = c(3/24, 4/24), 
                                                 fit_function = "mmer", var_estim_method = "REML"))

print(as.vector(c(REMLEfitsommer$error_variance, REMLEfitsommer$rand_eff_variance)))
print(norm_vec(c(REMLEfitsommer$error_variance, REMLEfitsommer$rand_eff_variance)))