diff --git a/R/basisfns.R b/R/basisfns.R
index 0e46e80..86b9000 100644
--- a/R/basisfns.R
+++ b/R/basisfns.R
@@ -26,7 +26,7 @@
#' to location and scale parameters.
#' @param df A data frame containing one row per basis function, typically for providing informative summaries.
#' @param regular logical indicating if the basis functions (of each resolution) are in a regular grid
-#' @details This constructor checks that all parameters are valid before constructing the basis functions.
+#' @details This constructor checks that all parameters are valid before constructing the basis functions.
#' The requirement that every function is encapsulated is tedious, but necessary for
#' FRK to work with a large range of basis functions in the future. Please see the example below which exemplifies
#' the process of constructing linear basis functions from scratch using this function.
@@ -62,11 +62,11 @@ Basis <- function(manifold, n, fn, pars, df, regular = FALSE) {
stop("manifold needs to be in the environment of every function")
if(!all(sapply(fn, function(f) names(formals(f)) == "s")))
stop("all functions need to take s as input argument")
-
+
# if(!is.logical(regular)) stop("regular should be of class logical")
## regular has been coded to 0 and 1. So, convert to numeric before constructing.
regular <- as.numeric(regular)
-
+
new("Basis", manifold = manifold, n = n, fn = fn, pars = pars, df = df, regular = regular)
}
@@ -98,23 +98,23 @@ Basis <- function(manifold, n, fn, pars, df, regular = FALSE) {
#' @export
local_basis <- function(manifold = sphere(), # default manifold is sphere
loc = matrix(c(1,0),nrow=1), # one centroid at (1,0)
- scale = 1,
- type = c("bisquare", "Gaussian", "exp", "Matern32"),
- res = 1,
+ scale = 1,
+ type = c("bisquare", "Gaussian", "exp", "Matern32"),
+ res = 1,
regular = FALSE) {
-
+
regular <- as.numeric(regular)
-
+
## Basic checks
if(!is.matrix(loc)) stop("loc needs to be a matrix")
if(!(dimensions(manifold) == ncol(loc))) stop("number of columns in loc needs to be the same as the number of manifold dimensions")
if(!(length(scale) == nrow(loc))) stop("need to have as many scale parameters as centroids")
type <- match.arg(type)
-
+
n <- nrow(loc) # n is the number of centroids
colnames(loc) <- c(outer("loc",1:ncol(loc), # label dimensions as loc1,loc2,...,locn
FUN = paste0))
-
+
fn <- pars <- list() # initialise lists of functions and parameters
for (i in 1:n) {
## Assign the functions as determined by user
@@ -127,15 +127,15 @@ local_basis <- function(manifold = sphere(), # default manifold is sphe
} else if (type=="Matern32") {
fn[[i]] <- .Matern32_wrapper(manifold,matrix(loc[i,],nrow=1),scale[i])
}
-
+
## Save the parameters to the parameter list (loc and scale)
pars[[i]] <- list(loc = matrix(loc[i,],nrow=1), scale=scale[i])
}
-
+
## Create a data frame which summarises info about the functions, and set the resolution.
df <- data.frame(loc, scale, res = res)
-
- ## For temporal basis functions, which are often constructed using
+
+ ## For temporal basis functions, which are often constructed using
## local_basis(), we'll set regular = TRUE if the following conditions hold:
if (is(manifold, "real_line") && length(unique(res)) == 1 && length(unique(scale)) == 1 && length(unique(diff(loc))) == 1) {
regular = TRUE
@@ -163,30 +163,30 @@ local_basis <- function(manifold = sphere(), # default manifold is sphe
#' @param tunit temporal unit, required when constructing a spatio-temporal basis. Should be the same as used for the BAUs. Can be "secs", "mins", "hours", "days", "years", etc.
#' @param ... unused
#' @details This function automatically places basis functions within the domain of interest. If the domain is a plane or the real line, then the object \code{data} is used to establish the domain boundary.
-#'
-#'Let \eqn{\phi(u)} denote the value of a basis function evaluated at \eqn{u = s - c},
-#'where \eqn{s} is a spatial coordinate and \eqn{c} is the basis-function centroid.
-#'The argument \code{type} can be either ``Gaussian'', in which case
-#'
+#'
+#'Let \eqn{\phi(u)} denote the value of a basis function evaluated at \eqn{u = s - c},
+#'where \eqn{s} is a spatial coordinate and \eqn{c} is the basis-function centroid.
+#'The argument \code{type} can be either ``Gaussian'', in which case
+#'
#'\ifelse{html}{{\out{ φ(u) = exp(-|u|²/2σ²),
}}}{\deqn{\phi(u) = \exp\left(-\frac{\|u \|^2}{2\sigma^2}\right),}{phi(u) = exp(-|u|^2/2sigma^2),}}
-#'
+#'
#'``bisquare'', in which case
-#'
+#'
#'\ifelse{html}{{\out{ φ(u) = (1 -(|u|/R)²)²,
}}}{\deqn{\phi(u) = \left(1- \left(\frac{\| u \|}{R}\right)^2\right)^2 I(\|u\| < R),}{\phi(u) = (1- (|u|/R)^2)^2 I(|u| < R),}}
-#'
+#'
#'``exp'', in which case
-#'
+#'
#'\ifelse{html}{{\out{ φ(u) = exp(-|u|/τ),
}}}{\deqn{\phi(u) = \exp\left(-\frac{\|u \|}{\tau}\right),}{phi(u) = exp(-|u|/tau),}}
-#'
+#'
#' or ``Matern32'', in which case
-#'
+#'
#' \ifelse{html}{{\out{ φ(u) = (1 + √3|u|/κ)exp(-√3|u|/κ),
}}}{\deqn{\phi(u) = \left(1 + \frac{\sqrt{3}\|u\|}{\kappa}\right)\exp\left(-\frac{\sqrt{3}\| u \|}{\kappa}\right),}{\phi(u) = (1 + \sqrt{3}|u|/\kappa)\exp(-\sqrt{3}|u|/\kappa),}}
-#'
+#'
#' where the parameters \eqn{\sigma, R, \tau} and \eqn{\kappa} are \code{scale} arguments.
#'
#' If the manifold is the real line, the basis functions are placed regularly inside the domain, and the number of basis functions at the coarsest resolution is dictated by the integer parameter \code{regular} which has to be greater than zero. On the real line, each subsequent resolution has twice as many basis functions. The scale of the basis function is set based on the minimum distance between the centre locations following placement. The scale is equal to the minimum distance if the type of basis function is Gaussian, exponential, or Matern32, and is equal to 1.5 times this value if the function is bisquare.
#'
-#' If the manifold is a plane, and \code{regular > 0}, then basis functions are placed regularly within the bounding box of \code{data}, with the smallest number of basis functions in each row or column equal to the value of \code{regular} in the coarsest resolution (note, this is just the smallest number of basis functions). Subsequent resolutions have twice the number of basis functions in each row or column. If \code{regular = 0}, then the function \code{fmesher::fm_nonconvex_hull_inla()} is used to construct a (non-convex) hull around the data. The buffer and smoothness of the hull is determined by the parameter \code{convex}. Once the domain boundary is found, \code{fmesher::fm_mesh_2d_inla()} is used to construct a triangular mesh such that the node vertices coincide with data locations, subject to some minimum and maximum triangular-side-length constraints. The result is a mesh that is dense in regions of high data density and not dense in regions of sparse data. Even basis functions are irregularly placed, the scale is taken to be a function of the minimum distance between basis function centres, as detailed above. This may be changed in a future revision of the package.
+#' If the manifold is a plane, and \code{regular > 0}, then basis functions are placed regularly within the bounding box of \code{data}, with the smallest number of basis functions in each row or column equal to three times the value of \code{regular} in the coarsest resolution. Subsequent resolutions have thrice the number of basis functions in each row or column. If \code{regular = 0}, then the function \code{fmesher::fm_nonconvex_hull_inla()} is used to construct a (non-convex) hull around the data. The buffer and smoothness of the hull is determined by the parameter \code{convex}. Once the domain boundary is found, \code{fmesher::fm_mesh_2d_inla()} is used to construct a triangular mesh such that the node vertices coincide with data locations, subject to some minimum and maximum triangular-side-length constraints. The result is a mesh that is dense in regions of high data density and not dense in regions of sparse data. Even basis functions are irregularly placed, the scale is taken to be a function of the minimum distance between basis function centres, as detailed above. This may be changed in a future revision of the package.
#'
#' If the manifold is the surface of a sphere, then basis functions are placed on the centroids of the discrete global grid (DGG), with the first basis resolution corresponding to the third resolution of the DGG (ISEA3H resolution 2, which yields 92 basis functions globally). It is not recommended to go above \code{nres == 3} (ISEA3H resolutions 2--4) for the whole sphere; \code{nres=3} yields a total of 1176 basis functions. Up to ISEA3H resolution 6 is available with \code{FRK}; for finer resolutions; please install \code{dggrids} from \code{https://github.com/andrewzm/dggrids} using \code{devtools}.
#'
@@ -230,16 +230,16 @@ auto_basis <- function(manifold = plane(),
buffer = 0,
tunit = NULL,
...) { # currently unused
-
+
## If called from FRK(), ... may contain regular and buffer.
if("buffer" %in% names(list(...))) buffer <- list(...)$buffer
if("regular" %in% names(list(...))) regular <- list(...)$regular
-
+
## If the data is spatio-temporal, construct spatio-temporal basis functions
- ## using the tensor product of spatial and temporal basis functions.
- ## Note that we put this here, because there is a call to auto_basis()
- ## when constructing the spatial basis, at which point the arguments will be
- ## checked.
+ ## using the tensor product of spatial and temporal basis functions.
+ ## Note that we put this here, because there is a call to auto_basis()
+ ## when constructing the spatial basis, at which point the arguments will be
+ ## checked.
if (is(manifold, "STplane") || is(manifold, "STsphere")) {
if(is.null(tunit)) {
stop("Please specify tunit (the temporal unit) if you want to construct a spatio-temporal basis. Note that this should be the same as the BAUs.")
@@ -250,18 +250,18 @@ auto_basis <- function(manifold = plane(),
} else if (is(manifold, "STsphere")) {
spatial_manifold <- sphere()
}
- B_spatial <- auto_basis(manifold = plane(), data = as(data, "Spatial"),
- regular = regular, nres = nres, prune = prune,
- max_basis = max_basis, subsamp = subsamp,
- type = type, isea3h_lo = isea3h_lo, bndary = bndary,
- scale_aperture = scale_aperture, verbose = verbose,
+ B_spatial <- auto_basis(manifold = plane(), data = as(data, "Spatial"),
+ regular = regular, nres = nres, prune = prune,
+ max_basis = max_basis, subsamp = subsamp,
+ type = type, isea3h_lo = isea3h_lo, bndary = bndary,
+ scale_aperture = scale_aperture, verbose = verbose,
buffer = buffer)
-
+
## Set up the temporal basis functions
end_loc <- max(data@endTime) - min(data@endTime)
- units(end_loc) <- "secs"
-
- ## We always compute the difference in terms of times. Then, we convert to
+ units(end_loc) <- "secs"
+
+ ## We always compute the difference in terms of times. Then, we convert to
## whatever unit was specified.
end_loc <- as.numeric(end_loc)
if (tunit == "secs") {
@@ -282,19 +282,19 @@ auto_basis <- function(manifold = plane(),
stop("tunit should be a one of: secs, mins, hours, days, weeeks, months, years")
}
end_loc <- ceiling(end_loc) + 1 # add 1 as a buffer
-
- B_temporal <- local_basis(manifold = real_line(),
- type = "Gaussian",
- loc = matrix(seq(0, end_loc, length.out = 10)),
- scale = rep(end_loc / 10, 10), regular = TRUE)
-
+
+ B_temporal <- local_basis(manifold = real_line(),
+ type = "Gaussian",
+ loc = matrix(seq(0, end_loc, length.out = 10)),
+ scale = rep(end_loc / 10, 10), regular = TRUE)
+
## Spatio-temporal basis functions generated with a tensor product
- B <- TensorP(B_spatial, B_temporal)
-
+ B <- TensorP(B_spatial, B_temporal)
+
return(B)
}
-
-
+
+
## Basic checks
type <- match.arg(type)
if(!is(manifold,"manifold"))
@@ -317,8 +317,8 @@ auto_basis <- function(manifold = plane(),
stop("buffer needs to be between 0 and 0.5")
if (buffer != 0 & regular != 1)
stop("A buffer of basis functions along the boundary can only be computed if regular == 1")
-
-
+
+
## If the user has specified a maximum number of basis functions then
## we need to add resolutions iteratively and stop when exceed the maximum
## of basis functions
@@ -340,18 +340,18 @@ auto_basis <- function(manifold = plane(),
bndary = bndary,
scale_aperture = scale_aperture,
verbose = 0) # force verbose to 0 for this procedure of finding the number of basis functions
-
+
tot_basis <- nbasis(G) # record the number of basis functions
}
nres <- nres - 1 # nres resolutions was too much, deduct by one
}
-
+
## Now call the local function with checked parameters
return(.auto_basis(manifold = manifold, data = data, regular = regular, nres = nres,
prune = prune, subsamp = subsamp, type = type, isea3h_lo = isea3h_lo,
- bndary = bndary, scale_aperture = scale_aperture, verbose = verbose,
+ bndary = bndary, scale_aperture = scale_aperture, verbose = verbose,
buffer = buffer))
-
+
}
@@ -365,15 +365,15 @@ auto_basis <- function(manifold = plane(),
isea3h_lo = 2,
bndary = NULL,
scale_aperture = ifelse(is(manifold,"sphere"),1,1.25),
- verbose = 0L,
+ verbose = 0L,
buffer = 0) {
-
-
+
+
type <- match.arg(type) # match type
m <- manifold # abbreviate for convenience
isea3h <- centroid <- res <- NULL # (suppress warnings, these are loaded from data)
coords <- coordinates(data) # data coordinates or centroids
-
+
## Irregular basis function placement can only proceed with fmesher. Throw an error
## if fmesher is not installed
if(is(m,"plane") & !regular & is.null(bndary)) {
@@ -381,17 +381,17 @@ auto_basis <- function(manifold = plane(),
stop("For irregularly-placed basis-function generation fmesher needs to be installed
for constructing basis function centres.")
}
-
+
## Subsamp can be used to make the placement/pruning process more efficient with big data.
## If subsamp < number of data points then sample the number of data randomly.
if(nrow(coords)>subsamp) {
coords <- coords[sample(nrow(coords),size=subsamp,replace=FALSE),]
}
-
+
## Find the x and y extent of the (subsampled) data
xrange <- range(coords[,1])
yrange <- range(coords[,2])
-
+
## If we are on the plane and want an irregular function placement, then call fmesher
## and find a nonconvex hull in which the basis functions can be enclosed. If
## a boundary is supplied as a matrix, convert it a segment object.
@@ -404,7 +404,7 @@ auto_basis <- function(manifold = plane(),
bndary_seg <- fmesher::fm_segm(bndary)
}
}
-
+
## If we want basis functions regularly placed on the plane then
## first calculate the aspect ratio (extent of y / extent of x)
## Then set the number of basis functions in y = regular if
@@ -421,13 +421,13 @@ auto_basis <- function(manifold = plane(),
ny <- nx * asp_ratio
}
}
-
+
## Initialise the variables loc and scale which we will fill in later (centroids and scales)
loc <- scale <- NULL
-
+
## Initialise G which will contain a list of basis function objects by resolution (1 resolution per list item)
G <- list()
-
+
## If we are on the sphere, then load the dggrids. This will be either
## from FRK if isea3h_lo + nres - 1 <= 6, or else from the dggrids package if
## sea3h_lo + nres - 1 > 6. See load_dggrids() for details
@@ -435,11 +435,11 @@ auto_basis <- function(manifold = plane(),
isea3h <- load_dggrids(res = isea3h_lo + nres - 1) %>%
dplyr::filter(res >= isea3h_lo) # keep only the resolutions we need
}
-
-
+
+
## Find possible grid points for basis locations in each resolution
for(i in 1:nres) {
-
+
## If we are on the plane and we do not want a regular set
if(is(m,"plane") & !regular) {
## Generate mesh using fmesher and use the triangle vertices as notes
@@ -450,10 +450,10 @@ auto_basis <- function(manifold = plane(),
boundary = list(bndary_seg), # boundary
max.edge = max(diff(xrange),diff(yrange))/(2*2.5^(i-1)),
cutoff = max(diff(xrange),diff(yrange))/(3*2.5^(i-1)))$loc[,1:2]
-
+
## Sometimes fmesher returns overlapping points, therefore find the unique set of locations
this_res_locs <- unique(this_res_locs)
-
+
## If we want a regular set of basis functions
} else if(is(m,"plane") & regular) {
## Make a grid that gets more dense as i increases
@@ -472,13 +472,13 @@ auto_basis <- function(manifold = plane(),
this_res_locs <- filter(isea3h,centroid==1,res==(i + isea3h_lo - 1))[c("lon","lat")] %>%
as.matrix() # and convert to matrix
}
-
+
## Setting the scales and Refinement stage
## Set scales: e.g., set to 1.5x the distance to nearest basis if bisquare
## Refine/Prune: Remove basis which are not influenced by data and re-find the scales
for(j in 1:2) { ## Need to go over this twice for refinement when pruning
## Note, prune will be removed in future, and basis functions
-
+
if(nrow(this_res_locs) == 1) { # if we only have one basis function
this_res_scales <-max(diff(xrange),diff(yrange))/2 # set the "base" scale parameter to span most of the range
} else {
@@ -493,7 +493,7 @@ auto_basis <- function(manifold = plane(),
diffx <- abs(xgrid[2] - xgrid[1])
diffy <- abs(ygrid[2] - ygrid[1])
this_res_scales <- rep(min(diffx, diffy), nrow(this_res_locs)) * scale_aperture
- } else {
+ } else {
batch_size <- 5000
batching <- cut(1:nrow(this_res_locs),
breaks = seq(0,nrow(this_res_locs) + batch_size,
@@ -510,9 +510,9 @@ auto_basis <- function(manifold = plane(),
this_res_scales <- unlist(list_scales)
}
}
-
+
}
-
+
## If we have more than one basis at this resolution (could be 0 because of pruning)
if(nrow(this_res_locs) >0)
## The following code can be used to verify that the following basis functions have a similar shape:
@@ -531,7 +531,7 @@ auto_basis <- function(manifold = plane(),
loc=this_res_locs,
scale=ifelse(type=="bisquare",1.5,0.7)*this_res_scales,
type=type, regular = regular)
-
+
## Now refine/prune these basis functions [deprecated]
if(prune > 0 & nrow(this_res_locs)>0) {
## Only in the first step
@@ -541,12 +541,12 @@ auto_basis <- function(manifold = plane(),
## at the data locations and summing over the function values. If the sum is less
## then 'prune' we remove the basis function
rm_idx <- which(colSums(eval_basis(this_res_basis,coords)) < prune)
-
+
## Throw a warning if all basis functions at a given resolution have been removed
if(length(rm_idx) == length(this_res_scales))
warning("prune is too large -- all functions at a resolution removed.
Consider also removing number of resolutions.")
-
+
## If there are basis functions to remove,then remove them and reset scales in j==2
if(length(rm_idx) > 0)
this_res_locs <- this_res_locs[-rm_idx,,drop=FALSE]
@@ -556,11 +556,11 @@ auto_basis <- function(manifold = plane(),
break
}
}
-
+
## Print the number of basis functions at each resolution if verbose
if(verbose)
cat("Number of basis functions at resolution",i,"=",nrow(this_res_locs))
-
+
## Now that we actually have the centroids and scales we can construct the basis functions for this
## resolution. If all the basis functions have been removed at this resolution don't do anything
if(nrow(this_res_locs) > 0) {
@@ -571,83 +571,83 @@ auto_basis <- function(manifold = plane(),
G[[i]]$res=i # put resolution in the basis function data frame
}
}
-
+
## Finally concatenate all the basis functions together using the S4 method concat
G_basis <- Reduce("concat", G)
-
- ## Add a buffer if requested (note that this is only allowed if regular == 1,
- ## which has already been checked at this stage of the function, i.e., if
- ## the user specifies a non-zero buffer with regular == 0, an error is
+
+ ## Add a buffer if requested (note that this is only allowed if regular == 1,
+ ## which has already been checked at this stage of the function, i.e., if
+ ## the user specifies a non-zero buffer with regular == 0, an error is
## thrown upon entry of auto_basis)
- if (buffer != 0) G_basis <- .buffer(G_basis, buffer = buffer, type = type)
-
+ if (buffer != 0) G_basis <- .buffer(G_basis, buffer = buffer, type = type)
+
return(G_basis)
}
## Function to add a buffer of basis functions along the boundary.
-## buffer indicates the percentage increase in the number of basis functions
+## buffer indicates the percentage increase in the number of basis functions
## in each dimension.
.buffer <- function(basis, buffer, type) {
-
+
L <- length(unique(basis@df$res)) # number of resolutions of basis functions
dimens <- basis@manifold@measure@dim # dimension of the manifold
dim_names <- c(paste0("loc",1:(dimens))) # names of each dimension
basis_locs_list <- list()
scale <- res <- c()
-
+
for (k in 1:L) { # for each resolution
-
+
## Extract basis function dataframe for the kth resolution
basis_df <- basis@df[basis@df$res == k, ]
-
+
## unique locations for all dimensions
## NB: need to consider rectangular domains too
- unique_locs <- lapply(basis_df[, dim_names, drop = FALSE], function(x) sort(unique(x)))
-
- ## Compute the distance between locations
+ unique_locs <- lapply(basis_df[, dim_names, drop = FALSE], function(x) sort(unique(x)))
+
+ ## Compute the distance between locations
## (assumes equally spaced basis functions in each dimension, i.e., regular = TRUE)
dist_btwn_locs <- sapply(unique_locs, function(x) diff(x)[1])
-
+
## Compute number of extra terms to add in each direction of each dimension
- ## (divide buffer by 2 because we add locations in both directions of a
+ ## (divide buffer by 2 because we add locations in both directions of a
## given dimension)
n_add <- sapply(unique_locs, function(x) ceiling(buffer/2 * length(x)))
-
- ## Compute the new unique locations of basis functions for each dimension.
- ## First, define a function to add n_add elements separated by a units to
+
+ ## Compute the new unique locations of basis functions for each dimension.
+ ## First, define a function to add n_add elements separated by a units to
## the head and tail of a vector x.
.add_head_tail <- function(x, n_add, a) {
- c(min(x) - (n_add:1) * a,
- x,
+ c(min(x) - (n_add:1) * a,
+ x,
max(x) + (1:n_add) * a)
}
-
+
## Now generate new unique locations in each dimension
tmp <- names(unique_locs)
names(tmp) <- names(unique_locs)
- unique_locs_new <- lapply(tmp,
- function(id, x, n, a) .add_head_tail(x[[id]], n[id], a[id]),
+ unique_locs_new <- lapply(tmp,
+ function(id, x, n, a) .add_head_tail(x[[id]], n[id], a[id]),
x = unique_locs, n = n_add, a = dist_btwn_locs)
-
- ## Now create new basis function locations, and append to previously
- ## computed basis function locations. Conveniently, expand.grid() allows
+
+ ## Now create new basis function locations, and append to previously
+ ## computed basis function locations. Conveniently, expand.grid() allows
## lists to be passed as arguments.
basis_locs_list[[k]] <- as.matrix(expand.grid(unique_locs_new))
-
+
scale <- c(scale, rep(basis_df$scale[1], nrow(basis_locs_list[[k]])))
-
+
res <- c(res, rep(k, nrow(basis_locs_list[[k]])))
}
-
+
## Convert to matrix
basis_locs <- do.call(rbind, basis_locs_list)
-
- ## Convert to a basis
- basis_buffer <- local_basis(manifold = basis@manifold,
- loc = basis_locs,
- scale = scale,
- type = type,
- res = res,
+
+ ## Convert to a basis
+ basis_buffer <- local_basis(manifold = basis@manifold,
+ loc = basis_locs,
+ scale = scale,
+ type = type,
+ res = res,
regular = basis@regular)
return(basis_buffer)
}
@@ -662,18 +662,18 @@ setMethod("TensorP",signature(Basis1="Basis",Basis2="Basis"),function(Basis1,Bas
## (This is typically the case, as time usually has one resolution of basis functions in FRK)
if(nres(Basis2) > 1)
stop("Only one basis can be multiresolution (Basis 1)")
-
+
## Extract data frames and dimensions
df1 <- data.frame(Basis1)
df2 <- data.frame(Basis2)
n1 <- dimensions(Basis1)
n2 <- dimensions(Basis2)
-
+
## Create long data frame with all possible centroid combinations
## (Kronecker of space and time)
expand.grid(df1[,1:n1,drop=FALSE],
df2[,1:n2,drop=FALSE])
-
+
## We adopt a space-first approach, where the spatial index changes quicker than the
## temporal index. So we repeat the centroids of Basis1 for as many time points as
## we have, and we repeat the time points so they look like 1111111,222222,33333 etc.
@@ -682,24 +682,24 @@ setMethod("TensorP",signature(Basis1="Basis",Basis2="Basis"),function(Basis1,Bas
df <- cbind(df1[rep(1:nrow(df1),times = nrow(df2)),1:n1,drop=FALSE],
df2[rep(1:nrow(df2),each = nrow(df1)),1:n2,drop=FALSE],
df1[rep(1:nrow(df1),times = nrow(df2)),"res",drop=FALSE])
-
+
## Change the names of the data frame to what is standard in this package
names(df) <- c(paste0("loc",1:(n1 + n2)),"res")
-
+
## Check number of basis functions and throw warning if necessary
if((nbasis_tot <- nbasis(Basis1) * nbasis(Basis2)) > 5000)
warning("Tensor product has resulted in more than 5000 functions.
- If you intend to use the EM algorithm during model fitting (i.e., method = 'EM'),
+ If you intend to use the EM algorithm during model fitting (i.e., method = 'EM'),
please reduce the number of spatial or temporal basis functions.")
-
+
regular <- Basis1@regular && Basis2@regular
-
+
## Create new Tensor Basis function from this information
new("TensorP_Basis",
Basis1 = Basis1,
Basis2 = Basis2,
n = nbasis_tot,
- df = df,
+ df = df,
regular = regular)
})
@@ -708,22 +708,22 @@ setMethod("TensorP",signature(Basis1="Basis",Basis2="Basis"),function(Basis1,Bas
#' @aliases eval_basis,Basis-matrix-method
setMethod("eval_basis",signature(basis="Basis",s="matrix"),
function(basis,s){
-
+
space_dim <- dimensions(basis) # spatial dimensions
n <- nrow(s) # number of points over which to evaluate basis functions
-
+
batching=cut(1:n, # create batches of computation (in batches of 10000)
breaks = seq(0, # break up into batches ot 10000
n+10000,
by=10000),
labels=F) # do not assign labels to batches
-
+
if(opts_FRK$get("parallel") > 1L) { # if we have a parallel backend when compute in parallel
-
+
clusterExport(opts_FRK$get("cl"), # export variables to the cluster
c("batching","basis","s","space_dim"),
envir=environment())
-
+
## Use parLapply to compute batches in parallel. The drop = FALSE in the end is required
## for when we have one spatial dimension
pnt_eval_list <- parLapply(opts_FRK$get("cl"),1:max(unique(batching)),
@@ -742,7 +742,7 @@ setMethod("eval_basis",signature(basis="Basis",s="matrix"),
s[idx,1:space_dim,drop = FALSE]))
})
}
-
+
## Finally concatenate all the bits together using rbind
do.call(rbind,pnt_eval_list)
})
@@ -751,7 +751,7 @@ setMethod("eval_basis",signature(basis="Basis",s="matrix"),
#' @aliases eval_basis,Basis-SpatialPointsDataFrame-method
setMethod("eval_basis",signature(basis="Basis",s="SpatialPointsDataFrame"),
function(basis,s){
-
+
## Now just evaluate the basis functions at the coordinates of the SpatialPoints
eval_basis(basis=basis,
s = coordinates(s))
@@ -762,25 +762,25 @@ setMethod("eval_basis",signature(basis="Basis",s="SpatialPointsDataFrame"),
#' @aliases eval_basis,Basis-SpatialPolygonsDataFrame-method
setMethod("eval_basis",signature(basis="Basis",s="SpatialPolygonsDataFrame"),
function(basis,s){
-
+
## Inform user this might take a while
cat("Averaging over polygons...\n")
-
-
+
+
## If we have a parallel backend use parLapply
if(opts_FRK$get("parallel") > 1L) {
-
+
## Export variavles we need to the cluster
clusterExport(opts_FRK$get("cl"),
c("basis","s"),envir=environment())
-
+
## Compute averaging over footprints in parallel
X <- parLapply(opts_FRK$get("cl"),1:length(s), function(i) {
samps <- .samps_in_polygon(basis,s,i) # sample 1000 times (fixed) uniformly in polygon
colSums(.point_eval_fn(basis@fn,samps))/nrow(samps) # This is the averaging
})
clusterEvalQ(opts_FRK$get("cl"), {gc()}) # clear the cluster memory
-
+
} else {
## Same as above but serially
X <- lapply(1:length(s), function(i) {
@@ -788,10 +788,10 @@ setMethod("eval_basis",signature(basis="Basis",s="SpatialPolygonsDataFrame"),
colSums(.point_eval_fn(basis@fn,samps))/nrow(samps)
})
}
-
+
X <- Reduce("rbind",X) # join the rows together
as(X,"Matrix") # coerce to Matrix if not already Matrix
-
+
})
#' @rdname eval_basis
@@ -814,7 +814,7 @@ setMethod("eval_basis",signature(basis="TensorP_Basis",s="matrix"),
n1 <- dimensions(basis@Basis1) # spatial dimensions
S1 <- eval_basis(basis@Basis1,s[,1:n1,drop=FALSE]) # evaluate at spatial locations
S2 <- eval_basis(basis@Basis2,s[,-(1:n1),drop=FALSE]) # evaluate at temporal locations
-
+
## Order: Construct S matrix over space and time
## This is ordered as first space then time, therefore we take S2[,1]*S1 as our first block
## Then S2[,2]*S1 as our second block etc.
@@ -845,11 +845,11 @@ setMethod("eval_basis",signature(basis="TensorP_Basis",s = "STFDF"),function(bas
first data point.")
tlocs <- matrix(s$t) # all time points
nt <- length(time(s)) # number of unique time points
-
+
S1 <- eval_basis(basis@Basis1,s[,1]) # evaluate over space (just take first time point)
S1 <- do.call("rbind",lapply(1:nt,function(x) S1)) # now just repeat that for the nt time points
S2 <- eval_basis(basis@Basis2,tlocs[,,drop=FALSE]) # evaluate over time (all time points)
-
+
## As in previous functions we compute S with space running fastest
XX <- lapply(1:ncol(S2),function(i) (S2[,i] * S1))
S <- quickcbind(XX) # a quick cbind method using sparse matrix construction
@@ -867,12 +867,12 @@ setMethod("remove_basis",signature(Basis="Basis"),function(Basis,rmidx) {
if(!all(rmidx %in% 1:ntot))
stop("Please make sure indices are numeric and within
1 and the number of basis functions.")
-
+
if (length(rmidx) == 0) {
cat("length(rmidx) == 0; not removing any basis functions\n")
return(Basis)
}
-
+
Basis_df <- data.frame(Basis) # extract data frame
Basis@fn <- Basis@fn[-rmidx] # remove functions
Basis@pars <- Basis@pars[-rmidx] # remove parameters
@@ -885,22 +885,22 @@ setMethod("remove_basis",signature(Basis="Basis"),function(Basis,rmidx) {
#' @rdname combine_basis
#' @aliases combine_basis,Basis-method
setMethod("combine_basis",signature(Basis_list="list"),function(Basis_list) {
-
+
if (!is.list(Basis_list))
stop("Basis_list should be a list of objects of class Basis")
-
+
if(!all(sapply(Basis_list, function(x) class(x) == "Basis")))
stop("Basis_list should be a list of objects of class Basis")
-
+
## Check that only one manifold is present:
if (sapply(Basis_list, manifold) %>% unique %>% length != 1) {
- stop("Multiple manifolds detected: Please only combine basis functions that
+ stop("Multiple manifolds detected: Please only combine basis functions that
are defined on the same manifold.")
} else {
manifold <- Basis_list[[1]]@manifold
}
-
- ## Find the number of unique resolutions in each basis object:
+
+ ## Find the number of unique resolutions in each basis object:
num_resolutions_in_each_basis <- length(unique(sapply(Basis_list, function(x) length(unique(x@df$res)))))
if (num_resolutions_in_each_basis != 1) {
stop("Currently only allow each element of Basis_list to contain basis functions at one resolution")
@@ -912,40 +912,40 @@ setMethod("combine_basis",signature(Basis_list="list"),function(Basis_list) {
df <- rbind(df, Basis_list[[i]]@df)
}
}
-
+
n <- sum(sapply(Basis_list, nbasis))
-
+
## Easiest just to combine the list slots using a for loop
fn <- pars <- list()
for (i in 1:length(Basis_list)) {
fn <- c(fn, Basis_list[[i]]@fn)
pars <- c(pars, Basis_list[[i]]@pars)
}
-
+
regular <- sapply(Basis_list, function(x) x@regular) %>% as.logical %>% all
-
-
- Basis(manifold = manifold,
- n = n,
- fn = fn,
- pars = pars,
- df = df,
+
+
+ Basis(manifold = manifold,
+ n = n,
+ fn = fn,
+ pars = pars,
+ df = df,
regular = regular)
})
#' @rdname remove_basis
#' @aliases remove_basis,Basis-method
setMethod("remove_basis", signature(Basis="Basis", rmidx = "SpatialPolygons"), function(Basis,rmidx) {
-
+
if (!is(Basis@manifold, "plane")) stop("remove_basis with SpatialPolygons only works on the plane")
-
- ## Using remove_basis()
+
+ ## Using remove_basis()
sp_df <- Basis@df
coordinates(sp_df) <- ~ loc1 + loc2
spatial_object <- rmidx
rmidx <- which(is.na(over(sp_df, spatial_object)))
Basis <- remove_basis(Basis, rmidx)
-
+
return(Basis)
})
@@ -1058,12 +1058,12 @@ setMethod("count_res",signature="TensorP_Basis", # Returns count by resolution f
res <- NULL # suppress bindings (it's in the data frame)
c1 <- count(data.frame(.Object@Basis1),res) # count spatial by resolution
c2 <- count(data.frame(.Object@Basis2),res) # count temporal by resolution
-
+
## In the below we define resolutions with the spatial resolution moving fastest
## So for example if we have resolutions 1,2, and 3 for space and 1, and 2 for time,
## then we will have resolutions 1,2,3,4,5,6, where 1,2,3 correspond to space 1,2,3 and
## time 1, and 4,5,6 correspond to space 1,2,3 and time 2
-
+
c_all <- NULL # initialise NULL
max_res_c1 <- c1$res[1] - 1 # initialise at one less than coarsest resolution
for( i in 1:nrow(c2)) { # for each temporal resolution (should be only one)
@@ -1133,7 +1133,7 @@ setMethod("BuildD",signature(G="Basis"),function(G){
nres <- nres(G) # number of resoluations
m <- manifold(G) # basis manifold
n <- dimensions(G) # manifld dimensions
-
+
## Since the first columns of df are ALWAYS the centroid coordinates
## we just take the first n columns when computing the distances
D_basis = lapply(1:nres,function(i) {
@@ -1168,31 +1168,31 @@ setMethod("BuildD",signature(G="TensorP_Basis"),function(G){
## else sample on the sphere by first drawing a box around the polygon (on the sphere)
## sampling in the box, and then keeping only those samples inside the polygon
} else if(is(basis@manifold,"sphere")){
-
+
## Find coordinates
coords <- data.frame(slot(s@polygons[[i]]@Polygons[[1]],"coords"))
-
+
## Find lat/lon range
rangelon <- range(coords$lon) + 180
rangelat <- range(coords$lat) + 90
-
+
## Find limits in transformed space
rangeu <- rangelon/360
rangev <- (cos(rangelat*2*pi/360) + 1)/2
-
+
## Sample in transformed space
samps <- cbind(runif(nMC,min=min(rangeu), max=max(rangeu)),
runif(nMC,min=min(rangev), max=max(rangev)))
-
+
## Re-transform back to spherical coordinates
samps[,1] <- 360*samps[,1] - 180
samps[,2] <- acos(2*samps[,2] -1) * 360 / (2*pi) - 90
-
+
## Find which points are in polygon (usually approx. 70%)
pip <- over(SpatialPoints(samps), # pip = points in polygon
SpatialPolygons(list(s@polygons[[i]]),1L))
samps <- samps[which(pip == 1),] # keep those samples for which pip == TRUE
-
+
}
samps
}
@@ -1259,7 +1259,7 @@ setMethod("concat",signature = "Basis",function(...) {
if(!(length(unique(sapply(sapply(l,manifold),type))) == 1))
stop("Basis need to be on the same manifold")
G <- l[[1]] # initialise first basis
-
+
## We are going to use internals (slots) since this is not exported
## This is more direct and safe in this case
for (i in 2:length(l)) { # add on the other basis functions