DE.heat.FEM.time {fdaPDE}

Spatio-temporal density initialization

Description

This function implements two methods for the density initialization procedure.

Usage

DE.heat.FEM.time(data, data_time, FEMbasis, mesh_time, lambda=NULL,
                            lambda_time=NULL, heatStep=0.1, heatIter=10, 
                            init="Heat", nFolds=5, search="tree", 
                            isTimeDiscrete=FALSE, flagMass=FALSE, flagLumped=FALSE)

Arguments

data

A matrix of dimensions #observations-by-ndim. Data are locations: each row corresponds to one point, the first column corresponds to the x-coordinates, the second column corresponds to the y-coordinates and, if ndim=3, the third column corresponds to the z-coordinates.

data_time

A vector of dimensions #observations. The i-th datum is the time instant during which the i-th location is observed (according to the order in which data are provided).

FEMbasis

A FEMbasis object describing the Finite Element basis, as created by create.FEM.basis.

mesh_time

A vector containing the b-splines knots for separable smoothing. It is the time mesh of the considered time domain (interval). Its nodes are in increasing order.

lambda

A scalar or vector of smoothing parameters in space. Default is NULL. It is useful only if init='Heat'.

lambda_time

A scalar or vector of smoothing parameters in time. Default is NULL. It is useful only if init='Heat'.

heatStep

A real specifying the time step for the discretized heat diffusion process.

heatIter

An integer specifying the number of iterations to perform the discretized heat diffusion process.

init

A string specifying the initialization procedure. It can be either 'Heat' or 'CV'.

nFolds

An integer specifying the number of folds used in the cross validation technique. It is useful only for the case init = 'CV'.

search

A flag to decide the search algorithm type (tree or naive or walking search algorithm).

isTimeDiscrete

A boolean specifying the time data type: TRUE for discrete (with many duplicates) time data; FALSE for continuous time data. Default is FALSE.

flagMass

A boolean specifying whether to consider full mass matrices (TRUE) or identity mass matrices (FALSE) for the computation of space and time penalty matrices. Default is FALSE.

flagLumped

A boolean specifying whether to perform mass lumping. This numerical technique presents computational advantages during the procedure involving a mass matrix inversion for the computation of the space penalty matrix. Default is FALSE. N.B. We suggest to put it TRUE in case of a large spatial domain or in case of a dense/refined spatial mesh.

Value

If init = 'Heat' it returns a matrix in which each column contains the initial vector for each possible pair (lambda, lambda_time). If init = 'CV' it returns the initial vector associated to the unique pair (lambda, lambda_time) given.

Examples

library(fdaPDE)

## Create a 2D mesh over a squared domain
Xbound <- seq(-3, 3, length.out = 10)
Ybound <- seq(-3, 3, length.out = 10)
grid_XY <- expand.grid(Xbound, Ybound)
Bounds <- grid_XY[(grid_XY$Var1 %in% c(-3, 3)) | (grid_XY$Var2 %in% c(-3, 3)), ]
mesh <- create.mesh.2D(nodes = Bounds, order = 1)
mesh <- refine.mesh.2D(mesh, maximum_area = 0.25, minimum_angle = 20)
FEMbasis <- create.FEM.basis(mesh)

## Create a 1D time mesh over a (non-negative) interval
mesh_time <- seq(0, 1, length.out=3)

## Generate data
n <- 50
set.seed(10)
x <- rnorm(n,0,2)
y <- rnorm(n,0,2)
locations <- cbind(x,y)
times <- runif(n,0,1)
data <- cbind(locations, times)

t <- 0.5 # time instant in which to evaluate the solution

plot(mesh)
sample <- data[abs(data[,3]-t)<0.05,1:2]
points(sample, col="red", pch=19, cex=1, main=paste('Sample | ', t-0.05,' < t < ', t+0.05))

## Density initialization
lambda = 0.1
lambda_time <- 0.001
sol = DE.heat.FEM.time(data = locations, data_time = times, FEMbasis = FEMbasis,
                       mesh_time = mesh_time, lambda = lambda, lambda_time = lambda_time,
                       heatStep=0.1, heatIter=10, init="Heat")

## Visualization

n = 100
X <- seq(-3, 3, length.out = n)
Y <- seq(-3, 3, length.out = n)
grid <- expand.grid(X, Y)

FEMfunction = FEM.time(sol$f_init[,1,1], mesh_time, FEMbasis, FLAG_PARABOLIC = FALSE)
evaluation <- eval.FEM.time(FEM.time = FEMfunction, locations = grid, time.instants = t)
image2D(x = X, y = Y, z = matrix(evaluation, n, n), col = heat.colors(100),
        xlab = "", ylab = "", contour = list(drawlabels = FALSE),
        main = paste("Initial guess at t = ", t), zlim=c(0,0.2), asp = 1)