<commit_before>#+setup, include=FALSE, cache=FALSE
library(knitr)
opts_chunk$set(fig.path='simulator-figures/', fig.align='center', fig.show='hold')
options(replace.assign=TRUE,width=80)
Sys.setlocale("LC_TIME", "C") #Needed for identical()
Sys.setlocale("LC_COLLATE", "C")

#' ## Load Packages

set.seed(4253)
library(plyr)
library(maptools)
library(raster)
library(rgdal)
print(sessionInfo())

#' ## Load data

#' County land areas

rd <- read.delim('2013_Gaz_counties_national.txt')
rd <- rd[, c('GEOID', 'ALAND')]

#' Hog farm counts by county

tmpf <- function(){
    data(county.fips, package='maps')
  cens <- read.csv('table-12-hogs-and-pigs-by-county.csv',strip.white=TRUE, na.strings='(D)',
                   stringsAsFactors=FALSE)
# The introduction to the reports says '-' represents 0 and (D) means deleted for privacy
  tmpf <- function(x) {
      if(is.character(x)){
          ret <- ifelse(x=='-', 0, x)
          type.convert(ret)
      }else{
          x
      }
  }
  cens <- lapply(as.list(cens), tmpf)
  cens <- data.frame(cens)
  test <- cens$STCOFIPS %in% county.fips$fips
  cens <- cens[test,]
  ind <- grep('farms, 2007)$', cens$ITEM)
  cens <- cens[ind,]
  cens
}
countyData <- tmpf()

#' ERS farm resource regions

regs <- read.csv('reglink.csv', skip=2,
                 colClasses=c(NA, NA, 'NULL', 'NULL'))

#' #' Get shipment flows data

tmpf <- function(){
    fp <- file.path('shipment-flows-origins-on-rows-dests-on-columns.csv')
    ep <- read.csv(fp, row.names=1)
    ep <- t(ep)
    data.matrix(ep)
}
flowMat <- tmpf()

#' Farms and sales by operation types

tmpf <- function(){
    censPath <- file.path('table26-2002.csv')
    cens <- read.csv(censPath,strip.white=TRUE, na.strings='(D)',
                     stringsAsFactors=FALSE)
    # The introduction to the reports says '-' represents 0 and (D) means deleted for privacy
    tmpf <- function(x) {
        if(is.character(x)){
            ret <- ifelse(x=='-', 0, x)
            type.convert(ret)
        }else{
            x
        }
    }
    cens <- lapply(as.list(cens), tmpf)
    cens <- data.frame(cens)
    test <- cens$GEO %in% state.name
    stateCounts <- cens[test,]

    key <- match(stateCounts$GEO, state.name)
    stateCounts$abb <- state.abb[key]
    stateCounts
}
farmsSales <- tmpf()

#' Farms and inventory

tmpf <- function(){
    censPath <- file.path('table19-2002.csv')
    cens <- read.csv(censPath,strip.white=TRUE, na.strings='(D)',
                     stringsAsFactors=FALSE)
    # The introduction to the reports says '-' represents 0 and (D) means deleted for privacy
    tmpf <- function(x) {
        if(is.character(x)){
            ret <- ifelse(x=='-', 0, x)
            type.convert(ret)
        }else{
            x
        }
    }
    cens <- lapply(as.list(cens), tmpf)
    cens <- data.frame(cens)
    test <- cens$GEO %in% state.name
    stateCounts <- cens[test,]

    key <- match(stateCounts$GEO, state.name)
    rownames(stateCounts) <- state.abb[key]
    stateCounts
}
fiAll <- tmpf()

balanceSheet <- read.csv('state-hogBalanceSheetDec2000Dec2001.csv',
                         na.strings='-', row.names=1, colClasses=c(state2='NULL'))


ea.proj <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs" ## NAD83 Lambert Azimuthal Equal Area

## cb.proj <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs" ## from .proj file
## cb <- readShapeSpatial("cb_2014_us_county_500k/cb_2014_us_county_500k")
##                         proj4string=CRS(cb.proj))
cb <- readOGR('cb_2014_us_county_500k', 'cb_2014_us_county_500k')
cb.ea <- spTransform(cb, ea.proj)

tmpf <- function(geo=cb.ea){
  sf <- countyData$STFIPS
  cf <- countyData$COFIPS
  sf <- formatC(sf, flag="0", format="d", width=2)
  cf <- formatC(cf, flag="0", format="d", width=3)
  test1 <- geo@data[, 'STATEFP'] %in% sf
  test2 <- geo@data[, 'COUNTYFP'] %in% cf
  test <- test1 & test2
  geo[test, ]
}
cb2 <- tmpf()

plotAllCountiesWithFarmData <- function(){
    plot(cb2)
}

