SMC fails with single observation counties

[christopherkenny]

This is a MWE of a bug where SMC fails when a county only has a single precinct.

library(redist)
#> 
#> Attaching package: 'redist'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(tidyverse)

bb <- st_sfc(st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0)))))
grid <- st_make_grid(bb, n = 10)
tb <- st_as_sf(grid)
tb <- tb %>% rename(geometry = x)
tb <- tb %>% mutate(pop = 1, counties = row_number())
box_map <- redist_map(tb, ndists = 4, pop_tol = 0.1, total_pop = pop)
#> Warning in redist_map(tb, ndists = 4, pop_tol = 0.1, total_pop = pop): Missing
#> CRS, assuming NAD83 (4269).
#> Projecting to CRS 3857

set.seed(1)
test <- redist_smc(box_map, 10, counties = counties)
#> SEQUENTIAL MONTE CARLO
#> Sampling 10 100-unit maps with 4 districts and population between 22.5 and 27.5.
#> Ensuring no more than 3 splits of the 100 administrative units.
#> Making split 1 of 3
#> Error in smc_plans(nsims, adj, counties, pop, ndists, pop_bounds[2], pop_bounds[1], : vector::_M_range_check: __n (which is 0) >= this->size() (which is 0)

^{Created on 2021-06-13 by the reprex package (v2.0.0)}

[christopherkenny]

Temporary workaround:

add_ghost_precincts <- function(map, counties){
  cty <- rlang::eval_tidy(rlang::enquo(counties), map)
  
  map_sub <- map %>% 
    slice(which(cty %in% names(which(table(cty) == 1)))) %>% 
    select(-.data$adj) %>% as.data.frame() %>% st_as_sf() %>% st_as_sf()

  if(nrow(map_sub) == 0){
    return(map)
  }
  
  map_sub <- st_centroid(map_sub)

  map_sub <- map_sub %>% mutate(across(starts_with(c('pop', 'vap', 'cvap')), function(x){0}))
  cty_sub <- rlang::eval_tidy(rlang::enquo(counties), map_sub)

  adj <- get_adj(map)
  adj_sub <- vector(mode = 'list', length = nrow(map_sub))
  adj_new <- c(adj, adj_sub)
  
  
  for(i in 1:nrow(map_sub)){
    adj_new <- geomander::add_edge(adj_new, i + nrow(map), which(cty_sub[i] == cty))
  }
  
  return(redist_map(bind_rows(map, map_sub), 
                    pop_tol = get_pop_tol(map), 
                    adj = adj_new, ndists = attr(map, 'ndists'),
                    total_pop = attr(map, 'pop_col')))
}

This way, you can run following the above example:

box_map_ghost <- box_map %>% add_ghost_precincts(counties)