SPAM Model

engine.r

# Run model simulation with given values for: 
# initial cull, maintenance cull, cell target density, allocated budget,
# area to cull, and "control area" (area for which proportional population is calculated)
# (last two are kept separate for the spatial optimisation, where we may choose to cull a subset
# of the target control area for efficiency)
run.sim = function(curr.init.cull, curr.maint.cull, curr.target.density, budget, C.mask, ctrl.mask)
{
  Nx = ncol(K); Ny = nrow(K); # gauge size of map
  IC = K # note - initial population is set to carrying capacity

  cull.cells = sum(C.mask > 0, na.rm=T) # number of cells to cull
  disp.av = (disp.max + disp.min) / 2 # values for STAR's probability of dispersal algorithm
  disp.range = disp.max - disp.min

  N.years = duration
  N.seasons = 2 # note - not yet allowing this to change
  pop = IC # pop is the current population map, starting with initial condition
  pop.store = vector('list', 1 + N.years * N.seasons) # create a list to store the entire population time series
  pop.store[[1]] = IC # first map is initial condition
  
  total.pop = rep(NA, 1 + N.years * N.seasons) # vector to store total population size

  total.pop[1] = sum(pop, na.rm=T) # our version - count ALL animals
#  total.pop[1] = sum(pop * (R==1), na.rm=T) # STAR version - only count animals within national park (R is generated from "regions.csv", see code at end)


  total.cost = rep(NA, N.years * N.seasons) # create variables for management stuff - we don't include the initial condition as we do for pop size
  P.N0 = rep(NA, N.years * N.seasons)
  P.N0.spec = rep(NA, N.years * N.seasons)
  total.cull = rep(NA, N.years * N.seasons)

  for (t in 1:N.years)
  {
    for (t2 in 1:N.seasons) # cycle over time
    {
      P = disp.av + disp.range / (0.9 * K) * (pop - 0.55 * K) # calculate probability of dispersal for each cell (using STAR algorithm)

      E = as.matrix(pop * P / 8 * (1 - 0.5 * EL)) # emigration rate for each cell
      Epad = matrix(0, Ny+2, Nx+2) # pad matrix out with a border of zeroes (for later calculations)
      Epad[1 + (1:Ny), 1 + (1:Nx)] = E
      Epad[is.na(Epad)] = 0 # note - this method means that dispersal rate is less in cells on the border - this may not necessarily be the case

      I = matrix(0, Ny, Nx) # immigration rate
      for (i in 0:2) # calculate immigration rate based on all surrounding cells (note - mathematically questionable to include diagonal cells for each timestep)
      {
        for (j in 0:2)
        {
          if (!(i==1 & j==1)) # except itself
            I = I + Epad[i + (1:Ny), j + (1:Nx)] 
        }
      }

#      I = I * (R==1) # STAR - only allow immigration into cells inside the park

      if (t2 == 1) # if we're in the right season for culling (note - culling only in "dry" season, not applicable for other cases)
      {
        if (t == 1) # if in first timestep, do initial cull, otherwise do maintenance cull
          C = pop * curr.init.cull
        else
          C = pop * curr.maint.cull
        C[pop <= (curr.target.density * cell.size)] = 0 # don't cull any cells that are already at/below the target density
#        C[20,6] = 0 # STAR - arbitrarily removed cell (error)
        C = C * C.mask # mask out unculled areas
      }
      else
        C = pop * 0                                                 

      CC = round(((cost.int*(pop/cell.size)^cost.slope)*C*helicopter.cost)*overhead,0) # work out cost of culling each cell
      CB = ((CC+1)/((budget/duration)/(cull.cells)))/PR  # work out cost-benefit value for each cell (note - based on *allocated budget*, not actual budget which is not yet known)
      
      if (t2 == 1 & t == 1)  # add cost benefit scores to total.CB
        total.CB = CB
      else
        total.CB = total.CB + CB

      pop = pop * exp(r.m * (1 - (pop / K)^theta)) - C - (E*8 - I) # update current population based on growth, culling, emigration and immigration

      pop[pop < (min.d * cell.size)] =  (min.d * cell.size) # any cells that go below a minimum threshold density are held at that threshold

      pop.store[[t * N.seasons + t2 - 1]] = pop # store population

      total.pop[t * N.seasons + t2 - 1] = sum(pop, na.rm=T) # calculate and store current total population
#      total.pop[t * N.seasons + t2 - 1] = sum(pop * (R==1), na.rm=T) # STAR version - store population within park only

      total.cost[t * N.seasons + t2 - 2] = sum(CC, na.rm=T) # calculate current total cost and number of animals culled
      total.cull[t * N.seasons + t2 - 2] = sum(C, na.rm=T)

      # STAR version - wrong cells selected for cost *AND* cull
#      total.cost[t * N.seasons + t2 - 2] = sum(CC[-22, -c(1,14)], na.rm=T) - sum(quantile(CC, c(0.4,0.6,0.8), na.rm=T))
#      total.cull[t * N.seasons + t2 - 2] = sum(C, na.rm=T) - sum(quantile(C, c(0,0.25,0.75), na.rm=T))

      # work out P.N0 (proportion of initial population) in control region only
      P.N0[t * N.seasons + t2 - 2] = sum(pop[ctrl.mask == 1], na.rm=T) / sum(K[ctrl.mask == 1], na.rm=T)
#      P.N0[t * N.seasons + t2 - 2] = mean(as.double(as.matrix(pop/K)[ctrl.mask == 1]), na.rm=T) # STAR version - calculates a "mean proportion by cell" rather than an overall proportion

    }
  }
  list(total.cost = total.cost, total.pop = total.pop, pop.ts = pop.store, P.N0 = P.N0, cost.benefit = total.CB, total.cull = total.cull)
}             


# global constant variables for use in optimisation
td = c(0.75, 0.5, 0.25, 0.1, 0.01) # different levels of target density to try
probs = c(0.01, 0.1, 0.25, 0.5, 0.75, 1) # different quantiles for spatial budget optimisation


# non-spatial budget (based on but using a different method to STAR version - didn't bother to replicate theirs entirely)
# note - didn't bother to put in progress bar
NS.B = function(budget, C.mask)
{
  # work out what happens if we ramp up culling for each target density level
  top = rep(NA, 5)
  for (i in 1:5) {res=run.sim(0.99, 0.99, td[i], budget, C.mask, C.mask); top[i] = sum(res$total.cost)}
  # alternatively, work out what happens if we do the bare minimum
  bottom = rep(NA, 5)
  for (i in 1:5) {res=run.sim(0.01, 0.01, td[i], budget, C.mask, C.mask); bottom[i] = sum(res$total.cost)}
  
  easy = which(top < budget)   # work out which target density levels still satisfies the budget even at the highest cull rate
  hard = which(bottom < budget) # work out which target density levels are achievable within budget doing the bare minimum
  
  if (length(easy) > 0) # if there's at least one case where we can get a 99% cull within budget, we're done
  {
    m = max(easy) # pick the hardest case where we can do this (smallest target density)
#    m = min(easy) # easiest case (even though we could still go under budget with more)  - STAR version
    cull = 0.99
  }
  else if (length(hard) == 0) # if we can't go under budget no matter what we do
  {
    m = 5
    cull = NA
  }
  else # if we can go under budget, but can't do it at the maximum 99% cull rate
  {
    m = min(which(top > budget & bottom < budget)) # look for easiest case (largest cell target density) where we can cull within budget - this is to maximise the cull rate we can get (this should always be m=1, or 0.75 density)
    sim = function(x) sum(run.sim(x, x, td[m], budget, C.mask, C.mask)$total.cost) - budget # make a function to tell us how far we are over budget for a particular cull rate x
    U = uniroot(sim, c(0.01, 0.99)) # find where we are pretty much exactly on budget between the min and max cull rates
    cull = floor(U$root * 100)/100 # round down the cull rate so that we are strictly within budget (note - this won't work if the required budget is INCREASING with decreasing cull rate, shouldn't be an issue though)
  }
  
  list(target.density = td[m], cull.rate = cull) # return the target density and cull rate
}


