################################################### ### Preliminaries ################################################### options(useFancyQuotes = FALSE, warn = -1) RECMAP_ENABLE_PARALLEL <- FALSE message("This manuscript requires the R packages: recmap (>=0.5.34), maps, noncensus.") message(paste("recmap version", packageVersion("recmap"), "is used.")) stopifnot(packageVersion("recmap") >= "0.5.34") stopifnot(packageVersion("noncensus") >= "0.1") stopifnot(packageVersion("maps") >= "2.0") library("recmap") library("maps") library("noncensus") library("colorspace") library("xtable") data("jss2711", package = "recmap") set.seed(1) ################################################### ### 1. Introduction ################################################### ## US Election 2004 - Figure 1 op <- par(mar = c(0, 0, 0, 0)) recmap:::.draw_recmap_us_state_ev() par(op) ################################################### ### 3. The package usage ################################################### ## 3.1. Input US.map <- data.frame(x = state.center$x, y = state.center$y, dx = sqrt(state.area) / 2 / (0.7 * 60 * cos(state.center$y * pi / 180)), dy = sqrt(state.area) / 2 / (0.7 * 60), z = sqrt(state.area), name = state.name) head(US.map) attr(US.map, "Map.name") <- "U.S." attr(US.map, "Map.area") <- "area" ## 3.2. Run library("recmap") ## Generate the rectangular statistical cartogram. US.cartogram <- recmap(US.map) ## 3.3. Output head(US.cartogram) op <- par(mfrow = c(1, 2), mar = c(0, 0, 0, 0)) plot.recmap(US.map, col.text = "darkred") plot(US.cartogram, col.text = "darkred") summary(US.cartogram) X <- checkerboard(8) all.equal(X, as.recmap(as.SpatialPolygonsDataFrame(X))) ################################################### ### 4. Implementation ################################################### ## Figure 3 (right) place_rectangle <- function(alpha = 0.0, a.x = 0, a.y = 0, a.dx = 1, a.dy = 1, b.dx = 2, b.dy = 2, text = "", ...) { rect(a.x - a.dx - b.dx, a.y - a.dy - b.dy, a.x + a.dx + b.dx, a.y + a.dy + b.dy, lwd = 3, border = "grey") rect(a.x - a.dx, a.y - a.dy, a.x + a.dx, a.y + a.dy) c <- recmap:::place_rectangle(a.x, a.y, a.dx, a.dy, b.dx, b.dy, alpha) rect(c$x - c$dx, c$y - c$dy, c$x + c$dx, c$y + c$dy, ...) text(c$x, c$y - c$dy, text, pos = 3) } op <- par(mfrow = c(3, 3), mar = c(0, 0, 3, 0)) beta <- seq(0, 0.35 * pi, length = 9) * rep(c(1, -1), 5)[0:9] rv <- lapply(1:length(beta), function(id) { beta0 <- beta[id] plot(0, 0, xlim = c(-5, 5), ylim = c(-5, 5), axes = FALSE, type = "n", main = substitute(paste(beta, "=", b), list(b = round(beta0, 2))), xlab = "", ylab = "") place_rectangle(alpha = 0, col = "#44444444") text(0, 0, "A") if (round(beta0, 1) != 1.1) { cc <- "#AA444444" } else { cc <- "#44AA4444" } place_rectangle(alpha = (pi/4 + beta0), col = cc, text = "B") }) op <- par(mar = c(4, 4, 1.5, 0.25)) plot(0, 0, type = "n", xlim = c(-7, 7), ylim = c(-7, 7), main = expression(f[pl]:bold(R)^6 %*% group("[", list(-pi, pi), "]") %->% bold(R)^2), asp = 1, xlab = "x return value", ylab = "y return value", axes = FALSE) n <- 90 cm <- rainbow_hcl(n, alpha = 0.5) alpha <- seq(-pi, pi, length = n) r <- lapply(1:n, function(idx) { place_rectangle(alpha[idx], col = cm[idx], border = "#88888877") }) legend("bottomleft", as.character(round(seq(-pi, pi, length = 9), 2)), bty = "n", fill = rainbow_hcl(9, alpha = 0.5), border = "white", col = "white", title = expression(paste(alpha, " argument", sep = " ")), cex = 0.66, horiz = FALSE, bg = "white") text(0, 0, "A") par(op) op <- par(mar = c(4, 4, 1, 1), mfrow = c(1, 3)) boxplot(number ~ sqrt(size), axes = FALSE, data = mbb_check, log = "y", cex = 0.75, subset = algorithm == "list", col = "red", at = 2:20 - 0.2, boxwex = 0.25) abline(v = sqrt(50), col = "lightgray", lwd = 3) boxplot(number ~ sqrt(size), data = mbb_check, log = "y", subset = algorithm == "multiset", at = 2:20 + 0.2, cex = 0.75, ylab = "number of MBB intersection calls", xlab = "number of map regions", boxwex = 0.25, add = TRUE, axes = FALSE) axis(2) axis(1, c(5, sqrt(50), 10, 15, 20), c("5x5", "US", "10x10", "15x15", "20x20")) box() legend("bottomright", c("C++ STL list", "C++ STL multiset"), col = c("red", "black"), pch = 16, cex = 1.0) ################################################################################ S <- aggregate(time_in_secs ~ size * sysname * algorithm, data = mbb_check, FUN = median) plot(time_in_secs ~ sqrt(size), data = S, axes = FALSE, log = "y", xlab = "number of map regions", ylab = "median aggregated process time [in seconds]", type = "n") abline(v = sqrt(50), col = "lightgray", lwd = 3) S.f <- S[S$sysname == "Linux" & S$algorithm == "list", c("size", "time_in_secs")] points(sqrt(S.f$size), S.f$time_in_secs, pch = "L", type = "b", col = "red", cex = 0.5) S.f <- S[S$sysname == "Linux" & S$algorithm == "multiset", c("size", "time_in_secs")] points(sqrt(S.f$size), S.f$time_in_secs, pch = "L", type = "b", cex = 0.5) S.f <- S[S$sysname == "Darwin" & S$algorithm == "list", c("size", "time_in_secs")] points(sqrt(S.f$size), S.f$time_in_secs, pch = "A", type = "b", col = "red", cex = 0.5) S.f <- S[S$sysname == "Darwin" & S$algorithm == "multiset", c("size", "time_in_secs")] points(sqrt(S.f$size), S.f$time_in_secs, pch = "A", type = "b", cex = 0.5) S.f <- S[S$sysname == "Windows" & S$algorithm == "list", c("size", "time_in_secs")] points(sqrt(S.f$size), S.f$time_in_secs, pch = "W", type = "b", col = "red", cex = 0.5) S.f <- S[S$sysname == "Windows" & S$algorithm == "multiset", c("size", "time_in_secs")] points(sqrt(S.f$size), S.f$time_in_secs, pch = "W", type = "b", cex = 0.5) legend("bottomright", c("C++ STL list", "C++ STL multiset", "-A- MacOSX 10.13.4 / 2.9 GHz Intel Core i7", "-L- Debian deb9.3 / 2.9 GHz Intel Core i7 (VM)", "-W- Windows 10 / 2.9 GHz Intel Core i7 (VM)"), col = c("red", "black", "white", "white", "white"), pch = c(16, 16, 0, 0, 0), cex = 1.0) axis(2) axis(1, c(5, sqrt(50), 10, 15, 20), c("5x5", "US", "10x10", "15x15", "20x20")) box() ################################################################################ plot(1 / time_in_secs ~ sqrt(size), type = "n", log = "y", xlab = "number of map regions", ylab = "number of rectangular cartograms generated per second", data = mbb_check, axes = FALSE, subset = algorithm == "multiset" & sysname == "Darwin") abline(v = sqrt(50), col = "lightgray", lwd = 3) boxplot(1 / time_in_secs ~ sqrt(size), data = mbb_check, cex = 0.75, add = TRUE, axes = FALSE, boxwex = 0.5, subset = algorithm == "multiset" & sysname == "Darwin") axis(2) axis(1, c(5, sqrt(50), 10, 15, 20), c("5x5", "US", "10x10", "15x15", "20x20")) box() legend("topright", c("C++ STL multiset / MacOSX 10.13.4 / 2.9 GHz Intel Core i7"), col = "black", pch = 16, cex = 1.0) ################################################### ### 5. Choose a metaheuristic ################################################### recmap:::.recmap.fitness Checkerboard <- checkerboard(8) summary(Checkerboard) ## 5.1. Greedy randomized adaptive search procedures res.GA <- recmapGA(Checkerboard, popSize = 50, run = 300, maxiter = 300, seed = 3) summary(res.GA$Cartogram) plot(res.GA$Cartogram, col = c("white", "white", "white", "black")[res.GA$Cartogram$z]) op <- par(mar = c(0, 0, 0, 0), mfrow = c(1, 3), bg = "azure") plot(cmp_GA_GRASP$GRASP$Map, border = "black", col = c("white", "white", "white", "black")[cmp_GA_GRASP$GRASP$Map$z]) plot(cmp_GA_GRASP$GRASP$Cartogram, border = "black", col = c("white", "white", "white", "black")[cmp_GA_GRASP$GRASP$Cartogram$z]) plot(cmp_GA_GRASP$GA$Cartogram, border = "black", col = c("white", "white", "white", "black")[cmp_GA_GRASP$GA$Cartogram$z]) op <- par(mar = c(0, 0, 0, 0), mfrow = c(1, 1), bg = "azure") bestsolution <- cmp_GA_GRASP$GA$GA@solution[1, ] Cartogram.Checkerboard <- recmap(Checkerboard[bestsolution, ]) idx <- order(Cartogram.Checkerboard$dfs.num) plot(Cartogram.Checkerboard, border = "black", col = c("white", "white", "white", "black")[Cartogram.Checkerboard$z]) ## draw placement order lines(Cartogram.Checkerboard$x[idx], Cartogram.Checkerboard$y[idx], col = rgb(1, 0, 0, alpha = 0.3), lwd = 4, cex = 0.5) text(Cartogram.Checkerboard$x[idx], Cartogram.Checkerboard$y[idx], 1:length(idx), pos = 1, col = rgb(1, 0, 0, alpha = 0.7)) points(Cartogram.Checkerboard$x[idx[1]], Cartogram.Checkerboard$y[idx[1]], lwd = 5, col = "red") text(Cartogram.Checkerboard$x[idx[1]], Cartogram.Checkerboard$y[idx[1]], "start", col = "red", pos = 3, cex = 0.75) points(Cartogram.Checkerboard$x[idx[length(idx)]], Cartogram.Checkerboard$y[idx[length(idx)]], cex = 1.25, lwd = 2, col = "red", pch = 5) ################################################### ### cmp_GA_GRASP ################################################### op <- par(mar = c(4, 4, 1.5, 0.5), mfrow = c(1, 1), bg = "white") plot(best ~ elapsedtime, data = cmp_GA_GRASP$cmp, type = "n", ylab = "best fitness value", xlab = "elapsed time [in seconds]") abline(v = 60, col = "lightgrey", lwd = 2) lines(cmp_GA_GRASP$cmp[cmp_GA_GRASP$cmp$algorithm == "GRASP", c("elapsedtime", "best")], type = "b", col = "grey", pch = 16) lines(cmp_GA_GRASP$cmp[cmp_GA_GRASP$cmp$algorithm == "GA", c("elapsedtime", "best")], type = "b", pch = 16) legend("bottomright", c("Evolutionary based Genetic Algorithm (GA)", "Greedy Randomized Adaptive Search Procedures (GRASP)"), col = c("black", "grey"), pch = 16, cex = 0.5) ################################################### ### GAseed ################################################### res <- lapply(c(1, 1, 2, 2, 3, 3), function(seed) { set.seed(seed) res <- recmapGA(Map = checkerboard(4), pmutation = 0.25) res$seed <- seed res }) op <- par(mfcol = c(2, 4), bg = "azure", mar = c(5, 0.5, 0.5, 0.5)) x <- recmap(checkerboard(4)) p <- paste(" = (", paste(1:length(x$z), collapse = ", "), ")", sep = "") plot(x, sub = substitute(paste(Pi["forward"], p), list(p = p)), col = c("white", "white", "white", "black")[x$z]) x <- recmap(checkerboard(4)[rev(1:16), ]) p <- paste(" = (", paste(rev(1:length(x$z)), collapse = ", "), ")", sep = "") plot(x, sub = substitute(paste(Pi[reverse], p), list(p = p)), col = c("white", "white", "white", "black")[x$z]) rv <- lapply(res, function(x) { p <- paste(" = (", paste(x$GA@solution[1, ], collapse = ", "), ")", sep = "") plot(x$Cartogram, col = c("white", "white", "white", "black")[x$Cartogram$z], sub = substitute(paste(Pi[seed], perm), list(perm = p, seed = x$seed))) }) fitness.weighted <- function(idxOrder, Map, ...) { Cartogram <- recmap(Map[idxOrder, ]) if (sum(Cartogram$topology.error == 100) > 0) { return(0) } S <- summary(Cartogram) dT <- max(Cartogram$topology.error) dR <- S["relative position error", ] dE <- (100 - S["screen filling [in %]", ]) / 100 ## weighting the objectives 1 / (c(0.2, 0.6, 0.2) %*% c(dT, dR, dE)) } US.map.best <- recmapGA(Map = US.map, fitness = fitness.weighted, maxiter = 100, maxFitness = 100, popSize = 50, keepBest = TRUE, pmutation = 0.35, seed = 1, parallel = FALSE) fitness <- unlist(lapply(US.map.best$GA@bestSol, function(x) { fitness.weighted(idxOrder = x[1, ], Map = US.map) })) op <- par(mfrow = c(1, 2)) plot(US.map.best$GA) GA.generation <- c(1, 6, 12, 25, 50, 100) abline(v = GA.generation, col = "darkred") image(t(sapply(US.map.best$GA@bestSol, function(x) { x[1, ] })), col = grey.colors(50), las = 2, axes = FALSE, xlab = "Generation", ylab = "Index order") axis(1, c(0.2, 0.4, 0.6, 0.8, 1), seq(20, 100, by = 20)) axis(2, c(0.2, 0.4, 0.6, 0.8, 1), seq(10, 50, by = 10)) abline(v = GA.generation / 100, col = "darkred") ################################################### ### bestSolution ################################################### par(mfrow = c(3, 2), mar = c(0, 0, 0, 0)) rv <- lapply(GA.generation, function(idx) { p <- US.map.best$GA@bestSol[[idx]][1, ] fitness <- fitness.weighted(idxOrder = p, Map = US.map) C <- recmap(US.map[p, ]) plot(C, col.text = "darkred") S <- summary(C) legend("bottomright", title = paste("Generation", idx), c(paste("dT =", max(C$topology.error)), paste("dR =", round(S[4, ], 2)), paste("dE =", round(S[5, ])))) }) recmap_state_x77 <- function(input, Map = US.map, DF = state.x77, cm = heat_hcl(10)) { # Join map and data.frame Map <- cbind(Map, DF, match(Map$name, row.names(DF))) attr(Map, "Map.name") <- "U.S." attr(Map, "Map.area") <- input$area # Filter Map <- Map[!Map$name %in% c("Hawaii", "Alaska"), ] # Set attribute for desired area Map$z <- Map[, input$area] res <- recmapGA(Map = Map, popSize = 300, maxiter = 30, run = 10) # Set attribute for the coloring S <- Map[res$GA@solution[1, ], input$color] col.idx <- round((length(cm) - 1) * (S - min(S)) / (max(S) - min(S))) + 1 # Have fun plot(res$Cartogram, col = cm[col.idx], col.text = "black") legend("bottomleft", c(paste("area:", input$area), paste("color:", input$color)), cex = 1.5) res } op <- par(mfrow = c(4, 1), mar = rep(0.25, 4), bg = "white") set.seed(1) recmapGA.x77 <- lapply(list(list(color = "Area", area = "Population"), list(color = "HS Grad", area = "Murder"), list(color = "HS Grad", area = "Income"), list(color = "Life Exp", area = "Illiteracy")), recmap_state_x77) par(op) ################################################### ### x77GA ################################################### op <- par(mar = c(5, 5, 3, 3), mfrow = c(4, 1)) res <- lapply(recmapGA.x77, function(x) { plot(x$GA) }) par(op) ################################################### ### us_counties ################################################### data("counties", package = "noncensus") get_county_mbb <- function(state = "colorado") { MBB <- lapply(map("county", state, plot = FALSE)$names, function(x) { r <- map("county", x, plot = FALSE) dx <- 0.5 * (r$range[2] - r$range[1]) dy <- 0.5 * (r$range[4] - r$range[3]) x <- r$range[1] + dx y <- r$range[3] + dy data.frame(polyname = r$name, x = x, y = y, dx = dx, dy = dy) }) MBB <- do.call("rbind", MBB) MBB <- merge(MBB, county.fips, by = "polyname") MBB$fips <- as.integer(MBB$fips) P <- data.frame(fips = paste(counties$state_fips, counties$county_fips, sep = ""), z = counties$population, name = counties$county_name) P$fips <- as.integer(levels(P$fips))[P$fips] M <- merge(MBB, P, by = "fips") class(M) <- c("recmap", "data.frame") M } op <- par(mfrow = c(5, 2), mar = c(0, 0, 0, 0), bg = "white") recmapGA.county <- lapply(c("california", "colorado", "florida", "new jersey", "new york"), function(state, seed = 1) { M <- get_county_mbb(state) attr(M, "Map.name") <- paste("U.S.", state) attr(M, "Map.area") <- "population" plot(M, col.text = "darkred") M$dx <- 1.25 * M$dx rv <- recmapGA(fitness = function(idxOrder, Map, ...) { Cartogram <- recmap(Map[idxOrder, ]) if (sum(Cartogram$topology.error == 100) > 0) { return(0) } 1 / sum(Cartogram$relpos.error) }, Map = M, popSize = nrow(M) * 5, maxiter = 200, run = 20, seed = seed, pmutation = 0.25) plot(rv$Cartogram, col.text = "darkred") rv }) ################################################### ### us_countiesGA ################################################### op <- par(mar = c(5, 5, 3, 3), mfrow = c(5, 1)) res <- lapply(recmapGA.county, function(x) { plot(x$GA) }) par(op) ################################################### ### shiny ################################################### ## library("shiny") ## recmap_shiny <- system.file("shiny-examples", package = "recmap") ## shiny::runApp(recmap_shiny, display.mode = "normal") ################################################### ### Switzerland_population ################################################### op <- par(mar = c(0, 0, 0, 0), mfrow = c(1, 1)) C <- Switzerland$Cartogram plot(C) idx <- rev(order(C$z))[1:50] text(C$x[idx], C$y[idx], C$name[idx], col = "red", cex = C$dx[idx] / strwidth(as.character(C$name[idx]))) par(op) ################################################### ### SBB ################################################### op <- par(mfrow = c(2, 1), mar = c(0, 0, 0, 0)) SBB.Map <- SBB$Map S <- SBB$Map S <- S[!is.na(S$x), ] S$dx <- 0.1 S$dy <- 0.1 S$z <- S$DTV S$name <- S$Bahnhof_Haltestelle plot.recmap(S) idx <- rev(order(S$z))[1:10] text(S$x[idx], S$y[idx], S$name[idx], col = "red", cex = 0.7) idx <- rev(order(S$z))[11:30] text(S$x[idx], S$y[idx], S$name[idx], col = "red", cex = 0.5) Cartogram.SBB <- recmap(SBB$Map) plot(Cartogram.SBB) idx <- rev(order(Cartogram.SBB$z))[1:40] text(Cartogram.SBB$x[idx], Cartogram.SBB$y[idx], Cartogram.SBB$name[idx], col = "red", cex = 1.25 * Cartogram.SBB$dx[idx] / strwidth(Cartogram.SBB$name[idx])) par(op) ################################################### ### fitnessSBB ################################################### fitness.SBB <- function(idxOrder, Map, ...) { Cartogram <- recmap(Map[idxOrder, ]) if (sum(Cartogram$topology.error == 100) > 1) { return(0) } 1 / sum(Cartogram$z / (sqrt(sum(Cartogram$z^2))) * Cartogram$relpos.error) } ################################################### ### UK_Cartogram ################################################### # construct the cartogram; # note that the index order is already optimized. add_NI <- function(S, DF) { NI.Electorate <- 1260955 NI.Pct_Leave <- 44.22 NI.Pct_Turnout <- 62.69 NI.Pct_Rejected <- 0.05 NI <- data.frame(x = min(S$x) - 3 * max(S$dx), y = min(S$y) + 0.7 * (max(S$y) - min(S$dy)), dx = 76524, dy = 76524, z = 1260955, name = "Northern\nIreland") rownames(NI) <- nrow(S) + 1 rbind(S[, c("x", "y", "dx", "dy", "z", "name")], NI) } UK.recmap <- recmap(UK$Map) # plot the cartogram; 50%-50% regions should have a white color! op <- par(mar = rep(0, 4), mfrow = c(1, 1)) UK_NI.recmap <- add_NI(UK.recmap) stopifnot((UK_NI.recmap[1, "dx"] * UK_NI.recmap[1, "dy"] / UK_NI.recmap[1, "z"]) - (UK_NI.recmap[371, "dx"] * UK_NI.recmap[371, "dy"] / UK_NI.recmap[371, "z"]) < 0.1) plot(UK_NI.recmap, col = diverge_hcl(100, alpha = 1.0)[round(c(UK$Map$Pct_Leave, 44.22))]) # add region names as label; should be readable by zoom into the pdf text(UK_NI.recmap$x, UK_NI.recmap$y, UK_NI.recmap$name, cex = UK_NI.recmap$dx / strwidth(UK_NI.recmap$name)) ################################################### ### UK_Map ################################################### DF <- data.frame(Pct_Leave = UK$Map$Pct_Leave, row.names = UK$Map$name) spplot(as.SpatialPolygonsDataFrame(UK$Map, DF), col.regions = diverge_hcl(16, alpha = 0.5), main = "Input England/Wales/Scottland") ################################################### ### UK_spplot ################################################### DF <- rbind( data.frame(Pct_Leave = UK$Map$Pct_Leave, Pct_Turnout = UK$Map$Pct_Turnout, Pct_Rejected = UK$Map$Pct_Rejected, row.names = UK$Map$name), data.frame(Pct_Leave = 44.22, Pct_Turnout = 62.69, Pct_Rejected = 0.05, row.names = "Northern\nIreland")) UK.sp <- as.SpatialPolygonsDataFrame(add_NI(UK.recmap), DF) summary(UK.sp) spplot(UK.sp, col.regions = diverge_hcl(19)[1:16], layout = c(3, 1)) ################################################### ### Table 2 ################################################### SUMMARY.USX77 <- do.call("rbind", lapply(recmapGA.x77, function(x) { x$Summary })) SUMMARY.USCOUNTY <- do.call("rbind", lapply(recmapGA.county, function(x) { x$Summary })) SUMMARY.CB8 <- recmapGA(checkerboard(8), seed = 1, maxiter = 10)$Summary SUMMARY.US <- US.map.best$Summary UK$Summary$Map.name <- "UK" UK$Summary$Map.area <- "number of electorates" SBB$Summary$Map.name <- "SBB" SBB$Summary$Map.area <- "passagier frequency" Switzerland$Summary$Map.name <- "CH" Switzerland$Summary$Map.area <- "population" SUMMARY <- do.call("rbind", list(SUMMARY.USX77, SUMMARY.USCOUNTY, SUMMARY.CB8, SUMMARY.US, Switzerland$Summary, SBB$Summary, UK$Summary)) xt <- xtable(SUMMARY[order(SUMMARY["Map.number.regions"]), ], caption = paste0("\\label{table:reproduced:overview}The spreadsheet ", "provides a summary of the statistical rectangular cartograms drawn ", "in Figures \\ref{figure:cmp_GA_GRASP}, \\ref{figure:bestSolution}, ", "\\ref{figure:x77}, \\ref{figure:us_census}, ", "\\ref{figure:Switzerland:population}, \\ref{figure:sbb}, and ", "\\ref{fig:brexit} ordered by the number of map regions. All listed ", "rectangular cartograms were processed on Intel(R) Core(TM) i5-2500 ", "CPU @ 3.30GHz having four cores running Debian 9 GNU/Linux.")) print(xt, floating.environment = "sidewaystable", rotate.colnames = TRUE, include.rownames = FALSE) ################################################### ### session info ################################################### toLatex(sessionInfo())