# R/tbmain.R In tbart: Teitz and Bart's p-Median Algorithm

#### Documented in allocateallocationseuc.distsmink.distsstar.diagramtb

```.complement <- function(ivec, imax) {
result <- rep(TRUE,imax)
result[ivec] <- FALSE
return(which(result))
}

# Implimentation of TB algorithm without 'forced' supply points

.tb <- function(d,guess,verbose=FALSE) {
config <- guess
n <- ncol(d)
repeat {
old.config <- config
config <- .bestswap(d,config,.complement(config,n))
if (verbose) {
cat("Configuration: ",config)
score <- .dtotal(d,config)
cat("  Score:",score,"\n")
}
if (all(old.config==config)) break
}
return(config)
}

# Implimentation of TB algorithm with 'forced' supply points

.tb2 <- function(d,guess,n_force, verbose=FALSE) {
config <- guess
n <- ncol(d)
repeat {
old.config <- config
if (missing(n_force)) {
config <- .bestswap(d,config,.complement(config,n))
} else {
config <- .bestswap2(d,config,.complement(config,n),n_force)
}
if (verbose) {
cat("Configuration: ",config)
score <- .dtotal(d,config)
cat("  Score:",score,"\n")
}
if (all(old.config==config)) break
}
return(config)
}

##' Teitz-Bart algorithm applied to a 'raw' distance matrix
##'
##' @param  d - A distance matrix (not necessarily Euclidean)
##' @param  guess - a guess at the set of p points constituting the \eqn{p}-median
##' @param  verbose - if TRUE print out each swap in the algorithm (default is FALSE)
##' @return Set of point indices for \eqn{p}-median (may be local optimum)
##'
##' @examples
##' x1 <- rnorm(100)
##' y1 <- rnorm(100)
##' d <- as.matrix(dist(cbind(x1,y1)))
##' tb.raw(d,c(1,2))
##' @export
##'

tb.raw <- .tb

##' Euclidean distances from a Spatial* or Spatial*DataFrame object
##'
##' @param  swdf1 - First Spatial*DataFrame object
##' @param  swdf2 - Second Spatial*DataFrame object (if omitted,  defaults to the same value as \code{swdf1})
##' @param  scale - allows re-scaling eg: value of 1000 means distances in km if coordinates of \code{swdf1}/\code{swdf2} in meters.
##' @return Distance matrix (if \code{swdf1} or \code{swdf2} not SpatialPoints*, distances are based on points obtained from \code{coordinates} function)
##'
##' @examples
##' data(meuse)
##' coordinates(meuse) <- ~x+y
##' euc.dists(meuse,scale=1000)
##'
##' @export
##'

euc.dists <- function(swdf1,swdf2,scale) {
if (missing(swdf2)) swdf2 <- swdf1
xy1 <- coordinates(swdf1)
xy2 <- coordinates(swdf2)
if (!missing(scale)) {
xy1 <- xy1/scale
xy2 <- xy2/scale
}
return(.dmat(xy1[,1],xy2[,1],xy1[,2],xy2[,2]))
}

##' Minkowski distances from a Spatial* or Spatial*DataFrame object
##'
##' @param  swdf1 - First Spatial*DataFrame object
##' @param  swdf2 - Second Spatial*DataFrame object (if omitted,  defaults to the same value as \code{swdf1})
##' @param  pwr - Minkowski exponent
##' @param  scale - allows re-scaling eg: value of 1000 means distances in km if coordinates of \code{swdf1}/\code{swdf2} in meters.
##' @param  weight - weight for each element in swdf1 (the demand locations)
##' @return Distance matrix (if \code{swdf1} or \code{swdf2} not SpatialPoints*,
##'   distances are based on points obtained from \code{coordinates} function)
##'
##' @examples
##' data(meuse)
##' coordinates(meuse) <- ~x+y
##' d1 <- mink.dists(meuse,pwr=1,scale=1000)   # Taxicab metric
##' d2 <- mink.dists(meuse,pwr=Inf,scale=1000) # Works for limiting case
##'
##' @export
##'

mink.dists <- function(swdf1,swdf2,pwr,scale,weight) {
if (missing(swdf2)) swdf2 <- swdf1
xy1 <- coordinates(swdf1)
xy2 <- coordinates(swdf2)
if (!missing(scale)) {
xy1 <- xy1/scale
xy2 <- xy2/scale
}
temp <- .dmatex(xy1[,1],xy2[,1],xy1[,2],xy2[,2],pwr)
if (missing(weight)) return(temp)
sweep(temp, 1, weight, "/")
}

##' Teitz-Bart algorithm applied to Spatial* and Spatial*DataFrame objects
##'
##' This reports the \eqn{p}-median set
##'
##' @param  swdf1 - first Spatial* or Spatial*DataFrame objects - the 'demand' set
##' @param  swdf2 - second Spatial* or Spatial*DataFrame objects - the 'supply' set (if omitted,  defaults to the same value as \code{swdf1})
##' @param  p - either a guess at the initial \eqn{p}-median set of a single integer indicating the size of the set (which is then chosen randomly)
##' @param  metric - the distance matrix (defaults to Euclidean computed via \code{euc.dists(swdf1,swdf2)} if not supplied)
##' @param  verbose - if TRUE print out each swap in the algorithm (default is FALSE)
##' @return Set of point indices for \eqn{p}-median (may be local optimum)
##'
##' @examples
##' data(meuse)
##' coordinates(meuse) <- ~x+y
##' tb(meuse,p=5)
##' @export
##'

tb <- function(swdf1,swdf2,p,metric,verbose=FALSE) {
if (missing(swdf2))  swdf2 <- swdf1
n.choices <- nrow(coordinates(swdf2))
if (length(p) == 1) p <- sample(n.choices,p)
if (missing(metric)) metric <- euc.dists(swdf1,swdf2)
result <- .tb(metric,p,verbose)
return(result)
}

##' Creates the lines for a 'star diagram'
##'
##' @param  swdf1 - first Spatial* or Spatial*DataFrame objects
##' @param  swdf2 - second Spatial* or Spatial*DataFrame objects (if omitted,  defaults to the same value as \code{swdf1})
##' @param  alloc - a list saying which coordinate in swdf2 is allocated to each point in swdf1 (if ommitted, looks for \code{allocation} column in \code{swdf1})
##'
##' @examples
##' data(meuse)
##' coordinates(meuse) <- ~x+y
##' allocations.list <- allocate(meuse,p=5)
##' star.lines <- star.diagram(meuse,alloc=allocations.list)
##' plot(star.lines)
##'
##' # Acquire allocations from swdf1
##' require(GISTools)
##' set.seed(461976) # Reproducibility
##' data(georgia)
##' georgia3 <- allocations(georgia2,p=8)
##' plot(georgia3,border='grey')
##'
##'
##' @export
##'

star.diagram <- function(swdf1,swdf2,alloc) {
if (missing(swdf2))  swdf2 <- swdf1
if (missing(alloc)) {
alloc <- swdf1\$allocation
}
co1 <- coordinates(swdf1)
co2 <- coordinates(swdf2)
result <- vector(nrow(co1),mode='list')
for (i in 1:nrow(co1))  result[[i]] <- Lines(list(Line(cbind(c(co1[i,1],co2[alloc[i],1]),c(co1[i,2],co2[alloc[i],2])))),ID=sprintf("Star%d",i))
sl <- SpatialLines(result)
sldf <- SpatialLinesDataFrame(sl,data.frame(allocate=alloc),match.ID=FALSE)
return(sldf)}

##' Teitz-Bart algorithm applied to Spatial* and Spatial*DataFrame objects
##'
##' Return demand Spatial*Dataframe with new columns giving allocation id and distance to supply point
##'
##' @param  swdf1 - first Spatial* or Spatial*DataFrame objects
##' @param  swdf2 - second Spatial* or Spatial*DataFrame objects (if omitted,  defaults to the same value as \code{swdf1})
##' @param  force - list of supply points or logical vector with length the same as the number of supply points that are forced to be used - eg e
##' @param  p - either a guess at the initial \eqn{p}-median set of a single integer indicating the size of the set (which is then chosen randomly)
##' @param  metric - the distance matrix (defaults to Euclidean computed via \code{euc.dists(swdf1,swdf2)} if not supplied)
##' @param  verbose - if TRUE print out each swap in the algorithm (default is FALSE)
##' @return Copy of swdf1 with extra data columns called \code{allocation} and \code{allocdist}
##' with indices for each element from the \eqn{p}-median set
##'
##' @examples
##'
##' require(RColorBrewer)
##' require(GISTools)
##' data(georgia)
##' georgia3 <- allocations(georgia2,p=5,force=c(1,120,44))
##' col.index <- match(georgia3\$allocation,unique(georgia3\$allocation))
##' col.alloc <- brewer.pal(5,'Accent')[col.index]
##' par(mfrow=c(1,2))
##' plot(georgia3,col=col.alloc)
##' choropleth(georgia3,georgia3\$allocdist)
##'
##'
##' # Use in conjunction with rgeos
##' require(rgeos)
##' require(GISTools)
##' georgia3 <- allocations(georgia2,p=5,force=c(1,120,44))
##' georgia4 <- gUnaryUnion(georgia3,georgia3\$allocation)
##' plot(georgia4)
##'
##' @export
##'

allocations <-function(swdf1,swdf2,force,p,metric,verbose=FALSE) {
if (missing(swdf2))  swdf2 <- swdf1
n.choices <- nrow(coordinates(swdf2))
if (missing(metric)) metric <- euc.dists(swdf1,swdf2)
if (missing(force)) {
nni <- allocate(swdf1,swdf2,p=p,metric=metric,verbose=verbose)
} else {
nni <- allocate(swdf1,swdf2,force=force,p=p,metric=metric,verbose=verbose)
}
n.demand <- nrow(coordinates(swdf1))
picker <- cbind(1:n.demand,nni)
dist <- metric[picker]
swdf1_df <- data.frame(swdf1)
swdf1_df <- cbind(swdf1_df,data.frame(allocation=nni,allocdist=dist))
if (class(geometry(swdf1)) == "SpatialPolygons") return(SpatialPolygonsDataFrame(swdf1,swdf1_df))
if (class(geometry(swdf1)) == "SpatialPoints") return(SpatialPointsDataFrame(swdf1,swdf1_df))
return(SpatialLinesDataFrame(swdf1,swdf1_df))
}

##' Teitz-Bart algorithm applied to Spatial* and Spatial*DataFrame objects
##'
##' This function returns the allocations for each demand point - in terms of the index number of
##' the record in \code{swdf2} assigned as the supply point.  This version is useful as part of
##' code inside other functions
##'
##' @param  swdf1 - first Spatial* or Spatial*DataFrame objects
##' @param  swdf2 - second Spatial* or Spatial*DataFrame objects (if omitted,  defaults to the same value as \code{swdf1})
##' @param  force - list of supply points or logical vector with length the same as the number of supply points that are forced to be used - eg existing outlets
##' @param  p - either a guess at the initial \eqn{p}-median set of a single integer indicating the size of the set (which is then chosen randomly)
##' @param  metric - the distance matrix (defaults to Euclidean computed via \code{euc.dists(swdf1,swdf2)} if not supplied)
##' @param  verbose - if TRUE print out each swap in the algorithm (default is FALSE)
##' @return List of nearest neigbour indices for each element from the \eqn{p}-median set
##'
##' @examples
##' data(meuse)
##' coordinates(meuse) <- ~x+y
##' allocate(meuse,p=5)
##'
##'
##'
##' require(RColorBrewer)
##' require(GISTools)
##' data(georgia)
##' allocations.list <- allocate(georgia2,p=5)
##' zones <- gUnaryUnion(georgia2,allocations.list)
##' plot(zones,col=brewer.pal(5,"Accent"))
##' points(coordinates(georgia2)[allocations.list,],pch=16,cex=2,col=rgb(1,0.5,0.5,0.1))
##'
##' @export
##'

allocate <-function(swdf1,swdf2,force,p,metric,verbose=FALSE) {
if (missing(swdf2))  swdf2 <- swdf1
n.choices <- nrow(coordinates(swdf2))
if (length(p) == 1) {
n.subset <- p
p <- sample(n.choices,p)
} else {
n.subset <- length(p)
}
if (missing(metric)) metric <- euc.dists(swdf1,swdf2)
if (missing(force)) {
indices <- .tb(metric,p,verbose)
} else {
if (is.logical(force)) force <- which(force)
p <- c(force,setdiff(p,force))[1:n.subset]
indices <- .tb2(metric,p,as.integer(length(force)), verbose)
}
nni <- .rviss(metric,indices)
return(nni)
}
```

## Try the tbart package in your browser

Any scripts or data that you put into this service are public.

tbart documentation built on May 2, 2019, 2:49 p.m.