# non-spatial control density (very similar to non-spatial budget, but instead looking for 1) *minimum* cull rate and 2) *maximum* target density)
# note - didn't bother to put in progress bar
NS.D = function(pn0, C.mask)
{
  top = rep(NA, 5)
  for (i in 1:5) {res=run.sim(0.99, 0.99, td[i], budget, C.mask, C.mask); top[i] = res$P.N0[length(res$P.N0) - 1]}
  bottom = rep(NA, 5)
  for (i in 1:5) {res=run.sim(0.01, 0.01, td[i], budget, C.mask, C.mask); bottom[i] = res$P.N0[length(res$P.N0) - 1]}
  
  hard = which(top < pn0)
  easy = which(bottom < pn0)
  
  if (length(easy) > 0)
  {
    m = min(easy) # easiest case
    cull = 0.99
  }
  else if (length(hard) == 0)
  {
    m = 1
    cull = NA
  }
  else # look for easiest case (lowest target density -> lowest cull rate required) where we're under target
  {
    m = max(which(top < pn0 & bottom > pn0)) # should pretty much always be m=5 where it exists
    sim = function(x) sum(run.sim(x, x, td[m], budget, C.mask, C.mask)$P.N0[length(res$P.N0) - 1]) - pn0
    U = uniroot(sim, c(0.01, 0.99))
    cull = ceiling(U$root * 100)/100
  }
  
  list(target.density = td[m], cull.rate = cull)
}



#spatial budget
S.B = function(budget, ctrl.mask)
{
  pb <- winProgressBar("0% done", "0% done", 0, 100, 0) # draw progress bar
  
  #start with 50% cull rate + target
  best = run.sim(0.5, 0.5, 0.5, budget, ctrl.mask, ctrl.mask) 
  # work out cost-benefit percentiles
  CB = best$cost.benefit  
  Q = quantile(best$cost.benefit, probs = probs , na.rm=T)

  for (t in 1:5) # iterate five times (note - no measure of convergence used in STAR, may not converge properly or waste time when already converged)
  {
#    print(paste('Iteration:', t)) # debugging output

    # initialise "best" variables
    best.run = 0 
    best.dens = 1 # set to a value that will be beatable by anything useful
    best.cullrate = 0
    best.td = 0
    
    for (i in 1:6) # run models culling for best X% of full culling cost-benefit map
    {
#      print(sprintf('Test case: %d%%', probs[i]*100)) # debugging output
      W = (best$cost.benefit <= Q[i]) + 0 # set W to be the map
      test = NS.B(budget, W) # run non-spatial budget optimiser on the map
      out = run.sim(test$cull.rate, test$cull.rate, test$target.density, budget, W, ctrl.mask) # run a sim of the result of the optimisation
      dred = out$P.N0[duration * 2 - 1] # store the current density reduction to compare to the best

#      cat(sprintf('Results:\nDensity = %.2f%%\nCull rate = %d%%\nTarget density = %.2f\n\n',  # debugging output
#        dred*100, round(test$cull.rate*100), test$target.density))

      if (dred < best.dens) # if it's the best, store its information
      {
        best.run = i
        best.dens = dred
        best.cullrate = test$cull.rate
        best.td = test$target.density
      }

      pc.done = round((6*(t-1) + i)*(10/3)) # calculate how far through we are
      info = sprintf("%d%% done", pc.done) # show current progress on progress bar
      setWinProgressBar(pb, pc.done, info, info)      
    }
    
    if (t < 5) # if we're not finished yet, run the best of the six quantiles with full culling to get a new cost-benefit map under those conditions
    {
      best = run.sim(best.cullrate, best.cullrate, best.td, budget, ctrl.mask, ctrl.mask)
      Q = quantile(best$cost.benefit, probs = probs , na.rm=T)      
    }
    else # otherwise work out the *final* culling map and simulation results based on the best run
    {
      W = (best$cost.benefit <= Q[best.run]) + 0
      final = run.sim(best.cullrate, best.cullrate, best.td, budget, W, ctrl.mask)
    }
    
#    print(sprintf('Best run: %d%%', probs[best.run]*100))  # debugging output

    CB = best$cost.benefit  # store cost-benefit map for next run
  }
  
  close(pb) # close progress bar
  list(test = list(target.density = best.td, cull.rate = best.cullrate), out = final, culling.area = W) # return best optimisation results, best simulation, and culling map
}










#stop('Functions loaded')


# Example and messy code used to test STAR

#source("C:/Users/nick/Desktop/alps stuff/newstar.r")
#setwd('C:/Users/nick/Desktop/alps stuff/')

#IC = read.csv('ic.csv', header=FALSE)
#image(t(as.matrix(IC))[,22:1], asp=1)

## STAR maps
#K = read.csv('cc.csv', header=FALSE)
#EL = read.csv('elevation.csv', header=FALSE)
#R = read.csv('regions.csv', header=FALSE)
#PR = read.csv('priority.csv', header=FALSE)

## Alps maps (convert csv file into raster)
#K = read.csv('horse.k.csv', header=FALSE)
#min.east = 365000
#min.north = 5795000
#test = raster(as.matrix(K), min.east - 5e3, min.east + 34e4 - 5e3, min.north - 5e3, min.north + 30e4 - 5e3, crs = "+proj=utm +zone=55 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
#writeRaster(test, filename = 'K.asc')
#EL = K * 0
#R = K * 0 + 1
#PR = R