getSamp <- function(spdf=cb2, cfps, sfps, n, plot.samp=FALSE){
    if(n == 0) {
        samp <- NULL
    } else {
        cfps <- formatC(cfps, flag="0", format="d", width=3)
        sfps <- formatC(sfps, flag="0", format="d", width=2)
        ind <- which(spdf@data[, 'COUNTYFP']==cfps & spdf@data[, 'STATEFP']==sfps)
        cty.poly <- spdf[ind, ]
        samp <- spsample(cty.poly, type='random', n=n, iter=10)
        if(plot.samp){
            plot(cty.poly)
            points(samp)
        }
    }
    samp
}

testSamplingOfMultiPolygonCounties <- function(){
   getSamp(cfps=21, sfps=25, n=100, plot.samp=TRUE)
}

testSamplingOfCountiesWithHoles <- function(){
   getSamp(cfps=15, sfps=51, n=100, plot.samp=TRUE)
}

coord.samps <- mapply(getSamp, cfps=countyData$COFIPS, sfps=countyData$STFIPS,
                      n=countyData$DATA, SIMPLIFY=FALSE)

r <- raster(cb2)
res(r) <- 1600 * 10
##adj <- adjacent(r, 1:ncell(r), directions=8) ## this could be quickly calculated as needed

tmpf <- function(xy) {
  if(is.null(xy)) {
    NULL
  } else {
    cellFromXY(object=r, xy=xy)
  }
}
cell.samps <- lapply(coord.samps, tmpf)

## plot cell counts
tab <- table(unlist(cell.samps))
r[as.integer(names(tab))] <- tab
plot(r)

## convert fips to abbreviation
stfips.rle <- list(lengths=countyData$DATA, values=countyData$STFIPS)
data(state.fips, package='maps')
key <- match(countyData$STFIPS, state.fips$fips)
countyData$abb <- as.character(state.fips$abb[key])

## create agent data structures

adf <- data.frame(cell=unlist(cell.samps),
                  abb=rep(countyData$abb, times=countyData$DATA),
                  infection.time=NA,
                  recovery.time=NA)
sp.nbs <- vector('list', length=nrow(adf))

## initialize infections

cases <- which(adf$cell %in% c(27799))

get.nbs <- function(cell, rast=r) {
  nb <- adjacent(x=rast, cell, directions=8)
  nb.cells <- nb[, 'to']
  which(adf$cell %in% c(cell, nb.cells))
}
sp.nbs[cases] <- lapply(adf$cell[cases], get.nbs)
adf$infection.time[cases] <- 0

nsteps <- 2
step <- 1
tprob <- 0.01
rprob <- 0.5

while(step < nsteps){
  new.cases <- repeat.cases <- integer(0)
  for (case in cases){
    contacts <- sp.nbs[[case]]
    rand <- runif(n=length(contacts))
    test <- rand < tprob
    if(any(test)){
      new.cases <- c(new.cases, contacts[test])
    }
    rand2 <- runif(n=1)
    if(rand2 > rprob){
      repeat.cases <- c(repeat.cases, case)
    } else {
      adf$recovery.time[case] <- step
    }
  }
  sp.nbs[new.cases] <- lapply(adf$cell[new.cases], get.nbs)
  adf$infection.time[new.cases] <- step
  cases <- c(new.cases, repeat.cases)
  cat('step: ', step, '\n')
  step <- step + 1
}




#'
#' ## Data preparation

#' Derive state-level summaries of county-level data

