halton.indices: Halton indices

Description

Compute and attach "inverse" or indices of the Halton sequence to points. Points can be an arbitrary set or a Halton lattice.

Usage

halton.indices(
  x,
  J = NULL,
  hl.bbox,
  bases = c(2, 3),
  index.name = "HaltonIndex",
  use.CRT = FALSE
)

Arguments

Either a data frame or a SpatialPoints* object. Suitable input objects are the output of functions halton.lattice (a data frame) and halton.lattice.polygon (a SpatialPointsDataFrame object).

If x is a data frame, it must either contain the names of coordinates columns as attribute "coordnames", or coordinates must be the first D columns of the data frame. I.e., coordinates are either x[,attr(x,"coordnames")] or x[,1: length(bases)].

This function works for dimensions >2 if x is a data.frame. SpatialPoints* objects are not defined for D>2.

A vector of length D containing base powers. J determines the size and shape of the smallest Halton boxes in D space. There are bases[i]^J[i] boxes over the i-th dimension of x's bounding box. Total number Halton boxes is prod(bases^J). The size of each box in the i-th dimension is delta[i]/ (bases[i]^J[i]), where delta[i] is the extent of x's bounding box along the i-th dimension. If J is NULL (the default), approximately length(x) boxes will be chosen (approx. one point per box) and boxes will be as square as possible.

hl.bbox

DX2 matrix containing bounding box of the full set of Halton boxes. First column of this matrix is the lower-left coordinate (i.e., minimums) of the bounding box. Second column is the upper-right coordinate (i.e., maximums) of the bounding box. For example, if D = 2, hl.bbox = matrix( c(min(x), min(y), max(x), max(y)),2). If hl.bbox is missing (the default), the bounding box of x is used, but expanded on the top and right by 1 percent to include any points exactly on the top and right boundaries. If hl.bbox is supplied, keep in mind that all point outside the box, or on the maximums (i.e., hl.bbox[,2]), will not be assigned Halton indices.

bases

A vector of length D containing Halton bases. These must be co-prime.

index.name

A character string giving the name of the column in the output data frame or SpatialPoints object to contain the Halton indices. This name is saved as an attribute attached to the output object.

use.CRT

A logical values specifying whether to invert the Halton sequence using the Chinese Remainder Theorem (CRT). The other method (use.CRT == FALSE) is a direct method, and is very fast, but requires multiple huge vectors be allocated (size of vectors is prod{bases^J}, see Details). As the number of points grows, eventually the direct method will not be able to allocate sufficient memory (tips: Make sure to run 64-bit R, and try increasing memory limit with memory.limit). The CRT method, while much (much) slower, does not require as much memory, and should eventually complete a much larger problem. Patience is required if your problem is big enough to require use.CRT == TRUE.

Value

If x is a data frame, x is returned with an addition column. The additional column is named index.name and stores the index of the Halton box containing the point represented on that line of x. If x is a SpatialPoints* object, a SpatialPointsDataFrame is returned containing the points in x. The attributes of the returned object have an additional column, the index of the Halton box containing the point. Name of the attribute is index.name. If multiple points fall in the same Halton box, their Halton indices are identical.

Details

Halton indices are the arguments to the Halton sequence. This routine is the inverse function for the Halton sequence. Given a point in D space, this routine computes the index (a non-negative integer) of the Halton sequence which maps to the Halton region containing the point.

For example, in 1D, with bases == 2, J == 3, and hl.bbox= matrix(c(0,1),1), all points in the interval [0,1/8) have Halton index equal to 0, all point in [1/8,2/8) have Halton index 4, points in [2/8,3/8) have index 2, etc. To check, note that the Halton sequence maps x (mod 8) = 4 to the interval [1/8,2/8), x (mod 8) = 2 are mapped to [2/8,3/8), etc. (i.e., check range(halton(200)[ (0:199)%% 8 == 4]) and range(halton(200)[ (0:199)%% 8 == 2]) )

Examples

Run this code

# NOT RUN {
# The following is equivalent to hip.point(WA.cities,25,J=c(3,2))
#
# Add tiny amount to right and top of bounding box because Halton boxes are 
# closed on the left and bottom.  This includes points exactly on the bounding lines.
bb <- bbox(WA.cities) + c(0,0,1,1) 

# Compute Halton indices
frame <- halton.indices( WA.cities, J=c(3,2), hl.bbox=bb  )

# Construct Halton frame
frame <- halton.frame( frame )

# Draw HAL sample
n <- 25
N.frame <- nrow(frame)
m <- floor(runif(1, 0, N.frame)) # Integer 0,...,N.frame-1
ind <- (((0:(n-1))+m) %% N.frame ) + 1  # Cycle around frame if necessary
samp <- frame[ind,]  # draw sample


# }