#park.mask = (R == 1) + 0
#park.mask[is.na(park.mask)] = 0 
#res = run.sim(0.4, 0.2, 0.25, 4.8e6, park.mask, park.mask)
#sum(res$total.cull)
#
#C.mask = (R > 0) + 0
#C.mask[is.na(C.max.mask)] = 0 
#out = S.B(4.8e6, C.mask)
#  
#
#
#
#
#
#ii = seq(0, 1, 0.01)
#tc = matrix(NA, 101, 5)
#pn0 = matrix(NA, 101, 5)
#
#for (j in 1:5)
#{
#  print(j)
#  for (i in 1:101) 
#  {
#    if (i%%10 == 0) print(i)
#    res = run.sim(ii[i], ii[i], td[j])
#    tc[i,j] = sum(res$total.cost)
#    pn0[i,j] = res$P.N0[duration * 2 - 1]    
#  }
#}
#
## Plot for budget
#par(lwd=2)
#palette(c('red', 'orange', 'yellow', 'green', 'blue'))
#plot(0:100, tc[,1]/1e6, xlab='Percentage culled', ylab='Cost ($mil)', type='o', 
#col=1, pch=19, ylim=c(1e5, max(tc))/1e6, log='y', axes=F)
#for (j in 2:5)
#  lines(0:100, tc[,j]/1e6, type='o', col=j, pch=19)
#abline(h=20, v=47, lwd=2, lty=2)
#axis(1)
#axis(2, at=c(0.1, 0.5, 1, 2, 5, 10, 20, 50, 100, 200))
#box(bty='L')
#legend(x = 80, y=2, legend=td, col=1:5, lwd=2, pch=19, title = 'Target density')
#
## Plot for relative density
#par(lwd=2)
#palette(c('red', 'orange', 'yellow', 'green', 'blue'))
#plot(0:100, pn0[,1], xlab='Percentage culled', ylab='Control density', type='o', 
#col=1, pch=19, ylim=c(0, 0.1))
#for (j in 2:5)
#  lines(0:100, pn0[,j], type='o', col=j, pch=19)
#abline(h=c(0.05, 0.055), v=91, lwd=2, lty=2)
#box(bty='L')
#legend(x = 20, y=0.04, legend=td, col=1:5, lwd=2, pch=19, title = 'Target density')

gui.r

fileChoose <- function(action="print", text = "Select a file...", type="open", ...) {
   gfile(text=text, type=type, ..., action = action, handler =function(h,...) {
     if(h$action != ""){
       print(h$action)
       do.call(h$action, list(h$file))
     }
   })}


fileSaveChoose <- function(action="print", text = "Select a file...", type="save", ...) {
  gfile(text=text, type=type, ..., action = action, handler =function(h,...) {
    if(h$action != ""){
      print(h$action)
      do.call(h$action, list(h$file))
   }
})
}

## Simple confirmation dialog
confirmDialog <- function(message, handler=NULL) {
   window <- gwindow("Confirm")
   group <- ggroup(container = window)
   gimage("info", dirname="stock", size="dialog", container=group)
  
   ## A group for the message and buttons
   inner.group <- ggroup(horizontal=FALSE, container = group)
   glabel(message, container=inner.group, expand=TRUE)
  
   ## A group to organize the buttons
   button.group <- ggroup(container = inner.group)
   ## Push buttons to right
   addSpring(button.group)
   gbutton("OK", handler = function(h,...) dispose(window), container=button.group)
   return()
}


# Define top menu list
spam.menu <- list()
spam.menu$File$Exit$handler = function(h,...) gtkMainQuit()