tmpf <- function() {
    ctyTots <- ddply(countyData, 'STCOFIPS', summarize, totalFarms=sum(DATA), STFIPS=STFIPS[1],
                     smallFarms = sum(DATA[
                         "Inventory \\ Total hogs and pigs \\ Farms by inventory \\ 1 to 24 (farms, 2007)" == ITEM]))
    ctyTots$nonSmallFarms <- with(ctyTots, totalFarms - smallFarms)
    mg <- merge(ctyTots, rd, by.x='STCOFIPS', by.y='GEOID')
    mg <- merge(regs, mg, by.x='Fips', by.y='STCOFIPS')
    mg$kmsq <- mg$ALAND/1e6<commit_msg>Simulation faster and producing multiple time series<commit_after>#+setup, include=FALSE, cache=FALSE
library(knitr)
opts_chunk$set(fig.path='simulator-figures/', fig.align='center', fig.show='hold')
options(replace.assign=TRUE,width=80)
Sys.setlocale("LC_TIME", "C") #Needed for identical()
Sys.setlocale("LC_COLLATE", "C")

#' ## Load Packages

set.seed(4253)
library(plyr)
library(maptools)
library(raster)
library(rgdal)
print(sessionInfo())

#' ## Load data

#' County land areas

rd <- read.delim('2013_Gaz_counties_national.txt')
rd <- rd[, c('GEOID', 'ALAND')]

#' Hog farm counts by county

tmpf <- function(){
    data(county.fips, package='maps')
  cens <- read.csv('table-12-hogs-and-pigs-by-county.csv',strip.white=TRUE, na.strings='(D)',
                   stringsAsFactors=FALSE)
# The introduction to the reports says '-' represents 0 and (D) means deleted for privacy
  tmpf <- function(x) {
      if(is.character(x)){
          ret <- ifelse(x=='-', 0, x)
          type.convert(ret)
      }else{
          x
      }
  }
  cens <- lapply(as.list(cens), tmpf)
  cens <- data.frame(cens)
  test <- cens$STCOFIPS %in% county.fips$fips
  cens <- cens[test,]
  ind <- grep('farms, 2007)$', cens$ITEM)
  cens <- cens[ind,]
  cens
}
countyData <- tmpf()

#' ERS farm resource regions

regs <- read.csv('reglink.csv', skip=2,
                 colClasses=c(NA, NA, 'NULL', 'NULL'))

#' #' Get shipment flows data

tmpf <- function(){
    fp <- file.path('shipment-flows-origins-on-rows-dests-on-columns.csv')
    ep <- read.csv(fp, row.names=1)
    ep <- t(ep)
    data.matrix(ep)
}
flowMat <- tmpf()

#' Farms and sales by operation types

tmpf <- function(){
    censPath <- file.path('table26-2002.csv')
    cens <- read.csv(censPath,strip.white=TRUE, na.strings='(D)',
                     stringsAsFactors=FALSE)
    # The introduction to the reports says '-' represents 0 and (D) means deleted for privacy
    tmpf <- function(x) {
        if(is.character(x)){
            ret <- ifelse(x=='-', 0, x)
            type.convert(ret)
        }else{
            x
        }
    }
    cens <- lapply(as.list(cens), tmpf)
    cens <- data.frame(cens)
    test <- cens$GEO %in% state.name
    stateCounts <- cens[test,]

    key <- match(stateCounts$GEO, state.name)
    stateCounts$abb <- state.abb[key]
    stateCounts
}
farmsSales <- tmpf()

#' Farms and inventory

tmpf <- function(){
    censPath <- file.path('table19-2002.csv')
    cens <- read.csv(censPath,strip.white=TRUE, na.strings='(D)',
                     stringsAsFactors=FALSE)
    # The introduction to the reports says '-' represents 0 and (D) means deleted for privacy
    tmpf <- function(x) {
        if(is.character(x)){
            ret <- ifelse(x=='-', 0, x)
            type.convert(ret)
        }else{
            x
        }
    }
    cens <- lapply(as.list(cens), tmpf)
    cens <- data.frame(cens)
    test <- cens$GEO %in% state.name
    stateCounts <- cens[test,]

    key <- match(stateCounts$GEO, state.name)
    rownames(stateCounts) <- state.abb[key]
    stateCounts
}
fiAll <- tmpf()

balanceSheet <- read.csv('state-hogBalanceSheetDec2000Dec2001.csv',
                         na.strings='-', row.names=1, colClasses=c(state2='NULL'))


ea.proj <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs" ## NAD83 Lambert Azimuthal Equal Area

## cb.proj <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs" ## from .proj file
## cb <- readShapeSpatial("cb_2014_us_county_500k/cb_2014_us_county_500k")
##                         proj4string=CRS(cb.proj))
cb <- readOGR('cb_2014_us_county_500k', 'cb_2014_us_county_500k')
cb.ea <- spTransform(cb, ea.proj)

tmpf <- function(geo=cb.ea){
  sf <- countyData$STFIPS
  cf <- countyData$COFIPS
  sf <- formatC(sf, flag="0", format="d", width=2)
  cf <- formatC(cf, flag="0", format="d", width=3)
  test1 <- geo@data[, 'STATEFP'] %in% sf
  test2 <- geo@data[, 'COUNTYFP'] %in% cf
  test <- test1 & test2
  geo[test, ]
}
cb2 <- tmpf()

plotAllCountiesWithFarmData <- function(){
    plot(cb2)
}

getSamp <- function(spdf=cb2, cfps, sfps, n, plot.samp=FALSE){
    if(n == 0) {
        samp <- NULL
    } else {
        cfps <- formatC(cfps, flag="0", format="d", width=3)
        sfps <- formatC(sfps, flag="0", format="d", width=2)
        ind <- which(spdf@data[, 'COUNTYFP']==cfps & spdf@data[, 'STATEFP']==sfps)
        cty.poly <- spdf[ind, ]
        samp <- spsample(cty.poly, type='random', n=n, iter=10)
        if(plot.samp){
            plot(cty.poly)
            points(samp)
        }
    }
    samp
}

testSamplingOfMultiPolygonCounties <- function(){
   getSamp(cfps=21, sfps=25, n=100, plot.samp=TRUE)
}

testSamplingOfCountiesWithHoles <- function(){
   getSamp(cfps=15, sfps=51, n=100, plot.samp=TRUE)
}

coord.samps <- mapply(getSamp, cfps=countyData$COFIPS, sfps=countyData$STFIPS,
                      n=countyData$DATA, SIMPLIFY=FALSE)

r <- raster(cb2)
res(r) <- 1600 * 10
#adj <- adjacent(r, 1:ncell(r), directions=8) ## this could be quickly calculated as needed



tmpf <- function(xy) {
  if(is.null(xy)) {
    NULL
  } else {
    cellFromXY(object=r, xy=xy)
  }
}
cell.samps <- lapply(coord.samps, tmpf)

## plot cell counts
tab <- table(unlist(cell.samps))
r[as.integer(names(tab))] <- tab
plot(r)

## convert fips to abbreviation
stfips.rle <- list(lengths=countyData$DATA, values=countyData$STFIPS)
data(state.fips, package='maps')
key <- match(countyData$STFIPS, state.fips$fips)
countyData$abb <- as.character(state.fips$abb[key])

## create agent data structures

adf <- data.frame(cell=unlist(cell.samps),
                  abb=rep(countyData$abb, times=countyData$DATA),
                  infection.time=NA,
                  recovery.time=NA)
sp.nbs <- vector('list', length=nrow(adf))

occupied.cells <- unique(adf$cell)
cell2id <- lapply(occupied.cells, function(x) which(adf$cell == x))
names(cell2id) <- occupied.cells

## initialize infections

                                        #cases <- which(adf$cell %in% c(27799))
adf$infection.time <- NA
adf$recovery.time <- NA
cases <- sample.int(nrow(adf), 10)

get.nbs <- function(cell, rast=r) {
  nb <- adjacent(x=rast, cell, directions=8)
  nb.cells <- nb[, 'to']
  cell.names <- as.character(c(cell, nb.cells))
  unlist(cell2id[cell.names])
}
sp.nbs[cases] <- lapply(adf$cell[cases], get.nbs)
adf$infection.time[cases] <- 0

nsteps <- 20
step <- 1
tprob <- 0.01
rprob <- 0.5

while(step < nsteps){
  new.cases <- repeat.cases <- integer(0)
  for (case in cases){
    contacts <- sp.nbs[[case]]
    is.susceptible <- is.na(adf$infection.time[contacts])
    contacts <- contacts[is.susceptible]
    rand <- runif(n=length(contacts))
    test <- rand < tprob
    if(any(test)){
      new.cases <- c(new.cases, contacts[test])
    }
    rand2 <- runif(n=1)
    if(rand2 > rprob){
      repeat.cases <- c(repeat.cases, case)
    } else {
      adf$recovery.time[case] <- step
    }
  }
  print(system.time(sp.nbs[new.cases] <- lapply(adf$cell[new.cases], get.nbs)))
  adf$infection.time[new.cases] <- step
  cases <- c(new.cases, repeat.cases)
  cat('step: ', step, '\n')
  step <- step + 1
}

## tabulation of case counts

step <- seq(1, to=nsteps)

events.by.state <- function(x, what) {
  tapply(adf[[what]] < x, adf$abb, sum, na.rm=TRUE)
}
cum.infections <- sapply(step, events.by.state, what='infection.time')
cum.recoveries <- sapply(step, events.by.state, what='recovery.time')
no.infected <- cum.infections - cum.recoveries

#'
#' ## Data preparation

#' Derive state-level summaries of county-level data

tmpf <- function() {
    ctyTots <- ddply(countyData, 'STCOFIPS', summarize, totalFarms=sum(DATA), STFIPS=STFIPS[1],
                     smallFarms = sum(DATA[
                         "Inventory \\ Total hogs and pigs \\ Farms by inventory \\ 1 to 24 (farms, 2007)" == ITEM]))
    ctyTots$nonSmallFarms <- with(ctyTots, totalFarms - smallFarms)
    mg <- merge(ctyTots, rd, by.x='STCOFIPS', by.y='GEOID')
    mg <- merge(regs, mg, by.x='Fips', by.y='STCOFIPS')
    mg$kmsq <- mg$ALAND/1e6
{'fraction_non_alphanumeric': 0.17269940650883714, 'fraction_numerical': 0.023413552468531926, 'mean_word_length': 2.994269340974212, 'pattern_counts': {'":': 0, '<': 242, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1824004', 'n_tokens_mistral': 6305, 'n_tokens_neox': 5702, 'n_words': 1135}