R/redistMPI.R
In kosukeimai/redist: Simulation Methods for Legislative Redistricting
Defines functions
redist.mcmc.mpi
ecutsAppend
redist.combine.mpi
ecutsMPI
Documented in
redist.combine.mpi
redist.mcmc.mpi
##############################################
# Authors: Ben Fifield & Alex Tarr #
# Created: 2015/06/01 #
# Last Revision: N/A #
# Institution: Princeton University #
# Purpose: R wrapper to run swMH() code w/ #
# parallel tempering (mpi) #
##############################################
ecutsMPI <- function(procID = procID, params = params, adj = adj, total_pop = total_pop, init_plan = init_plan, swaps = swaps) {
## Load redist library
library(redist)
if (is.na(params$savename)) {
fname <- paste0("log", procID)
} else {
fname <- paste0("log", procID, "_", params$savename)
}
if (params$verbose) {
sink(fname)
## Initialize ##
divider <- c(paste(rep("=", 20), sep = "", collapse = ""), "\n")
cat("\n", append = TRUE)
cat(divider, append = TRUE)
cat("redist.mcmc.mpi(): Automated Redistricting Simulation Using
Markov Chain Monte Carlo w/ Parallel Tempering \n\n", append = TRUE)
}
## Extract variables
if (is.na(group_pop)) {
group_pop <- NULL
}
if (is.na(counties)) {
counties <- NULL
}
if (is.na(cities)) {
cities <- NULL
}
if (is.na(areasvec)) {
areasvec <- NULL
}
if (is.na(ssdmat)) {
ssdmat <- matrix(1, 2, 2)
}
if (is.na(borderlength_mat)) {
borderlength_mat <- NULL
}
if (is.na(params$adjswaps)) {
adjswaps <- NULL
} else {
adjswaps <- params$adjswaps
}
if (is.na(params$freq)) {
freq <- NULL
} else {
freq <- params$freq
}
if (is.na(params$constraint)) {
constraint <- NULL
} else {
constraint <- strsplit(params$constraint, ",")[[1]]
}
if (is.na(params$constraintweights)) {
constraintweights <- NULL
} else {
constraintweights <- as.numeric(strsplit(params$constraintweights, ",")[[1]])
}
if (is.na(params$nsims)) {
nsims <- NULL
} else {
nsims <- params$nsims
}
if (is.na(params$nloop)) {
nloop <- NULL
} else {
nloop <- params$nloop
}
if (is.na(params$eprob)) {
eprob <- NULL
} else {
eprob <- params$eprob
}
if (is.na(params$pop_tol)) {
pop_tol <- NULL
} else {
pop_tol <- params$pop_tol
}
if (is.na(params$lambda)) {
lambda <- NULL
} else {
lambda <- params$lambda
}
if (is.na(params$maxiterrsg)) {
maxiterrsg <- NULL
} else {
maxiterrsg <- params$maxiterrsg
}
if (is.na(params$contiguitymap)) {
contiguitymap <- NULL
} else {
contiguitymap <- params$contiguitymap
}
if (is.na(params$loopscompleted)) {
loopscompleted <- NULL
} else {
loopscompleted <- params$loopscompleted
}
if (is.na(params$rngseed)) {
rngseed <- NULL
} else {
rngseed <- params$rngseed
}
if (is.na(params$ndists)) {
ndists <- NULL
} else {
ndists <- params$ndists
}
if (is.na(params$savename)) {
savename <- NULL
} else {
savename <- params$savename
}
if (sum(is.na(init_plan)) == length(init_plan)) {
init_plan <- NULL
}
if (is.na(params$compactness_metric)) {
compactness_metric <- NULL
} else {
compactness_metric <- params$compactness_metric
}
nthin <- params$nthin
## Run redist preprocessing function
preprocout <- redist.preproc(adj = adj, total_pop = total_pop,
init_plan = init_plan, ndists = ndists,
pop_tol = pop_tol,
temper = FALSE,
betaseq = NULL, betaweights = NULL,
adjswaps = adjswaps,
maxiterrsg = maxiterrsg)
## Set betas - if tempering, modified later
beta <- params$beta
weightpop <- preprocout$params$weightpop
weightcompact <- preprocout$params$weightcompact
weightvra <- preprocout$params$weightvra
weightsimilar <- preprocout$params$weightsimilar
weightcountysplit <- preprocout$params$weightcountysplit
weighthinge <- preprocout$params$weighthinge
weightqps <- preprocout$params$weightqps
temper <- "parallel"
## Find procID involved in swaps (non-adjacent only)
if (!adjswaps) {
swapIts <- which(swaps == procID, arr.ind = TRUE)[, 2]
}
## Set seed before first iteration of algorithm if provided by user
if (!is.null(rngseed) & is.numeric(rngseed)) {
set.seed(rngseed)
}
## Get starting loop value
loopstart <- loopscompleted + 1
for (i in loopstart:nloop) {
if (adjswaps) {
nsimsAdj <- rep(freq, nsims/freq)
}
else {
## Construct adjusted "nsims" vector
tempIts <- swapIts[swapIts <= nsims*i & swapIts > nsims*(i - 1)]
## Swap partners
partner <- swaps[, tempIts][swaps[, tempIts] != procID]
nsimsAdj <- c(tempIts, nsims*i) - c((i - 1)*nsims, tempIts)
nsimsAdj <- nsimsAdj[nsimsAdj > 0] # Corrects issue with swaps occurring on nsims*loop
}
## Get initial partition
if (i > loopstart) {
cds <- algout$partitions[, nsims]
if (temper == "parallel") {
beta <- algout$beta_sequence[nsims]
}
if (!is.null(rngseed) & is.numeric(rngseed)) {
set.seed(algout$randseed)
}
rm(list = "algout")
algout <- list()
} else {
## Reload the data if restarting
if (loopstart > 1) {
## Load the data
load(paste(savename, "_proc", procID, "_loop", i - 1, ".RData", sep = ""))
## Load the temperature adjacency matrix (need to specify WD)
load(paste0(savename, "_tempadjMat.RData"))
tempadj <- tempadjMat[nrow(tempadjMat), ]
## Load the swapping schedule (need to specify WD)
load(paste0(savename, "_swaps.RData"))
## Stop if number of simulations per loop is different
if (nsims != ncol(algout[[1]])) {
stop("Please specify the same number of simulations per
loop across all loops")
}
cds <- algout$partitions[, nsims]
if (temper == "parallel") {
beta <- algout$beta_sequence[nsims]
}
if (!is.null(rngseed) & is.numeric(rngseed)) {
set.seed(algout$randseed)
}
rm(list = "algout")
algout <- list()
} else {
cds <- preprocout$data$init_plan
## Initialize algout object (for use in ecutsAppend)
algout <- list()
## Temperature Adjacency Matrix
if (adjswaps) {
tempadjMat <- tempadj
}
}
}
#######################
## Run the algorithm ##
#######################
for (j in 1:length(nsimsAdj)) {
cat("Swap ", j, " out of ", length(nsimsAdj), " swaps.\n",
append = TRUE)
## Run algorithm
temp <- swMH(aList = preprocout$data$adjlist,
cdvec = cds,
popvec = preprocout$data$total_pop,
constraints = list(),
nsims = nsimsAdj[j],
eprob = eprob,
pct_dist_parity = preprocout$params$pctdistparity,
beta_sequence = preprocout$params$betaseq,
beta_weights = preprocout$params$betaweights,
lambda = lambda,
beta = beta,
adapt_beta = "none",
adjswap = preprocout$params$adjswaps,
exact_mh = 0,
adapt_eprob = 0,
adapt_lambda = 0)
## Combine data
algout <- ecutsAppend(algout, temp)
## Get temperature
beta <- temp$beta_sequence[nsimsAdj[j]]
## Get likelihood
like <- temp$energy_psi[nsimsAdj[j]]
## Use MPI to exchange swap information
##
## Tag guide: 1 -> likelihood sent
## 2 -> temperature sent
## 3 -> acceptance probability sent
##
if (adjswaps) {
## Determine which nodes are swapping
tempvra <- swaps[(i - 1)*nsims + j*freq]
## Get node indices
temps <- tempadj[tempvra:(tempvra + 1)]
## Communication step
if (procID %in% temps) {
## Determine partner
partner <- temps[procID != temps]
## Send commands (blocking)
Rmpi::mpi.send.Robj(like, dest = partner, tag = 1)
Rmpi::mpi.send.Robj(beta, dest = partner, tag = 2)
## Receive commands (blocking)
likePart <- Rmpi::mpi.recv.Robj(partner, tag = 1)
betaPart <- Rmpi::mpi.recv.Robj(partner, tag = 2)
## Higher ranked process communicates random
## draw to lower ranked process
if (partner < procID) {
accept <- runif(1)
Rmpi::mpi.send.Robj(accept, dest = partner, tag = 3)
} else {
accept <- Rmpi::mpi.recv.Robj(partner, tag = 3)
}
## Compute acceptance probability (for now, population only)
cat("Components of likelihood:\n", append = TRUE)
cat("like =", like, "\n", append = TRUE)
cat("beta =", beta, "\n", append = TRUE)
cat("betaPart =", betaPart, "\n", append = TRUE)
cat("likePart =", likePart, "\n", append = TRUE)
a_like <- -1*betaPart*like
b_like <- -1*beta*likePart
x_like <- -1*beta*like
y_like <- -1*betaPart*likePart
prob <- exp(a_like + b_like - x_like - y_like)
if (prob > accept) {
## Exchange temperature values
beta <- betaPart
## Adjust temperature adjacency list
tempadj[tempvra:(tempvra + 1)] <- tempadj[(tempvra + 1):tempvra]
## Send temperature adjacency list
if (procID == tempadj[tempvra + 1]) {
oProcs <- tempadj[!(tempadj %in% temps)]
for (k in 1:length(oProcs)) {
Rmpi::mpi.send.Robj(tempadj, dest = oProcs[k], tag = 4)
}
}
} else {
if (procID == tempadj[tempvra]) {
oProcs <- tempadj[!(tempadj %in% temps)]
for (k in 1:length(oProcs)) {
Rmpi::mpi.send.Robj(tempadj, dest = oProcs[k], tag = 4)
}
}
}
} else {
tempadj <- Rmpi::mpi.recv.Robj(tempadj[tempvra], tag = 4)
}
} else {
if (j != length(nsimsAdj) || length(nsimsAdj) == length(tempIts)) {
## Swap proposed
## Send commands (blocking)
Rmpi::mpi.send.Robj(like, dest = partner[j], tag = 1)
Rmpi::mpi.send.Robj(beta, dest = partner[j], tag = 2)
## Receive commands (blocking)
likePart <- Rmpi::mpi.recv.Robj(partner[j], tag = 1)
betaPart <- Rmpi::mpi.recv.Robj(partner[j], tag = 2)
## Higher ranked process communicates random
## draw to lower ranked process
if (partner[j] < procID) {
accept <- runif(1)
Rmpi::mpi.send.Robj(accept, dest = partner[j], tag = 3)
}
else {
accept <- Rmpi::mpi.recv.Robj(partner[j], tag = 3)
}
## Compute acceptance probability
a_like <- -1*betaPart*like
b_like <- -1*beta*likePart
x_like <- -1*beta*like
y_like <- -1*betaPart*likePart
prob <- exp(a_like + b_like - x_like - y_like)
if (prob > accept) {
## Exchange temperature values
beta <- betaPart
}
}
}
## Update inputs to swMH
cds <- temp$partitions[, nsimsAdj[j]]
## Update tempadjMat
if (adjswaps) {
## Update temperature adjacency matrix
tempadjMat <- rbind(tempadjMat, tempadj)
}
## End loop over j
}
class(algout) <- "redist"
## Save random number state if setting the seed
if (!is.null(rngseed)) {
algout$randseed <- .Random.seed[3]
}
## Save output
if (nloop > 1) {
save(algout, file = paste(savename, "_proc", procID, "_loop",
i, ".RData", sep = ""))
## Save temperature adjacency matrix
if (adjswaps) {
save(tempadjMat, file = paste0(savename, "_tempadjMat.RData"))
}
## Save swaps
save(swaps, file = paste0(savename, "_swaps.RData"))
} else if (!is.null(savename)) {
save(algout, file = paste(savename, "_chain",
algout$beta_sequence[nsims],
".RData", sep = ""))
}
## End loop over i
}
###############################
## Combine and save the data ##
###############################
if (nloop > 1) {
if (procID == tempadj[1]) {
redist.combine.mpi(savename = savename, nloop = nloop,
nthin = nthin, tempadj = tempadj)
}
} else if (!is.null(savename)) {
if (procID == tempadj[1]) {
save(algout, file = paste(savename, ".RData", sep = ""))
}
}
if (params$verbose) {
cat("\n", append = TRUE)
cat(divider, append = TRUE)
cat("redist.mcmc.mpi() simulations finished.\n", append = TRUE)
sink()
}
## End function
}
#' Combine successive runs of \code{redist.mcmc.mpi}
#'
#' \code{redist.combine.mpi} is used to combine successive runs of
#' \code{redist.mcmc.mpi} into a single data object
#'
#' @usage redist.combine.mpi(savename, nloop, nthin, tempadj)
#'
#' @param savename The name (without the loop or \code{.RData} suffix)
#' of the saved simulations.
#' @param nloop The number of loops being combined.
#' @param nthin How much to thin the simulations being combined.
#' @param tempadj The temperature adjacency object saved by
#' \code{redist.mcmc.mpi}.
#'
#' @details This function allows users to combine multiple successive runs of
#' \code{redist.mcmc.mpi} into a single \code{redist} object for analysis.
#'
#' @return \code{redist.combine.mpi} returns an object of class "redist".
#' The object \code{redist} is a list that contains the following components (the
#' inclusion of some components is dependent on whether tempering
#' techniques are used):
#' \item{plans}{Matrix of congressional district assignments generated by the
#' algorithm. Each row corresponds to a geographic unit, and each column
#' corresponds to a simulation.}
#' \item{distance_parity}{Vector containing the maximum distance from parity for
#' a particular simulated redistricting plan.}
#' \item{mhdecisions}{A vector specifying whether a proposed redistricting plan
#' was accepted (1) or rejected (0) in a given iteration.}
#' \item{mhprob}{A vector containing the Metropolis-Hastings acceptance
#' probability for each iteration of the algorithm.}
#' \item{pparam}{A vector containing the draw of the \code{p} parameter for each
#' simulation, which dictates the number of swaps attempted.}
#' \item{constraint_pop}{A vector containing the value of the population
#' constraint for each accepted redistricting plan.}
#' \item{constraint_compact}{A vector containing the value of the compactness
#' constraint for each accepted redistricting plan.}
#' \item{constraint_vra}{A vector containing the value of the
#' vra constraint for each accepted redistricting plan.}
#' \item{constraint_similar}{A vector containing the value of the similarity
#' constraint for each accepted redistricting plan.}
#' \item{constraint_qps}{A vector containing the value of the
#' QPS constraint for each accepted redistricting plan.}
#' \item{beta_sequence}{A vector containing the value of beta for each iteration
#' of the algorithm. Returned when tempering is being used.}
#' \item{mhdecisions_beta}{A vector specifying whether a proposed beta value was
#' accepted (1) or rejected (0) in a given iteration of the algorithm. Returned
#' when tempering is being used.}
#' \item{mhprob_beta}{A vector containing the Metropolis-Hastings acceptance
#' probability for each iteration of the algorithm. Returned when tempering
#' is being used.}
#'
#' @references Fifield, Benjamin, Michael Higgins, Kosuke Imai and Alexander Tarr.
#' (2016) "A New Automated Redistricting Simulator Using Markov Chain Monte
#' Carlo." Working Paper. Available at
#' \url{http://imai.princeton.edu/research/files/redist.pdf}.
#'
#' @examples
#' \dontrun{
#' # Cannot run on machines without Rmpi
#' data(fl25)
#' data(fl25_enum)
#' data(fl25_adj)
#'
#' ## Code to run the simulations in Figure 4 in Fifield, Higgins, Imai and
#' ## Tarr (2015)
#'
#' ## Get an initial partition
#' init_plan <- fl25_enum$plans[, 5118]
#'
#' ## Run the algorithm
#' redist.mcmc.mpi(adj = fl25_adj, total_pop = fl25$pop,
#' init_plan = init_plan, nsims = 10000, nloops = 2, savename = "test")
#' out <- redist.combine.mpi(savename = "test", nloop = 2,
#' nthin = 10, tempadj = tempAdjMat)
#' }
#' @concept post
#' @export
redist.combine.mpi <- function(savename, nloop, nthin, tempadj) {
##############################
## Set up container objects ##
##############################
load(paste(savename, "_proc", tempadj[1], "_loop1.RData", sep = ""))
names_obj <- names(algout)
## Create containers
nr <- nrow(algout$partitions)
nc <- ncol(algout$partitions)
partitions <- matrix(NA, nrow = nr, ncol = (nc*nloop/nthin))
veclist <- vector(mode = "list", length = length(algout) - 1)
for (i in 1:length(veclist)) {
veclist[[i]] <- rep(NA, (nc*nloop/nthin))
}
## Indices for thinning
indthin <- which((1:nc) %% nthin == 0)
####################################
## Combine data in multiple loops ##
####################################
for (i in 1:nloop) {
## Load data
load(paste(savename, "_proc", tempadj[1], "_loop", i, ".RData", sep = ""))
ind <- ((i - 1)*(nc/nthin) + 1):(i*(nc/nthin))
## Store objects together
for (j in 1:length(algout)) {
if (j == 1) {
partitions[1:nr, ind] <- algout$partitions[, indthin]
} else {
veclist[[j - 1]][ind] <- algout[[j]][indthin]
}
}
}
#################################
## Store data in algout object ##
#################################
algout <- vector(mode = "list", length = length(algout))
for (i in 1:length(algout)) {
if (i == 1) {
algout[[i]] <- partitions
} else {
algout[[i]] <- veclist[[i - 1]]
}
}
names(algout) <- names_obj
#########################
## Set class of object ##
#########################
class(algout) <- "redist"
#################
## Save object ##
#################
save(algout, file = paste(savename, ".RData", sep = ""))
}
ecutsAppend <- function(algout, ndata) {
if (length(algout) == 0) {
algout <- ndata
} else {
names_obj <- names(algout)
for (i in 1:length(algout)) {
if (i == 1) {
algout[[i]] <- cbind(algout[[i]], ndata[[i]])
} else {
algout[[i]] <- c(algout[[i]], ndata[[i]])
}
}
names(algout) <- names_obj
}
algout
}
#' MCMC Redistricting Simulator using MPI
#'
#' \code{redist.mcmc.mpi} is used to simulate Congressional redistricting
#' plans using Markov Chain Monte Carlo methods.
#'
#'
#' @param adj An adjacency matrix, list, or object of class
#' "SpatialPolygonsDataFrame."
#' @param total_pop A vector containing the populations of each geographic
#' unit.
#' @param nsims The number of simulations run before a save point.
#' @param ndists The number of congressional districts. The default is
#' \code{NULL}.
#' @param init_plan A vector containing the congressional district labels
#' of each geographic unit. The default is \code{NULL}. If not provided, random
#' and contiguous congressional district assignments will be generated using
#' \code{redist.rsg}.
#' @param loopscompleted Number of save points reached by the
#' algorithm. The default is \code{0}.
#' @param nloop The total number of save points for the algorithm. The
#' default is \code{1}. Note that the total number of simulations run
#' will be \code{nsims} * \code{nloop}.
#' @param nthin The amount by which to thin the Markov Chain. The default
#' is \code{1}.
#' @param eprob The probability of keeping an edge connected. The default
#' is \code{0.05}.
#' @param lambda The parameter determining the number of swaps to attempt
#' each iteration of the algorithm. The number of swaps each iteration is
#' equal to Pois(\code{lambda}) + 1. The default is \code{0}.
#' @param pop_tol The strength of the hard population
#' constraint. \code{pop_tol} = 0.05 means that any proposed swap that
#' brings a district more than 5\% away from population parity will be
#' rejected. The default is \code{NULL}.
#' @param group_pop A vector of populations for some sub-group of
#' interest. The default is \code{NULL}.
#' @param areasvec A vector of precinct areas for discrete Polsby-Popper.
#' The default is \code{NULL}.
#' @param counties A vector of county membership assignments. The default is \code{NULL}.
#' @param borderlength_mat A matrix of border length distances, where
#' the first two columns are the indices of precincts sharing a border and
#' the third column is its distance. Default is \code{NULL}.
#' @param ssdmat A matrix of squared distances between geographic
#' units. The default is \code{NULL}.
#' @param compactness_metric The compactness metric to use when constraining on
#' compactness. Default is \code{fryer-holden}, the other implemented option
#' is \code{polsby-popper}.
#' @param rngseed Allows the user to set the seed for the
#' simulations. Default is \code{NULL}.
#' @param constraint Which constraint to apply. Accepts any combination of \code{compact},
#' \code{vra}, \code{population}, \code{similarity}, or \code{none}
#' (no constraint applied). The default is NULL.
#' @param constraintweights The weights to apply to each constraint. Should be a vector
#' the same length as constraint. Default is NULL.
#' @param betaseq Sequence of beta values for tempering. The default is
#' \code{powerlaw} (see Fifield et. al (2015) for details).
#' @param betaseqlength Length of beta sequence desired for
#' tempering. The default is \code{10}.
#' @param adjswaps Flag to restrict swaps of beta so that only
#' values adjacent to current constraint are proposed. The default is
#' \code{TRUE}.
#' @param freq Frequency of between-chain swaps. Default to once every 100
#' iterations
#' @param savename Filename to save simulations. Default is \code{NULL}.
#' @param maxiterrsg Maximum number of iterations for random seed-and-grow
#' algorithm to generate starting values. Default is 5000.
#' @param verbose Whether to print initialization statement. Default is
#' \code{TRUE}.
#' @param cities integer vector of cities for QPS constraint.
#'
#'
#' @details This function allows users to simulate redistricting plans
#' using Markov Chain Monte Carlo methods. Several constraints
#' corresponding to substantive requirements in the redistricting process
#' are implemented, including population parity and geographic
#' compactness. In addition, the function includes multiple-swap and
#' parallel tempering functionality in MPI to improve the mixing of the Markov
#' Chain.
#'
#' @return \code{redist.mcmc.mpi} returns an object of class "redist". The object
#' \code{redist} is a list that contains the following components (the
#' inclusion of some components is dependent on whether tempering
#' techniques are used):
#' \item{partitions}{Matrix of congressional district assignments generated by the
#' algorithm. Each row corresponds to a geographic unit, and each column
#' corresponds to a simulation.}
#' \item{distance_parity}{Vector containing the maximum distance from parity for
#' a particular simulated redistricting plan.}
#' \item{mhdecisions}{A vector specifying whether a proposed redistricting plan
#' was accepted (1) or rejected (0) in a given iteration.}
#' \item{mhprob}{A vector containing the Metropolis-Hastings acceptance
#' probability for each iteration of the algorithm.}
#' \item{pparam}{A vector containing the draw of the \code{p} parameter for each
#' simulation, which dictates the number of swaps attempted.}
#' \item{constraint_pop}{A vector containing the value of the population
#' constraint for each accepted redistricting plan.}
#' \item{constraint_compact}{A vector containing the value of the compactness
#' constraint for each accepted redistricting plan.}
#' \item{constraint_vra}{A vector containing the value of the
#' vra constraint for each accepted redistricting plan.}
#' \item{constraint_similar}{A vector containing the value of the similarity
#' constraint for each accepted redistricting plan.}
#' \item{beta_sequence}{A vector containing the value of beta for each iteration
#' of the algorithm. Returned when tempering is being used.}
#' \item{mhdecisions_beta}{A vector specifying whether a proposed beta value was
#' accepted (1) or rejected (0) in a given iteration of the algorithm. Returned
#' when tempering is being used.}
#' \item{mhprob_beta}{A vector containing the Metropolis-Hastings acceptance
#' probability for each iteration of the algorithm. Returned when tempering
#' is being used.}
#'
#' @references Fifield, Benjamin, Michael Higgins, Kosuke Imai and Alexander
#' Tarr. (2016) "A New Automated Redistricting Simulator Using Markov Chain Monte
#' Carlo." Working Paper. Available at
#' \url{http://imai.princeton.edu/research/files/redist.pdf}.
#'
#' @concept simulate
#' @examples
#' \dontrun{
#' # Cannot run on machines without Rmpi
#' data(fl25)
#' data(fl25_enum)
#' data(fl25_adj)
#'
#' ## Code to run the simulations in Figure 4 in Fifield, Higgins, Imai and
#' ## Tarr (2015)
#'
#' ## Get an initial partition
#' init_plan <- fl25_enum$plans[, 5118]
#'
#' ## Run the algorithm
#' redist.mcmc.mpi(adj = fl25_adj, total_pop = fl25$pop,
#' init_plan = init_plan, nsims = 10000, savename = "test")
#' }
#' @export
redist.mcmc.mpi <- function(adj, total_pop, nsims, ndists = NA,
init_plan = NULL,
loopscompleted = 0, nloop = 1, nthin = 1,
eprob = 0.05, lambda = 0,
pop_tol = NA, group_pop = NA, areasvec = NA,
counties = NA, borderlength_mat = NA,
ssdmat = NA, compactness_metric = "fryer-holden", rngseed = NA,
constraint = NA, constraintweights = NA,
betaseq = "powerlaw", betaseqlength = 10,
adjswaps = TRUE, freq = 100, savename = NA,
maxiterrsg = 5000, verbose = FALSE,
cities = NULL) {
contiguitymap <- "rooks"
## Check if Rmpi library is installed
if (!requireNamespace("Rmpi", quietly = TRUE)) {
stop("You must install package 'Rmpi' to use this function. Please install it if you wish to continue."
, call. = FALSE)
}
## ## Load Rmpi library
## if (!is.loaded("mpi_initialize")) {
## library("Rmpi")
## }
##########################
## Is anything missing? ##
##########################
if (missing(adj)) {
stop("Please supply adjacency matrix or list")
}
if (missing(total_pop)) {
stop("Please supply vector of geographic unit populations")
}
if (missing(nsims)) {
stop("Please supply number of simulations to run algorithm")
}
if (is.na(ndists) & is.null(init_plan)) {
stop("Please provide either the desired number of congressional districts
or an initial set of congressional district assignments")
}
if (nloop > 1 & missing(savename)) {
stop("Please supply save directory if saving simulations at checkpoints")
}
###################
## Preprocessing ##
###################
## Augment init_plan if necessary
nrow.init <- ifelse(is.null(init_plan), 0, nrow(init_plan))
ncol.init <- ifelse(is.null(init_plan), ndists, ncol(init_plan))
if (nrow.init < betaseqlength) {
init_plan <- rbind(init_plan, matrix(NA, betaseqlength - nrow.init, ncol.init))
}
## Generate temperature sequence (power law)
temp <- rep(NA, betaseqlength)
for (i in 1:betaseqlength) {
temp[i] <- 0.1^((i - 1)/(betaseqlength - 1)) - .1
}
beta <- temp/0.9
target.beta <- beta[1]
## Generate swapping sequence
if (adjswaps) {
swaps <- matrix(NA, 1, nsims*(nloop - loopscompleted))
## partner <- matrix(NA,1,nits)
for (i in 1:length(swaps)) {
if (i %% freq == 0) {
swaps[i] <- sample(1:(betaseqlength - 1), size = 1)
}
## Initial temperature adjacency
tempadj <- 1:betaseqlength
}
} else {
swaps <- matrix(NA, 2, nsims*(nloop - loopscompleted))
for (i in 1:ncol(swaps)) {
if (i %% freq == 0) {
swaps[, i] <- sample(1:betaseqlength, size = 2)
}
}
}
## Create parameters list to distribute across nodes
params <- expand.grid(nsims = nsims, nloop = nloop, eprob = eprob,
ndists = ndists, lambda = lambda, pop_tol = pop_tol,
beta = beta, target.beta = target.beta,
constraint = paste(constraint, collapse = ","),
constraintweights = paste(constraintweights, collapse = ","),
compactness_metric = compactness_metric,
betaseqlength = betaseqlength, adjswaps = adjswaps,
nthin = nthin, freq = freq, maxiterrsg = maxiterrsg,
contiguitymap = contiguitymap, verbose = verbose,
loopscompleted = loopscompleted, rngseed = rngseed,
savename = savename, stringsAsFactors = FALSE)
##################
## Spawn Slaves ##
##################
## Note this will not work on Windows platform
Rmpi::mpi.spawn.Rslaves(nslaves = betaseqlength)
## Get processor ID for each slave
Rmpi::mpi.bcast.cmd(procID <- Rmpi::mpi.comm.rank())
#########################
## Send Data to Slaves ##
#########################
## Swapping Schedule
Rmpi::mpi.bcast.Robj2slave(swaps)
## Temperature adjacency
if (adjswaps) {
Rmpi::mpi.bcast.Robj2slave(tempadj)
}
## Adjacency Object
Rmpi::mpi.bcast.Robj2slave(adj)
## Population Vector
Rmpi::mpi.bcast.Robj2slave(total_pop)
## Initial Plans
init_plan <- split(init_plan, f = 1:nrow(init_plan))
Rmpi::mpi.scatter.Robj2slave(init_plan)
## Group population vector
Rmpi::mpi.bcast.Robj2slave(group_pop)
## Areas vector
Rmpi::mpi.bcast.Robj2slave(areasvec)
## County memberhsip vector
Rmpi::mpi.bcast.Robj2slave(counties)
## Border distance matrix
Rmpi::mpi.bcast.Robj2slave(borderlength_mat)
## Squared-distance matrix
Rmpi::mpi.bcast.Robj2slave(ssdmat)
## Parameters List
params <- split(params, f = 1:nrow(params))
Rmpi::mpi.scatter.Robj2slave(params)
## Send ecutsMPI function to slaves
Rmpi::mpi.bcast.Robj2slave(ecutsMPI)
## Send ecutsAppend function to slaves
Rmpi::mpi.bcast.Robj2slave(ecutsAppend)
## Execute ecutsMPI program on each slave
Rmpi::mpi.bcast.cmd(ecutsMPI(procID, params, adj, total_pop, init_plan, swaps))
## Close slaves
Rmpi::mpi.close.Rslaves()
## Terminate MPI processes and close R
Rmpi::mpi.quit()
}
kosukeimai/redist documentation built on March 28, 2024, 7:36 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions redist.mcmc.mpi ecutsAppend redist.combine.mpi ecutsMPI
Documented in redist.combine.mpi redist.mcmc.mpi
##############################################
# Authors: Ben Fifield & Alex Tarr #
# Created: 2015/06/01 #
# Last Revision: N/A #
# Institution: Princeton University #
# Purpose: R wrapper to run swMH() code w/ #
# parallel tempering (mpi) #
##############################################
ecutsMPI <- function(procID = procID, params = params, adj = adj, total_pop = total_pop, init_plan = init_plan, swaps = swaps) {
## Load redist library
library(redist)
if (is.na(params$savename)) {
fname <- paste0("log", procID)
} else {
fname <- paste0("log", procID, "_", params$savename)
}
if (params$verbose) {
sink(fname)
## Initialize ##
divider <- c(paste(rep("=", 20), sep = "", collapse = ""), "\n")
cat("\n", append = TRUE)
cat(divider, append = TRUE)
cat("redist.mcmc.mpi(): Automated Redistricting Simulation Using
Markov Chain Monte Carlo w/ Parallel Tempering \n\n", append = TRUE)
}
## Extract variables
if (is.na(group_pop)) {
group_pop <- NULL
}
if (is.na(counties)) {
counties <- NULL
}
if (is.na(cities)) {
cities <- NULL
}
if (is.na(areasvec)) {
areasvec <- NULL
}
if (is.na(ssdmat)) {
ssdmat <- matrix(1, 2, 2)
}
if (is.na(borderlength_mat)) {
borderlength_mat <- NULL
}
if (is.na(params$adjswaps)) {
adjswaps <- NULL
} else {
adjswaps <- params$adjswaps
}
if (is.na(params$freq)) {
freq <- NULL
} else {
freq <- params$freq
}
if (is.na(params$constraint)) {
constraint <- NULL
} else {
constraint <- strsplit(params$constraint, ",")[[1]]
}
if (is.na(params$constraintweights)) {
constraintweights <- NULL
} else {
constraintweights <- as.numeric(strsplit(params$constraintweights, ",")[[1]])
}
if (is.na(params$nsims)) {
nsims <- NULL
} else {
nsims <- params$nsims
}
if (is.na(params$nloop)) {
nloop <- NULL
} else {
nloop <- params$nloop
}
if (is.na(params$eprob)) {
eprob <- NULL
} else {
eprob <- params$eprob
}
if (is.na(params$pop_tol)) {
pop_tol <- NULL
} else {
pop_tol <- params$pop_tol
}
if (is.na(params$lambda)) {
lambda <- NULL
} else {
lambda <- params$lambda
}
if (is.na(params$maxiterrsg)) {
maxiterrsg <- NULL
} else {
maxiterrsg <- params$maxiterrsg
}
if (is.na(params$contiguitymap)) {
contiguitymap <- NULL
} else {
contiguitymap <- params$contiguitymap
}
if (is.na(params$loopscompleted)) {
loopscompleted <- NULL
} else {
loopscompleted <- params$loopscompleted
}
if (is.na(params$rngseed)) {
rngseed <- NULL
} else {
rngseed <- params$rngseed
}
if (is.na(params$ndists)) {
ndists <- NULL
} else {
ndists <- params$ndists
}
if (is.na(params$savename)) {
savename <- NULL
} else {
savename <- params$savename
}
if (sum(is.na(init_plan)) == length(init_plan)) {
init_plan <- NULL
}
if (is.na(params$compactness_metric)) {
compactness_metric <- NULL
} else {
compactness_metric <- params$compactness_metric
}
nthin <- params$nthin
## Run redist preprocessing function
preprocout <- redist.preproc(adj = adj, total_pop = total_pop,
init_plan = init_plan, ndists = ndists,
pop_tol = pop_tol,
temper = FALSE,
betaseq = NULL, betaweights = NULL,
adjswaps = adjswaps,
maxiterrsg = maxiterrsg)
## Set betas - if tempering, modified later
beta <- params$beta
weightpop <- preprocout$params$weightpop
weightcompact <- preprocout$params$weightcompact
weightvra <- preprocout$params$weightvra
weightsimilar <- preprocout$params$weightsimilar
weightcountysplit <- preprocout$params$weightcountysplit
weighthinge <- preprocout$params$weighthinge
weightqps <- preprocout$params$weightqps
temper <- "parallel"
## Find procID involved in swaps (non-adjacent only)
if (!adjswaps) {
swapIts <- which(swaps == procID, arr.ind = TRUE)[, 2]
}
## Set seed before first iteration of algorithm if provided by user
if (!is.null(rngseed) & is.numeric(rngseed)) {
set.seed(rngseed)
}
## Get starting loop value
loopstart <- loopscompleted + 1
for (i in loopstart:nloop) {
if (adjswaps) {
nsimsAdj <- rep(freq, nsims/freq)
}
else {
## Construct adjusted "nsims" vector
tempIts <- swapIts[swapIts <= nsims*i & swapIts > nsims*(i - 1)]
## Swap partners
partner <- swaps[, tempIts][swaps[, tempIts] != procID]
nsimsAdj <- c(tempIts, nsims*i) - c((i - 1)*nsims, tempIts)
nsimsAdj <- nsimsAdj[nsimsAdj > 0] # Corrects issue with swaps occurring on nsims*loop
}
## Get initial partition
if (i > loopstart) {
cds <- algout$partitions[, nsims]
if (temper == "parallel") {
beta <- algout$beta_sequence[nsims]
}
if (!is.null(rngseed) & is.numeric(rngseed)) {
set.seed(algout$randseed)
}
rm(list = "algout")
algout <- list()
} else {
## Reload the data if restarting
if (loopstart > 1) {
## Load the data
load(paste(savename, "_proc", procID, "_loop", i - 1, ".RData", sep = ""))
## Load the temperature adjacency matrix (need to specify WD)
load(paste0(savename, "_tempadjMat.RData"))
tempadj <- tempadjMat[nrow(tempadjMat), ]
## Load the swapping schedule (need to specify WD)
load(paste0(savename, "_swaps.RData"))
## Stop if number of simulations per loop is different
if (nsims != ncol(algout[[1]])) {
stop("Please specify the same number of simulations per
loop across all loops")
}
cds <- algout$partitions[, nsims]
if (temper == "parallel") {
beta <- algout$beta_sequence[nsims]
}
if (!is.null(rngseed) & is.numeric(rngseed)) {
set.seed(algout$randseed)
}
rm(list = "algout")
algout <- list()
} else {
cds <- preprocout$data$init_plan
## Initialize algout object (for use in ecutsAppend)
algout <- list()
## Temperature Adjacency Matrix
if (adjswaps) {
tempadjMat <- tempadj
}
}
}
#######################
## Run the algorithm ##
#######################
for (j in 1:length(nsimsAdj)) {
cat("Swap ", j, " out of ", length(nsimsAdj), " swaps.\n",
append = TRUE)
## Run algorithm
temp <- swMH(aList = preprocout$data$adjlist,
cdvec = cds,
popvec = preprocout$data$total_pop,
constraints = list(),
nsims = nsimsAdj[j],
eprob = eprob,
pct_dist_parity = preprocout$params$pctdistparity,
beta_sequence = preprocout$params$betaseq,
beta_weights = preprocout$params$betaweights,
lambda = lambda,
beta = beta,
adapt_beta = "none",
adjswap = preprocout$params$adjswaps,
exact_mh = 0,
adapt_eprob = 0,
adapt_lambda = 0)
## Combine data
algout <- ecutsAppend(algout, temp)
## Get temperature
beta <- temp$beta_sequence[nsimsAdj[j]]
## Get likelihood
like <- temp$energy_psi[nsimsAdj[j]]
## Use MPI to exchange swap information
##
## Tag guide: 1 -> likelihood sent
## 2 -> temperature sent
## 3 -> acceptance probability sent
##
if (adjswaps) {
## Determine which nodes are swapping
tempvra <- swaps[(i - 1)*nsims + j*freq]
## Get node indices
temps <- tempadj[tempvra:(tempvra + 1)]
## Communication step
if (procID %in% temps) {
## Determine partner
partner <- temps[procID != temps]
## Send commands (blocking)
Rmpi::mpi.send.Robj(like, dest = partner, tag = 1)
Rmpi::mpi.send.Robj(beta, dest = partner, tag = 2)
## Receive commands (blocking)
likePart <- Rmpi::mpi.recv.Robj(partner, tag = 1)
betaPart <- Rmpi::mpi.recv.Robj(partner, tag = 2)
## Higher ranked process communicates random
## draw to lower ranked process
if (partner < procID) {
accept <- runif(1)
Rmpi::mpi.send.Robj(accept, dest = partner, tag = 3)
} else {
accept <- Rmpi::mpi.recv.Robj(partner, tag = 3)
}
## Compute acceptance probability (for now, population only)
cat("Components of likelihood:\n", append = TRUE)
cat("like =", like, "\n", append = TRUE)
cat("beta =", beta, "\n", append = TRUE)
cat("betaPart =", betaPart, "\n", append = TRUE)
cat("likePart =", likePart, "\n", append = TRUE)
a_like <- -1*betaPart*like
b_like <- -1*beta*likePart
x_like <- -1*beta*like
y_like <- -1*betaPart*likePart
prob <- exp(a_like + b_like - x_like - y_like)
if (prob > accept) {
## Exchange temperature values
beta <- betaPart
## Adjust temperature adjacency list
tempadj[tempvra:(tempvra + 1)] <- tempadj[(tempvra + 1):tempvra]
## Send temperature adjacency list
if (procID == tempadj[tempvra + 1]) {
oProcs <- tempadj[!(tempadj %in% temps)]
for (k in 1:length(oProcs)) {
Rmpi::mpi.send.Robj(tempadj, dest = oProcs[k], tag = 4)
}
}
} else {
if (procID == tempadj[tempvra]) {
oProcs <- tempadj[!(tempadj %in% temps)]
for (k in 1:length(oProcs)) {
Rmpi::mpi.send.Robj(tempadj, dest = oProcs[k], tag = 4)
}
}
}
} else {
tempadj <- Rmpi::mpi.recv.Robj(tempadj[tempvra], tag = 4)
}
} else {
if (j != length(nsimsAdj) || length(nsimsAdj) == length(tempIts)) {
## Swap proposed
## Send commands (blocking)
Rmpi::mpi.send.Robj(like, dest = partner[j], tag = 1)
Rmpi::mpi.send.Robj(beta, dest = partner[j], tag = 2)
## Receive commands (blocking)
likePart <- Rmpi::mpi.recv.Robj(partner[j], tag = 1)
betaPart <- Rmpi::mpi.recv.Robj(partner[j], tag = 2)
## Higher ranked process communicates random
## draw to lower ranked process
if (partner[j] < procID) {
accept <- runif(1)
Rmpi::mpi.send.Robj(accept, dest = partner[j], tag = 3)
}
else {
accept <- Rmpi::mpi.recv.Robj(partner[j], tag = 3)
}
## Compute acceptance probability
a_like <- -1*betaPart*like
b_like <- -1*beta*likePart
x_like <- -1*beta*like
y_like <- -1*betaPart*likePart
prob <- exp(a_like + b_like - x_like - y_like)
if (prob > accept) {
## Exchange temperature values
beta <- betaPart
}
}
}
## Update inputs to swMH
cds <- temp$partitions[, nsimsAdj[j]]
## Update tempadjMat
if (adjswaps) {
## Update temperature adjacency matrix
tempadjMat <- rbind(tempadjMat, tempadj)
}
## End loop over j
}
class(algout) <- "redist"
## Save random number state if setting the seed
if (!is.null(rngseed)) {
algout$randseed <- .Random.seed[3]
}
## Save output
if (nloop > 1) {
save(algout, file = paste(savename, "_proc", procID, "_loop",
i, ".RData", sep = ""))
## Save temperature adjacency matrix
if (adjswaps) {
save(tempadjMat, file = paste0(savename, "_tempadjMat.RData"))
}
## Save swaps
save(swaps, file = paste0(savename, "_swaps.RData"))
} else if (!is.null(savename)) {
save(algout, file = paste(savename, "_chain",
algout$beta_sequence[nsims],
".RData", sep = ""))
}
## End loop over i
}
###############################
## Combine and save the data ##
###############################
if (nloop > 1) {
if (procID == tempadj[1]) {
redist.combine.mpi(savename = savename, nloop = nloop,
nthin = nthin, tempadj = tempadj)
}
} else if (!is.null(savename)) {
if (procID == tempadj[1]) {
save(algout, file = paste(savename, ".RData", sep = ""))
}
}
if (params$verbose) {
cat("\n", append = TRUE)
cat(divider, append = TRUE)
cat("redist.mcmc.mpi() simulations finished.\n", append = TRUE)
sink()
}
## End function
}
#' Combine successive runs of \code{redist.mcmc.mpi}
#'
#' \code{redist.combine.mpi} is used to combine successive runs of
#' \code{redist.mcmc.mpi} into a single data object
#'
#' @usage redist.combine.mpi(savename, nloop, nthin, tempadj)
#'
#' @param savename The name (without the loop or \code{.RData} suffix)
#' of the saved simulations.
#' @param nloop The number of loops being combined.
#' @param nthin How much to thin the simulations being combined.
#' @param tempadj The temperature adjacency object saved by
#' \code{redist.mcmc.mpi}.
#'
#' @details This function allows users to combine multiple successive runs of
#' \code{redist.mcmc.mpi} into a single \code{redist} object for analysis.
#'
#' @return \code{redist.combine.mpi} returns an object of class "redist".
#' The object \code{redist} is a list that contains the following components (the
#' inclusion of some components is dependent on whether tempering
#' techniques are used):
#' \item{plans}{Matrix of congressional district assignments generated by the
#' algorithm. Each row corresponds to a geographic unit, and each column
#' corresponds to a simulation.}
#' \item{distance_parity}{Vector containing the maximum distance from parity for
#' a particular simulated redistricting plan.}
#' \item{mhdecisions}{A vector specifying whether a proposed redistricting plan
#' was accepted (1) or rejected (0) in a given iteration.}
#' \item{mhprob}{A vector containing the Metropolis-Hastings acceptance
#' probability for each iteration of the algorithm.}
#' \item{pparam}{A vector containing the draw of the \code{p} parameter for each
#' simulation, which dictates the number of swaps attempted.}
#' \item{constraint_pop}{A vector containing the value of the population
#' constraint for each accepted redistricting plan.}
#' \item{constraint_compact}{A vector containing the value of the compactness
#' constraint for each accepted redistricting plan.}
#' \item{constraint_vra}{A vector containing the value of the
#' vra constraint for each accepted redistricting plan.}
#' \item{constraint_similar}{A vector containing the value of the similarity
#' constraint for each accepted redistricting plan.}
#' \item{constraint_qps}{A vector containing the value of the
#' QPS constraint for each accepted redistricting plan.}
#' \item{beta_sequence}{A vector containing the value of beta for each iteration
#' of the algorithm. Returned when tempering is being used.}
#' \item{mhdecisions_beta}{A vector specifying whether a proposed beta value was
#' accepted (1) or rejected (0) in a given iteration of the algorithm. Returned
#' when tempering is being used.}
#' \item{mhprob_beta}{A vector containing the Metropolis-Hastings acceptance
#' probability for each iteration of the algorithm. Returned when tempering
#' is being used.}
#'
#' @references Fifield, Benjamin, Michael Higgins, Kosuke Imai and Alexander Tarr.
#' (2016) "A New Automated Redistricting Simulator Using Markov Chain Monte
#' Carlo." Working Paper. Available at
#' \url{http://imai.princeton.edu/research/files/redist.pdf}.
#'
#' @examples
#' \dontrun{
#' # Cannot run on machines without Rmpi
#' data(fl25)
#' data(fl25_enum)
#' data(fl25_adj)
#'
#' ## Code to run the simulations in Figure 4 in Fifield, Higgins, Imai and
#' ## Tarr (2015)
#'
#' ## Get an initial partition
#' init_plan <- fl25_enum$plans[, 5118]
#'
#' ## Run the algorithm
#' redist.mcmc.mpi(adj = fl25_adj, total_pop = fl25$pop,
#' init_plan = init_plan, nsims = 10000, nloops = 2, savename = "test")
#' out <- redist.combine.mpi(savename = "test", nloop = 2,
#' nthin = 10, tempadj = tempAdjMat)
#' }
#' @concept post
#' @export
redist.combine.mpi <- function(savename, nloop, nthin, tempadj) {
##############################
## Set up container objects ##
##############################
load(paste(savename, "_proc", tempadj[1], "_loop1.RData", sep = ""))
names_obj <- names(algout)
## Create containers
nr <- nrow(algout$partitions)
nc <- ncol(algout$partitions)
partitions <- matrix(NA, nrow = nr, ncol = (nc*nloop/nthin))
veclist <- vector(mode = "list", length = length(algout) - 1)
for (i in 1:length(veclist)) {
veclist[[i]] <- rep(NA, (nc*nloop/nthin))
}
## Indices for thinning
indthin <- which((1:nc) %% nthin == 0)
####################################
## Combine data in multiple loops ##
####################################
for (i in 1:nloop) {
## Load data
load(paste(savename, "_proc", tempadj[1], "_loop", i, ".RData", sep = ""))
ind <- ((i - 1)*(nc/nthin) + 1):(i*(nc/nthin))
## Store objects together
for (j in 1:length(algout)) {
if (j == 1) {
partitions[1:nr, ind] <- algout$partitions[, indthin]
} else {
veclist[[j - 1]][ind] <- algout[[j]][indthin]
}
}
}
#################################
## Store data in algout object ##
#################################
algout <- vector(mode = "list", length = length(algout))
for (i in 1:length(algout)) {
if (i == 1) {
algout[[i]] <- partitions
} else {
algout[[i]] <- veclist[[i - 1]]
}
}
names(algout) <- names_obj
#########################
## Set class of object ##
#########################
class(algout) <- "redist"
#################
## Save object ##
#################
save(algout, file = paste(savename, ".RData", sep = ""))
}
ecutsAppend <- function(algout, ndata) {
if (length(algout) == 0) {
algout <- ndata
} else {
names_obj <- names(algout)
for (i in 1:length(algout)) {
if (i == 1) {
algout[[i]] <- cbind(algout[[i]], ndata[[i]])
} else {
algout[[i]] <- c(algout[[i]], ndata[[i]])
}
}
names(algout) <- names_obj
}
algout
}
#' MCMC Redistricting Simulator using MPI
#'
#' \code{redist.mcmc.mpi} is used to simulate Congressional redistricting
#' plans using Markov Chain Monte Carlo methods.
#'
#'
#' @param adj An adjacency matrix, list, or object of class
#' "SpatialPolygonsDataFrame."
#' @param total_pop A vector containing the populations of each geographic
#' unit.
#' @param nsims The number of simulations run before a save point.
#' @param ndists The number of congressional districts. The default is
#' \code{NULL}.
#' @param init_plan A vector containing the congressional district labels
#' of each geographic unit. The default is \code{NULL}. If not provided, random
#' and contiguous congressional district assignments will be generated using
#' \code{redist.rsg}.
#' @param loopscompleted Number of save points reached by the
#' algorithm. The default is \code{0}.
#' @param nloop The total number of save points for the algorithm. The
#' default is \code{1}. Note that the total number of simulations run
#' will be \code{nsims} * \code{nloop}.
#' @param nthin The amount by which to thin the Markov Chain. The default
#' is \code{1}.
#' @param eprob The probability of keeping an edge connected. The default
#' is \code{0.05}.
#' @param lambda The parameter determining the number of swaps to attempt
#' each iteration of the algorithm. The number of swaps each iteration is
#' equal to Pois(\code{lambda}) + 1. The default is \code{0}.
#' @param pop_tol The strength of the hard population
#' constraint. \code{pop_tol} = 0.05 means that any proposed swap that
#' brings a district more than 5\% away from population parity will be
#' rejected. The default is \code{NULL}.
#' @param group_pop A vector of populations for some sub-group of
#' interest. The default is \code{NULL}.
#' @param areasvec A vector of precinct areas for discrete Polsby-Popper.
#' The default is \code{NULL}.
#' @param counties A vector of county membership assignments. The default is \code{NULL}.
#' @param borderlength_mat A matrix of border length distances, where
#' the first two columns are the indices of precincts sharing a border and
#' the third column is its distance. Default is \code{NULL}.
#' @param ssdmat A matrix of squared distances between geographic
#' units. The default is \code{NULL}.
#' @param compactness_metric The compactness metric to use when constraining on
#' compactness. Default is \code{fryer-holden}, the other implemented option
#' is \code{polsby-popper}.
#' @param rngseed Allows the user to set the seed for the
#' simulations. Default is \code{NULL}.
#' @param constraint Which constraint to apply. Accepts any combination of \code{compact},
#' \code{vra}, \code{population}, \code{similarity}, or \code{none}
#' (no constraint applied). The default is NULL.
#' @param constraintweights The weights to apply to each constraint. Should be a vector
#' the same length as constraint. Default is NULL.
#' @param betaseq Sequence of beta values for tempering. The default is
#' \code{powerlaw} (see Fifield et. al (2015) for details).
#' @param betaseqlength Length of beta sequence desired for
#' tempering. The default is \code{10}.
#' @param adjswaps Flag to restrict swaps of beta so that only
#' values adjacent to current constraint are proposed. The default is
#' \code{TRUE}.
#' @param freq Frequency of between-chain swaps. Default to once every 100
#' iterations
#' @param savename Filename to save simulations. Default is \code{NULL}.
#' @param maxiterrsg Maximum number of iterations for random seed-and-grow
#' algorithm to generate starting values. Default is 5000.
#' @param verbose Whether to print initialization statement. Default is
#' \code{TRUE}.
#' @param cities integer vector of cities for QPS constraint.
#'
#'
#' @details This function allows users to simulate redistricting plans
#' using Markov Chain Monte Carlo methods. Several constraints
#' corresponding to substantive requirements in the redistricting process
#' are implemented, including population parity and geographic
#' compactness. In addition, the function includes multiple-swap and
#' parallel tempering functionality in MPI to improve the mixing of the Markov
#' Chain.
#'
#' @return \code{redist.mcmc.mpi} returns an object of class "redist". The object
#' \code{redist} is a list that contains the following components (the
#' inclusion of some components is dependent on whether tempering
#' techniques are used):
#' \item{partitions}{Matrix of congressional district assignments generated by the
#' algorithm. Each row corresponds to a geographic unit, and each column
#' corresponds to a simulation.}
#' \item{distance_parity}{Vector containing the maximum distance from parity for
#' a particular simulated redistricting plan.}
#' \item{mhdecisions}{A vector specifying whether a proposed redistricting plan
#' was accepted (1) or rejected (0) in a given iteration.}
#' \item{mhprob}{A vector containing the Metropolis-Hastings acceptance
#' probability for each iteration of the algorithm.}
#' \item{pparam}{A vector containing the draw of the \code{p} parameter for each
#' simulation, which dictates the number of swaps attempted.}
#' \item{constraint_pop}{A vector containing the value of the population
#' constraint for each accepted redistricting plan.}
#' \item{constraint_compact}{A vector containing the value of the compactness
#' constraint for each accepted redistricting plan.}
#' \item{constraint_vra}{A vector containing the value of the
#' vra constraint for each accepted redistricting plan.}
#' \item{constraint_similar}{A vector containing the value of the similarity
#' constraint for each accepted redistricting plan.}
#' \item{beta_sequence}{A vector containing the value of beta for each iteration
#' of the algorithm. Returned when tempering is being used.}
#' \item{mhdecisions_beta}{A vector specifying whether a proposed beta value was
#' accepted (1) or rejected (0) in a given iteration of the algorithm. Returned
#' when tempering is being used.}
#' \item{mhprob_beta}{A vector containing the Metropolis-Hastings acceptance
#' probability for each iteration of the algorithm. Returned when tempering
#' is being used.}
#'
#' @references Fifield, Benjamin, Michael Higgins, Kosuke Imai and Alexander
#' Tarr. (2016) "A New Automated Redistricting Simulator Using Markov Chain Monte
#' Carlo." Working Paper. Available at
#' \url{http://imai.princeton.edu/research/files/redist.pdf}.
#'
#' @concept simulate
#' @examples
#' \dontrun{
#' # Cannot run on machines without Rmpi
#' data(fl25)
#' data(fl25_enum)
#' data(fl25_adj)
#'
#' ## Code to run the simulations in Figure 4 in Fifield, Higgins, Imai and
#' ## Tarr (2015)
#'
#' ## Get an initial partition
#' init_plan <- fl25_enum$plans[, 5118]
#'
#' ## Run the algorithm
#' redist.mcmc.mpi(adj = fl25_adj, total_pop = fl25$pop,
#' init_plan = init_plan, nsims = 10000, savename = "test")
#' }
#' @export
redist.mcmc.mpi <- function(adj, total_pop, nsims, ndists = NA,
init_plan = NULL,
loopscompleted = 0, nloop = 1, nthin = 1,
eprob = 0.05, lambda = 0,
pop_tol = NA, group_pop = NA, areasvec = NA,
counties = NA, borderlength_mat = NA,
ssdmat = NA, compactness_metric = "fryer-holden", rngseed = NA,
constraint = NA, constraintweights = NA,
betaseq = "powerlaw", betaseqlength = 10,
adjswaps = TRUE, freq = 100, savename = NA,
maxiterrsg = 5000, verbose = FALSE,
cities = NULL) {
contiguitymap <- "rooks"
## Check if Rmpi library is installed
if (!requireNamespace("Rmpi", quietly = TRUE)) {
stop("You must install package 'Rmpi' to use this function. Please install it if you wish to continue."
, call. = FALSE)
}
## ## Load Rmpi library
## if (!is.loaded("mpi_initialize")) {
## library("Rmpi")
## }
##########################
## Is anything missing? ##
##########################
if (missing(adj)) {
stop("Please supply adjacency matrix or list")
}
if (missing(total_pop)) {
stop("Please supply vector of geographic unit populations")
}
if (missing(nsims)) {
stop("Please supply number of simulations to run algorithm")
}
if (is.na(ndists) & is.null(init_plan)) {
stop("Please provide either the desired number of congressional districts
or an initial set of congressional district assignments")
}
if (nloop > 1 & missing(savename)) {
stop("Please supply save directory if saving simulations at checkpoints")
}
###################
## Preprocessing ##
###################
## Augment init_plan if necessary
nrow.init <- ifelse(is.null(init_plan), 0, nrow(init_plan))
ncol.init <- ifelse(is.null(init_plan), ndists, ncol(init_plan))
if (nrow.init < betaseqlength) {
init_plan <- rbind(init_plan, matrix(NA, betaseqlength - nrow.init, ncol.init))
}
## Generate temperature sequence (power law)
temp <- rep(NA, betaseqlength)
for (i in 1:betaseqlength) {
temp[i] <- 0.1^((i - 1)/(betaseqlength - 1)) - .1
}
beta <- temp/0.9
target.beta <- beta[1]
## Generate swapping sequence
if (adjswaps) {
swaps <- matrix(NA, 1, nsims*(nloop - loopscompleted))
## partner <- matrix(NA,1,nits)
for (i in 1:length(swaps)) {
if (i %% freq == 0) {
swaps[i] <- sample(1:(betaseqlength - 1), size = 1)
}
## Initial temperature adjacency
tempadj <- 1:betaseqlength
}
} else {
swaps <- matrix(NA, 2, nsims*(nloop - loopscompleted))
for (i in 1:ncol(swaps)) {
if (i %% freq == 0) {
swaps[, i] <- sample(1:betaseqlength, size = 2)
}
}
}
## Create parameters list to distribute across nodes
params <- expand.grid(nsims = nsims, nloop = nloop, eprob = eprob,
ndists = ndists, lambda = lambda, pop_tol = pop_tol,
beta = beta, target.beta = target.beta,
constraint = paste(constraint, collapse = ","),
constraintweights = paste(constraintweights, collapse = ","),
compactness_metric = compactness_metric,
betaseqlength = betaseqlength, adjswaps = adjswaps,
nthin = nthin, freq = freq, maxiterrsg = maxiterrsg,
contiguitymap = contiguitymap, verbose = verbose,
loopscompleted = loopscompleted, rngseed = rngseed,
savename = savename, stringsAsFactors = FALSE)
##################
## Spawn Slaves ##
##################
## Note this will not work on Windows platform
Rmpi::mpi.spawn.Rslaves(nslaves = betaseqlength)
## Get processor ID for each slave
Rmpi::mpi.bcast.cmd(procID <- Rmpi::mpi.comm.rank())
#########################
## Send Data to Slaves ##
#########################
## Swapping Schedule
Rmpi::mpi.bcast.Robj2slave(swaps)
## Temperature adjacency
if (adjswaps) {
Rmpi::mpi.bcast.Robj2slave(tempadj)
}
## Adjacency Object
Rmpi::mpi.bcast.Robj2slave(adj)
## Population Vector
Rmpi::mpi.bcast.Robj2slave(total_pop)
## Initial Plans
init_plan <- split(init_plan, f = 1:nrow(init_plan))
Rmpi::mpi.scatter.Robj2slave(init_plan)
## Group population vector
Rmpi::mpi.bcast.Robj2slave(group_pop)
## Areas vector
Rmpi::mpi.bcast.Robj2slave(areasvec)
## County memberhsip vector
Rmpi::mpi.bcast.Robj2slave(counties)
## Border distance matrix
Rmpi::mpi.bcast.Robj2slave(borderlength_mat)
## Squared-distance matrix
Rmpi::mpi.bcast.Robj2slave(ssdmat)
## Parameters List
params <- split(params, f = 1:nrow(params))
Rmpi::mpi.scatter.Robj2slave(params)
## Send ecutsMPI function to slaves
Rmpi::mpi.bcast.Robj2slave(ecutsMPI)
## Send ecutsAppend function to slaves
Rmpi::mpi.bcast.Robj2slave(ecutsAppend)
## Execute ecutsMPI program on each slave
Rmpi::mpi.bcast.cmd(ecutsMPI(procID, params, adj, total_pop, init_plan, swaps))
## Close slaves
Rmpi::mpi.close.Rslaves()
## Terminate MPI processes and close R
Rmpi::mpi.quit()
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.