# Hierarchial GUI definition
window <- gwindow("SPAM Model", visible=T, width=1200, height=900, handler=function(h,...){gtkMainQuit()})
  add.menu <- gmenu(spam.menu,container=window)
  left.pane <- ggroup(container=window,horizontal=F)
    input.pad <- gnotebook(container=left.pane,expand=T)
      file.group <- ggroup(container=input.pad,label="Rasters",horizontal=F)
        K.label <- glabel("",container=file.group)
        K.button <- gbutton("Load Carrying Capacity Raster", container=file.group,
                             handler=function(h,...){
                               K.file=fileChoose("print")
                               if(is.na(K.file)) return
                               K.rast <<- raster(K.file)
                               K<<- as.matrix(K.rast)
                               svalue(K.label)=K.file
                               visible(K.plot.raster)=T
                               plot(K.rast)
                             })
        EL.label <- glabel("",container=file.group)
        EL.button <- gbutton("Load Elevation Raster", container=file.group,
                    handler=function(h,...){
                      EL.file=fileChoose("print")
                      if(is.na(EL.file)) return
                      EL.rast <<- raster(EL.file)
                      EL <<- as.matrix(EL.rast)
                      svalue(EL.label)=EL.file
                      visible(EL.plot.raster)=T
                      plot(EL.rast)
                    })
        PR.label <- glabel("",container=file.group)
        PR.button <- gbutton("Load Priority Raster", container=file.group,
                     handler=function(h,...){
                       PR.file=fileChoose("print")
                       if(is.na(PR.file)) return
                       PR.rast <<- raster(PR.file)
                       PR<<- as.matrix(PR.rast)
                       svalue(PR.label)=PR.file
                       visible(PR.plot.raster)=T
                       plot(PR.rast)
                     })        
        C.label <- glabel("",container=file.group)
        C.button <- gbutton("Load Culling Mask Raster", container=file.group,
                     handler=function(h,...){
                       C.file=fileChoose("print")
                       if(is.na(C.file)) return
                       C.mask.rast <<- raster(C.file)
                       C.mask<<- as.matrix(C.mask.rast)
                       svalue(C.label)=C.file
                       visible(C.plot.raster)=T
                       plot(C.mask.rast)
                     })     
        raster.pad <- gnotebook(container=file.group,expand=T)
          K.plot.raster <- ggraphics(container=raster.pad,label="Carrying Capacity")
          EL.plot.raster <- ggraphics(container=raster.pad,label="Elevation")
          PR.plot.raster <- ggraphics(container=raster.pad,label="Priority")
          C.plot.raster <- ggraphics(container=raster.pad,label="Culling Mask")

      var.group <- ggroup(container=input.pad,label="Parameters",horizontal=F)

        budget.group=ggroup(container=var.group)
        budget.label=glabel("Maximum available budget (used for optimisation)",container=budget.group)
        addSpring(budget.group)
        budget.slider=gspinbutton(from=0,to=20000000,by=500000,value=budget,container=budget.group,
                              handler=function(h,...){
                                budget <<-svalue(budget.slider)
                              })
        area.target.group=ggroup(container=var.group)
        area.target.label=glabel("Target proportional population in region (used for optimisation)",container=area.target.group)
        addSpring(area.target.group)
        area.target.slider=gspinbutton(from=0,to=1,by=.05,value=area.target,container=area.target.group,
                          handler=function(h,...){
                            area.target <<-svalue(area.target.slider)
                          })
        elev.d.mod.group=ggroup(container=var.group)
        elev.d.mod.label=glabel("Dispersal rate in high elevation area (where elevation = 1 in raster)",container=elev.d.mod.group)
        addSpring(elev.d.mod.group)
        elev.d.mod.slider=gspinbutton(from=0,to=1,by=.05,value=elev.d.mod,container=elev.d.mod.group,
                               handler=function(h,...){
                                 elev.d.mod <<-svalue(elev.d.mod.slider)
                          })
        init.cull.group=ggroup(container=var.group)
        init.cull.label=glabel("Culling rate in first season",container=init.cull.group)
        addSpring(init.cull.group)
        init.cull.slider=gspinbutton(from=0,to=1,by=.05,value=init.cull,container=init.cull.group,
                              handler=function(h,...){
                                init.cull <<-svalue(init.cull.slider)
                          })
        maint.cull.group=ggroup(container=var.group)
        maint.cull.label=glabel("Culling rate after first season",container=maint.cull.group)
        addSpring(maint.cull.group)
        maint.cull.slider=gspinbutton(from=0,to=1,by=.05,value=maint.cull,container=maint.cull.group,
                             handler=function(h,...){
                               maint.cull <<-svalue(maint.cull.slider)
                          })
        target.density.group=ggroup(container=var.group)
        target.density.label=glabel("Species density at which culling stops in a cell",container=target.density.group)
        addSpring(target.density.group)
        target.density.slider=gspinbutton(from=0,to=1,by=.05,value=target.density,container=target.density.group,
                              handler=function(h,...){
                                target.density <<-svalue(target.density.slider)
                          })
        cell.size.group=ggroup(container=var.group)
        cell.size.label=glabel("Cell size (km^2)",container=cell.size.group)
        addSpring(cell.size.group)
        cell.size.slider=gspinbutton(from=0,to=1,by=.05,value=cell.size,container=cell.size.group,
                                  handler=function(h,...){
                                    cell.size <<-svalue(cell.size.slider)
                          })
        r.m.group=ggroup(container=var.group)
        r.m.label=glabel("Growth rate of species",container=r.m.group)
        addSpring(r.m.group)
        r.m.slider=gspinbutton(from=0,to=1,by=.05,value=r.m,container=r.m.group,
                             handler=function(h,...){
                               r.m <<-svalue(r.m.slider)
                          })
        theta.group=ggroup(container=var.group)
        theta.label=glabel("Shape parameter",container=theta.group)
        addSpring(theta.group)
        theta.slider=gspinbutton(from=0,to=2,by=.1,value=theta,container=theta.group,
                       handler=function(h,...){
                         theta <<-svalue(theta.slider)
                       })
        min.d.group=ggroup(container=var.group)
        min.d.label=glabel("Minimum possible population in a cell",container=min.d.group)
        addSpring(min.d.group)
        min.d.slider=gspinbutton(from=0,to=1,by=.01,value=min.d,container=min.d.group,
                         handler=function(h,...){
                           min.d <<-svalue(min.d.slider)
                         })
        disp.max.group=ggroup(container=var.group)
        disp.max.label=glabel("Maximum dispersal probability",container=disp.max.group)
        addSpring(disp.max.group)
        disp.max.slider=gspinbutton(from=0,to=1,by=.01,value=disp.max,container=disp.max.group,
                         handler=function(h,...){
                           disp.max <<-svalue(disp.max.slider)
                         })
        disp.min.group=ggroup(container=var.group)
        disp.min.label=glabel("Minimum dispersal probability",container=disp.min.group)
        addSpring(disp.min.group)
        disp.min.slider=gspinbutton(from=0,to=1,by=.01,value=disp.min,container=disp.min.group,
                            handler=function(h,...){
                              disp.min <<-svalue(disp.min.slider)
                        })
        cost.int.group=ggroup(container=var.group)
        cost.int.label=glabel("Intercept of cost per animal equation",container=cost.int.group)
        addSpring(cost.int.group)
        cost.int.slider=gspinbutton(from=0,to=1,by=.01,value=cost.int,container=cost.int.group,
                            handler=function(h,...){
                              cost.int <<-svalue(cost.int.slider)
                            })
        cost.slope.group=ggroup(container=var.group)
        cost.slope.label=glabel("Slope of cost per animal equation",container=cost.slope.group)
        addSpring(cost.slope.group)
        cost.slope.slider=gspinbutton(from=-3,to=3,by=.01,value=cost.slope,container=cost.slope.group,
                            handler=function(h,...){
                              cost.slope <<-svalue(cost.slope.slider)
                            })
        helicopter.cost.group=ggroup(container=var.group)
        helicopter.cost.label=glabel("Hourly cost of helicopter",container=helicopter.cost.group)
        addSpring(helicopter.cost.group)
        helicopter.cost.slider=gspinbutton(from=500,to=2000,by=100,value=helicopter.cost,container=helicopter.cost.group,
                              handler=function(h,...){
                                helicopter.cost <<-svalue(helicopter.cost.slider)
                              })
        overhead.group=ggroup(container=var.group)
        overhead.label=glabel("Overhead cost (so 18% default)",container=overhead.group)
        addSpring(overhead.group)
        overhead.slider=gspinbutton(from=1,to=2,by=.01,value=overhead,container=overhead.group,
                              handler=function(h,...){
                                overhead <<-svalue(overhead.slider)
                              })
        duration.group=ggroup(container=var.group)
        duration.label=glabel("Overhead cost (so 18% default)",container=duration.group)
        addSpring(duration.group)
        duration.slider=gspinbutton(from=1,to=20,by=1,value=duration,container=duration.group,
                            handler=function(h,...){
                              duration <<-svalue(duration.slider)
                            })
        Load.Params.button <- gbutton("Load Parameter File", container=var.group,
                                      handler=function(h,...){
                                        param.file=fileChoose("print")
                                        if(is.na(param.file)) return
                                        load_params(param.file,new=T)
                                      })
        Save.Params.button <- gbutton("Save Parameter File", container=var.group,
                              handler=function(h,...){
                                param.file=fileSaveChoose("print")
                                if(is.na(param.file)) return
                                save_params(param.file)
                              })

    run.group <- ggroup(container=left.pane)
      run.model.button <- gbutton("Run Model",container=run.group,
                           handler=function(h,...){
                             run.model()
                           })
      run.NSB.button <- gbutton("Run Non-Spatial Budget",container=run.group,
                            handler=function(h,...){
                              run.NSB()
                            })
      run.NSD.button <- gbutton("Run Non-Spatial Density",container=run.group,
                          handler=function(h,...){
                            run.NSD()
                          })
      run.SB.button <- gbutton("Run Spatial Budget",container=run.group,
                          handler=function(h,...){
                            run.SB()
                          })


  right.pane <- ggroup(container=window,horizontal=F)
    output.pad <- gnotebook(container=right.pane,expand=T)
      pop.plot.page <- ggraphics(container=output.pad,label="Population-Time Plot")
      cost.plot.page <- ggraphics(container=output.pad,label="Cost-Time Plot")
      pn0.plot.page <- ggraphics(container=output.pad,label="Propotional Population")
      pop.ts.group <- ggroup(container=output.pad,label="Pop TS",horizontal=F)
        pop.ts.page <- ggraphics(container=pop.ts.group,label="Population TS")
        anim.but.group <- ggroup(container=pop.ts.group)
        prev.but <- gbutton("Previous",container=anim.but.group,
                            handler=function(h,...){
                              curr.t <<- max(1, curr.t - 1) # set time back by 1 (unless it's already at the first frame)
                              draw.pop.ts(curr.t) # draw
                            })
        next.but <- gbutton("Next",container=anim.but.group,
                    handler=function(h,...){
                      curr.t <<- min(2 * duration + 1, curr.t + 1) # set time forward by 1 (unless it's already at the last frame, note assuming 2 per year)
                      draw.pop.ts(curr.t) # draw
                    })
      CB.page <- ggraphics(container=output.pad,label="CB Map")
      OCM.page <- ggraphics(container=output.pad,label="OCM Map")

    
