Example to run tripestimation

This example runs with the data set hoopoe1 included in the GeoLight package. The process for defining twilight segments is not included here, since this is planned to be run on a website UI in the future.

For that reason, the process that would otherwise look like this is avoided by directly loading a data set called d1 below.

 #NOTE: this is pseudocode, it's not for use
 library(GeoLight)
 data(hoopoe1)
 d <- hoopoe1
 d$segment <- defineSegments(d$datetime, d$light)
 #remove non-twilight data
 d1 <- d1[!is.na(d$segment), ]
 # ensure segment is a clean sequence of 1:ntwilights
 d1$segment <- unclass(factor(d1$segment))

Set up and run tripEstimation with data previously prepared.

Make the calibration function.

library(GeoLight)
library(trip)
library(tripEstimation)
library(maptools)


setwd("C:\\Users\\user\\Desktop\\MikeSumner\\runtripEstimation")
##source("BASsegments_m.r")


## example data is in GeoLight
data(hoopoe1)
d <- data.frame(id=1:nrow(hoopoe1),gmt=as.POSIXct(hoopoe1$datetime,"GMT"),
				light=hoopoe1$light)


## basic bounding box for entire study
lon.min <- -20
lon.max <- 55
lat.min <- -20
lat.max <- 50

## "start" is the known location, we use for calibration and
## for fixing a single coordinate in the model (we don't know the "end")
start <- c(7.1, 46.3)
## end <-

## subset of entire data set that we extract for calibration
d.calib <- d[d$gmt > as.POSIXct("2008-07-15 00:02:00") & d$gmt < as.POSIXct("2008-07-25 00:00:00", tz="GMT"),]

## apply all candidate twilights (in real terms this is an interactive
## web-site process that requires a fair bit of user inspection and
## thought)

calib <- mkCalibration(d.calib, known = start, elim = c(-10, 8), choose = FALSE)

## load the prepared example data set directly from the web
## NOTE: MDS doesn't know how to do this from https on GitHub
load(url("http://staff.acecrc.org.au/~mdsumner/tripEstimation/examples/runHoopoe1/d1.RData", open = "rb"))

names(d1)[1] <- "gmt"
d1 <- d1[!is.na(d1$segment), ]
d1$segment <- unclass(factor(d1$segment))

## can we assume that the "release" and/or "recapture" site/s are
## known, these are provided to the initialization process and simply
## assigned as fixed coordinates if known

fixed.release <- TRUE
fixed.recapture <- FALSE

## how many twilights plus fixed points?
m <- max(d1$segment) + fixed.release + fixed.recapture

Create a mask (only on land, but allowed to fly over Mediterranean).

## "topo" is a basic "xyz" list as used by ?image, with a topo$z matrix that
## encodes spatially TRUE/FALSE for land/sea (plus a bit of jigger)
data(wrld_simpl)
bb <- matrix(c(-20,55,-20,50),byrow=T,ncol=2)
xx <- seq(bb[1,1] - 1, bb[1,2]+ 1, length = 200)
yy <- seq(bb[2,1] - 1, bb[2,2]+ 1, length = 250)
res <- overlay(SpatialPoints(expand.grid(x = xx, y = yy)), wrld_simpl)
zz <- !is.na(res)
topo <- list(x = xx, y = yy, z = matrix(zz, length(xx), length(yy)))
topo$z[which(xx>=-6 & xx<=10),which(yy <= 41 & yy >=30)] <- TRUE
topo$z[which(xx>=0 & xx<=45),which(yy <= 50 & yy >=29)] <- TRUE
rm(res)
##image(topo)

See Creating simple masks for more options.

Create the lookup function and other objects required by the model

## lookup function that encapsulates the mask object, this function
## takes a matrix of coordinates of [lon,lat] and returns a logical
## vector for each lon,lat pair of whether it is allowed by the mask

lookup <- mkLookup(topo, FALSE)


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # Proposal functions, these create candidate coordinates
# for individual iteration steps in the MCMC
## m twilights, 3 params for Xs (lon,lat,attenutation), 2 params for Zs

proposal.x <- norm.proposal(m, 3, c(0.5, 0.5, 0.05))
proposal.z <- norm.proposal(m-1, 2, c(1, 1))

speed.mu <- 50/3.6
speed.sigma <- 20

Create and begin run of the model .

## create the model
d.model <- solar.model(d1$segment,d1$gmt,d1$light,
                       proposal.x = proposal.x$proposal,
                       proposal.z = proposal.z$proposal,
                       fix.release = fixed.release,
                       fix.recapture = fixed.recapture,
                       mask.x = lookup,
                       mask.z = lookup,
                       calibration=calib,
                       behav = "speed",
                       behav.dist = "log",
                       light.sigma=3,
                       k.sigma = 4,
                       behav.mean = speed.mu,
                       behav.sd = speed.sigma,
                       ekstrom = c(-8, 5, 128))

## find some starting points that satisfy the land
xx <- seq(lon.min, lon.max, length = 30)
yy <- seq(lat.min, lat.max, length = 30)
plg <- position.logp(d.model, xx, yy, xrest = -0.5,
start = start, end = FALSE, prob = 0.9)


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# initialization

X <- plg$X

xfile <- "initx.bin"
zfile <- "initz.bin"

ch <- metropolis0(d.model,iters=10,thin=1,start.x = X,
                  start.z = (X[-m,1:2]+X[-1,1:2])/2)

chain.write(xfile, ch$x)
chain.write(zfile, ch$z)

while(!all(ch$mask.x) | !all(ch$mask.z)) {
    ch <- metropolis0(d.model, iters=20,
    				  thin=10, start.x=ch$last.x,
				      start.z=ch$last.z)
    chain.write(xfile, ch$x, append = TRUE)
    chain.write(zfile, ch$z, append = TRUE)
}


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Estimation


xfile <- "burnx.bin"
zfile <- "burnz.bin"

ch <- metropolis(d.model,iters=100,thin=10,start.x = ch$last.x,
                  start.z = ch$last.z)

chain.write(xfile, ch$x)
chain.write(zfile, ch$z)


for (i in 1:3) {

    ch <- metropolis(d.model,iters=200,thin=10,start.x = ch$last.x,
                      start.z = ch$last.z)

    chain.write(xfile, ch$x, append = TRUE)
    chain.write(zfile, ch$z, append = TRUE)
    ##plot(wrld_simpl, xlim=c(lon.min, lon.max-30),
    	##	ylim=c(lat.min+20, lat.max), col = "grey")
    ##lines(ch$x[,,1], lwd = 2)
    ##lines(ch$last.x,col="red")
}


for (i in 1:3) {

    ch <- metropolis(d.model,iters=2000,thin=10,start.x = ch$last.x,
                      start.z = ch$last.z)

    chain.write(xfile, ch$x, append = TRUE)
    chain.write(zfile, ch$z, append = TRUE)
    proposal.x$tune(ch$x, 0.3)
    proposal.z$tune(ch$z, 0.3)
}



## START SAVING
xfile <- "X.bin"
zfile <- "Z.bin"
for (i in 1:20) {
    ch <- metropolis(d.model,iters=2000,thin=10,start.x = ch$last.x,
                     start.z = ch$last.z)
    xm <- apply(ch$z, 1:2, mean)
    plot(xm)
    map(add = TRUE)
    points(xm)
    proposal.x$tune(ch$x, 0.3)
    proposal.z$tune(ch$z, 0.3)
    chain.write(xfile, ch$x, append = i > 1)
    chain.write(zfile, ch$z, append = i > 1)
}