Nothing
R/smoothing_time.R
In fdaPDE: Physics-Informed Spatial and Functional Data Analysis
Defines functions
smooth.FEM.time
Documented in
smooth.FEM.time
#' @useDynLib fdaPDE
#' @import Matrix plot3D rgl
#' @importFrom grDevices heat.colors palette
#' @importFrom graphics plot segments points lines
#' @importFrom methods is
NULL
#' Space-time regression with differential regularization
#' @param locations A matrix where each row specifies the spatial coordinates \code{x} and \code{y} (and \code{z} if ndim=3) of the corresponding observations in the vector \code{observations}.
#' This parameter can be \code{NULL}. In this case, if also the incidence matrix is \code{NULL} the spatial coordinates are assumed to coincide with the nodes of the \code{mesh}.
#' @param time_locations A vector containing the times of the corresponding observations in the vector \code{observations}.
#' This parameter can be \code{NULL}. In this case the temporal locations are assumed to coincide with the nodes of the \code{time_mesh}.
#' @param observations A matrix of #locations x #time_locations with the observed data values over the spatio-temporal domain.
#' The spatial locations of the observations can be specified with the \code{locations} argument.
#' @param FEMbasis A \code{FEMbasis} object describing the Finite Element basis, as created by \code{\link{create.FEM.basis}}.
#' @param time_mesh A vector specifying the time mesh.
#' @param covariates A #observations-by-#covariates matrix where each row represents the covariates associated with the corresponding observed data value in \code{observations}.
#' @param PDE_parameters A list specifying the parameters of the PDE in the regularizing term. Default is NULL, i.e. regularization is by means of the Laplacian (stationary, isotropic case).
#' If the PDE is elliptic it must contain: \code{K}, a 2-by-2 matrix of diffusion coefficients. This induces an anisotropic
#' smoothing with a preferential direction that corresponds to the first eigenvector of the diffusion matrix K; \code{b}, a vector of length 2 of advection coefficients. This induces a
#' smoothing only in the direction specified by the vector \code{b}; \code{c}, a scalar reaction coefficient. \code{c} induces a shrinkage of the surface to zero
#' If the PDE is space-varying it must contain: \code{K}, a function that for each spatial location in the spatial domain
#' (indicated by the vector of the 2 spatial coordinates) returns a 2-by-2 matrix of diffusion coefficients. This induces an anisotropic
#' smoothing with a local preferential direction that corresponds to the first eigenvector of the diffusion matrix K.The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be
#' an array with dimensions 2-by-2-by-#points.\code{b}, a function that for each spatial location in the spatial domain returns
#' a vector of length 2 of transport coefficients. This induces a local smoothing only in the direction specified by the vector \code{b}. The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be
#' a matrix with dimensions 2-by-#points; \code{c}, a function that for each spatial location in the spatial domain returns a scalar reaction coefficient.
#' \code{c} induces a shrinkage of the surface to zero. The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be
#' a vector with length #points; \code{u}, a function that for each spatial location in the spatial domain returns a scalar reaction coefficient.
#' \code{u} induces a reaction effect. The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be
#' a vector with length #points.
#' For 2.5D and 3D only the Laplacian is available (\code{PDE_parameters=NULL})
#' @param BC A list with two vectors:
#' \code{BC_indices}, a vector with the indices in \code{nodes} of boundary nodes where a Dirichlet Boundary Condition should be applied;
#' \code{BC_values}, a vector with the values that the spatial field must take at the nodes indicated in \code{BC_indices}.
#' @param incidence_matrix A #regions-by-#triangles/tetrahedrons matrix where the element (i,j) equals 1 if the j-th triangle/tetrahedron is in the i-th region and 0 otherwise.
#' This is only for areal data. In case of pointwise data, this parameter is set to \code{NULL}.
#' @param areal.data.avg Boolean. It involves the computation of Areal Data. If \code{TRUE} the areal data are averaged, otherwise not.
#' @param FLAG_MASS Boolean. This parameter is considered only for separable problems i.e. when \code{FLAG_PARABOLIC==FALSE}. If \code{TRUE} the mass matrix in space and in time are used, if \code{FALSE} they are substituted with proper identity matrices.
#' @param FLAG_PARABOLIC Boolean. If \code{TRUE} the parabolic problem problem is selected, if \code{FALSE} the separable one.
#' @param FLAG_ITERATIVE Boolean. If \code{TRUE} the iterative method is selected, if \code{FALSE} the monolithic one.
#' @param threshold This parameter is used for arresting algorithm iterations. Algorithm stops when two successive iterations lead to improvement in penalized log-likelihood smaller than threshold.
#' Default value \code{threshold = 10^(-4)}.
#' @param max.steps This parameter is used to limit the maximum number of iteration.
#' Default value \code{max.steps=50}.
#' @param IC Initial condition needed in case of parabolic problem i.e. when \code{FLAG_PARABOLIC==TRUE}.
#' If \code{FLAG_PARABOLIC==FALSE} this parameter is ignored. If \code{FLAG_PARABOLIC=TRUE} and \code{IC=NULL} it is necessary to provide
#' also data at the initial time. IC will be estimated from them.
#' @param search a flag to decide the search algorithm type (tree or naive or walking search algorithm).
#' @param bary.locations A list with three vectors:
#' \code{locations}, location points which are same as the given locations options. (checks whether both locations are the same);
#' \code{element ids}, a vector of element id of the points from the mesh where they are located;
#' \code{barycenters}, a vector of barycenter of points from the located element.
#' @param family This parameter specify the distribution within exponential family used for GLM model.
#' The following distribution are implemented: "binomial", "exponential", "gamma", "poisson", "gaussian", "invgaussian".
#' The default link function for binomial is \code{logit} if you want either \code{probit} or \code{clogloc} set \code{family = "probit"}, \code{family = "cloglog"}.
#' @param mu0 This parameter is a vector that set the starting point for FPIRLS algorithm. It represent an initial guess of the location parameter.
#' Default is set to observation for non binary distribution while equal to \code{0.5(observations + 0.5)} for binary data.
#' @param scale.param Dispersion parameter of the chosen distribution. This is only required for "gamma", "gaussian", "invgaussian".
#' User may specify the parameter as a positive real number. If the parameter is not supplied, it is estimated from data according to Wilhelm Sangalli 2016.
#' @param threshold.FPIRLS This parameter is used for arresting algorithm iterations. Algorithm stops when two successive iterations lead to improvement in penalized log-likelihood smaller than threshold.FPIRLS.
#' Default value \code{threshold.FPIRLS = 0.0002020}.
#' @param max.steps.FPIRLS This parameter is used to limit the maximum number of iteration.
#' Default value \code{max.steps.FPIRLS=15}.
#' @param lambda.selection.criterion This parameter is used to select the optimization method related to smoothing parameter \code{lambda}.
#' The following methods are implemented: 'grid', further optimization methods are yet to come.
#' The 'grid' is a pure evaluation method, therefore a vector of \code{lambda} testing penalizations must be provided.
#' Default value \code{lambda.selection.criterion='grid'}
#' @param DOF.evaluation This parameter is used to identify if and how degrees of freedom computation has to be performed.
#' The following possibilities are allowed: NULL, 'exact' and 'stochastic'
#' In the former case no degree of freedom is computed, while the other two methods enable computation.
#' Stochastic computation of DOFs may be slightly less accurate than its deterministic counterpart, but is highly suggested for meshes of more than 5000 nodes, being fairly less time consuming.
#' Default value \code{DOF.evaluation=NULL}
#' @param lambda.selection.lossfunction This parameter is used to understand if some loss function has to be evaluated.
#' The following possibilities are allowed: NULL and 'GCV' (generalized cross validation)
#' The former case is that of \code{lambda.selection.criterion='grid'} pure evaluation, while the second can be employed for optimization methods.
#' Default value \code{lambda.selection.lossfunction=NULL}
#' @param lambdaS A scalar or vector of spatial smoothing parameters.
#' @param lambdaT A scalar or vector of temporal smoothing parameters.
#' @param DOF.stochastic.realizations This parameter is considered only when \code{DOF.evaluation = 'stochastic'}.
#' It is a positive integer that represents the number of uniform random variables used in stochastic GCV computation.
#' Default value \code{DOF.stochastic.realizations=100}.
#' @param DOF.stochastic.seed This parameter is considered only when \code{DOF.evaluation = 'stochastic'}.
#' It is a positive integer that represents user defined seed employed in stochastic GCV computation.
#' Default value \code{DOF.stochastic.seed=0}.
#' @param DOF.matrix Matrix of degrees of freedom. This parameter can be used if the DOF.matrix corresponding to \code{lambdaS} and \code{lambdaT} is available from precedent computation. This allows to save time
#' since the computation of the DOFs is the most expensive part of GCV.
#' @param GCV.inflation.factor Tuning parameter used for the estimation of GCV. Default value \code{GCV.inflation.factor = 1.0}.
#' It is advised to set it grather than 1 to avoid overfitting.
#' @param lambda.optimization.tolerance Tolerance parameter, a double between 0 and 1 that fixes how much precision is required by the optimization method: the smaller the parameter, the higher the accuracy.
#' Used only if \code{lambda.selection.criterion='newton'} or \code{lambda.selection.criterion='newton_fd'}, thus ot implemented yet.
#' Default value \code{lambda.optimization.tolerance=0.05}.
#' @param inference.data.object.time An \code{\link{inferenceDataObjectTime}} that stores all the information regarding inference over the linear and nonlinear parameters of the model. This parameter needs to be
#' consistent with \code{covariates} and mesh dimension number, otherwise will be discarded. If set and well defined, the function will have in output the inference results. It is suggested to create this object via \code{\link{inferenceDataObjectTimeBuilder}} function, so that the object is guaranteed to be well defined.
#' @return A list with the following variables:
#' \describe{
#' \item{\code{fit.FEM.time}}{A \code{FEM.time} object that represents the fitted spatio-temporal field.}
#' \item{\code{PDEmisfit.FEM.time}}{A \code{FEM.time} object that represents the misfit of the penalized PDE.}
#' \item{\code{beta}}{If \code{covariates} is not \code{NULL}, a matrix with number of rows equal to the number of covariates and number of columns equal to length of lambda. The \code{j}th column represents the vector of regression coefficients when
#' the smoothing parameter is equal to \code{lambda[j]}.}
#' \item{\code{edf}}{If GCV is \code{TRUE}, a scalar or matrix with the trace of the smoothing matrix for each combination of the smoothing parameters specified in \code{lambdaS} and \code{lambdaT}.}
#' \item{\code{stderr}}{If GCV is \code{TRUE}, a scalar or matrix with the estimate of the standard deviation of the error for each combination of the smoothing parameters specified in \code{lambdaS} and \code{lambdaT}.}
#' \item{\code{GCV}}{If GCV is \code{TRUE}, a scalar or matrix with the value of the GCV criterion for each combination of the smoothing parameters specified in \code{lambdaS} and \code{lambdaT}.}
#' \item{\code{bestlambda}}{If GCV is \code{TRUE}, a 2-elements vector with the indices of smoothing parameters returning the lowest GCV}
#' \item{\code{ICestimated}}{If FLAG_PARABOLIC is \code{TRUE} and IC is \code{NULL}, a list containing a \code{FEM} object with the initial conditions, the value of the smoothing parameter lambda returning the lowest GCV and, in presence of covariates, the estimated beta coefficients}
#' \item{\code{bary.locations}}{A barycenter information of the given locations if the locations are not mesh nodes.}
#' \item{\code{inference}}{A list set only if a well defined \code{\link{inferenceDataObjectTime}} is passed as parameter to the function; contains all inference outputs required:
#' \describe{
#' \item{\code{p_values}}{list of lists set only if at least one p-value is required; contains the p-values divided by implementation:
#' \describe{
#' \item{\code{wald}}{list containing all the Wald p-values required, in the same order of the \code{type} list in \code{inference.data.object.time}.
#' If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of \code{coeff} matrix in \code{inference.data.object.time}.
#' }
#' \item{\code{speckman}}{list containing all the Speckman p-values required, in the same order of the \code{type} list in \code{inference.data.object.time}.
#' If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of \code{coeff} matrix in \code{inference.data.object.time}.
#' }
#' \item{\code{eigen_sign_flip}}{list containing all the Eigen-Sign-Flip p-values required, in the same order of the \code{type} list in \code{inference.data.object.time}.
#' If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of \code{coeff} matrix in \code{inference.data.object.time}.
#' }
#' }
#' }
#' \item{\code{CI}}{list of lists set only if at least one confidence interval is required; contains the confidence intervals divided by implementation:
#' \describe{
#' \item{\code{wald}}{list containing all the Wald confidence intervals required, in the same order of the \code{type} list in \code{inference.data.object.time}.
#' Each item is a matrix with 3 columns and p rows, p being the number of rows of \code{coeff} matrix in \code{inference.data.object.time}; each row is the CI for the corresponding row of \code{coeff} matrix.
#' }
#' \item{\code{speckman}}{list containing all the Speckman confidence intervals required, in the same order of the \code{type} list in \code{inference.data.object.time}.
#' Each item is a matrix with 3 columns and p rows, p being the number of rows of \code{coeff} matrix in \code{inference.data.object.time}; each row is the CI for the corresponding row of \code{coeff} matrix.
#' }
#' }
#' }
#' }
#' }
#' }
#' @description Space-time regression with differential regularization. Space-varying covariates can be included in the model. The technique accurately handle data distributed over irregularly shaped domains. Moreover, various conditions can be imposed at the domain boundaries.
#' @usage smooth.FEM.time(locations = NULL, time_locations = NULL, observations, FEMbasis,
#' time_mesh=NULL, covariates = NULL, PDE_parameters = NULL, BC = NULL,
#' incidence_matrix = NULL, areal.data.avg = TRUE,
#' FLAG_MASS = FALSE, FLAG_PARABOLIC = FALSE, FLAG_ITERATIVE = FALSE,
#' threshold = 10^(-4), max.steps = 50, IC = NULL,
#' search = "tree", bary.locations = NULL,
#' family = "gaussian", mu0 = NULL, scale.param = NULL,
#' threshold.FPIRLS = 0.0002020, max.steps.FPIRLS = 15,
#' lambda.selection.criterion = "grid", DOF.evaluation = NULL,
#' lambda.selection.lossfunction = NULL, lambdaS = NULL, lambdaT = NULL,
#' DOF.stochastic.realizations = 100, DOF.stochastic.seed = 0,
#' DOF.matrix = NULL, GCV.inflation.factor = 1, lambda.optimization.tolerance = 0.05,
#' inference.data.object.time=NULL)
#' @export
#' @references #' @references Arnone, E., Azzimonti, L., Nobile, F., & Sangalli, L. M. (2019). Modeling
#' spatially dependent functional data via regression with differential regularization.
#' Journal of Multivariate Analysis, 170, 275-295.
#' Bernardi, M. S., Sangalli, L. M., Mazza, G., & Ramsay, J. O. (2017). A penalized
#' regression model for spatial functional data with application to the analysis of the
#' production of waste in Venice province.
#' Stochastic Environmental Research and Risk Assessment, 31(1), 23-38.
#' Federico Ferraccioli, Laura M. Sangalli & Livio Finos (2022). Some first inferential tools for spatial regression
#' with differential regularization. Journal of Multivariate Analysis, 189, 104866.
#' @examples
#' library(fdaPDE)
#'
#' data(horseshoe2D)
#' boundary_nodes = horseshoe2D$boundary_nodes
#' boundary_segments = horseshoe2D$boundary_segments
#' locations = horseshoe2D$locations
#' time_locations = seq(0,1,length.out = 5)
#'
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#'
#' space_time_locations = cbind(rep(time_locations,each=nrow(mesh$nodes)),
#' rep(mesh$nodes[,1],5),rep(mesh$nodes[,2],5))
#'
#' FEMbasis = create.FEM.basis(mesh)
#' lambdaS = 10^-1
#' lambdaT = 10^-1
#' data = fs.test(space_time_locations[,2],
#' space_time_locations[,3])*cos(pi*space_time_locations[,1]) +
#' rnorm(nrow(space_time_locations), sd = 0.5)
#' data = matrix(data, nrow = nrow(mesh$nodes), ncol = length(time_locations), byrow = TRUE)
#'
#' solution = smooth.FEM.time(observations = data, time_locations = time_locations,
#' FEMbasis = FEMbasis, lambdaS = lambdaS, lambdaT = lambdaT)
#' plot(solution$fit.FEM)
smooth.FEM.time<-function(locations = NULL, time_locations = NULL, observations, FEMbasis, time_mesh=NULL,
covariates = NULL, PDE_parameters = NULL, BC = NULL,
incidence_matrix = NULL, areal.data.avg = TRUE,
FLAG_MASS = FALSE, FLAG_PARABOLIC = FALSE,FLAG_ITERATIVE = FALSE, threshold = 10^(-4), max.steps = 50, IC = NULL,
search = "tree", bary.locations = NULL,
family = "gaussian", mu0 = NULL, scale.param = NULL,
threshold.FPIRLS = 0.0002020, max.steps.FPIRLS = 15,
lambda.selection.criterion = "grid", DOF.evaluation = NULL, lambda.selection.lossfunction = NULL,
lambdaS = NULL, lambdaT = NULL, DOF.stochastic.realizations = 100, DOF.stochastic.seed = 0, DOF.matrix = NULL, GCV.inflation.factor = 1, lambda.optimization.tolerance = 0.05,
inference.data.object.time = NULL)
{
if(is(FEMbasis$mesh, "mesh.2D"))
{
ndim = 2
mydim = 2
}else if(is(FEMbasis$mesh, "mesh.2.5D"))
{
ndim = 3
mydim = 2
}else if(is(FEMbasis$mesh, "mesh.3D"))
{
ndim = 3
mydim = 3
}else if(is(FEMbasis$mesh, "mesh.1.5D"))
{
ndim = 2
mydim = 1
}else
{
stop('Unknown mesh class')
}
##################### Checking parameters, sizes and conversion ################################
# Preliminary consistency of optimization parameters
if(lambda.selection.criterion == "grid")
{
optim = 0
}else if(lambda.selection.criterion == "newton")
{
optim = 1
}else if(lambda.selection.criterion == "newton_fd")
{
optim = 2
}else
{
stop("'lambda.selection.criterion' must belong to the following list: 'none', 'grid', 'newton', 'newton_fd'.")
}
if(is.null(DOF.evaluation))
{
optim = c(optim,0)
}else if(DOF.evaluation == 'stochastic')
{
optim = c(optim,1)
}else if(DOF.evaluation == 'exact')
{
optim = c(optim,2)
}else
{
stop("'DOF.evaluation' must be NULL, 'stochastic' or 'exact'.")
}
if(is.null(lambda.selection.lossfunction))
{
optim = c(optim,0)
}else if(lambda.selection.lossfunction == 'GCV')
{
optim = c(optim,1)
}else
{
stop("'lambda.selection.lossfunction' has to be 'GCV'.")
}
if(any(lambdaS<=0) || any(lambdaT<=0))
stop("'lambda' can not be less than or equal to 0")
if(optim[2]!=0 & optim[3]!=1)
{
warning("Dof are computed, setting 'lambda.selection.lossfunction' to 'GCV'")
optim[3]=1
}
if(optim[1]==1 & optim[2]!=2)
{
warning("This method needs evaluate DOF in an 'exact' way, selecting 'DOF.evaluation'='exact'")
optim[2] = 2
}
if(!is.null(BC) & optim[1]==1)
{
warning("'newton' 'lambda.selection.criterion' can't be performed with non-NULL boundary conditions, using 'newton_fd' instead")
optim[1] = 2
}
if((optim[1]==2 & optim[2]==0) || (optim[1]==0 & optim[2]==0 & optim[3]==1 & is.null(DOF.matrix)))
{
warning("This method needs evaluate DOF, selecting 'DOF.evaluation'='stochastic'")
optim[2] = 1
}
if(optim[1]!=0 & optim[3]==0)
{
warning("An optimized method needs a loss function to perform the evaluation, selecting 'lambda.selection.lossfunction' as 'GCV'")
optim[3] = 1
}
if(optim[1]>0 & (is.null(lambdaS) || is.null(lambdaT)))
{
warning("No initial point is given: automatic initialization of Newton method")
}
# Search algorithm
if(search=="naive"){
search=1
}else if(search=="tree"){
search=2
}else if(search=="walking" & is(FEMbasis$mesh, "mesh.2.5D")){
stop("walking search is not available for mesh class mesh.2.5D.")
}else if(search=="walking" & is(FEMbasis$mesh, "mesh.1.5D")){
stop("walking search is not available for mesh class mesh.1.5D.")
}else if(search=="walking" & !is(FEMbasis$mesh, "mesh.2.5D") & !is(FEMbasis$mesh, "mesh.1.5D")){
search=3
}else{
stop("'search' must must belong to the following list: 'naive', 'tree' or 'walking'.")
}
# If locations is null but bary.locations is not null, use the locations in bary.locations
if(is.null(locations) & !is.null(bary.locations))
{
locations = bary.locations$locations
locations = as.matrix(locations)
}
## Converting to format for internal usage
if(!is.null(locations))
locations = as.matrix(locations)
if(!is.null(time_locations))
time_locations = as.matrix(time_locations)
if(!is.null(time_mesh))
time_mesh = as.matrix(time_mesh)
observations = as.matrix(observations)
if(!is.null(covariates))
covariates = as.matrix(covariates)
if(!is.null(DOF.matrix))
DOF.matrix = as.matrix(DOF.matrix)
if(!is.null(incidence_matrix))
incidence_matrix = as.matrix(incidence_matrix)
if(!is.null(IC))
IC = as.matrix(IC)
if(!is.null(BC))
{
BC$BC_indices = as.matrix(BC$BC_indices)
BC$BC_values = as.matrix(BC$BC_values)
}
if(!is.null(lambdaS))
lambdaS = as.matrix(lambdaS)
if(!is.null(lambdaT))
lambdaT = as.matrix(lambdaT)
space_varying = checkSmoothingParameters_time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh = time_mesh,
covariates = covariates, PDE_parameters = PDE_parameters, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold, max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance)
# only if inference is required
if(!is.null(inference.data.object.time)){
# Check that GCV is set for inference
if(inference.data.object.time@definition==1 && is.null(lambda.selection.lossfunction) && (!is.numeric(lambdaS) && !is.null(lambdaT)) && (nrow(lambdaS)!=1 || ncol(lambdaS)!=1 || nrow(lambdaT)!=1 || ncol(lambdaT)!=1)){
warning("Inference is not defined when lambda grid is provided without GCV")
inference.data.object.time=new("inferenceDataObjectTime", test = as.integer(0), interval =as.integer(0), type = as.integer(0), component = as.integer(0), exact = as.integer(0), dim = as.integer(0), n_cov = as.integer(0),
locations = matrix(data=0, nrow = 1 ,ncol = 1), locations_indices = as.integer(0), locations_are_nodes = as.integer(0), time_locations = 0, coeff = matrix(data=0, nrow = 1 ,ncol = 1), beta0 = -1, f0 = function(){},
f0_eval = -1, f_var = as.integer(0), quantile = -1, alpha = 0, n_flip = as.integer(1000), tol_fspai = -1, definition=as.integer(0))
}
if(inference.data.object.time@definition==1 && FLAG_ITERATIVE == T){
warning("Inference is not provided when iterative method is selected")
inference.data.object.time=new("inferenceDataObjectTime", test = as.integer(0), interval =as.integer(0), type = as.integer(0), component = as.integer(0), exact = as.integer(0), dim = as.integer(0), n_cov = as.integer(0),
locations = matrix(data=0, nrow = 1 ,ncol = 1), locations_indices = as.integer(0), locations_are_nodes = as.integer(0), time_locations = 0, coeff = matrix(data=0, nrow = 1 ,ncol = 1), beta0 = -1, f0 = function(){},
f0_eval = -1, f_var = as.integer(0), quantile = -1, alpha = 0, n_flip = as.integer(1000), tol_fspai = -1, definition=as.integer(0))
}
}
# If I have PDE non-sv case I need (constant) matrices as parameters
if(!is.null(PDE_parameters) & space_varying==FALSE)
{
PDE_parameters$K = as.matrix(PDE_parameters$K)
PDE_parameters$b = as.matrix(PDE_parameters$b)
PDE_parameters$c = as.matrix(PDE_parameters$c)
}
checkSmoothingParametersSize_time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh = time_mesh,
covariates = covariates, PDE_parameters = PDE_parameters, incidence_matrix = incidence_matrix,
BC = BC, space_varying = space_varying, ndim = ndim, mydim = mydim,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC, IC = IC,
lambdaS = lambdaS, lambdaT = lambdaT, DOF.matrix = DOF.matrix)
# Further check
dim_obs_inf <- dim(observations)[1]
observations<-as.vector(observations)
if(is.null(time_locations))
{
if(FLAG_PARABOLIC && !is.null(IC))
time_locations <- time_mesh[2:length(time_mesh)]
else
time_locations <- time_mesh
}
if(is.null(time_mesh))
{
if(FLAG_PARABOLIC && !is.null(IC))
time_mesh <- rbind(2*time_locations[1]-time_locations[2],time_locations)
else
time_mesh<-time_locations
}
# Checking inference data: after time_locations definition
# Most of the checks have already been carried out by inferenceDataObjectBuilderTime function
if(!is.null(locations))
inference.data.object.time <- checkInferenceParametersTime(inference.data.object.time,ncol(covariates),time_locations,locations,FEMbasis$mesh$nodes, FLAG_PARABOLIC, is.null(IC)) #checking inference data consistency, constructing default object in NULL case
else
inference.data.object.time <- checkInferenceParametersTime(inference.data.object.time,ncol(covariates),time_locations,FEMbasis$mesh$nodes[1:dim_obs_inf,],FEMbasis$mesh$nodes, FLAG_PARABOLIC, is.null(IC)) #checking inference data consistency, constructing default object in NULL case
# Check whether the locations coincide with the mesh nodes (should be put after all the validations)
if (!is.null(locations))
{
if(dim(locations)[1]==dim(FEMbasis$mesh$nodes)[1] & dim(locations)[2]==dim(FEMbasis$mesh$nodes)[2])
{
sum1=0
sum2=0
for (i in 1:nrow(locations))
{
sum1 = sum1 + abs(locations[i,1]-FEMbasis$mesh$nodes[i,1])
sum2 = sum2 + abs(locations[i,2]-FEMbasis$mesh$nodes[i,2])
}
if (sum1==0 & sum2==0)
{
message("No search algorithm is used because the locations coincide with the nodes.")
locations = NULL #In principle, R uses pass-by-value semantics in its function calls. So put ouside of checkSmoothingParameters function.
}
}
}
# OPTIMIZATION NOT IMPLEMENTED FOR GAM
if (family != 'gaussian' & optim[1] != 0)
stop("'lambda.selection.criterion' = 'grid' is the only method implemented for GAM problems")
# FAMILY CHECK
family_admit = c("binomial", "exponential", "gamma", "poisson", "gaussian", "Gaussian")
if (sum(family == family_admit) == 0)
stop("'family' parameter required.\nCheck if it is one of the following: binomial, exponential, gamma, poisson, gaussian")
################## End checking parameters, sizes and conversion #############################
if (family == "gaussian")
{
if(is(FEMbasis$mesh, "mesh.2D") & is.null(PDE_parameters))
{
bigsol = NULL
bigsol = CPP_smooth.FEM.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold , max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.2D") & !is.null(PDE_parameters) & space_varying==FALSE)
{
bigsol = NULL
bigsol = CPP_smooth.FEM.PDE.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, PDE_parameters=PDE_parameters, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold, max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.2D") & !is.null(PDE_parameters) & space_varying==TRUE)
{
bigsol = NULL
bigsol = CPP_smooth.FEM.PDE.sv.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, PDE_parameters=PDE_parameters, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold, max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.2.5D"))
{
bigsol = NULL
bigsol = CPP_smooth.manifold.FEM.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold , max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.3D") & is.null(PDE_parameters))
{
bigsol = NULL
bigsol = CPP_smooth.volume.FEM.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold , max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.3D") & !is.null(PDE_parameters) & space_varying==FALSE)
{
bigsol = NULL
bigsol = CPP_smooth.volume.FEM.PDE.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, PDE_parameters=PDE_parameters, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold, max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.3D") & !is.null(PDE_parameters) & space_varying==TRUE)
{
bigsol = NULL
bigsol = CPP_smooth.volume.FEM.PDE.sv.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, PDE_parameters=PDE_parameters, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold, max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.1.5D"))
{
bigsol = NULL
bigsol = CPP_smooth.graph.FEM.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold , max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}
# ---------- Solution -----------
N = nrow(FEMbasis$mesh$nodes)
M = ifelse(FLAG_PARABOLIC,length(time_mesh)-1,length(time_mesh) + 2);
dim_1 = ifelse(optim[1]==0 & is.null(DOF.matrix) & optim[3]==0, length(lambdaS), 1)
dim_2 = ifelse(optim[1]==0 & is.null(DOF.matrix) & optim[3]==0, length(lambdaT), 1)
if(is.null(IC) && FLAG_PARABOLIC)
IC = bigsol[[27]]$coeff
if(FLAG_PARABOLIC)
{
f = array(dim=c(length(IC)+M*N,dim_1,dim_2))
g = array(dim=c(length(IC)+M*N,dim_1,dim_2))
for (i in 1:dim_1)
for (j in 1:dim_2)
{
f[,i,j] = c(IC,bigsol[[1]][1:(N*M),i+(j-1)*dim_1])
g[,i,j] = c(rep(0,length(IC)),bigsol[[1]][(N*M+1):(2*N*M),i+(j-1)*dim_1])
}
}
else
{
f = array(data=bigsol[[1]][1:(N*M),],dim = c(N*M,dim_1,dim_2))
g = array(data=bigsol[[1]][(N*M+1):(2*N*M),],dim = c(N*M,dim_1,dim_2))
}
GCV_ = bigsol[[3]]
if(!is.null(covariates))
{
if(optim[1]==0 & is.null(DOF.matrix) & optim[3]==0)
beta = array(data=bigsol[[5]],dim=c(ncol(covariates),length(lambdaS),length(lambdaT)))
else
beta = array(data=bigsol[[5]],dim=c(ncol(covariates),1,1))
}
else
beta = NULL
if(all(is.na(bigsol[[27]])))
ICestimated = NULL
else
ICestimated = list(IC.FEM=bigsol[[27]],bestlambdaindex=bigsol[[28]],bestlambda=bigsol[[29]],beta=bigsol[[30]])
bestlambda = bigsol[[4]]+1
if(optim[1]!=0) # newton or newton_fd
{
bestlambda = max(bestlambda[1], bestlambda[2])
}
if(optim[1]==0)
{
if(bestlambda[1] == 1 || bestlambda[1] == length(lambdaS))
warning("Your optimal 'GCV' is on the border of lambdaS sequence")
if(bestlambda[2] == 1 || bestlambda[2] == length(lambdaT))
warning("Your optimal 'GCV' is on the border of lambdaT sequence")
}
if (is.null(lambda.selection.lossfunction))
sd = -1
else
sd = sqrt(bigsol[[15]])
solution = list(
f = f,
g = g,
z_hat = bigsol[[13]],
rmse = bigsol[[14]],
estimated_sd=sd
)
term = bigsol[[19]]
ot = bigsol[[20]]
if(term == 1)
termination = "reached tolerance"
else if(term == 2)
termination = "reached max number iterations"
else
termination = "uninformative"
if(ot == 0)
optimization_type = "full optimization"
else if(ot == 1)
optimization_type = "full DOF grid"
else
optimization_type = "uninformative"
optimization = list(
lambda_solution = bigsol[[16]],
lambda_position = bestlambda,
optimization_details = list(
iterations = bigsol[[18]],
termination = termination,
optimization_type = optimization_type),
dof = bigsol[[2]],
lambda_vector = matrix(c(bigsol[[21]], bigsol[[22]]), nrow=length(bigsol[[21]]), ncol=2),
GCV_vector = GCV_
)
time = bigsol[[23]]
# Save information of Tree Mesh
tree_mesh = list(
treelev = bigsol[[6]][1],
header_orig= bigsol[[7]],
header_scale = bigsol[[8]],
node_id = bigsol[[9]][,1],
node_left_child = bigsol[[9]][,2],
node_right_child = bigsol[[9]][,3],
node_box= bigsol[[10]])
# Reconstruct FEMbasis with tree mesh
mesh.class= class(FEMbasis$mesh)
if (is.null(FEMbasis$mesh$treelev)) { #if doesn't exist the tree information
FEMbasis$mesh = append(FEMbasis$mesh, tree_mesh)
} #if already exist the tree information, don't append
class(FEMbasis$mesh) = mesh.class
# Save information of Barycenter
if (is.null(bary.locations)) {
bary.locations = list(locations=locations, element_ids = bigsol[[11]], barycenters = bigsol[[12]])
}
class(bary.locations) = "bary.locations"
# Make FEM.time objects
fit.FEM.time = FEM.time(f, time_mesh, FEMbasis, FLAG_PARABOLIC)
PDEmisfit.FEM.time = FEM.time(g, time_mesh, FEMbasis, FLAG_PARABOLIC)
# Save statistics and intervals
if(inference.data.object.time@definition==1){
inference = {}
n_loc_inference = ifelse(FLAG_PARABOLIC==TRUE,dim(inference.data.object.time@locations)[1]*(length(inference.data.object.time@time_locations)-1),dim(inference.data.object.time@locations)[1]*length(inference.data.object.time@time_locations))
confidence_intervals = matrix(data = bigsol[[25]], nrow = 2*3*length(inference.data.object.time@type), ncol = max(dim(inference.data.object.time@coeff)[1], n_loc_inference))
p_val = matrix(data = bigsol[[24]], nrow = dim(inference.data.object.time@coeff)[1]+1, ncol = length(inference.data.object.time@type))
for(i in 1:length(inference.data.object.time@type)){ # each element is a different inferential setting
if(inference.data.object.time@interval[i]!=0){ # Intervals requested by this setting, adding them to the right implementation position
if(inference.data.object.time@component[i]!=2){ # intervals for beta were requested, just for safety, actually it cannot be equal to 2
p = dim(inference.data.object.time@coeff)[1]
ci_beta=t(confidence_intervals[(3*(2*i-2)+1):(3*(2*i-2)+3),1:p])
if(inference.data.object.time@type[i]==1){
inference$beta$CI$wald[[length(inference$beta$CI$wald)+1]] = ci_beta
inference$beta$CI$wald=as.list(inference$beta$CI$wald)
}
else if(inference.data.object.time@type[i]==2){
inference$beta$CI$speckman[[length(inference$beta$CI$speckman)+1]] = ci_beta
inference$beta$CI$speckman=as.list(inference$beta$CI$speckman)
}
else if(inference.data.object.time@type[i]==3){
if(ci_beta[2]> 10^20){
warning("ESF CI bisection algorithm did not converge, returning NA")
for(h in 1:nrow(ci_beta))
for(k in 1:ncol(ci_beta))
ci_beta[h,k]=NA
}
inference$beta$CI$eigen_sign_flip[[length(inference$beta$CI$eigen_sign_flip)+1]] = ci_beta
inference$beta$CI$eigen_sign_flip=as.list(inference$beta$CI$eigen_sign_flip)
}
else if(inference.data.object.time@type[i]==4){
if(ci_beta[2]> 10^20){
warning("Enhanced ESF CI bisection algorithm did not converge, returning NA")
for(h in 1:nrow(ci_beta))
for(k in 1:ncol(ci_beta))
ci_beta[h,k]=NA
}
inference$beta$CI$enh_eigen_sign_flip[[length(inference$beta$CI$enh_eigen_sign_flip)+1]] = ci_beta
inference$beta$CI$enh_eigen_sign_flip=as.list(inference$beta$CI$enh_eigen_sign_flip)
}
}
if(inference.data.object.time@component[i]!=1){ # intervals for f were requested
n_loc = dim(inference.data.object.time@locations)[1]*length(inference.data.object.time@time_locations)
ci_f=t(confidence_intervals[(3*(2*i-1)+1):(3*(2*i-1)+3),])
if(inference.data.object.time@type[i]==1){ # wald confidence intervals for f
inference$f$CI$wald[[length(inference$f$CI$wald)+1]] = ci_f
inference$f$CI$wald=as.list(inference$f$CI$wald)
}
}
}
if(inference.data.object.time@test[i]!=0){ # Test requested by this setting, adding them to the right implementation position
statistics=p_val[,i]
if(inference.data.object.time@component[i]!=2){ # test on beta was requested, just for safety, it cannot be otherwise
beta_statistics = statistics[1:dim(inference.data.object.time@coeff)[1]]
p_values = numeric()
if(inference.data.object.time@type[i]==3 || inference.data.object.time@type[i]==4){ # eigen-sign-flip p-value is already computed in cpp code
if(inference.data.object.time@test[i]==1){
# one-at-the-time tests
p_values = beta_statistics
}
else if(inference.data.object.time@test[i]==2){
# simultaneous test
p_values = beta_statistics[1]
}
}else{
# Compute p-values
if(inference.data.object.time@test[i]==1){ # Wald and Speckman return statistics and needs computation of p-values (no internal use of distributions quantiles)
# one-at-the-time-tests
p_values = numeric(length(beta_statistics))
for(l in 1:length(beta_statistics)){
p_values[l] = 2*pnorm(-abs(beta_statistics[l]))
}
}
else if(inference.data.object.time@test[i]==2){
# simultaneous tests
p = dim(inference.data.object.time@coeff)[1]
p_values = 1-pchisq(beta_statistics[1], p)
}
}
# add p-values in the right position
if(inference.data.object.time@type[i]==1){
inference$beta$p_values$wald[[length(inference$beta$p_values$wald)+1]] = p_values
inference$beta$p_values$wald=as.list(inference$beta$p_values$wald)
}
else if(inference.data.object.time@type[i]==2){
inference$beta$p_values$speckman[[length(inference$beta$p_values$speckman)+1]] = p_values
inference$beta$p_values$speckman=as.list(inference$beta$p_values$speckman)
}
else if(inference.data.object.time@type[i]==3){
inference$beta$p_values$eigen_sign_flip[[length(inference$beta$p_values$eigen_sign_flip)+1]] = p_values
inference$beta$p_values$eigen_sign_flip=as.list(inference$beta$p_values$eigen_sign_flip)
}
else if(inference.data.object.time@type[i]==4){
inference$beta$p_values$enh_eigen_sign_flip[[length(inference$beta$p_values$enh_eigen_sign_flip)+1]] = p_values
inference$beta$p_values$enh_eigen_sign_flip=as.list(inference$beta$p_values$enh_eigen_sign_flip)
}
}
if(inference.data.object.time@component[i]!=1){ # test on f was requested
p_value = statistics[length(statistics)]
# add p-value in the right position
if(inference.data.object.time@type[i]==1){
inference$f$p_values$wald[[length(inference$f$p_values$wald)+1]] = p_value
inference$f$p_values$wald=as.list(inference$f$p_values$wald)
}
}
}
}
if(inference.data.object.time@f_var==1){
f_variances = matrix(data = bigsol[[26]], nrow = length(solution$z_hat), ncol = 1)
inference$f_var = f_variances
}
# Prepare return list
reslist = list(fit.FEM.time = fit.FEM.time, PDEmisfit.FEM.time = PDEmisfit.FEM.time, solution = solution,
optimization = optimization, beta = beta, time = time, ICestimated=ICestimated, bary.locations = bary.locations, inference=inference)
}else{
# Prepare return list
reslist = list(fit.FEM.time = fit.FEM.time, PDEmisfit.FEM.time = PDEmisfit.FEM.time, solution = solution,
optimization = optimization, beta = beta, time = time, ICestimated=ICestimated, bary.locations = bary.locations)
}
return(reslist)
}
else
{
#Iterative solver not implemented for GAM
if(FLAG_ITERATIVE)
warning("Iterative solver not implemented for GAM, switching to block solver")
FLAG_ITERATIVE = FALSE
threshold = 10^(-4)
max.steps = 50
#----------------------------------------------------#
############# GAMs: FPIRLS algorithm #################
#----------------------------------------------------#
checkGAMParameters(observations = observations, max.steps.FPIRLS = max.steps.FPIRLS,
mu0 = mu0, scale.param = scale.param,
threshold.FPIRLS = threshold.FPIRLS, family = family)
if (is(FEMbasis$mesh, "mesh.2D") & is.null(PDE_parameters)) {
bigsol = NULL
bigsol = CPP_smooth.GAM.FEM.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates, ndim = ndim,
mydim = mydim, BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if (is(FEMbasis$mesh, "mesh.2D") & !is.null(PDE_parameters) & space_varying == FALSE) {
bigsol = NULL
bigsol = CPP_smooth.GAM.FEM.PDE.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates,
PDE_parameters = PDE_parameters, ndim = ndim, mydim = mydim,
BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE = FLAG_ITERATIVE,
threshold = threshold, max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if (is(FEMbasis$mesh, "mesh.2D") & !is.null(PDE_parameters) & space_varying == TRUE) {
bigsol = NULL
bigsol = CPP_smooth.GAM.FEM.PDE.sv.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates,
PDE_parameters = PDE_parameters, ndim = ndim, mydim = mydim,
BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE = FLAG_ITERATIVE,
threshold = threshold, max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
} else if( is(FEMbasis$mesh, "mesh.2.5D") & is.null(PDE_parameters)){
bigsol = NULL
bigsol = CPP_smooth.manifold.GAM.FEM.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates, ndim = ndim,
mydim = mydim, BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if(is(FEMbasis$mesh, "mesh.3D") & is.null(PDE_parameters)){
bigsol = NULL
bigsol = CPP_smooth.volume.GAM.FEM.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates, ndim = ndim,
mydim = mydim, BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if(is(FEMbasis$mesh, "mesh.3D") & !is.null(PDE_parameters) & space_varying == FALSE){
bigsol = NULL
bigsol = CPP_smooth.volume.GAM.FEM.PDE.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates,
PDE_parameters = PDE_parameters, ndim = ndim, mydim = mydim,
BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if(is(FEMbasis$mesh, "mesh.3D") & !is.null(PDE_parameters) & space_varying == TRUE){
bigsol = NULL
bigsol = CPP_smooth.volume.GAM.FEM.PDE.sv.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates,
PDE_parameters = PDE_parameters, ndim = ndim, mydim = mydim,
BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if(is(FEMbasis$mesh, "mesh.1.5D") & is.null(PDE_parameters)){
bigsol = NULL
bigsol = CPP_smooth.graph.GAM.FEM.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates, ndim = ndim,
mydim = mydim, BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else{
stop("Not implemented for !is.null(PDE_parameters). Try Laplacian regularization.")
}
ICindx = 16
N = nrow(FEMbasis$mesh$nodes)
M = ifelse(FLAG_PARABOLIC, length(time_mesh) - 1, length(time_mesh) + 2)
if (is.null(IC) && FLAG_PARABOLIC)
IC = bigsol[[ICindx]]$coeff[, bigsol[[ICindx + 1]]]
if (FLAG_PARABOLIC) {
f = array(dim = c(length(IC) + M * N, length(lambdaS), length(lambdaT)))
for (i in 1:length(lambdaS))
for (j in 1:length(lambdaT))
f[, i, j] = c(IC, bigsol[[1]][1:(N * M), i + (j - 1) * length(lambdaS)])
}else
f = array(data = bigsol[[1]][1:(N * M), ],
dim = c(N * M, length(lambdaS), length(lambdaT)))
if (FLAG_PARABOLIC) {
g = array(dim = c(length(IC) + M * N, length(lambdaS), length(lambdaT)))
for (i in 1:length(lambdaS))
for (j in 1:length(lambdaT))
g[, i, j] = c(rep(0, length(IC)), bigsol[[1]][(N * M + 1):(2 * N * M),
i + (j - 1) * length(lambdaS)])
}else
g = array(data = bigsol[[1]][(N * M + 1):(2 * N * M), ],
dim = c(N * M, length(lambdaS), length(lambdaT)))
dof = bigsol[[2]]
GCV_ = bigsol[[3]]
bestlambda = bigsol[[4]] + 1
if (!is.null(covariates))
beta = array(data = bigsol[[5]], dim = c(ncol(covariates), length(lambdaS), length(lambdaT)))
else
beta = NULL
if (all(is.na(bigsol[[ICindx]])))
ICestimated = NULL
else
ICestimated = list(IC.FEM = bigsol[[ICindx]], bestlambdaindex = bigsol[[ICindx + 1]],
bestlambda = bigsol[[ICindx + 2]], beta = bigsol[[ICindx + 3]])
# Save information of Tree Mesh
tree_mesh = list(treelev = bigsol[[6]][1], header_orig = bigsol[[7]],
header_scale = bigsol[[8]], node_id = bigsol[[9]][, 1],
node_left_child = bigsol[[9]][, 2],
node_right_child = bigsol[[9]][, 3], node_box = bigsol[[10]])
# Reconstruct FEMbasis with tree mesh
mesh.class = class(FEMbasis$mesh)
if (is.null(FEMbasis$mesh$treelev)) { # if doesn't exist the tree
# information
FEMbasis$mesh = append(FEMbasis$mesh, tree_mesh)
} # if already exist the tree information, don't append
class(FEMbasis$mesh) = mesh.class
# Save information of Barycenter
if (is.null(bary.locations)) {
bary.locations = list(locations = locations, element_ids = bigsol[[11]],
barycenters = bigsol[[12]])
}
class(bary.locations) = "bary.locations"
# Make FEM.time objects
fit.FEM.time = FEM.time(f, time_mesh, FEMbasis, FLAG_PARABOLIC)
PDEmisfit.FEM.time = FEM.time(g, time_mesh, FEMbasis, FLAG_PARABOLIC)
# Prepare return list
reslist = NULL
nm <- length(observations)
fn.eval = array(dim = c(nm, length(lambdaS), length(lambdaT)))
if(FLAG_PARABOLIC)
{
if(!is.null(bigsol[[ICindx]]))
ICfn.eval = bigsol[[ICindx+4]][, bigsol[[ICindx + 1]]]
else
ICfn.eval = NULL
nm <- nm - length(ICfn.eval)
for (i in 1:length(lambdaS))
for (j in 1:length(lambdaT))
fn.eval[, i, j] = c(ICfn.eval, bigsol[[13]][1:nm, i + (j - 1) * length(lambdaS)])
}else{
for (i in 1:length(lambdaS))
for (j in 1:length(lambdaT))
fn.eval[, i, j] = bigsol[[13]][1:nm, i + (j - 1) * length(lambdaS)]
}
J_minima = bigsol[[14]]
variance.est = bigsol[[15]]
if(!is.numeric(variance.est[1])){
variance.est <- NULL
} else if (variance.est[1] < 0) variance.est = NULL
if (!is.null(lambda.selection.lossfunction)) {
stderr = sqrt(GCV_ * (sum(!is.na(observations)) - dof) / sum(!is.na(observations)))
reslist = list(fit.FEM.time = fit.FEM.time,
PDEmisfit.FEM.time = PDEmisfit.FEM.time, beta = beta,
edf = dof, GCV = GCV_, stderr = stderr,
bestlambda = bestlambda, ICestimated = ICestimated,
bary.locations = bary.locations, fn.eval = fn.eval,
J_minima = J_minima, variance.est = variance.est)
}else{
reslist = list(fit.FEM.time = fit.FEM.time,
PDEmisfit.FEM.time = PDEmisfit.FEM.time, beta = beta,
ICestimated = ICestimated, bary.locations = bary.locations,
fn.eval = fn.eval, J_minima = J_minima,
variance.est = variance.est)
}
return(reslist)
}
}
Try the fdaPDE package in your browser
Any scripts or data that you put into this service are public.
fdaPDE documentation built on Sept. 17, 2024, 5:09 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions smooth.FEM.time
Documented in smooth.FEM.time
#' @useDynLib fdaPDE
#' @import Matrix plot3D rgl
#' @importFrom grDevices heat.colors palette
#' @importFrom graphics plot segments points lines
#' @importFrom methods is
NULL
#' Space-time regression with differential regularization
#' @param locations A matrix where each row specifies the spatial coordinates \code{x} and \code{y} (and \code{z} if ndim=3) of the corresponding observations in the vector \code{observations}.
#' This parameter can be \code{NULL}. In this case, if also the incidence matrix is \code{NULL} the spatial coordinates are assumed to coincide with the nodes of the \code{mesh}.
#' @param time_locations A vector containing the times of the corresponding observations in the vector \code{observations}.
#' This parameter can be \code{NULL}. In this case the temporal locations are assumed to coincide with the nodes of the \code{time_mesh}.
#' @param observations A matrix of #locations x #time_locations with the observed data values over the spatio-temporal domain.
#' The spatial locations of the observations can be specified with the \code{locations} argument.
#' @param FEMbasis A \code{FEMbasis} object describing the Finite Element basis, as created by \code{\link{create.FEM.basis}}.
#' @param time_mesh A vector specifying the time mesh.
#' @param covariates A #observations-by-#covariates matrix where each row represents the covariates associated with the corresponding observed data value in \code{observations}.
#' @param PDE_parameters A list specifying the parameters of the PDE in the regularizing term. Default is NULL, i.e. regularization is by means of the Laplacian (stationary, isotropic case).
#' If the PDE is elliptic it must contain: \code{K}, a 2-by-2 matrix of diffusion coefficients. This induces an anisotropic
#' smoothing with a preferential direction that corresponds to the first eigenvector of the diffusion matrix K; \code{b}, a vector of length 2 of advection coefficients. This induces a
#' smoothing only in the direction specified by the vector \code{b}; \code{c}, a scalar reaction coefficient. \code{c} induces a shrinkage of the surface to zero
#' If the PDE is space-varying it must contain: \code{K}, a function that for each spatial location in the spatial domain
#' (indicated by the vector of the 2 spatial coordinates) returns a 2-by-2 matrix of diffusion coefficients. This induces an anisotropic
#' smoothing with a local preferential direction that corresponds to the first eigenvector of the diffusion matrix K.The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be
#' an array with dimensions 2-by-2-by-#points.\code{b}, a function that for each spatial location in the spatial domain returns
#' a vector of length 2 of transport coefficients. This induces a local smoothing only in the direction specified by the vector \code{b}. The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be
#' a matrix with dimensions 2-by-#points; \code{c}, a function that for each spatial location in the spatial domain returns a scalar reaction coefficient.
#' \code{c} induces a shrinkage of the surface to zero. The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be
#' a vector with length #points; \code{u}, a function that for each spatial location in the spatial domain returns a scalar reaction coefficient.
#' \code{u} induces a reaction effect. The function must support recycling for efficiency reasons, thus if the input parameter is a #point-by-2 matrix, the output should be
#' a vector with length #points.
#' For 2.5D and 3D only the Laplacian is available (\code{PDE_parameters=NULL})
#' @param BC A list with two vectors:
#' \code{BC_indices}, a vector with the indices in \code{nodes} of boundary nodes where a Dirichlet Boundary Condition should be applied;
#' \code{BC_values}, a vector with the values that the spatial field must take at the nodes indicated in \code{BC_indices}.
#' @param incidence_matrix A #regions-by-#triangles/tetrahedrons matrix where the element (i,j) equals 1 if the j-th triangle/tetrahedron is in the i-th region and 0 otherwise.
#' This is only for areal data. In case of pointwise data, this parameter is set to \code{NULL}.
#' @param areal.data.avg Boolean. It involves the computation of Areal Data. If \code{TRUE} the areal data are averaged, otherwise not.
#' @param FLAG_MASS Boolean. This parameter is considered only for separable problems i.e. when \code{FLAG_PARABOLIC==FALSE}. If \code{TRUE} the mass matrix in space and in time are used, if \code{FALSE} they are substituted with proper identity matrices.
#' @param FLAG_PARABOLIC Boolean. If \code{TRUE} the parabolic problem problem is selected, if \code{FALSE} the separable one.
#' @param FLAG_ITERATIVE Boolean. If \code{TRUE} the iterative method is selected, if \code{FALSE} the monolithic one.
#' @param threshold This parameter is used for arresting algorithm iterations. Algorithm stops when two successive iterations lead to improvement in penalized log-likelihood smaller than threshold.
#' Default value \code{threshold = 10^(-4)}.
#' @param max.steps This parameter is used to limit the maximum number of iteration.
#' Default value \code{max.steps=50}.
#' @param IC Initial condition needed in case of parabolic problem i.e. when \code{FLAG_PARABOLIC==TRUE}.
#' If \code{FLAG_PARABOLIC==FALSE} this parameter is ignored. If \code{FLAG_PARABOLIC=TRUE} and \code{IC=NULL} it is necessary to provide
#' also data at the initial time. IC will be estimated from them.
#' @param search a flag to decide the search algorithm type (tree or naive or walking search algorithm).
#' @param bary.locations A list with three vectors:
#' \code{locations}, location points which are same as the given locations options. (checks whether both locations are the same);
#' \code{element ids}, a vector of element id of the points from the mesh where they are located;
#' \code{barycenters}, a vector of barycenter of points from the located element.
#' @param family This parameter specify the distribution within exponential family used for GLM model.
#' The following distribution are implemented: "binomial", "exponential", "gamma", "poisson", "gaussian", "invgaussian".
#' The default link function for binomial is \code{logit} if you want either \code{probit} or \code{clogloc} set \code{family = "probit"}, \code{family = "cloglog"}.
#' @param mu0 This parameter is a vector that set the starting point for FPIRLS algorithm. It represent an initial guess of the location parameter.
#' Default is set to observation for non binary distribution while equal to \code{0.5(observations + 0.5)} for binary data.
#' @param scale.param Dispersion parameter of the chosen distribution. This is only required for "gamma", "gaussian", "invgaussian".
#' User may specify the parameter as a positive real number. If the parameter is not supplied, it is estimated from data according to Wilhelm Sangalli 2016.
#' @param threshold.FPIRLS This parameter is used for arresting algorithm iterations. Algorithm stops when two successive iterations lead to improvement in penalized log-likelihood smaller than threshold.FPIRLS.
#' Default value \code{threshold.FPIRLS = 0.0002020}.
#' @param max.steps.FPIRLS This parameter is used to limit the maximum number of iteration.
#' Default value \code{max.steps.FPIRLS=15}.
#' @param lambda.selection.criterion This parameter is used to select the optimization method related to smoothing parameter \code{lambda}.
#' The following methods are implemented: 'grid', further optimization methods are yet to come.
#' The 'grid' is a pure evaluation method, therefore a vector of \code{lambda} testing penalizations must be provided.
#' Default value \code{lambda.selection.criterion='grid'}
#' @param DOF.evaluation This parameter is used to identify if and how degrees of freedom computation has to be performed.
#' The following possibilities are allowed: NULL, 'exact' and 'stochastic'
#' In the former case no degree of freedom is computed, while the other two methods enable computation.
#' Stochastic computation of DOFs may be slightly less accurate than its deterministic counterpart, but is highly suggested for meshes of more than 5000 nodes, being fairly less time consuming.
#' Default value \code{DOF.evaluation=NULL}
#' @param lambda.selection.lossfunction This parameter is used to understand if some loss function has to be evaluated.
#' The following possibilities are allowed: NULL and 'GCV' (generalized cross validation)
#' The former case is that of \code{lambda.selection.criterion='grid'} pure evaluation, while the second can be employed for optimization methods.
#' Default value \code{lambda.selection.lossfunction=NULL}
#' @param lambdaS A scalar or vector of spatial smoothing parameters.
#' @param lambdaT A scalar or vector of temporal smoothing parameters.
#' @param DOF.stochastic.realizations This parameter is considered only when \code{DOF.evaluation = 'stochastic'}.
#' It is a positive integer that represents the number of uniform random variables used in stochastic GCV computation.
#' Default value \code{DOF.stochastic.realizations=100}.
#' @param DOF.stochastic.seed This parameter is considered only when \code{DOF.evaluation = 'stochastic'}.
#' It is a positive integer that represents user defined seed employed in stochastic GCV computation.
#' Default value \code{DOF.stochastic.seed=0}.
#' @param DOF.matrix Matrix of degrees of freedom. This parameter can be used if the DOF.matrix corresponding to \code{lambdaS} and \code{lambdaT} is available from precedent computation. This allows to save time
#' since the computation of the DOFs is the most expensive part of GCV.
#' @param GCV.inflation.factor Tuning parameter used for the estimation of GCV. Default value \code{GCV.inflation.factor = 1.0}.
#' It is advised to set it grather than 1 to avoid overfitting.
#' @param lambda.optimization.tolerance Tolerance parameter, a double between 0 and 1 that fixes how much precision is required by the optimization method: the smaller the parameter, the higher the accuracy.
#' Used only if \code{lambda.selection.criterion='newton'} or \code{lambda.selection.criterion='newton_fd'}, thus ot implemented yet.
#' Default value \code{lambda.optimization.tolerance=0.05}.
#' @param inference.data.object.time An \code{\link{inferenceDataObjectTime}} that stores all the information regarding inference over the linear and nonlinear parameters of the model. This parameter needs to be
#' consistent with \code{covariates} and mesh dimension number, otherwise will be discarded. If set and well defined, the function will have in output the inference results. It is suggested to create this object via \code{\link{inferenceDataObjectTimeBuilder}} function, so that the object is guaranteed to be well defined.
#' @return A list with the following variables:
#' \describe{
#' \item{\code{fit.FEM.time}}{A \code{FEM.time} object that represents the fitted spatio-temporal field.}
#' \item{\code{PDEmisfit.FEM.time}}{A \code{FEM.time} object that represents the misfit of the penalized PDE.}
#' \item{\code{beta}}{If \code{covariates} is not \code{NULL}, a matrix with number of rows equal to the number of covariates and number of columns equal to length of lambda. The \code{j}th column represents the vector of regression coefficients when
#' the smoothing parameter is equal to \code{lambda[j]}.}
#' \item{\code{edf}}{If GCV is \code{TRUE}, a scalar or matrix with the trace of the smoothing matrix for each combination of the smoothing parameters specified in \code{lambdaS} and \code{lambdaT}.}
#' \item{\code{stderr}}{If GCV is \code{TRUE}, a scalar or matrix with the estimate of the standard deviation of the error for each combination of the smoothing parameters specified in \code{lambdaS} and \code{lambdaT}.}
#' \item{\code{GCV}}{If GCV is \code{TRUE}, a scalar or matrix with the value of the GCV criterion for each combination of the smoothing parameters specified in \code{lambdaS} and \code{lambdaT}.}
#' \item{\code{bestlambda}}{If GCV is \code{TRUE}, a 2-elements vector with the indices of smoothing parameters returning the lowest GCV}
#' \item{\code{ICestimated}}{If FLAG_PARABOLIC is \code{TRUE} and IC is \code{NULL}, a list containing a \code{FEM} object with the initial conditions, the value of the smoothing parameter lambda returning the lowest GCV and, in presence of covariates, the estimated beta coefficients}
#' \item{\code{bary.locations}}{A barycenter information of the given locations if the locations are not mesh nodes.}
#' \item{\code{inference}}{A list set only if a well defined \code{\link{inferenceDataObjectTime}} is passed as parameter to the function; contains all inference outputs required:
#' \describe{
#' \item{\code{p_values}}{list of lists set only if at least one p-value is required; contains the p-values divided by implementation:
#' \describe{
#' \item{\code{wald}}{list containing all the Wald p-values required, in the same order of the \code{type} list in \code{inference.data.object.time}.
#' If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of \code{coeff} matrix in \code{inference.data.object.time}.
#' }
#' \item{\code{speckman}}{list containing all the Speckman p-values required, in the same order of the \code{type} list in \code{inference.data.object.time}.
#' If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of \code{coeff} matrix in \code{inference.data.object.time}.
#' }
#' \item{\code{eigen_sign_flip}}{list containing all the Eigen-Sign-Flip p-values required, in the same order of the \code{type} list in \code{inference.data.object.time}.
#' If one-at-the-time tests are required, the corresponding item is a vector of p values ordered as the rows of \code{coeff} matrix in \code{inference.data.object.time}.
#' }
#' }
#' }
#' \item{\code{CI}}{list of lists set only if at least one confidence interval is required; contains the confidence intervals divided by implementation:
#' \describe{
#' \item{\code{wald}}{list containing all the Wald confidence intervals required, in the same order of the \code{type} list in \code{inference.data.object.time}.
#' Each item is a matrix with 3 columns and p rows, p being the number of rows of \code{coeff} matrix in \code{inference.data.object.time}; each row is the CI for the corresponding row of \code{coeff} matrix.
#' }
#' \item{\code{speckman}}{list containing all the Speckman confidence intervals required, in the same order of the \code{type} list in \code{inference.data.object.time}.
#' Each item is a matrix with 3 columns and p rows, p being the number of rows of \code{coeff} matrix in \code{inference.data.object.time}; each row is the CI for the corresponding row of \code{coeff} matrix.
#' }
#' }
#' }
#' }
#' }
#' }
#' @description Space-time regression with differential regularization. Space-varying covariates can be included in the model. The technique accurately handle data distributed over irregularly shaped domains. Moreover, various conditions can be imposed at the domain boundaries.
#' @usage smooth.FEM.time(locations = NULL, time_locations = NULL, observations, FEMbasis,
#' time_mesh=NULL, covariates = NULL, PDE_parameters = NULL, BC = NULL,
#' incidence_matrix = NULL, areal.data.avg = TRUE,
#' FLAG_MASS = FALSE, FLAG_PARABOLIC = FALSE, FLAG_ITERATIVE = FALSE,
#' threshold = 10^(-4), max.steps = 50, IC = NULL,
#' search = "tree", bary.locations = NULL,
#' family = "gaussian", mu0 = NULL, scale.param = NULL,
#' threshold.FPIRLS = 0.0002020, max.steps.FPIRLS = 15,
#' lambda.selection.criterion = "grid", DOF.evaluation = NULL,
#' lambda.selection.lossfunction = NULL, lambdaS = NULL, lambdaT = NULL,
#' DOF.stochastic.realizations = 100, DOF.stochastic.seed = 0,
#' DOF.matrix = NULL, GCV.inflation.factor = 1, lambda.optimization.tolerance = 0.05,
#' inference.data.object.time=NULL)
#' @export
#' @references #' @references Arnone, E., Azzimonti, L., Nobile, F., & Sangalli, L. M. (2019). Modeling
#' spatially dependent functional data via regression with differential regularization.
#' Journal of Multivariate Analysis, 170, 275-295.
#' Bernardi, M. S., Sangalli, L. M., Mazza, G., & Ramsay, J. O. (2017). A penalized
#' regression model for spatial functional data with application to the analysis of the
#' production of waste in Venice province.
#' Stochastic Environmental Research and Risk Assessment, 31(1), 23-38.
#' Federico Ferraccioli, Laura M. Sangalli & Livio Finos (2022). Some first inferential tools for spatial regression
#' with differential regularization. Journal of Multivariate Analysis, 189, 104866.
#' @examples
#' library(fdaPDE)
#'
#' data(horseshoe2D)
#' boundary_nodes = horseshoe2D$boundary_nodes
#' boundary_segments = horseshoe2D$boundary_segments
#' locations = horseshoe2D$locations
#' time_locations = seq(0,1,length.out = 5)
#'
#' mesh = create.mesh.2D(nodes = rbind(boundary_nodes, locations), segments = boundary_segments)
#'
#' space_time_locations = cbind(rep(time_locations,each=nrow(mesh$nodes)),
#' rep(mesh$nodes[,1],5),rep(mesh$nodes[,2],5))
#'
#' FEMbasis = create.FEM.basis(mesh)
#' lambdaS = 10^-1
#' lambdaT = 10^-1
#' data = fs.test(space_time_locations[,2],
#' space_time_locations[,3])*cos(pi*space_time_locations[,1]) +
#' rnorm(nrow(space_time_locations), sd = 0.5)
#' data = matrix(data, nrow = nrow(mesh$nodes), ncol = length(time_locations), byrow = TRUE)
#'
#' solution = smooth.FEM.time(observations = data, time_locations = time_locations,
#' FEMbasis = FEMbasis, lambdaS = lambdaS, lambdaT = lambdaT)
#' plot(solution$fit.FEM)
smooth.FEM.time<-function(locations = NULL, time_locations = NULL, observations, FEMbasis, time_mesh=NULL,
covariates = NULL, PDE_parameters = NULL, BC = NULL,
incidence_matrix = NULL, areal.data.avg = TRUE,
FLAG_MASS = FALSE, FLAG_PARABOLIC = FALSE,FLAG_ITERATIVE = FALSE, threshold = 10^(-4), max.steps = 50, IC = NULL,
search = "tree", bary.locations = NULL,
family = "gaussian", mu0 = NULL, scale.param = NULL,
threshold.FPIRLS = 0.0002020, max.steps.FPIRLS = 15,
lambda.selection.criterion = "grid", DOF.evaluation = NULL, lambda.selection.lossfunction = NULL,
lambdaS = NULL, lambdaT = NULL, DOF.stochastic.realizations = 100, DOF.stochastic.seed = 0, DOF.matrix = NULL, GCV.inflation.factor = 1, lambda.optimization.tolerance = 0.05,
inference.data.object.time = NULL)
{
if(is(FEMbasis$mesh, "mesh.2D"))
{
ndim = 2
mydim = 2
}else if(is(FEMbasis$mesh, "mesh.2.5D"))
{
ndim = 3
mydim = 2
}else if(is(FEMbasis$mesh, "mesh.3D"))
{
ndim = 3
mydim = 3
}else if(is(FEMbasis$mesh, "mesh.1.5D"))
{
ndim = 2
mydim = 1
}else
{
stop('Unknown mesh class')
}
##################### Checking parameters, sizes and conversion ################################
# Preliminary consistency of optimization parameters
if(lambda.selection.criterion == "grid")
{
optim = 0
}else if(lambda.selection.criterion == "newton")
{
optim = 1
}else if(lambda.selection.criterion == "newton_fd")
{
optim = 2
}else
{
stop("'lambda.selection.criterion' must belong to the following list: 'none', 'grid', 'newton', 'newton_fd'.")
}
if(is.null(DOF.evaluation))
{
optim = c(optim,0)
}else if(DOF.evaluation == 'stochastic')
{
optim = c(optim,1)
}else if(DOF.evaluation == 'exact')
{
optim = c(optim,2)
}else
{
stop("'DOF.evaluation' must be NULL, 'stochastic' or 'exact'.")
}
if(is.null(lambda.selection.lossfunction))
{
optim = c(optim,0)
}else if(lambda.selection.lossfunction == 'GCV')
{
optim = c(optim,1)
}else
{
stop("'lambda.selection.lossfunction' has to be 'GCV'.")
}
if(any(lambdaS<=0) || any(lambdaT<=0))
stop("'lambda' can not be less than or equal to 0")
if(optim[2]!=0 & optim[3]!=1)
{
warning("Dof are computed, setting 'lambda.selection.lossfunction' to 'GCV'")
optim[3]=1
}
if(optim[1]==1 & optim[2]!=2)
{
warning("This method needs evaluate DOF in an 'exact' way, selecting 'DOF.evaluation'='exact'")
optim[2] = 2
}
if(!is.null(BC) & optim[1]==1)
{
warning("'newton' 'lambda.selection.criterion' can't be performed with non-NULL boundary conditions, using 'newton_fd' instead")
optim[1] = 2
}
if((optim[1]==2 & optim[2]==0) || (optim[1]==0 & optim[2]==0 & optim[3]==1 & is.null(DOF.matrix)))
{
warning("This method needs evaluate DOF, selecting 'DOF.evaluation'='stochastic'")
optim[2] = 1
}
if(optim[1]!=0 & optim[3]==0)
{
warning("An optimized method needs a loss function to perform the evaluation, selecting 'lambda.selection.lossfunction' as 'GCV'")
optim[3] = 1
}
if(optim[1]>0 & (is.null(lambdaS) || is.null(lambdaT)))
{
warning("No initial point is given: automatic initialization of Newton method")
}
# Search algorithm
if(search=="naive"){
search=1
}else if(search=="tree"){
search=2
}else if(search=="walking" & is(FEMbasis$mesh, "mesh.2.5D")){
stop("walking search is not available for mesh class mesh.2.5D.")
}else if(search=="walking" & is(FEMbasis$mesh, "mesh.1.5D")){
stop("walking search is not available for mesh class mesh.1.5D.")
}else if(search=="walking" & !is(FEMbasis$mesh, "mesh.2.5D") & !is(FEMbasis$mesh, "mesh.1.5D")){
search=3
}else{
stop("'search' must must belong to the following list: 'naive', 'tree' or 'walking'.")
}
# If locations is null but bary.locations is not null, use the locations in bary.locations
if(is.null(locations) & !is.null(bary.locations))
{
locations = bary.locations$locations
locations = as.matrix(locations)
}
## Converting to format for internal usage
if(!is.null(locations))
locations = as.matrix(locations)
if(!is.null(time_locations))
time_locations = as.matrix(time_locations)
if(!is.null(time_mesh))
time_mesh = as.matrix(time_mesh)
observations = as.matrix(observations)
if(!is.null(covariates))
covariates = as.matrix(covariates)
if(!is.null(DOF.matrix))
DOF.matrix = as.matrix(DOF.matrix)
if(!is.null(incidence_matrix))
incidence_matrix = as.matrix(incidence_matrix)
if(!is.null(IC))
IC = as.matrix(IC)
if(!is.null(BC))
{
BC$BC_indices = as.matrix(BC$BC_indices)
BC$BC_values = as.matrix(BC$BC_values)
}
if(!is.null(lambdaS))
lambdaS = as.matrix(lambdaS)
if(!is.null(lambdaT))
lambdaT = as.matrix(lambdaT)
space_varying = checkSmoothingParameters_time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh = time_mesh,
covariates = covariates, PDE_parameters = PDE_parameters, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold, max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance)
# only if inference is required
if(!is.null(inference.data.object.time)){
# Check that GCV is set for inference
if(inference.data.object.time@definition==1 && is.null(lambda.selection.lossfunction) && (!is.numeric(lambdaS) && !is.null(lambdaT)) && (nrow(lambdaS)!=1 || ncol(lambdaS)!=1 || nrow(lambdaT)!=1 || ncol(lambdaT)!=1)){
warning("Inference is not defined when lambda grid is provided without GCV")
inference.data.object.time=new("inferenceDataObjectTime", test = as.integer(0), interval =as.integer(0), type = as.integer(0), component = as.integer(0), exact = as.integer(0), dim = as.integer(0), n_cov = as.integer(0),
locations = matrix(data=0, nrow = 1 ,ncol = 1), locations_indices = as.integer(0), locations_are_nodes = as.integer(0), time_locations = 0, coeff = matrix(data=0, nrow = 1 ,ncol = 1), beta0 = -1, f0 = function(){},
f0_eval = -1, f_var = as.integer(0), quantile = -1, alpha = 0, n_flip = as.integer(1000), tol_fspai = -1, definition=as.integer(0))
}
if(inference.data.object.time@definition==1 && FLAG_ITERATIVE == T){
warning("Inference is not provided when iterative method is selected")
inference.data.object.time=new("inferenceDataObjectTime", test = as.integer(0), interval =as.integer(0), type = as.integer(0), component = as.integer(0), exact = as.integer(0), dim = as.integer(0), n_cov = as.integer(0),
locations = matrix(data=0, nrow = 1 ,ncol = 1), locations_indices = as.integer(0), locations_are_nodes = as.integer(0), time_locations = 0, coeff = matrix(data=0, nrow = 1 ,ncol = 1), beta0 = -1, f0 = function(){},
f0_eval = -1, f_var = as.integer(0), quantile = -1, alpha = 0, n_flip = as.integer(1000), tol_fspai = -1, definition=as.integer(0))
}
}
# If I have PDE non-sv case I need (constant) matrices as parameters
if(!is.null(PDE_parameters) & space_varying==FALSE)
{
PDE_parameters$K = as.matrix(PDE_parameters$K)
PDE_parameters$b = as.matrix(PDE_parameters$b)
PDE_parameters$c = as.matrix(PDE_parameters$c)
}
checkSmoothingParametersSize_time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh = time_mesh,
covariates = covariates, PDE_parameters = PDE_parameters, incidence_matrix = incidence_matrix,
BC = BC, space_varying = space_varying, ndim = ndim, mydim = mydim,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC, IC = IC,
lambdaS = lambdaS, lambdaT = lambdaT, DOF.matrix = DOF.matrix)
# Further check
dim_obs_inf <- dim(observations)[1]
observations<-as.vector(observations)
if(is.null(time_locations))
{
if(FLAG_PARABOLIC && !is.null(IC))
time_locations <- time_mesh[2:length(time_mesh)]
else
time_locations <- time_mesh
}
if(is.null(time_mesh))
{
if(FLAG_PARABOLIC && !is.null(IC))
time_mesh <- rbind(2*time_locations[1]-time_locations[2],time_locations)
else
time_mesh<-time_locations
}
# Checking inference data: after time_locations definition
# Most of the checks have already been carried out by inferenceDataObjectBuilderTime function
if(!is.null(locations))
inference.data.object.time <- checkInferenceParametersTime(inference.data.object.time,ncol(covariates),time_locations,locations,FEMbasis$mesh$nodes, FLAG_PARABOLIC, is.null(IC)) #checking inference data consistency, constructing default object in NULL case
else
inference.data.object.time <- checkInferenceParametersTime(inference.data.object.time,ncol(covariates),time_locations,FEMbasis$mesh$nodes[1:dim_obs_inf,],FEMbasis$mesh$nodes, FLAG_PARABOLIC, is.null(IC)) #checking inference data consistency, constructing default object in NULL case
# Check whether the locations coincide with the mesh nodes (should be put after all the validations)
if (!is.null(locations))
{
if(dim(locations)[1]==dim(FEMbasis$mesh$nodes)[1] & dim(locations)[2]==dim(FEMbasis$mesh$nodes)[2])
{
sum1=0
sum2=0
for (i in 1:nrow(locations))
{
sum1 = sum1 + abs(locations[i,1]-FEMbasis$mesh$nodes[i,1])
sum2 = sum2 + abs(locations[i,2]-FEMbasis$mesh$nodes[i,2])
}
if (sum1==0 & sum2==0)
{
message("No search algorithm is used because the locations coincide with the nodes.")
locations = NULL #In principle, R uses pass-by-value semantics in its function calls. So put ouside of checkSmoothingParameters function.
}
}
}
# OPTIMIZATION NOT IMPLEMENTED FOR GAM
if (family != 'gaussian' & optim[1] != 0)
stop("'lambda.selection.criterion' = 'grid' is the only method implemented for GAM problems")
# FAMILY CHECK
family_admit = c("binomial", "exponential", "gamma", "poisson", "gaussian", "Gaussian")
if (sum(family == family_admit) == 0)
stop("'family' parameter required.\nCheck if it is one of the following: binomial, exponential, gamma, poisson, gaussian")
################## End checking parameters, sizes and conversion #############################
if (family == "gaussian")
{
if(is(FEMbasis$mesh, "mesh.2D") & is.null(PDE_parameters))
{
bigsol = NULL
bigsol = CPP_smooth.FEM.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold , max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.2D") & !is.null(PDE_parameters) & space_varying==FALSE)
{
bigsol = NULL
bigsol = CPP_smooth.FEM.PDE.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, PDE_parameters=PDE_parameters, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold, max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.2D") & !is.null(PDE_parameters) & space_varying==TRUE)
{
bigsol = NULL
bigsol = CPP_smooth.FEM.PDE.sv.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, PDE_parameters=PDE_parameters, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold, max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.2.5D"))
{
bigsol = NULL
bigsol = CPP_smooth.manifold.FEM.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold , max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.3D") & is.null(PDE_parameters))
{
bigsol = NULL
bigsol = CPP_smooth.volume.FEM.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold , max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.3D") & !is.null(PDE_parameters) & space_varying==FALSE)
{
bigsol = NULL
bigsol = CPP_smooth.volume.FEM.PDE.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, PDE_parameters=PDE_parameters, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold, max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.3D") & !is.null(PDE_parameters) & space_varying==TRUE)
{
bigsol = NULL
bigsol = CPP_smooth.volume.FEM.PDE.sv.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, PDE_parameters=PDE_parameters, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold, max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}else if(is(FEMbasis$mesh, "mesh.1.5D"))
{
bigsol = NULL
bigsol = CPP_smooth.graph.FEM.time(locations = locations, time_locations = time_locations, observations = observations, FEMbasis = FEMbasis, time_mesh=time_mesh,
covariates = covariates, ndim = ndim, mydim = mydim, BC = BC,
incidence_matrix = incidence_matrix, areal.data.avg = areal.data.avg,
FLAG_MASS = FLAG_MASS, FLAG_PARABOLIC = FLAG_PARABOLIC,FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold , max.steps = max.steps, IC = IC,
search = search, bary.locations = bary.locations,
optim = optim, lambdaS = lambdaS, lambdaT = lambdaT, DOF.stochastic.realizations = DOF.stochastic.realizations, DOF.stochastic.seed = DOF.stochastic.seed, DOF.matrix = DOF.matrix, GCV.inflation.factor = GCV.inflation.factor, lambda.optimization.tolerance = lambda.optimization.tolerance,
inference.data.object = inference.data.object.time)
}
# ---------- Solution -----------
N = nrow(FEMbasis$mesh$nodes)
M = ifelse(FLAG_PARABOLIC,length(time_mesh)-1,length(time_mesh) + 2);
dim_1 = ifelse(optim[1]==0 & is.null(DOF.matrix) & optim[3]==0, length(lambdaS), 1)
dim_2 = ifelse(optim[1]==0 & is.null(DOF.matrix) & optim[3]==0, length(lambdaT), 1)
if(is.null(IC) && FLAG_PARABOLIC)
IC = bigsol[[27]]$coeff
if(FLAG_PARABOLIC)
{
f = array(dim=c(length(IC)+M*N,dim_1,dim_2))
g = array(dim=c(length(IC)+M*N,dim_1,dim_2))
for (i in 1:dim_1)
for (j in 1:dim_2)
{
f[,i,j] = c(IC,bigsol[[1]][1:(N*M),i+(j-1)*dim_1])
g[,i,j] = c(rep(0,length(IC)),bigsol[[1]][(N*M+1):(2*N*M),i+(j-1)*dim_1])
}
}
else
{
f = array(data=bigsol[[1]][1:(N*M),],dim = c(N*M,dim_1,dim_2))
g = array(data=bigsol[[1]][(N*M+1):(2*N*M),],dim = c(N*M,dim_1,dim_2))
}
GCV_ = bigsol[[3]]
if(!is.null(covariates))
{
if(optim[1]==0 & is.null(DOF.matrix) & optim[3]==0)
beta = array(data=bigsol[[5]],dim=c(ncol(covariates),length(lambdaS),length(lambdaT)))
else
beta = array(data=bigsol[[5]],dim=c(ncol(covariates),1,1))
}
else
beta = NULL
if(all(is.na(bigsol[[27]])))
ICestimated = NULL
else
ICestimated = list(IC.FEM=bigsol[[27]],bestlambdaindex=bigsol[[28]],bestlambda=bigsol[[29]],beta=bigsol[[30]])
bestlambda = bigsol[[4]]+1
if(optim[1]!=0) # newton or newton_fd
{
bestlambda = max(bestlambda[1], bestlambda[2])
}
if(optim[1]==0)
{
if(bestlambda[1] == 1 || bestlambda[1] == length(lambdaS))
warning("Your optimal 'GCV' is on the border of lambdaS sequence")
if(bestlambda[2] == 1 || bestlambda[2] == length(lambdaT))
warning("Your optimal 'GCV' is on the border of lambdaT sequence")
}
if (is.null(lambda.selection.lossfunction))
sd = -1
else
sd = sqrt(bigsol[[15]])
solution = list(
f = f,
g = g,
z_hat = bigsol[[13]],
rmse = bigsol[[14]],
estimated_sd=sd
)
term = bigsol[[19]]
ot = bigsol[[20]]
if(term == 1)
termination = "reached tolerance"
else if(term == 2)
termination = "reached max number iterations"
else
termination = "uninformative"
if(ot == 0)
optimization_type = "full optimization"
else if(ot == 1)
optimization_type = "full DOF grid"
else
optimization_type = "uninformative"
optimization = list(
lambda_solution = bigsol[[16]],
lambda_position = bestlambda,
optimization_details = list(
iterations = bigsol[[18]],
termination = termination,
optimization_type = optimization_type),
dof = bigsol[[2]],
lambda_vector = matrix(c(bigsol[[21]], bigsol[[22]]), nrow=length(bigsol[[21]]), ncol=2),
GCV_vector = GCV_
)
time = bigsol[[23]]
# Save information of Tree Mesh
tree_mesh = list(
treelev = bigsol[[6]][1],
header_orig= bigsol[[7]],
header_scale = bigsol[[8]],
node_id = bigsol[[9]][,1],
node_left_child = bigsol[[9]][,2],
node_right_child = bigsol[[9]][,3],
node_box= bigsol[[10]])
# Reconstruct FEMbasis with tree mesh
mesh.class= class(FEMbasis$mesh)
if (is.null(FEMbasis$mesh$treelev)) { #if doesn't exist the tree information
FEMbasis$mesh = append(FEMbasis$mesh, tree_mesh)
} #if already exist the tree information, don't append
class(FEMbasis$mesh) = mesh.class
# Save information of Barycenter
if (is.null(bary.locations)) {
bary.locations = list(locations=locations, element_ids = bigsol[[11]], barycenters = bigsol[[12]])
}
class(bary.locations) = "bary.locations"
# Make FEM.time objects
fit.FEM.time = FEM.time(f, time_mesh, FEMbasis, FLAG_PARABOLIC)
PDEmisfit.FEM.time = FEM.time(g, time_mesh, FEMbasis, FLAG_PARABOLIC)
# Save statistics and intervals
if(inference.data.object.time@definition==1){
inference = {}
n_loc_inference = ifelse(FLAG_PARABOLIC==TRUE,dim(inference.data.object.time@locations)[1]*(length(inference.data.object.time@time_locations)-1),dim(inference.data.object.time@locations)[1]*length(inference.data.object.time@time_locations))
confidence_intervals = matrix(data = bigsol[[25]], nrow = 2*3*length(inference.data.object.time@type), ncol = max(dim(inference.data.object.time@coeff)[1], n_loc_inference))
p_val = matrix(data = bigsol[[24]], nrow = dim(inference.data.object.time@coeff)[1]+1, ncol = length(inference.data.object.time@type))
for(i in 1:length(inference.data.object.time@type)){ # each element is a different inferential setting
if(inference.data.object.time@interval[i]!=0){ # Intervals requested by this setting, adding them to the right implementation position
if(inference.data.object.time@component[i]!=2){ # intervals for beta were requested, just for safety, actually it cannot be equal to 2
p = dim(inference.data.object.time@coeff)[1]
ci_beta=t(confidence_intervals[(3*(2*i-2)+1):(3*(2*i-2)+3),1:p])
if(inference.data.object.time@type[i]==1){
inference$beta$CI$wald[[length(inference$beta$CI$wald)+1]] = ci_beta
inference$beta$CI$wald=as.list(inference$beta$CI$wald)
}
else if(inference.data.object.time@type[i]==2){
inference$beta$CI$speckman[[length(inference$beta$CI$speckman)+1]] = ci_beta
inference$beta$CI$speckman=as.list(inference$beta$CI$speckman)
}
else if(inference.data.object.time@type[i]==3){
if(ci_beta[2]> 10^20){
warning("ESF CI bisection algorithm did not converge, returning NA")
for(h in 1:nrow(ci_beta))
for(k in 1:ncol(ci_beta))
ci_beta[h,k]=NA
}
inference$beta$CI$eigen_sign_flip[[length(inference$beta$CI$eigen_sign_flip)+1]] = ci_beta
inference$beta$CI$eigen_sign_flip=as.list(inference$beta$CI$eigen_sign_flip)
}
else if(inference.data.object.time@type[i]==4){
if(ci_beta[2]> 10^20){
warning("Enhanced ESF CI bisection algorithm did not converge, returning NA")
for(h in 1:nrow(ci_beta))
for(k in 1:ncol(ci_beta))
ci_beta[h,k]=NA
}
inference$beta$CI$enh_eigen_sign_flip[[length(inference$beta$CI$enh_eigen_sign_flip)+1]] = ci_beta
inference$beta$CI$enh_eigen_sign_flip=as.list(inference$beta$CI$enh_eigen_sign_flip)
}
}
if(inference.data.object.time@component[i]!=1){ # intervals for f were requested
n_loc = dim(inference.data.object.time@locations)[1]*length(inference.data.object.time@time_locations)
ci_f=t(confidence_intervals[(3*(2*i-1)+1):(3*(2*i-1)+3),])
if(inference.data.object.time@type[i]==1){ # wald confidence intervals for f
inference$f$CI$wald[[length(inference$f$CI$wald)+1]] = ci_f
inference$f$CI$wald=as.list(inference$f$CI$wald)
}
}
}
if(inference.data.object.time@test[i]!=0){ # Test requested by this setting, adding them to the right implementation position
statistics=p_val[,i]
if(inference.data.object.time@component[i]!=2){ # test on beta was requested, just for safety, it cannot be otherwise
beta_statistics = statistics[1:dim(inference.data.object.time@coeff)[1]]
p_values = numeric()
if(inference.data.object.time@type[i]==3 || inference.data.object.time@type[i]==4){ # eigen-sign-flip p-value is already computed in cpp code
if(inference.data.object.time@test[i]==1){
# one-at-the-time tests
p_values = beta_statistics
}
else if(inference.data.object.time@test[i]==2){
# simultaneous test
p_values = beta_statistics[1]
}
}else{
# Compute p-values
if(inference.data.object.time@test[i]==1){ # Wald and Speckman return statistics and needs computation of p-values (no internal use of distributions quantiles)
# one-at-the-time-tests
p_values = numeric(length(beta_statistics))
for(l in 1:length(beta_statistics)){
p_values[l] = 2*pnorm(-abs(beta_statistics[l]))
}
}
else if(inference.data.object.time@test[i]==2){
# simultaneous tests
p = dim(inference.data.object.time@coeff)[1]
p_values = 1-pchisq(beta_statistics[1], p)
}
}
# add p-values in the right position
if(inference.data.object.time@type[i]==1){
inference$beta$p_values$wald[[length(inference$beta$p_values$wald)+1]] = p_values
inference$beta$p_values$wald=as.list(inference$beta$p_values$wald)
}
else if(inference.data.object.time@type[i]==2){
inference$beta$p_values$speckman[[length(inference$beta$p_values$speckman)+1]] = p_values
inference$beta$p_values$speckman=as.list(inference$beta$p_values$speckman)
}
else if(inference.data.object.time@type[i]==3){
inference$beta$p_values$eigen_sign_flip[[length(inference$beta$p_values$eigen_sign_flip)+1]] = p_values
inference$beta$p_values$eigen_sign_flip=as.list(inference$beta$p_values$eigen_sign_flip)
}
else if(inference.data.object.time@type[i]==4){
inference$beta$p_values$enh_eigen_sign_flip[[length(inference$beta$p_values$enh_eigen_sign_flip)+1]] = p_values
inference$beta$p_values$enh_eigen_sign_flip=as.list(inference$beta$p_values$enh_eigen_sign_flip)
}
}
if(inference.data.object.time@component[i]!=1){ # test on f was requested
p_value = statistics[length(statistics)]
# add p-value in the right position
if(inference.data.object.time@type[i]==1){
inference$f$p_values$wald[[length(inference$f$p_values$wald)+1]] = p_value
inference$f$p_values$wald=as.list(inference$f$p_values$wald)
}
}
}
}
if(inference.data.object.time@f_var==1){
f_variances = matrix(data = bigsol[[26]], nrow = length(solution$z_hat), ncol = 1)
inference$f_var = f_variances
}
# Prepare return list
reslist = list(fit.FEM.time = fit.FEM.time, PDEmisfit.FEM.time = PDEmisfit.FEM.time, solution = solution,
optimization = optimization, beta = beta, time = time, ICestimated=ICestimated, bary.locations = bary.locations, inference=inference)
}else{
# Prepare return list
reslist = list(fit.FEM.time = fit.FEM.time, PDEmisfit.FEM.time = PDEmisfit.FEM.time, solution = solution,
optimization = optimization, beta = beta, time = time, ICestimated=ICestimated, bary.locations = bary.locations)
}
return(reslist)
}
else
{
#Iterative solver not implemented for GAM
if(FLAG_ITERATIVE)
warning("Iterative solver not implemented for GAM, switching to block solver")
FLAG_ITERATIVE = FALSE
threshold = 10^(-4)
max.steps = 50
#----------------------------------------------------#
############# GAMs: FPIRLS algorithm #################
#----------------------------------------------------#
checkGAMParameters(observations = observations, max.steps.FPIRLS = max.steps.FPIRLS,
mu0 = mu0, scale.param = scale.param,
threshold.FPIRLS = threshold.FPIRLS, family = family)
if (is(FEMbasis$mesh, "mesh.2D") & is.null(PDE_parameters)) {
bigsol = NULL
bigsol = CPP_smooth.GAM.FEM.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates, ndim = ndim,
mydim = mydim, BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if (is(FEMbasis$mesh, "mesh.2D") & !is.null(PDE_parameters) & space_varying == FALSE) {
bigsol = NULL
bigsol = CPP_smooth.GAM.FEM.PDE.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates,
PDE_parameters = PDE_parameters, ndim = ndim, mydim = mydim,
BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE = FLAG_ITERATIVE,
threshold = threshold, max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if (is(FEMbasis$mesh, "mesh.2D") & !is.null(PDE_parameters) & space_varying == TRUE) {
bigsol = NULL
bigsol = CPP_smooth.GAM.FEM.PDE.sv.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates,
PDE_parameters = PDE_parameters, ndim = ndim, mydim = mydim,
BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE = FLAG_ITERATIVE,
threshold = threshold, max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
} else if( is(FEMbasis$mesh, "mesh.2.5D") & is.null(PDE_parameters)){
bigsol = NULL
bigsol = CPP_smooth.manifold.GAM.FEM.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates, ndim = ndim,
mydim = mydim, BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if(is(FEMbasis$mesh, "mesh.3D") & is.null(PDE_parameters)){
bigsol = NULL
bigsol = CPP_smooth.volume.GAM.FEM.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates, ndim = ndim,
mydim = mydim, BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if(is(FEMbasis$mesh, "mesh.3D") & !is.null(PDE_parameters) & space_varying == FALSE){
bigsol = NULL
bigsol = CPP_smooth.volume.GAM.FEM.PDE.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates,
PDE_parameters = PDE_parameters, ndim = ndim, mydim = mydim,
BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if(is(FEMbasis$mesh, "mesh.3D") & !is.null(PDE_parameters) & space_varying == TRUE){
bigsol = NULL
bigsol = CPP_smooth.volume.GAM.FEM.PDE.sv.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates,
PDE_parameters = PDE_parameters, ndim = ndim, mydim = mydim,
BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else if(is(FEMbasis$mesh, "mesh.1.5D") & is.null(PDE_parameters)){
bigsol = NULL
bigsol = CPP_smooth.graph.GAM.FEM.time(locations = locations, time_locations = time_locations,
observations = observations, FEMbasis = FEMbasis,
time_mesh = time_mesh, covariates = covariates, ndim = ndim,
mydim = mydim, BC = BC, incidence_matrix = incidence_matrix,
areal.data.avg = areal.data.avg, FLAG_MASS = FLAG_MASS,
FLAG_PARABOLIC = FLAG_PARABOLIC, FLAG_ITERATIVE=FLAG_ITERATIVE, threshold = threshold,
max.steps = max.steps, IC = IC, FAMILY = family,
mu0 = mu0, max.steps.FPIRLS = max.steps.FPIRLS,
scale.param = scale.param, threshold.FPIRLS = threshold.FPIRLS,
search = search, bary.locations = bary.locations, optim = optim,
lambdaS = lambdaS, lambdaT = lambdaT,
DOF.stochastic.realizations = DOF.stochastic.realizations,
DOF.stochastic.seed = DOF.stochastic.seed,
DOF.matrix = DOF.matrix,
GCV.inflation.factor = GCV.inflation.factor,
lambda.optimization.tolerance = lambda.optimization.tolerance)
}else{
stop("Not implemented for !is.null(PDE_parameters). Try Laplacian regularization.")
}
ICindx = 16
N = nrow(FEMbasis$mesh$nodes)
M = ifelse(FLAG_PARABOLIC, length(time_mesh) - 1, length(time_mesh) + 2)
if (is.null(IC) && FLAG_PARABOLIC)
IC = bigsol[[ICindx]]$coeff[, bigsol[[ICindx + 1]]]
if (FLAG_PARABOLIC) {
f = array(dim = c(length(IC) + M * N, length(lambdaS), length(lambdaT)))
for (i in 1:length(lambdaS))
for (j in 1:length(lambdaT))
f[, i, j] = c(IC, bigsol[[1]][1:(N * M), i + (j - 1) * length(lambdaS)])
}else
f = array(data = bigsol[[1]][1:(N * M), ],
dim = c(N * M, length(lambdaS), length(lambdaT)))
if (FLAG_PARABOLIC) {
g = array(dim = c(length(IC) + M * N, length(lambdaS), length(lambdaT)))
for (i in 1:length(lambdaS))
for (j in 1:length(lambdaT))
g[, i, j] = c(rep(0, length(IC)), bigsol[[1]][(N * M + 1):(2 * N * M),
i + (j - 1) * length(lambdaS)])
}else
g = array(data = bigsol[[1]][(N * M + 1):(2 * N * M), ],
dim = c(N * M, length(lambdaS), length(lambdaT)))
dof = bigsol[[2]]
GCV_ = bigsol[[3]]
bestlambda = bigsol[[4]] + 1
if (!is.null(covariates))
beta = array(data = bigsol[[5]], dim = c(ncol(covariates), length(lambdaS), length(lambdaT)))
else
beta = NULL
if (all(is.na(bigsol[[ICindx]])))
ICestimated = NULL
else
ICestimated = list(IC.FEM = bigsol[[ICindx]], bestlambdaindex = bigsol[[ICindx + 1]],
bestlambda = bigsol[[ICindx + 2]], beta = bigsol[[ICindx + 3]])
# Save information of Tree Mesh
tree_mesh = list(treelev = bigsol[[6]][1], header_orig = bigsol[[7]],
header_scale = bigsol[[8]], node_id = bigsol[[9]][, 1],
node_left_child = bigsol[[9]][, 2],
node_right_child = bigsol[[9]][, 3], node_box = bigsol[[10]])
# Reconstruct FEMbasis with tree mesh
mesh.class = class(FEMbasis$mesh)
if (is.null(FEMbasis$mesh$treelev)) { # if doesn't exist the tree
# information
FEMbasis$mesh = append(FEMbasis$mesh, tree_mesh)
} # if already exist the tree information, don't append
class(FEMbasis$mesh) = mesh.class
# Save information of Barycenter
if (is.null(bary.locations)) {
bary.locations = list(locations = locations, element_ids = bigsol[[11]],
barycenters = bigsol[[12]])
}
class(bary.locations) = "bary.locations"
# Make FEM.time objects
fit.FEM.time = FEM.time(f, time_mesh, FEMbasis, FLAG_PARABOLIC)
PDEmisfit.FEM.time = FEM.time(g, time_mesh, FEMbasis, FLAG_PARABOLIC)
# Prepare return list
reslist = NULL
nm <- length(observations)
fn.eval = array(dim = c(nm, length(lambdaS), length(lambdaT)))
if(FLAG_PARABOLIC)
{
if(!is.null(bigsol[[ICindx]]))
ICfn.eval = bigsol[[ICindx+4]][, bigsol[[ICindx + 1]]]
else
ICfn.eval = NULL
nm <- nm - length(ICfn.eval)
for (i in 1:length(lambdaS))
for (j in 1:length(lambdaT))
fn.eval[, i, j] = c(ICfn.eval, bigsol[[13]][1:nm, i + (j - 1) * length(lambdaS)])
}else{
for (i in 1:length(lambdaS))
for (j in 1:length(lambdaT))
fn.eval[, i, j] = bigsol[[13]][1:nm, i + (j - 1) * length(lambdaS)]
}
J_minima = bigsol[[14]]
variance.est = bigsol[[15]]
if(!is.numeric(variance.est[1])){
variance.est <- NULL
} else if (variance.est[1] < 0) variance.est = NULL
if (!is.null(lambda.selection.lossfunction)) {
stderr = sqrt(GCV_ * (sum(!is.na(observations)) - dof) / sum(!is.na(observations)))
reslist = list(fit.FEM.time = fit.FEM.time,
PDEmisfit.FEM.time = PDEmisfit.FEM.time, beta = beta,
edf = dof, GCV = GCV_, stderr = stderr,
bestlambda = bestlambda, ICestimated = ICestimated,
bary.locations = bary.locations, fn.eval = fn.eval,
J_minima = J_minima, variance.est = variance.est)
}else{
reslist = list(fit.FEM.time = fit.FEM.time,
PDEmisfit.FEM.time = PDEmisfit.FEM.time, beta = beta,
ICestimated = ICestimated, bary.locations = bary.locations,
fn.eval = fn.eval, J_minima = J_minima,
variance.est = variance.est)
}
return(reslist)
}
}
Try the fdaPDE package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.