text.1 <- gtext("",container=right.pane,expand=T,font.attr=c(family="monospace"))

init.r

#!/usr/bin/Rscript

## Initiate spam model

## Load required libraries
options(guiToolkit="RGtk2")
#options(guiToolkit="tcltk")

packages = c("gWidgets","gWidgetsRGtk2","raster","fields","RGtk2","cairoDevice")

for(the.package in packages){
  if (!require(the.package,character.only=T)) {
    inst.text=paste("install.packages(",the.package,")",sep="")
    eval(parse(text=inst.text))
  }
  if (!require(the.package,character.only=T)) {
    print(paste("Could not install package '",the.packge,"', please contact your sysadmin.",sep=""))
    return() 
  }
}

#require(gWidgets)
#require(gWidgetsRGtk2)
#require(raster)
#require(fields)
#require(RGtk2)
#require(cairoDevice)


## Load default parameters, or a specified file.  If a new file is being loaded (ie. after interface
## definition), then update slider contents to reflect new file

load_params=function(the_file="params.csv",new_file=F){
  DATA=read.csv(the_file,header=F)
  eval(parse(text=paste(DATA$V1,"<<-",DATA$V2,sep="")))
  if(new_file){
    for(the_row in seq_along(DATA$V1)){
      to_run=paste("svalue(",DATA$V1[the_row],".slider)=",DATA$V2[the_row],sep="")
      eval(parse(text=to_run),envir=.GlobalEnv)
    }
  }
}

## Save parameters to a new file

save_params=function(the_file="params.csv"){
  param_list=c('budget','area.target','elev.d.mod','init.cull','maint.cull','target.density','cell.size','r.m','theta','min.d','disp.max','disp.min','cost.int','cost.slope','helicopter.cost','overhead','duration')
  params = c(budget,area.target,elev.d.mod,init.cull,maint.cull,target.density,cell.size,r.m,theta,min.d,disp.max,disp.min,cost.int,cost.slope,helicopter.cost,overhead,duration)
  save_frame=data.frame(param_list,params,stringsAsFactors=F)
  write.table(save_frame,the_file,col.names=F,row.names=F,sep=",",quote=F)
}
## Load functions

load_params()

source("spam_ui.r")
source("engine.r")
source("gui.r")

#pause_for_input <- function()
#{
#  message("Press ENTER to continue")
#  invisible(scan(n = 0, quiet = TRUE))
#}

#pause_for_input()
#while(readline()!="x"){}
gtkMain()

params.csv

spam_ui.r

# clear all variables
#rm(list=ls(all=TRUE))


# Load parameters from a .csv file into variable space
#load.params = function()
#{
  #update.graphics() # clears graphics window
  #cat(rep('\n',64)) # clears console window

#  try.error = TRUE; exit = FALSE # to start while loop

#  while (try.error & !exit) # while we keep getting an error and the user hasn't chosen to exit
#  {
#    filename = try(file.choose(), silent = TRUE) # choose filename using a dialog
#    if (inherits(filename, 'try-error')) # if the user has chosen to exit
#      exit = TRUE
#    else
#    {
#      ext = extension(filename)
#      if (ext == '.csv')  # does the filename have the right extension?
#      {
#        DATA = try(read.csv(filename, header = FALSE), silent = TRUE) # try to read it
#        try.error = inherits(DATA, 'try-error') # if it fails, we have an error
#      }
#      else   
#        try.error = TRUE # if wrong extension, we have an error
    
#      if (try.error) winDialog(NULL, "File not compatible! Please select file in either .csv or raster format") # complain to user about error
#    }
#  }
  
#  if (!exit) # if the user hasn't exited, continue processing
#  {
    # creates an R file (tmp.r) based on the csv
    # first column is parameter names, second column is parameter values (anything else won't be read by this script)
    # NOTE: security risk - user could conceivably blow things up with incorrect and/or malicious entries in csv file
 #   e = try(cat(sprintf('%s = %f\n', as.character(DATA$V1), DATA$V2), file='tmp.r'), silent = TRUE) 
#    error = FALSE    
#    if (inherits(e, 'try.error')) # if for some reason R can't write the R file
#      error = TRUE
#    else
#    {
#      e = try(source('tmp.r'), silent = TRUE) # try to run the R file, thus generating variables with the names and values given
#      if (inherits(e, 'try.error'))
#        error = TRUE
#    }
    
 #   if (error) # if it didn't work for either reason
#      winDialog(NULL, "Error loading parameters file. Check file for errors")
#    else
#    {
#      winDialog(NULL, 'Data successfully loaded (we hope...)')
#      PARAMS.FILE <<- filename # generate a global variable of the filename for use in edit.params()
#    }
#  }
#  else # if the user exited
#  {
#    winDialog(NULL, 'Data not loaded')
#  }

#  cat(rep('\n',64)) # clears console window
#  NULL # return nothing
#}


# Edits an already loaded parameters file using Notepad
edit.params = function()
{
  shell(paste('notepad.exe "', PARAMS.FILE, '"', sep=''))  # Runs Notepad and uses a global variable for the filename - R will wait until Notepad is closed before continuing
  DATA = try(read.csv(PARAMS.FILE, header = FALSE), silent = TRUE)  # Read the newly-edited file into memory

  error = inherits(DATA, 'try-error')  # error checking - file read error
  if (error) 
    winDialog(NULL, "File read error")
  else
  {
    e = try(cat(sprintf('%s = %f\n', as.character(DATA$V1), DATA$V2), file='tmp.r'), silent = TRUE) # create R file - see load.params() for details
    
    if (inherits(e, 'try.error'))
      error = TRUE
    else
    {
      e = try(source('tmp.r'), silent = TRUE)
      if (inherits(e, 'try.error'))
        error = TRUE
    }
  }
  
  if (error)
    winDialog(NULL, "Error loading parameters file. Check file for errors")
  else
  {
    winDialog(NULL, 'Data successfully loaded (we hope...)')
  }  
    
}


# small routine to plot a gridded raster on the graphics window
# (Note that resizing window puts the raster out of line with the grid. 
# I suspect this is a problem with the "image" routine as the grid lines 
# are still in the right places [can confirm by running axis(1) and axis(2) 
# but the image itself is not)
plot.raster = function(DATA, ...) 
{
    image.plot(seq(0, 0.9, l=ncol(DATA)+1), seq(0, 0.7, l=nrow(DATA)+1), t(DATA)[, nrow(DATA):1], add=TRUE, col = rainbow(1024)[1:350], ...)
    for (i in 0:ncol(DATA)) lines(0.9*rep(i/ncol(DATA), 2), c(0, 0.7))
    for (i in 0:nrow(DATA)) lines(c(0,0.9), 0.7*rep(i/nrow(DATA), 2))
}


