# R/generation_kernel_periodic.R In ardeeshany/FLAME: FLAME for high-dimensional Function-on-Scalar regression problems

#### Documented in generation_kernel_periodic

#' Definition of the eigenfunctions and eigenvalues of the periodic kernel
#'
#' @description Given the period \code{p} and the smoothing parameter \eqn{\sigma}, it returns the
#' evaluation of the eigenfunctions of the kernel on the grid \code{domain}
#' and the correspondent eigenvalues.
#'
#' @param period scalar. Period of the kernel.
#' @param parameter scalar. Parameter to tune the smoothness level of the kernel.
#' The \eqn{\sigma} parameter introduced in Details.
#' @param domain vector. \code{m}-length vector for the abscissa grid
#' of the kernel.
#' @param thres scalar. Threshold for the identification of the significant
#' eigenvalues of the kernel. The number of significant eigennvalues \code{J} is
#' the minimum \eqn{J} s.t.
#' \deqn{
#' \sum_{j = 1}^J \theta_j \geq \textrm{thres} \sum_{j = 1}^{\infty} \theta_j.
#' } Default is 0.99.
#' @param return.derivatives bool. If \code{TRUE} the function returns the matrix of
#' the derivatives of the selected eigenfunctions evaluated on the time domain.
#'  Default is \code{FALSE}.
#'
#' @return list containing
#'  \itemize{
#'  \item \code{eigenvect} \code{m} \eqn{\times} \code{J} matrix of
#'  the eigenfunctions of the kernel evaluated on the \code{domain}.
#'  \item \code{eigenval} \code{J}-length  vector of the
#'  eigenvalues of the kernel
#'  \item \code{derivatives}. if \code{return.derivatives = TRUE}.
#'  \code{m-1} \eqn{\times} \code{J} matrix of the derivatives of
#'  the eigenfunctions.
#'}
#'
#' @details The periodic kernel of period \eqn{p} defined in this function is
#' \deqn{
#' K(x, y) = \sigma^2 \exp \left\{  -2/\sigma \sin^2\left(\frac{\pi | x - y|}{p} \right)\right\}.
#' }
#' where \eqn{\sigma} is the smoothing parameter.
#' @export
#'
#' @examples
#' param_kernel <- 8
#' T_domain <- seq(0, 1, length = 50)
#' kernel_here <- generation_kernel_periodic(period = 1/2,
#'                                           parameter = param_kernel,
#'                                           domain = T_domain,
#'                                           thres = 1-10^{-16},
#'                                           return.derivatives = TRUE)
#' names(kernel_here)
#'
generation_kernel_periodic <- function(period = NULL, parameter = NULL, domain, thres = 0.99,
return.derivatives = FALSE)
{

## [email protected]: Shouldn't it be M_integ <- (length(domain)-1)/diff(range(domain))?
M_integ <- length(domain)/diff(range(domain))

if (length(period) != 1 )
{
stop ("please provide the period of the kernel.") # period is as proportion of the domain
}

if (length(parameter) != 1 )
{
stop ("please provide the parameters of the kernel. See the help page
of the function for the parameter definition")
}

dist_periodic <- function(x, y, sigma, p)
{
sigma^2* exp(- (2*sin(pi*abs(x-y)/(p))^2)/(sigma))
}

## @A: difference between apply, lapply and sapply:
## apply  - When you want to apply a function to the rows or columns of a matrix
## lapply - When you want to apply a function to each element of a list in turn and get a list back.
## sapply - When you want to apply a function to each element of a list in turn, but you want a vector back, rather than a list.

generate_matrix <- function(domain, sigma, p){
## @A: it means we apply the function(x) on each element of "domain" but the output is a vector.
sapply(domain, function(x){
sapply(domain, function(y) {dist_periodic(x,y, sigma,p)})
})
}

## @A: the following is another way to define generate_matrix but takes more time
## generate_matrix=function(domain,sigma,p){
#      G2=matrix(NA,length(domain),length(domain))
#          for(i in 1:length(domain)){
#          for(j in 1:length(domain)){
#            G2[i,j]=  dist_periodic(domain[i],domain[j], sigma,p)
#          }
#      }
# return(G2)
#       }

kernel_matrix <- generate_matrix(domain, sigma = parameter, p = period)

kernel_def <- eigen(kernel_matrix)

## [email protected]: Again, shouldn't it be kernel_def$vectors <- kernel_def$vectors
kernel_def$values <- kernel_def$values/M_integ
kernel_def$vectors <- kernel_def$vectors*sqrt(M_integ)

# isolate the number of significant eigenvalues
num_eigen <- which((cumsum(kernel_def$values)/sum(kernel_def$values))>thres)[1]
eigen_chosen <- 1:num_eigen

# defintion of the eigenvalues and eigenvectros of the kernel

autoval <- kernel_def$values[eigen_chosen] autovett <- kernel_def$vectors[,eigen_chosen]

if (return.derivatives == TRUE)
{
diff_autovett <- apply(autovett, 2, diff)/(domain[2]-domain[1])
return(list(eigenvect = autovett, eigenval = autoval,
derivatives = diff_autovett))
}else{
return(list(eigenvect = autovett, eigenval = autoval))
}

}

ardeeshany/FLAME documentation built on Sept. 25, 2017, 9:55 a.m.