## Read in the data table:
Z=read.table("AnnArbor.txt",sep=",", row.names=NULL, header=T)

## Create a sequence of the lower edge of the temp intervals
range.bottom = seq(from=-20, to=70, by = 10)

mean.vec = NULL
var.vec = NULL
n.obs.vec = NULL

for ( lower in range.bottom )
{
  upper = lower + 9

  ## get row indices corresponding to days
  ## with minimum temperatures in the range [lower, upper]
  indices = which( (Z[,4] >= lower) & (Z[,4] <= upper) )
  
  ## keep track of the number of observations in this interval
  n.obs.vec = c(n.obs.vec, length(indices) )

  ## compute the mean, and variance of the max-seal-level-pressure
  ## then add it to the end of the respective vectors.
  mean.vec = c(mean.vec, mean(Z[indices, 11]) )
  var.vec = c(var.vec, var(Z[indices, 11]))
}

## we now have mean.vec and var.vec, 
## we can look at them
mean.vec
var.vec
n.obs.vec



## Infection/Outbreak

people = array(2, 100)
people[1] = 1

infected = which(people == 1)

p.death = 0.2
for (j in infected) 
{
  if (runif(1) < p.death) { people[j] = 3 }
}

susceptible = which(people == 2)
to.be.infected = sample(susceptible)[1]
people[to.be.infected] = 1

#####################

nrep = 1e3
p.recover = 0.1
p.death = 0.2
p.infect = 0.2

num.died = NULL

for (k in 1:nrep) 
{
  ## Initialize the population.
  
  ## 1 = Infected
  ## 2 = Susceptible 
  ## 3 = Dead
 
  people = array(2, 100)
  people[1] = 1
 
  ## Count the weeks.
  n = 0

  ## Iterate until nobody is infected.
  while (any(people==1)) 
  {
    n = n+1
   
    ## Figure out which infected people die
    infected = which(people == 1)
    for (j in infected) 
    {
      if (runif(1) < p.death) { people[j] = 3 }
    }
       
    ## Figure out which remaining infected people infect others, 
    ## and if they recover
    infected = which(people == 1)
    for (j in infected) 
    {
      if (runif(1) < p.infect) 
      {
        susceptible = which(people == 2)
        if ( sum(people == 2) > 0 )
        {
          to.be.infected = sample(susceptible)[1]
          people[to.be.infected] = 1
        }
      }
      if (runif(1) < p.recover) { people[j] = 2 }
    }             
  }
  
  ## find the number of people that died during this outbreak
  num.died[k] = sum(people == 3)       
}