# prompts the user to load a raster (or .csv) file and returns it (similar to load.params())
load.raster = function(message) 
{
  #update.graphics()
  #cat(rep('\n',64))
  try.error = TRUE
  exit = FALSE

  while (try.error & !exit)
  {
    filename = try(file.choose(), silent = TRUE)
    if (inherits(filename, 'try-error'))
      exit = TRUE
    else
    {
      ext = extension(filename)
      if (ext == '.csv')  
        DATA = try(read.csv(filename, header = FALSE), silent = TRUE) 
      else   
        DATA = try(as.matrix(raster(filename)), silent = TRUE)
      try.error = inherits(DATA, 'try-error')
      if (try.error) winDialog(NULL, "File not compatible! Please select file in either .csv or raster format")
    }
  }
  
  if (!exit)
  {
    winDialog(NULL, 'Data successfully loaded (we hope...)')
    cat(rep('\n',64))
    plot.raster(DATA) # draws loaded raster to screen
    DATA
  }
  else
  {
    winDialog(NULL, 'Data not loaded')
    cat(rep('\n',64))
    NULL
  }
}


# updates graphics information and menus based on current variables
update.status = function()   
{
  lights = rep(NA, 5) # creates the five "lights" showing which required information has been loaded (1 = not loaded = red, 2 = loaded = green)
  env = ls('.GlobalEnv') # loads the names of all variables currently in R's workspace
  
  if ("out" %in% env) # "out" contains the current model run - if it exists, we can enable the results menu
    action = 'enable'
  else
    action = 'disable'

  # enable or disable all Results menu items (these have already been created in open.graphics())
  results.items = c('Total population plot', 'Total cost plot', 'Reduction in control region', 'Population time series', 'Cost-benefit value map', 'Optimal culling map')
  for (i in 1:5)
    winMenuAddItem('$Graph2Main/Results', results.items[i], action) 
  
  if ("optim.cull" %in% env) # "optim.cull" contains the optimal culling map from the spatial optimisation - if it exists, we can enable the result separately
    action = 'enable'
  else
    action = 'disable'
  winMenuAddItem('$Graph2Main/Results', results.items[6], action) 


  variables = c('K', 'EL', 'PR', 'C.mask', 'duration')
  for (i in 1:5)
    lights[i] = ((variables[i] %in% env) + 1) # if variable is in memory, light the corresponding button green

  # if ALL variables are in memory, enable Model menu items
  if (all(lights == 2)) 
    action = 'enable'
  else
    action = 'disable'
  run.menu = c('$Graph2Main/Run', rep('$Graph2Main/Run/Optimisation', 3))
  run.items = c('Model', 'Non-spatial budget', 'Non-spatial density', 'Spatial budget')
  for (i in 1:4)
    winMenuAddItem(run.menu[i], run.items[i], action)

  # enable or disable editing each file in the Edit menu based on whether it has been loaded 
  edit.items = c('Carrying capacity', 'Elevation', 'Priority', 'Culling', 'Parameters')
  for (i in 1:5)
  {
    if (lights[i] == 2) 
      winMenuAddItem('$Graph2Main/Edit', edit.items[i] , 'enable')
    else
      winMenuAddItem('$Graph2Main/Edit', edit.items[i], 'disable')
  }    

  # If the model has not been run yet, render the lights (otherwise the space will be used for model results text info)
  if (!("out" %in% env))
  {
    palette(c('red', 'green'))
    points(rep(0, 5), 0.75 + 0.05 * (4:0), col=lights, pch=19, cex=1.5)    
  }
}


# Update graphics information based on whether model has been run (clears screen)
update.graphics = function()
{
  #dev.set(2) # make sure we are using graphics window 2
  #plot(c(0,1), c(0,1), col='white', axes=FALSE, xlab='', ylab='') # clear screen
  env = ls('.GlobalEnv') # loads all variable names in R workspace

  if ("out" %in% env) # if model has been run
  {
    if ("test" %in% env) # if an optimisation has been done
    {
      insert(text.1,sprintf('Total cost: $%d\nAnimals culled: %d\nTotal population post-cull: %d\nProportion remaining in control region: %.1f%%\n\nCull rate: %d%%\nTarget density: %.2f / sq. km', 
        sum(out$total.cost), round(sum(out$total.cull)), round(out$total.pop[length(out$total.pop)]), out$P.N0[length(out$P.N0) - 1]*100, round(test$cull.rate*100), test$target.density))
    }
    else # if a model but no optimisation
    {
      insert(text.1,sprintf('Total cost: $%d\nAnimals culled: %d\nTotal population post-cull: %d\nProportion remaining in control region: %.1f%%', 
        sum(out$total.cost), round(sum(out$total.cull)), round(out$total.pop[length(out$total.pop)]), out$P.N0[length(out$P.N0) - 1]*100))
    }
  }
  #else # draw labels for "lights"
  #{
  #  insert(text.1,'Required files:', pos=4)
  #  insert(text.1,'Carrying capacity', pos=4)
  #  insert(text.1,'Elevation', pos=4)
  #  insert(text.1,'Priority', pos=4 )
  #  insert(text.1,'Culling', pos=4)
  #  insert(text.1,'Parameters')      
  #}

  # updates other graphics information
 # update.status()
}


# initialises graphics - loads all menus
#open.graphics = function()
#{
  # Close all windows then open the one we want
  #graphics.off()
  #dev.new(2)

  # tricky business (see e.g. http://tolstoy.newcastle.edu.au/R/e3/help/07/11/3254.html)
  #winMenuAdd('$Graph2Main/File')
  #winMenuAdd('$Graph2Main/History')
  #winMenuAdd('$Graph2Main/Resize')
  #winMenuDel('$Graph2Main/Resize')
  #winMenuDel('$Graph2Main/History')
  #winMenuDel('$Graph2Main/File')

  # Load menus and function calls
  #winMenuAddItem('$Graph2Main/Load', 'Carrying capacity', 'K = load.raster("Please select file containing species\' carrying capacity"); update.status()')
  #winMenuAddItem('$Graph2Main/Load', 'Elevation', 'EL = load.raster("Please select file containing elevation data"); update.status()')
  #winMenuAddItem('$Graph2Main/Load', 'Priority', 'PR = load.raster("Please select file containing priority information"); update.status()')
  #winMenuAddItem('$Graph2Main/Load', 'Culling', 'C.mask = load.raster("Please select file containing culling information"); update.status()')
  #winMenuAddItem('$Graph2Main/Load', 'Parameters', 'load.params(); update.status()')

  #winMenuAddItem('$Graph2Main/Edit', 'Carrying capacity', 'K = edit.raster(K, 0)')
  #winMenuAddItem('$Graph2Main/Edit', 'Elevation', 'EL = edit.raster(EL, 0, 1)')
  #winMenuAddItem('$Graph2Main/Edit', 'Priority', 'PR = edit.raster(PR, 0)')
  #winMenuAddItem('$Graph2Main/Edit', 'Culling', 'C.mask = edit.raster(C.mask, 0, 1)')
  #winMenuAddItem('$Graph2Main/Edit', 'Parameters', 'edit.params()')

  #winMenuAddItem('$Graph2Main/Run', 'Model', 'run.model()')
  #winMenuAddItem('$Graph2Main/Run/Optimisation', 'Non-spatial budget', 'run.NSB()')
  #winMenuAddItem('$Graph2Main/Run/Optimisation', 'Non-spatial density', 'run.NSD()')
  #winMenuAddItem('$Graph2Main/Run/Optimisation', 'Spatial budget', 'run.SB()')  

  #winMenuAddItem('$Graph2Main/Results', 'Total population plot', 'draw.pop.plot()')
  #winMenuAddItem('$Graph2Main/Results', 'Total cost plot', 'draw.cost.plot()')
  #winMenuAddItem('$Graph2Main/Results', 'Reduction in control region', 'draw.pn0.plot()')
  #winMenuAddItem('$Graph2Main/Results', '-', '')
  #winMenuAddItem('$Graph2Main/Results', 'Population time series', 'draw.pop.anim()')
  #winMenuAddItem('$Graph2Main/Results', 'Cost-benefit value map', 'draw.cb()')
  #winMenuAddItem('$Graph2Main/Results', 'Optimal culling map', 'draw.ocm()')

  # Update all graphics information (this also calls update.status())
  #update.graphics()

  #cat(rep('\n',64))
