halton.indices.vector: Halton indices for an entire vector of coordinates

Description

Computes Halton indices of an entire vector of points by matching them with a vector of the Halton sequence. This function is relatively fast, but can only handle reasonably sized vectors.

Usage

halton.indices.vector(
  hl.coords,
  n.boxes,
  D = 2,
  b = c(2, 3),
  delta = c(1, 1),
  ll.corner = c(0, 0)
)

Arguments

hl.coords

nXD vector of coordinates for points

n.boxes

DX1 vector containing number of Halton boxes in each dimension.

Number of dimensions

DX1 vector of bases to use for each dimension

delta

DX1 vector of study area extents in each dimension. Study area is delta[i] units wide in dimension i.

ll.corner

DX1 vector containing minimum coordinates in all dimensions.

Value

A nX1 vector of Halton indices corresponding to points in hl.coords.

Details

The Halton sequence maps the non-negative integers (the Halton indices) to D-space. This routine does the inverse. Given a point in D-space and a grid of Halton boxes, the point's Halton index is any integer N which gets mapped to the Halton box containing the point. (i.e., any integer in the set ${x:N = x mod C}$, where $C$ = prod(n.boxes)).

This routine uses the Halton sequence and modular arithmetic to find Halton indices. This means several vectors of size nrow(hl.coords) must be created. Depending on memory, this approach fails for a sufficiently large number of points. When this routine fails, see the slower halton.indices.CRT, which computes indices by solving the Chinese Remainder Theorem.

Examples

Run this code

# NOT RUN {
# Compute Halton box index for one value
pt <- data.frame(x=0.43, y=0.64)
n.boxes <- c(16,9) 
halton.indices.vector(pt, n.boxes) # should equal 70

# The following should also equal 70
pt <- data.frame(x=143, y=164)
halton.indices.vector(pt, n.boxes, delta=c(100,100), ll.corner=c(100,100)) 

# Plot Halton boxes and indices to check
b <- c(2,3)
J <- c(4,2)  # or, J <- log(n.boxes) / log(b) # = (log base 2 of 16, log base 3 of 9)
hl.ind <- halton( prod(n.boxes), 2,0 )
plot(c(0,1),c(0,1),type="n")
for( i in J[1]:1) abline(v=(0:b[1]^i)/b[1]^i, lwd=J[1]+1-i, col=i)
for( i in J[2]:1) abline(h=(0:b[2]^i)/b[2]^i, lwd=J[2]+1-i, col=i)
for( i in 1:prod(n.boxes)){
  box.center <- (floor(n.boxes*hl.ind[i,]+.Machine$double.eps*10) + 1-.5)/n.boxes
  text(box.center[1],box.center[2], i-1, adj=.5)  
}
points(pt$x, pt$y, col=6, pch=16, cex=2)


# }