halton.indices.CRT: Halton indices by the Chinese Remainder Theorem (CRT)

Description

Computes Halton indices of D-dimensional points by solving the Chinese Remainder Theorem. This function is slightly slower than halton.indices.vector, but it works for large problems.

Usage

halton.indices.CRT(
  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. No points can be outside the bounding box or exactly on the right or top boundary. See Details.

n.boxes

DX1 vector containing number of Halton boxes in each dimension.

Number of dimensions

DX1 vector of bases to use in the Halton sequence.

delta

DX1 vector of study area bounding box extents in each dimension. Study area is bounded by a cube in D space, which is delta[i] units wide in dimension i. Area of bounding cube is prod{delta} units to the D power.

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 solves the Chinese Remainder Theorem to find Halton indices. This routine loops over the points in hl.coords, and as such minimizes memory usage but sacrifices speed. For small problems, see halton.indices.vector, which computes indices by actually placing points in Halton boxes to find their indices.

No point can be less than it's corresponding ll.corner. No point can be equal to or greater than it's corresponding ll.corner + delta.

Note: n.boxes is checked for compatibility with b. That is, log(n.boxes, b) must all be integers.

Examples

Run this code

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

# Plot Halton boxes and indices to check.  
# pt should plot in box labeled 70
b <- c(2,3)
J <- log(n.boxes,b)  # J must be integers
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)

# Longer vector
tmp <- data.frame(x=(0:100)/101,y=.2)
n.boxes <- c(16,9)
tmp.crt <- halton.indices.CRT(tmp, n.boxes)

# }