#}


# Runs a basic STAR-style simulation based on the current parameters and variables
run.model = function()
{
  out <<- run.sim(init.cull, maint.cull, target.density, budget, C.mask, C.mask) # outputs model result to "out"
  draw.pop.plot()
  draw.cost.plot()
  draw.pn0.plot()
  curr.t <<- 1 
  draw.pop.ts(curr.t)
  draw.cb()
  #draw.ocm()
  #winDialog(NULL, 'Model done!')
  update.graphics() # update graphics (in particular, output some text model results to screen)
  #cat(rep('\n',64))    
}


# Runs a non-spatial budget optimisation based on the current parameters and variables
run.NSB = function()
{
  test <<- NS.B(budget, C.mask) # outputs results of optimisation to "test"
  if (is.na(test$cull.rate)) # if the budget is too small to give a result
  {
    confirmDialog('Budget too small!')
    rm(test, envir = globalenv())
  }
  else # run model based on found optimal parameters
  {
    out <<- run.sim(test$cull.rate, test$cull.rate, test$target.density, budget, C.mask, C.mask) 
    confirmDialog( 'Optimisation done!')
  }
  update.graphics() # update graphics (in particular, output some text model + optim results to screen)
  draw.pop.plot()
  draw.cost.plot()
  draw.pn0.plot() 
  curr.t <<- 1 
  draw.pop.ts(curr.t)
  draw.cb()
  #cat(rep('\n',64))    
}


# Runs a non-spatial control-area density optimisation based on the current parameters and variables
run.NSD = function()
{
  test <<- NS.D(area.target, C.mask) # outputs results of optimisation to "test"
  if (is.na(test$cull.rate)) # if the required density is too small to give a result
  {
    confirmDialog('Can\'t hit target density even at 99% cull rate!')
    rm(test, envir = globalenv())
  }
  else # run model based on found optimal parameters
  {
    out <<- run.sim(test$cull.rate, test$cull.rate, test$target.density, budget, C.mask, C.mask) 
    confirmDialog('Optimisation done!')
  }
  update.graphics() # update graphics (in particular, output some text model + optim results to screen)
  draw.pop.plot()
  draw.cost.plot()
  draw.pn0.plot()
  curr.t <<- 1 
  draw.pop.ts(curr.t)
  draw.cb()
  #cat(rep('\n',64))    
}


# Runs a spatial budget optimisation based on the current parameters and variables
run.SB = function()
{
  test2 <<- S.B(budget, C.mask) # outputs results of optimisation to "test"
  if (is.na(test2$test$cull.rate)) # if the budget is too small to give a result
  {
    confirmDialog('Budget too small!')
    rm(test2, envir = globalenv())
  }
  else # run model based on found optimal parameters
  {
    optim.cull <<- test2$culling.area
    out <<- test2$out
    test <<- test2$test
    confirmDialog('Optimisation done!')
  }
  draw.ocm()
  draw.pop.plot()
  draw.cost.plot()
  draw.pn0.plot() 
  curr.t <<- 1 
  draw.pop.ts(curr.t)
  draw.cb()
  update.graphics() # update graphics (in particular, output some text model + optim results to screen)
  #cat(rep('\n',64))    
}


# Draws total population against time (note that we assume 2 time periods per year as per STAR)
draw.pop.plot = function()
{
  #update.graphics()
  visible(pop.plot.page)=T
  plot.new()
  x.v = seq(0.1, 1, l = length(out$total.pop)) # define plot area
  y.v = 0.1 + out$total.pop / max(out$total.pop) * 0.5
  
  points(x.v, y.v, pch=19, ylim=c(0, 0.6))

  v = pretty(seq(0, duration, 0.5)) # create axes independently (based on actual data and not graphics coordinates)
  axis(1, at = v * 0.9 / max(v) + 0.1, labels = v, pos=0.1)

  v = pretty(c(0, max(out$total.pop)))
  axis(2, at = 0.1 + v * 0.5 / max(out$total.pop), labels = v, pos=0.1)

  text(0.6, 0.0, 'Years')  # create axis labels
  text(0.0, 0.4, 'Population', srt=90)
  #cat(rep('\n',64))    
}


# Draws total cost against time (note that currently culling is set to only happen once per year, i.e. every second time period)
draw.cost.plot = function()
{
  #update.graphics()
  visible(cost.plot.page)=T
  plot.new()
  x.v = seq(0.1, 1, l = length(out$total.cost)) # define plot area
  y.v = 0.1 + out$total.cost / max(out$total.cost) * 0.5
  points(x.v, y.v, pch=19, ylim=c(0, 0.6))

  v = pretty(seq(0.5, duration, 0.5)) # create axes independently (based on actual data and not graphics coordinates)
  axis(1, at = v * 0.9 / max(v) + 0.1, labels = v, pos=0.1)

  v = pretty(c(0, max(out$total.cost)))
  axis(2, at = 0.1 + v * 0.5 / max(out$total.cost), labels = v / 1e6, pos=0.1)

  text(0.6, 0.0, 'Years')    # create axis labels
  text(0.0, 0.4, 'Cost (million $)', srt=90)
  #cat(rep('\n',64))    
}


# Draws proportional population in control area (against carrying capacity - note that currently carrying capacity is the initial condition)
draw.pn0.plot = function()
{
  #update.graphics()
  visible(pn0.plot.page)=T
  plot.new()
  P.N0 = c(1, out$P.N0) # set original proportion as 1 (before any culling or movement)
  x.v = seq(0.1, 1, l = length(P.N0))
  y.v = 0.1 + P.N0 / max(P.N0) * 0.5
  points(x.v, y.v, pch=19, ylim=c(0, 0.6))

  v = pretty(seq(0, duration, 0.5)) # create axes independently (based on actual data and not graphics coordinates)
  axis(1, at = v * 0.9 / max(v) + 0.1, labels = v, pos=0.1)

  v = pretty(c(0, max(P.N0)))
  axis(2, at = 0.1 + v * 0.5 / max(P.N0), labels = v, pos=0.1)

  text(0.6, 0.0, 'Years') # create axis labels 
  text(0.0, 0.4, 'Proportion surviving', srt=90)
  #cat(rep('\n',64))
}


# Draw a single frame of the population map time series
draw.pop.ts = function(t)
{
  #update.graphics() # clear screen
  visible(pop.ts.page)=T
  plot.new()
  #lines(0.4 + 0.1 * c(0,0,1,1,0), 0.8 + 0.1 * c(0,1,1,0,0)) # draw "Back" button
  #lines(0.8 + 0.1 * c(0,0,1,1,0), 0.8 + 0.1 * c(0,1,1,0,0)) # "Forward" button
  #lines(0.6 + 0.1 * c(0,0,1,1,0), 0.9 + 0.1 * c(0,1,1,0,0)) # "Stop" button
  #polygon(0.4 + 0.1 * c(1,3,3,1)/4, 0.8 + 0.1 * c(2,1,3,2)/4, col='grey') # coloured triangles
  #polygon(0.8 + 0.1 * c(3,1,1,3)/4, 0.8 + 0.1 * c(2,1,3,2)/4, col='grey')  
  #text(0.65, 0.95, 'Stop', col='red') # "Stop" text
  text(0.65, 0.85, sprintf('Time\n%.1f', (t-1)/2)  ) # Current time frame info (note assuming 2 frames per year)

  DATA = out$pop.ts[[t]] # load and plot required time frame
  plot.raster(DATA, zlim = c(0, range(out$pop.ts, na.rm=TRUE)[2]))
  #cat(rep('\n',64))    
}


# What to do if ANY mouse button is pressed when viewing population map time series
mousedown.1 = function(buttons, x, y) 
{
  x = grconvertX(x, "ndc", "user") # change coordinates from window to axis coords
  y = grconvertY(y, "ndc", "user")

  if (x > 0.6 & x < 0.7 & y > 0.9 & y < 1) # if mouse is inside "Stop" button
  {
    return(invisible(1)) # stop event handling
  }
  else if (y > 0.8 & y < 0.9) 
  {
    if (x > 0.4 & x < 0.5) # if mouse is inside "Back" button
    {
      curr.t <<- max(1, curr.t - 1) # set time back by 1 (unless it's already at the first frame)
      draw.pop.ts(curr.t) # draw
    }
    else if (x > 0.8 & x < 0.9) # if mouse is inside "Forward" button
    {
      curr.t <<- min(2 * duration + 1, curr.t + 1) # set time forward by 1 (unless it's already at the last frame, note assuming 2 per year)
      draw.pop.ts(curr.t) # draw
    }
  }
  NULL
}


# Handle population map time series
draw.pop.anim = function()
{
  
  #update.graphics()  
  curr.t <<- 1 # start at first frame
  draw.pop.ts(curr.t)
  
  winDialog(NULL, 'Press "Stop" button when done, before selecting other menu options') # User message - selecting any other option won't work while event handling is on
  
  setGraphicsEventHandlers(onMouseDown = mousedown.1) # tell R how to handle mouse clicks
  getGraphicsEvent()  # start event handling so user can examine time series
  update.graphics()  # when user is done, clear screen
  
  cat(rep('\n',64))  
}


# Draw cost-benefit map
draw.cb = function()
{
  #update.graphics()
  visible(CB.page)=T
  plot.new()
  DATA = out$cost.benefit
  plot.raster(DATA)
  
  #cat(rep('\n',64))  
}


# Draw optimal culling map
draw.ocm = function()
{
  #update.graphics()
  visible(OCM.page)=T
  plot.new()
  DATA = test2$culling.area
  plot.raster(DATA)
  
  #cat(rep('\n',64))  
}


# What to do if a mouse button is pressed when raster editing
mousedown.2 = function(buttons, x, y) 
{
  if (MOUSEDOWN) # If a mouse button is currently being pressed
  {
    x = grconvertX(x, "ndc", "user")  # change coordinates from window to axis coords
    y = grconvertY(y, "ndc", "user")
  
    L.x = nrow(DATA.OUT) # work out how big the current raster is in cells
    L.y = ncol(DATA.OUT)
  
    if (x > 0.6 & x < 0.7 & y > 0.9 & y < 1) # if the user has pressed the "Save" button
    {
      SAVE <<- TRUE
      return(invisible(1))
    }
    else if (x > 0.7 & x < 0.8 & y > 0.9 & y < 1) # if the user has pressed the "Cancel" button
    {
      SAVE <<- FALSE
      return(invisible(1))
    }
    else if (x > 0 & x < 0.9 & y > 0 & y < 0.7) # if the user has clicked within the raster
    {
      j = 1 + floor(x * L.y /0.9) # work out what cell we're in
      i = 1 + floor(y * L.x /0.7)
      DATA.OUT[L.x + 1 - i, j] <<- min(max(MIN, DATA.OUT[L.x + 1 - i, j] + (buttons - 1)), MAX) # as long as within MIN-MAX range: subract 1 if left button clicked, add 1 if right clicked(do nothing if middle mouse)
      update.graphics() # clear screen
      lines(0.6 + 0.1 * c(0,0,1,1,0), 0.9 + 0.1 * c(0,1,1,0,0)) # redraw "Save" button
      text(0.65, 0.95, 'Save', col='dark green')
      lines(0.7 + 0.1 * c(0,0,1,1,0), 0.9 + 0.1 * c(0,1,1,0,0)) # redraw "Cancel" button
      text(0.75, 0.95, 'Cancel', col='dark red')      
      plot.raster(DATA.OUT) # draw current raster
    }
  }
  NULL
}


# Handle raster editing
edit.raster = function(DATA, min.val = -Inf, max.val = Inf)
{
  update.graphics()  
  lines(0.6 + 0.1 * c(0,0,1,1,0), 0.9 + 0.1 * c(0,1,1,0,0)) # draw "Save" button
  text(0.65, 0.95, 'Save', col='dark green')
  lines(0.7 + 0.1 * c(0,0,1,1,0), 0.9 + 0.1 * c(0,1,1,0,0)) # draw "Cancel" button
  text(0.75, 0.95, 'Cancel', col='dark red')
  
  plot.raster(DATA) # draw current raster (DATA)
  
  DATA.OUT <<- DATA # create global variable of current raster for access by event handling routine
  SAVE <<- TRUE # default to save
  MIN <<- min.val # create global variable of max and min values for event handling
  MAX <<- max.val
  MOUSEDOWN <<- FALSE # start with mouse unclicked
  
  winDialog(NULL, 'Left-click a cell to subtract 1 from value; right-click to add 1') # User instructions
  
  # If mouse clicked, let handler know; if mouse moved, continue changing cells that are dragged over; if mouse unclicked, let handler know to stop editing cells
  setGraphicsEventHandlers(onMouseDown = function(buttons, x, y) {MOUSEDOWN <<- TRUE; mousedown.2(buttons, x, y)}, onMouseMove = mousedown.2, onMouseUp = function(buttons, x, y) {MOUSEDOWN <<- FALSE; NULL})
  getGraphicsEvent() # start event handling so user can edit raster
  update.graphics() # clear screen when done
  cat(rep('\n',64))    

  if (SAVE) # if raster is saved, overwrite original raster with edited raster
    DATA.OUT
  else # otherwise, just return the original raster
    DATA
    
  # note that the raster file itself is not edited (unlike the parameters csv file when parameter editing)
}


# once all functions are loaded into memory, start!
#open.graphics()