This R Markdown Notebook accompanies two SEEC Toolbox Seminars given by Res Altwegg at the University of Cape Town, one on 25 August 2016 and the other one on 31 May 2018.

Occupancy models are becoming increasingly popular and these two seminars were intended to introduce the basic ideas of the single-season occupancy model. I show how to fit these models in a fequentist way using the R package `unmarked’, and in a Bayesian way using JAGS.

Occupancy can be a useful state variable in many situations. E.g. if we want to know where a species occurs and where it doesn’t occur, but we are less worried about density of individuals in places where a species does occur. For example, metapopulation ecology often deals with patch occupancy as a state variable; some of the IUCN re-listing criteria are based on occupancy (the range of a species, and area of occupancy); and managers may need to know which species occur in which of the sites they are managing.

One way to think of occupancy is to imagine a set of \(n\) discrete sites. Occupancy \(\Psi\) is then the proportion of sites that are occupied \(\frac{occupied}{n}\). If we knew exactly which sites are occupied and not occupied, we could use \(\frac{occupied}{n}\) as an estimator of occupancy. And we could use methods like logistic regression, \(logit(\Psi) = f(covariates)\) to examine patterns in occupancy. The trouble is that we usually only know where a species was detected or not. If we visit every single site we are interested in and record whether we find the species, there is still a good chance that our information is incomplete because we missed the species in many places where it actually occurs. That is because the detection probability \(p\) is usually less than 1. The fraction of sites that are {} to be occupied is thus in reality \(\Psi \times p\), which is only clost to \(\Psi\) if \(p \sim 1\). And a logistic regression on these detections would actually model \(logit(\Psi \times p) = f(covariates)\). Any pattern found in this way could be the result of variation in occupancy or detection, or both. We are interested in occupancy and so need to find a way to separate the detection process from the occupancy process. For this, we need some kind of information on the observation process. One way to get at the observation process is by visiting sites repeatedly, each time recording whether the species was seen or not.

Simple occupancy model with constant \(\Psi\) and \(p\)

Here is a simulation to demonstrate these ideas and how a single-season occupancy model can be used to estimate true occupancy. I use an artifical landscape of \(20 \times 20\) grid cells, each of them occupied with probability \(\Psi = 0.5\). With each visit, we have a detection probability \(p = 0.3\) of finding the species if the grid cell is occupied.

set.seed(246)
occ.prob <- 0.5
det.prob <- 0.3
grid <- expand.grid(long=c(1:20),lat=c(1:20))
grid$occ <- rbinom(dim(grid)[1],1,prob = occ.prob) # occupancy
grid$col <- ifelse(grid$occ==1,'orange','grey')
grid$detected1 <- rbinom(dim(grid)[1],1,prob = grid$occ * det.prob) # detections on first survey
grid$detected2 <- rbinom(dim(grid)[1],1,prob = grid$occ * det.prob) # detections on second survey
grid$detected3 <- rbinom(dim(grid)[1],1,prob = grid$occ * det.prob) # detections on third survey

Let’s plot the simulated landscape. Truly occupied grid cells are orange; detections black. A grey cell means that the species was not detected.

op <- par(mfrow=c(2,2), mar=c(1,1,1,1))
plot(lat~long,col=col,pch= 15, data=grid,ann=F,axes=F,cex=1.5)
plot(lat~long,col='grey',pch= 15, data=grid,ann=F,axes=F,cex=1.5)
points(lat~long,col=detected1,pch= 15, data=grid[grid$detected1==1,],cex=1.5)
plot(lat~long,col='grey',pch= 15, data=grid,ann=F,axes=F,cex=1.5)
points(lat~long,col=detected2,pch= 15, data=grid[grid$detected2==1,],cex=1.5)
plot(lat~long,col='grey',pch= 15, data=grid,ann=F,axes=F,cex=1.5)
points(lat~long,col=detected3,pch= 15, data=grid[grid$detected3==1,],cex=1.5)
par(op)

For the analysis, we need the data in the format of survey histories:

y <- grid[,5:7]
head(y)
dim(y)
[1] 400   3

In how many grid cells did we detect the species at least once?

sum(apply(y,1,sum)>0)/400 
[1] 0.3775

This `naive’ estimate of occupancy is lower than the true occupancy probability, which – as we know in this case because we simulated the data – is 0.5. Let’s see whether an occupancy model retrieves a better estimate.

require(unmarked)
y <- grid[,5:7]
ex1 <- unmarkedFrameOccu(y = y)
summary(ex1)
unmarkedFrame Object

400 sites
Maximum number of observations per site: 3 
Mean number of observations per site: 3 
Sites with at least one detection: 151 

Tabulation of y observations:
  0   1 
991 209 
m1 <- occu(~1 ~1, data=ex1)
summary(m1)

Call:
occu(formula = ~1 ~ 1, data = ex1)

Occupancy (logit-scale):
 Estimate    SE    z P(>|z|)
    0.253 0.209 1.21   0.227

Detection (logit-scale):
 Estimate    SE     z  P(>|z|)
   -0.803 0.143 -5.62 1.88e-08

AIC: 1084.622 
Number of sites: 400
optim convergence code: 0
optim iterations: 15 
Bootstrap iterations: 0 
backTransform(m1, type="state")
Backtransformed linear combination(s) of Occupancy estimate(s)

 Estimate     SE LinComb (Intercept)
    0.563 0.0515   0.253           1

Transformation: logistic 
backTransform(m1, type="det")
Backtransformed linear combination(s) of Detection estimate(s)

 Estimate     SE LinComb (Intercept)
    0.309 0.0305  -0.803           1

Transformation: logistic

Yes, it does! Our estimated occupancy (0.452) and detection probabilities (0.328) are close (given the standard errors) to the values used in the simulation. But does our model adequately reflect the structure in our data? We can check that with a Goodness-of-fit test using the parametric bootstrap routing in `unmarked’:

# goodness of fit
chisq <- function(fm) {  # calculates chi-square value
  observed <- getY(fm@data)
  expected <- fitted(fm)
  sum((observed - expected)^2/expected)
  }
pb <- parboot(m1, statistic=chisq, nsim=100)
plot(pb, main = "")

Now, we can fit the same model in a Bayesian analysis using JAGS:

# code the model in BUGS
sink("occ.txt")
cat("
    model {
    # Ecological model for true occurrence
    for (i in 1:nsites) {
    z[i] ~ dbern(psi)
    # Observation model for replicated detection/nondetection observations
    p.eff[i] <- z[i] * p
    for (j in 1:nvisits) {
    y[i,j] ~ dbern(p.eff[i])
    } #j
    } #i
    # Priors
    psi ~ dunif(0, 1)
    p ~ dunif(0, 1)
    # Derived quantities
    occ.fs <- sum(z[])
    }
    ",fill = TRUE)
sink()
library(jagsUI)
occ.data <- list(y = y, nsites = nrow(y), nvisits = ncol(y))
# Initial values
zst <- apply(y, 1, max)
# Observed occurrence as starting values for z
inits <- function() list(z = zst)
# Parameters monitored
params <- c("psi", "p", "occ.fs")
# MCMC settings
nc <- 3   ;   ni <- 5000   ;   nb <- 2000   ;   nt <- 2
out.occ <- jags(occ.data, inits, params, "occ.txt", n.chains=nc, n.iter=ni, n.burn = nb, n.thin=nt)

Processing function input....... 

Converting data frame 'y' to matrix.

Done. 
 
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 1200
   Unobserved stochastic nodes: 402
   Total graph size: 2008

Initializing model

Adaptive phase..... 
Adaptive phase complete 
 

 Burn-in phase, 2000 iterations x 3 chains 
 

  |                                                        
  |                                                  |   0%
  |                                                        
  |*                                                 |   2%
  |                                                        
  |**                                                |   4%
  |                                                        
  |***                                               |   6%
  |                                                        
  |****                                              |   8%
  |                                                        
  |*****                                             |  10%
  |                                                        
  |******                                            |  12%
  |                                                        
  |*******                                           |  14%
  |                                                        
  |********                                          |  16%
  |                                                        
  |*********                                         |  18%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |***********                                       |  22%
  |                                                        
  |************                                      |  24%
  |                                                        
  |*************                                     |  26%
  |                                                        
  |**************                                    |  28%
  |                                                        
  |***************                                   |  30%
  |                                                        
  |****************                                  |  32%
  |                                                        
  |*****************                                 |  34%
  |                                                        
  |******************                                |  36%
  |                                                        
  |*******************                               |  38%
  |                                                        
  |********************                              |  40%
  |                                                        
  |*********************                             |  42%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |***********************                           |  46%
  |                                                        
  |************************                          |  48%
  |                                                        
  |*************************                         |  50%
  |                                                        
  |**************************                        |  52%
  |                                                        
  |***************************                       |  54%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |*****************************                     |  58%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |*******************************                   |  62%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |*********************************                 |  66%
  |                                                        
  |**********************************                |  68%
  |                                                        
  |***********************************               |  70%
  |                                                        
  |************************************              |  72%
  |                                                        
  |*************************************             |  74%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |***************************************           |  78%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |*****************************************         |  82%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |*******************************************       |  86%
  |                                                        
  |********************************************      |  88%
  |                                                        
  |*********************************************     |  90%
  |                                                        
  |**********************************************    |  92%
  |                                                        
  |***********************************************   |  94%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |************************************************* |  98%
  |                                                        
  |**************************************************| 100%

Sampling from joint posterior, 3000 iterations x 3 chains 
 

  |                                                        
  |                                                  |   0%
  |                                                        
  |*                                                 |   2%
  |                                                        
  |**                                                |   4%
  |                                                        
  |***                                               |   6%
  |                                                        
  |****                                              |   8%
  |                                                        
  |*****                                             |  10%
  |                                                        
  |******                                            |  12%
  |                                                        
  |*******                                           |  14%
  |                                                        
  |********                                          |  16%
  |                                                        
  |*********                                         |  18%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |***********                                       |  22%
  |                                                        
  |************                                      |  24%
  |                                                        
  |*************                                     |  26%
  |                                                        
  |**************                                    |  28%
  |                                                        
  |***************                                   |  30%
  |                                                        
  |****************                                  |  32%
  |                                                        
  |*****************                                 |  34%
  |                                                        
  |******************                                |  36%
  |                                                        
  |*******************                               |  38%
  |                                                        
  |********************                              |  40%
  |                                                        
  |*********************                             |  42%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |***********************                           |  46%
  |                                                        
  |************************                          |  48%
  |                                                        
  |*************************                         |  50%
  |                                                        
  |**************************                        |  52%
  |                                                        
  |***************************                       |  54%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |*****************************                     |  58%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |*******************************                   |  62%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |*********************************                 |  66%
  |                                                        
  |**********************************                |  68%
  |                                                        
  |***********************************               |  70%
  |                                                        
  |************************************              |  72%
  |                                                        
  |*************************************             |  74%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |***************************************           |  78%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |*****************************************         |  82%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |*******************************************       |  86%
  |                                                        
  |********************************************      |  88%
  |                                                        
  |*********************************************     |  90%
  |                                                        
  |**********************************************    |  92%
  |                                                        
  |***********************************************   |  94%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |************************************************* |  98%
  |                                                        
  |**************************************************| 100%

Calculating statistics....... 

Done. 
plot(out.occ)

print(out.occ$summary, digits=2)
           mean     sd   2.5%    25%    50%    75%  97.5% Rhat n.eff overlap0 f
psi        0.57  0.053   0.47   0.53   0.57   0.60   0.68    1  4500        0 1
p          0.31  0.030   0.25   0.29   0.31   0.33   0.37    1  4500        0 1
occ.fs   227.49 18.402 196.00 214.00 226.00 239.00 267.00    1  4500        0 1
deviance 839.95 39.978 765.70 812.55 838.13 866.15 921.63    1  4500        0 1

The Bayesian analysis returns estimates that are very similar to the frequentist estimates. This is as it should be since we used non-informative priors.

Occupancy model with site-level covariates

Ok, this was a very simple occupancy model, with constant occupancy and detection probabilities. Usually, we are interested in how occupancy varies among sites. I next simulate a data set where occupancy probabilities increase from left to right – or from west to east, if you wish. At the same time, detection probabilities decrease from west to east. So when we look at the raw detection data, it is not obvious that there are any trends in these data at all. The trends in the ecological and observation process cancel each other out. Let’s see whether we can disentangle the two effects using an occupancy model.

library(boot)
set.seed(246)
grid <- expand.grid(long=c(1:20),lat=c(1:20))
occ.prob <- inv.logit(-2+0.2*grid$long) # gradient in occupancy probability
det.prob <- inv.logit(2-0.2*grid$long) # opposite gradient in detection probability
grid$occ <- rbinom(dim(grid)[1],1,prob = occ.prob) # occupancy
grid$col <- ifelse(grid$occ==1,'orange','grey')
grid$detected1 <- rbinom(dim(grid)[1],1,prob = grid$occ * det.prob) # detections on first survey
grid$detected2 <- rbinom(dim(grid)[1],1,prob = grid$occ * det.prob) # detections on second survey
grid$detected3 <- rbinom(dim(grid)[1],1,prob = grid$occ * det.prob) # detections on third survey
# plot the true occupancy and detections during each of the three surveys
op <- par(mfrow=c(2,2), mar=c(1,1,1,1))
plot(lat~long,col=col,pch= 15, data=grid,ann=F,axes=F,cex=1.5)
plot(lat~long,col='grey',pch= 15, data=grid,ann=F,axes=F,cex=1.5)
points(lat~long,col=detected1,pch= 15, data=grid[grid$detected1==1,],cex=1.5)
plot(lat~long,col='grey',pch= 15, data=grid,ann=F,axes=F,cex=1.5)
points(lat~long,col=detected2,pch= 15, data=grid[grid$detected2==1,],cex=1.5)
plot(lat~long,col='grey',pch= 15, data=grid,ann=F,axes=F,cex=1.5)
points(lat~long,col=detected3,pch= 15, data=grid[grid$detected3==1,],cex=1.5)
par(op)
# the true gradients look like this
op <- par(mfrow=c(2,1), mar=c(1,4,1,1), oma=c(3,0,0,0))

plot(1:20, inv.logit(-2+0.2*1:20), ylim=c(0,1), type='l', axes=F, ylab='true occupancy', lwd=2, xlim=c(0,20))  # true occupancy
axis(2, at=c(0,0.5,1), las=1)
plot(1:20, inv.logit(2-0.2*1:20), ylim=c(0,1), type='l', axes=F, ylab='true detection', lwd=2,xlim=c(0,20))  # true detection probability
axis(2, at=c(0,0.5,1), las=1)
axis(1, at=c(0,10,20))
mtext("Longitude", 1, 1.5, outer=T)
par(op)

The figures above show the gradients in occupancy and detection. They run in opposite directions and therefore we can’t see a clear pattern in the raw data.

We analyse this data set using `unmarked’:

require(unmarked)
y <- grid[,5:7]
siteCovs <- data.frame(long=grid[,"long"])
ex1 <- unmarkedFrameOccu(y = y, siteCovs = siteCovs)
summary(ex1)
unmarkedFrame Object

400 sites
Maximum number of observations per site: 3 
Mean number of observations per site: 3 
Sites with at least one detection: 141 

Tabulation of y observations:
  0   1 
977 223 

Site-level covariates:
      long      
 Min.   : 1.00  
 1st Qu.: 5.75  
 Median :10.50  
 Mean   :10.50  
 3rd Qu.:15.25  
 Max.   :20.00  
m1 <- occu(~long ~long, data=ex1)
summary(m1)

Call:
occu(formula = ~long ~ long, data = ex1)

Occupancy (logit-scale):
            Estimate     SE     z  P(>|z|)
(Intercept)   -2.732 0.3872 -7.05 1.73e-12
long           0.282 0.0474  5.95 2.75e-09

Detection (logit-scale):
            Estimate     SE    z  P(>|z|)
(Intercept)    1.794 0.3147  5.7 1.19e-08
long          -0.182 0.0223 -8.2 2.44e-16

AIC: 992.3903 
Number of sites: 400
optim convergence code: 0
optim iterations: 52 
Bootstrap iterations: 0

This model correctly returns the gradients: a positive effect of longitude on occupancy and a negative one on detection. In fact, the estimated values are close to the ones we used in the simulations (\(logit(\Psi) = -2 + 0.2 \times longitude\) and \(logit(p) = 2 - 0.2 \times longitude\)).

But this model was of course the true generating model, which we knew in this case. Let’s fit some simpler models and see whether model selection picks the correct model:

m2 <- occu(~1 ~long, data=ex1)
summary(m2)

Call:
occu(formula = ~1 ~ long, data = ex1)

Occupancy (logit-scale):
            Estimate     SE     z  P(>|z|)
(Intercept)   -1.368 0.2796 -4.89 9.98e-07
long           0.105 0.0249  4.21 2.57e-05

Detection (logit-scale):
 Estimate    SE     z P(>|z|)
   -0.309 0.131 -2.36  0.0184

AIC: 1047.358 
Number of sites: 400
optim convergence code: 0
optim iterations: 26 
Bootstrap iterations: 0 
m3 <- occu(~1 ~1, data=ex1)
summary(m3)

Call:
occu(formula = ~1 ~ 1, data = ex1)

Occupancy (logit-scale):
 Estimate    SE     z P(>|z|)
   -0.269 0.137 -1.97  0.0487

Detection (logit-scale):
 Estimate    SE     z P(>|z|)
   -0.286 0.127 -2.24  0.0249

AIC: 1066.687 
Number of sites: 400
optim convergence code: 0
optim iterations: 23 
Bootstrap iterations: 0 
ml <- fitList("psi(long)p(long)" = m1, "psi(long)p(.)"=m2, "psi(.)p(.)" = m3)
modSel(ml)
                 nPars     AIC delta   AICwt cumltvWt
psi(long)p(long)     4  992.39  0.00 1.0e+00     1.00
psi(long)p(.)        3 1047.36 54.97 1.2e-12     1.00
psi(.)p(.)           2 1066.69 74.30 7.4e-17     1.00

It does. We have clear evidence for model m1 (including the effect of longitude on both \(\Psi\) and \(p\)) to be the best model in the set.

What do the fitted relationships look like? We use the model to predict occupancy and detection as a function of longitude.

backTransform(linearComb(m1, coefficients = c(1,0), type = "state"))  # occupancy at long = 0
Backtransformed linear combination(s) of Occupancy estimate(s)

 Estimate     SE LinComb (Intercept) long
   0.0611 0.0222   -2.73           1    0

Transformation: logistic 
backTransform(linearComb(m1, coefficients = c(1,20), type = "state"))  # occupancy at long = 20
Backtransformed linear combination(s) of Occupancy estimate(s)

 Estimate     SE LinComb (Intercept) long
    0.948 0.0306     2.9           1   20

Transformation: logistic 
newData <- data.frame(long=1:20)
predict(m1, type = "state", newdata = newData, appendData=TRUE)
pr.s <- predict(m1, type = "state", newdata = newData, appendData=TRUE)
pr.d <- predict(m1, type = "det", newdata = newData, appendData=TRUE)
op <- par(mfrow=c(2,1), mar=c(1,4,1,1), oma=c(3,0,0,0))
plot(1:20, inv.logit(-2+0.2*1:20), ylim=c(0,1), type='l', axes=F, ylab='true occupancy', lwd=2, xlim=c(0,20))  # true occupancy
points(Predicted ~ long, data=pr.s)
axis(2, at=c(0,0.5,1), las=1)
plot(1:20, inv.logit(2-0.2*1:20), ylim=c(0,1), type='l', axes=F, ylab='true detection', lwd=2, xlim=c(0,20))  # true detection probability
points(Predicted ~ long, data=pr.d)
axis(2, at=c(0,0.5,1), las=1)
axis(1, at=c(0,10,20))
mtext("Longitude", 1, 1.5, outer=T)
par(op)

We see that the fitted relationships (small circles) closely follow the true relationships (solid lines). The predicted occupancy probabilities for our simulated landscape correctly describe the gradient from west to east and look like this:

newData <- data.frame(long = grid$long)
grid$Predictedm1 <- predict(m1, type = "state", newdata = newData)$Predicted
op <- par(mfcol=c(1,2), mar=c(1,1,1,1))
plot(lat~long,col=col,pch= 15, data=grid,ann=F,axes=F,cex=1.5)
nc <- as.numeric(cut(grid$Predictedm1,20))
plot(lat~long,pch= 15, data=grid,ann=F,axes=F,cex=1.5,col=colorRampPalette(c("white", "darkorange3"))(20)[nc])

And now analyse the same model in JAGS:

sink("occ_cov.txt")
cat("
    model {
    # Ecological model for true occurrence
    for (i in 1:nsites) {
    z[i] ~ dbern(psi[i])
    
    # Observation model for replicated detection/nondetection observations
    p.eff[i] <- z[i] * p[i]
    for (j in 1:nvisits) {
    y[i,j] ~ dbern(p.eff[i])
    } #j
    } #i
    
    # covariates
    for (i in 1:nsites){
      logit(psi[i]) = beta.psi[1] + beta.psi[2] * long[i]
      logit(p[i]) = beta.p[1] + beta.p[2] * long[i]
    }
    # Priors
    for (b in 1:2){
      beta.psi[b] ~ dnorm(0, 0.01)
      beta.p[b] ~ dnorm(0, 0.01)
    }
    
    # Derived quantities
    occ.fs <- sum(z[])
    }
    ",fill = TRUE)
sink()
library(jagsUI)
occ.data <- list(y = grid[,5:7], nsites = nrow(y), nvisits = ncol(y), long=grid[,"long"])
# Initial values
zst <- apply(y, 1, max)
# Observed occurrence as starting values for z
inits <- function() list(z = zst, beta.psi=runif(2,-3,3), beta.p=runif(2,-3,3))
# Parameters monitored
params <- c("beta.psi", "beta.p", "occ.fs")
# MCMC settings
nc <- 3   ;   ni <- 5000   ;   nb <- 2000   ;   nt <- 2
out.occ.cov <- jags(occ.data, inits, params, "occ_cov.txt", n.chains=nc, n.iter=ni, n.burn = nb, n.thin=nt)

Processing function input....... 

Converting data frame 'y' to matrix.

Done. 
 
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 1200
   Unobserved stochastic nodes: 404
   Total graph size: 2530

Initializing model
Deleting model
Error in jags.model(file = model.file, data = data, inits = inits, n.chains = n.chains,  : 
  Error in node z[20]
Node inconsistent with parents

And plotting the fitted relationships:

That’s it.

An example with a real data set the SABAP2 data for bald ibis

Above, we looked at simulated data, which has the advantage that we know the truth and so we can easily check that our model does the correct thing. Usually, of course, we don’t know truth but are interested in extracting occupancy patterns from real data. Here is an example using data from the Second Southern African Bird Atlas Project. I chose the bald ibis as an example.

bi.m <- read.csv("baldibisdata.csv")
# start occupancy analysis
require(unmarked)
write.csv(bi.m[,c("Pentad","jday","Spp","lhours")],"bi.csv",row.names = F)
bi.umf <- csvToUMF("bi.csv",long=TRUE, type = "unmarkedFrameOccu")
summary(bi.umf)
unmarkedFrame Object

3220 sites
Maximum number of observations per site: 719 
Mean number of observations per site: 9.2 
Sites with at least one detection: 370 

Tabulation of y observations:
      0       1    <NA> 
  28597    1017 2285566 

Observation-level covariates:
     lhours          JulianDate     
 Min.   :0.7       Min.   :  1.0    
 1st Qu.:0.7       1st Qu.:181.0    
 Median :1.1       Median :356.0    
 Mean   :1.1       Mean   :358.6    
 3rd Qu.:1.4       3rd Qu.:537.0    
 Max.   :3.9       Max.   :731.0    
 NA's   :2285566   NA's   :2285566  
# have a look at the observation-level covariate
head(obsCovs(bi.umf,matrices=T)$lhours)
          [,1]      [,2]     [,3]     [,4]     [,5]      [,6]     [,7]     [,8]      [,9]    [,10]    [,11]     [,12]    [,13]     [,14]
[1,] 0.6931472        NA       NA       NA       NA        NA       NA       NA        NA       NA       NA        NA       NA        NA
     [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
[1,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
     [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58]
[1,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
     [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80]
[1,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
     [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100] [,101]
[1,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA     NA     NA
     [,102] [,103] [,104] [,105] [,106] [,107] [,108] [,109] [,110] [,111] [,112] [,113] [,114] [,115] [,116] [,117] [,118] [,119] [,120]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,121] [,122] [,123] [,124] [,125] [,126] [,127] [,128] [,129] [,130] [,131] [,132] [,133] [,134] [,135] [,136] [,137] [,138] [,139]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,140] [,141] [,142] [,143] [,144] [,145] [,146] [,147] [,148] [,149] [,150] [,151] [,152] [,153] [,154] [,155] [,156] [,157] [,158]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,159] [,160] [,161] [,162] [,163] [,164] [,165] [,166] [,167] [,168] [,169] [,170] [,171] [,172] [,173] [,174] [,175] [,176] [,177]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,178] [,179] [,180] [,181] [,182] [,183] [,184] [,185] [,186] [,187] [,188] [,189] [,190] [,191] [,192] [,193] [,194] [,195] [,196]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,197] [,198] [,199] [,200] [,201] [,202] [,203] [,204] [,205] [,206] [,207] [,208] [,209] [,210] [,211] [,212] [,213] [,214] [,215]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,216] [,217] [,218] [,219] [,220] [,221] [,222] [,223] [,224] [,225] [,226] [,227] [,228] [,229] [,230] [,231] [,232] [,233] [,234]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,235] [,236] [,237] [,238] [,239] [,240] [,241] [,242] [,243] [,244] [,245] [,246] [,247] [,248] [,249] [,250] [,251] [,252] [,253]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,254] [,255] [,256] [,257] [,258] [,259] [,260] [,261] [,262] [,263] [,264] [,265] [,266] [,267] [,268] [,269] [,270] [,271] [,272]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,273] [,274] [,275] [,276] [,277] [,278] [,279] [,280] [,281] [,282] [,283] [,284] [,285] [,286] [,287] [,288] [,289] [,290] [,291]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,292] [,293] [,294] [,295] [,296] [,297] [,298] [,299] [,300] [,301] [,302] [,303] [,304] [,305] [,306] [,307] [,308] [,309] [,310]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,311] [,312] [,313] [,314] [,315] [,316] [,317] [,318] [,319] [,320] [,321] [,322] [,323] [,324] [,325] [,326] [,327] [,328] [,329]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,330] [,331] [,332] [,333] [,334] [,335] [,336] [,337] [,338] [,339] [,340] [,341] [,342] [,343] [,344] [,345] [,346] [,347] [,348]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,349] [,350] [,351] [,352] [,353] [,354] [,355] [,356] [,357] [,358] [,359] [,360] [,361] [,362] [,363] [,364] [,365] [,366] [,367]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,368] [,369] [,370] [,371] [,372] [,373] [,374] [,375] [,376] [,377] [,378] [,379] [,380] [,381] [,382] [,383] [,384] [,385] [,386]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,387] [,388] [,389] [,390] [,391] [,392] [,393] [,394] [,395] [,396] [,397] [,398] [,399] [,400] [,401] [,402] [,403] [,404] [,405]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,406] [,407] [,408] [,409] [,410] [,411] [,412] [,413] [,414] [,415] [,416] [,417] [,418] [,419] [,420] [,421] [,422] [,423] [,424]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,425] [,426] [,427] [,428] [,429] [,430] [,431] [,432] [,433] [,434] [,435] [,436] [,437] [,438] [,439] [,440] [,441] [,442] [,443]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,444] [,445] [,446] [,447] [,448] [,449] [,450] [,451] [,452] [,453] [,454] [,455] [,456] [,457] [,458] [,459] [,460] [,461] [,462]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,463] [,464] [,465] [,466] [,467] [,468] [,469] [,470] [,471] [,472] [,473] [,474] [,475] [,476] [,477] [,478] [,479] [,480] [,481]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,482] [,483] [,484] [,485] [,486] [,487] [,488] [,489] [,490] [,491] [,492] [,493] [,494] [,495] [,496] [,497] [,498] [,499] [,500]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,501] [,502] [,503] [,504] [,505] [,506] [,507] [,508] [,509] [,510] [,511] [,512] [,513] [,514] [,515] [,516] [,517] [,518] [,519]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,520] [,521] [,522] [,523] [,524] [,525] [,526] [,527] [,528] [,529] [,530] [,531] [,532] [,533] [,534] [,535] [,536] [,537] [,538]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,539] [,540] [,541] [,542] [,543] [,544] [,545] [,546] [,547] [,548] [,549] [,550] [,551] [,552] [,553] [,554] [,555] [,556] [,557]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,558] [,559] [,560] [,561] [,562] [,563] [,564] [,565] [,566] [,567] [,568] [,569] [,570] [,571] [,572] [,573] [,574] [,575] [,576]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,577] [,578] [,579] [,580] [,581] [,582] [,583] [,584] [,585] [,586] [,587] [,588] [,589] [,590] [,591] [,592] [,593] [,594] [,595]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,596] [,597] [,598] [,599] [,600] [,601] [,602] [,603] [,604] [,605] [,606] [,607] [,608] [,609] [,610] [,611] [,612] [,613] [,614]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,615] [,616] [,617] [,618] [,619] [,620] [,621] [,622] [,623] [,624] [,625] [,626] [,627] [,628] [,629] [,630] [,631] [,632] [,633]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,634] [,635] [,636] [,637] [,638] [,639] [,640] [,641] [,642] [,643] [,644] [,645] [,646] [,647] [,648] [,649] [,650] [,651] [,652]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,653] [,654] [,655] [,656] [,657] [,658] [,659] [,660] [,661] [,662] [,663] [,664] [,665] [,666] [,667] [,668] [,669] [,670] [,671]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,672] [,673] [,674] [,675] [,676] [,677] [,678] [,679] [,680] [,681] [,682] [,683] [,684] [,685] [,686] [,687] [,688] [,689] [,690]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,691] [,692] [,693] [,694] [,695] [,696] [,697] [,698] [,699] [,700] [,701] [,702] [,703] [,704] [,705] [,706] [,707] [,708] [,709]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
     [,710] [,711] [,712] [,713] [,714] [,715] [,716] [,717] [,718] [,719]
[1,]     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
 [ reached getOption("max.print") -- omitted 5 rows ]
# add site-level covariates (for some reason, this didn't work automatically with 'csvToUMF')
siteCovs(bi.umf) <- data.frame(GDD0=aggregate(bi.m$GDD0,list(bi.m$Pentad),mean)$x,
                           GDD5=aggregate(bi.m$GDD5,list(bi.m$Pentad),mean)$x,
                           MTCO=aggregate(bi.m$MTCO,list(bi.m$Pentad),mean)$x,
                           MTWA=aggregate(bi.m$MTWA,list(bi.m$Pentad),mean)$x,
                           AET.PET=aggregate(bi.m$AET.PET,list(bi.m$Pentad),mean)$x,
                           Wet.Intensity=aggregate(bi.m$Wet.Intensity,list(bi.m$Pentad),mean)$x,
                           Dry.Intensity=aggregate(bi.m$Dry.Intensity,list(bi.m$Pentad),mean)$x)
summary(bi.umf)
unmarkedFrame Object

3220 sites
Maximum number of observations per site: 719 
Mean number of observations per site: 9.2 
Sites with at least one detection: 370 

Tabulation of y observations:
      0       1    <NA> 
  28597    1017 2285566 

Site-level covariates:
      GDD0           GDD5           MTCO            MTWA          AET.PET       Wet.Intensity    Dry.Intensity    
 Min.   :2254   Min.   : 820   Min.   : 0.70   Min.   :10.80   Min.   :0.2370   Min.   :  0.00   Min.   :-108.60  
 1st Qu.:5642   1st Qu.:3817   1st Qu.: 8.60   1st Qu.:20.10   1st Qu.:0.4430   1st Qu.: 19.20   1st Qu.: -58.90  
 Median :6186   Median :4361   Median :10.90   Median :21.60   Median :0.5370   Median : 37.60   Median : -48.50  
 Mean   :6343   Mean   :4520   Mean   :11.21   Mean   :21.88   Mean   :0.5606   Mean   : 46.42   Mean   : -45.39  
 3rd Qu.:7089   3rd Qu.:5264   3rd Qu.:13.50   3rd Qu.:23.60   3rd Qu.:0.6870   3rd Qu.: 68.30   3rd Qu.: -31.40  
 Max.   :8437   Max.   :6612   Max.   :18.20   Max.   :27.40   Max.   :0.9930   Max.   :183.10   Max.   :   0.00  

Observation-level covariates:
     lhours          JulianDate     
 Min.   :0.7       Min.   :  1.0    
 1st Qu.:0.7       1st Qu.:181.0    
 Median :1.1       Median :356.0    
 Mean   :1.1       Mean   :358.6    
 3rd Qu.:1.4       3rd Qu.:537.0    
 Max.   :3.9       Max.   :731.0    
 NA's   :2285566   NA's   :2285566  
m1 <- occu(~ lhours ~ 1, data = bi.umf)
summary(m1)

Call:
occu(formula = ~lhours ~ 1, data = bi.umf)

Occupancy (logit-scale):
 Estimate     SE     z  P(>|z|)
    -1.27 0.0643 -19.7 1.38e-86

Detection (logit-scale):
            Estimate     SE      z  P(>|z|)
(Intercept)  -1.3467 0.1012 -13.31 2.01e-40
lhours        0.0987 0.0846   1.17 2.43e-01

AIC: 6308.95 
Number of sites: 3220
optim convergence code: 0
optim iterations: 28 
Bootstrap iterations: 0 
m2 <- occu(~ lhours ~ MTCO + MTWA + AET.PET + Wet.Intensity, data = bi.umf)
summary(m2)

Call:
occu(formula = ~lhours ~ MTCO + MTWA + AET.PET + Wet.Intensity, 
    data = bi.umf)

Occupancy (logit-scale):
              Estimate     SE      z  P(>|z|)
(Intercept)     -2.637 2.6700 -0.988 3.23e-01
MTCO            -0.388 0.0932 -4.169 3.06e-05
MTWA            -0.316 0.1250 -2.532 1.14e-02
AET.PET         24.908 3.3814  7.366 1.75e-13
Wet.Intensity   -0.057 0.0120 -4.764 1.90e-06

Detection (logit-scale):
            Estimate     SE      z  P(>|z|)
(Intercept)   -1.479 0.1005 -14.71 5.43e-49
lhours         0.163 0.0829   1.96 4.97e-02

AIC: 5657.779 
Number of sites: 3220
optim convergence code: 0
optim iterations: 101 
Bootstrap iterations: 0 
# fitted relationships
# --------------------
op<-par(mfrow=c(2,2), mar=c(4,1,1,2),oma=c(1,5,1,0))
MTCO.p <- data.frame(MTCO = seq(min(bi.m$MTCO), max(bi.m$MTCO), length=100), MTWA = mean(bi.m$MTWA), AET.PET = mean(bi.m$AET.PET), Wet.Intensity = mean(bi.m$Wet.Intensity))
plot(MTCO.p$MTCO,predict(m2,type='state',newdata=MTCO.p)$Predicted, type='l', axes=F, xlab='MTCO',ylim=c(0,1),ylab='')
axis(1,seq(min(bi.m$MTCO), max(bi.m$MTCO), length=3))
axis(2,c(0,0.5,1),las=1)
MTWA.p <- data.frame(MTCO = mean(bi.m$MTCO), MTWA = seq(min(bi.m$MTWA), max(bi.m$MTWA), length=100), AET.PET = mean(bi.m$AET.PET), Wet.Intensity = mean(bi.m$Wet.Intensity))
plot(MTWA.p$MTWA,predict(m2,type='state',newdata=MTWA.p)$Predicted, type='l',axes=F, xlab='MTWA',ylim=c(0,1),ylab='')
axis(1, seq(min(bi.m$MTWA), max(bi.m$MTWA), length=3))
AET.PET.p <- data.frame(MTCO = mean(bi.m$MTCO), MTWA = mean(bi.m$MTWA), AET.PET = seq(min(bi.m$AET.PET), max(bi.m$AET.PET), length=100), Wet.Intensity = mean(bi.m$Wet.Intensity))
plot(AET.PET.p$AET.PET,predict(m2,type='state',newdata=AET.PET.p)$Predicted, type='l',axes=F, xlab='AET/PET',ylim=c(0,1),ylab='')
axis(1,seq(min(bi.m$AET.PET), max(bi.m$AET.PET), length=3))
axis(2,c(0,0.5,1),las=1)
Wet.Intensity.p <- data.frame(MTCO = mean(bi.m$MTCO), MTWA = mean(bi.m$MTWA), AET.PET = mean(bi.m$AET.PET), Wet.Intensity = seq(min(bi.m$Wet.Intensity), max(bi.m$Wet.Intensity), length=100))
plot(Wet.Intensity.p$Wet.Intensity,predict(m2,type='state',newdata=Wet.Intensity.p)$Predicted, type='l',axes=F, xlab='Wet.Intensity',ylim=c(0,1),ylab='')
axis(1,seq(min(bi.m$Wet.Intensity), max(bi.m$Wet.Intensity), length=3))
mtext("Occupancy probability",side=2,line=3,outer=T)
par(op)
# predictions for drawing the map
# -------------------------------
pred <- predict(m2, type='state', newData=siteCovs(bi.umf)[,3:6], appendData=F)
pred$Pentad <- aggregate(bi.m$GDD0,list(bi.m$Pentad),mean)$Group.1
# coordinates
pred$lat <- -(as.numeric(substr(pred$Pentad,1,2)) + (as.numeric(substr(pred$Pentad,3,4))+2.5)/60)
pred$long <- as.numeric(substr(pred$Pentad,6,7)) + (as.numeric(substr(pred$Pentad,8,9))+2.5)/60
# plot the bald ibis map
library(rgdal)
Loading required package: sp

Attaching package: ‘sp’

The following object is masked from ‘package:unmarked’:

    coordinates

rgdal: version: 1.2-16, (SVN revision 701)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
 Path to GDAL shared files: /usr/share/gdal/2.2
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.2-5 
ZAmap <- readOGR(dsn = getwd(), layer ="RSA_Provinces_WGS84") # map of South Africa
OGR data source with driver: ESRI Shapefile 
Source: "/home/res/Dropbox/res/centre for statistical ecology/SEEC toolbox seminars/2018_05_31 Bayesian Occupancy", layer: "RSA_Provinces_WGS84"
with 11 features
It has 10 fields
Integer64 fields read as strings:  hectares 
proj4string(ZAmap)
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
coordinates(pred)= ~lat+long
proj4string(pred) <-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
nc <- as.numeric(cut(pred$Predicted,20))
par(mar=c(1,1,1,1))

plot(lat~long, data=pred, pch=15, cex=0.6, ann=F, axes=F, xlim=c(15,33), ylim = c(-35, -21), col=colorRampPalette(c("grey", "darkorange3"))(20)[nc])
plot(ZAmap, add=T)
par(mar=c(0,0,0,0))

plot(lat~long, data=bi.m, pch=15, cex=0.6, ann=F, axes=F, xlim=c(15,33), ylim = c(-35, -21), col="grey")
points(lat~long, data=bi.m[bi.m$Spp==1,], pch=15, cex=0.6,col='black')
plot(ZAmap, add=T)

# code the model in BUGS
sink("occ_ibis.txt")
cat("
    model {
    # Ecological model for true occurrence
    for (i in 1:nsites) {
      z[i] ~ dbern(psi[i])
    
      # Observation model for replicated detection/nondetection observations
      for (j in 1:ncards[i]) {
        p.eff[i,j] <- z[i] * p[i,j]
        y[i,j] ~ dbern(p.eff[i,j])
      } #j
    } #i
    
    # covariates
    for (i in 1:nsites){
    logit(psi[i]) = beta.psi[1] + beta.psi[2] * MTCO[i] + beta.psi[3] * MTWA[i] + beta.psi[4] * AET.PET[i] + beta.psi[5] * Wet.Intensity[i]
      for (j in 1:ncards[i]){
        logit(p[i,j]) = beta.p[1] + beta.p[2] * lhours[i,j]
      } #j
    } #i
  
    # Priors
    for (b in 1:5){
      beta.psi[b] ~ dnorm(0, 0.01)
    }
    for (b in 1:2){
      beta.p[b] ~ dnorm(0, 0.01)
    }
    }
    ",fill = TRUE)
sink()
library(jagsUI)
y <- bi.umf@y
dim(y)  # 3220 sites, max 719 cards
[1] 3220  719
ncards <- as.numeric(table(bi.m$Pentad))
MTCO=scale(aggregate(bi.m$MTCO,list(bi.m$Pentad),mean)$x)
MTWA=scale(aggregate(bi.m$MTWA,list(bi.m$Pentad),mean)$x)
AET.PET=scale(aggregate(bi.m$AET.PET,list(bi.m$Pentad),mean)$x)
Wet.Intensity=scale(aggregate(bi.m$Wet.Intensity,list(bi.m$Pentad),mean)$x)
occ.data <- list(y = y, nsites = nrow(y), ncards = ncards,
                 lhours=matrix(bi.umf@obsCovs$lhours, nrow=nrow(y), byrow=T),
                 MTCO=as.numeric(MTCO),
                 MTWA=as.numeric(MTWA),
                 AET.PET=as.numeric(AET.PET),
                 Wet.Intensity=as.numeric(Wet.Intensity))
# Initial values
zst <- apply(y, 1, function(x){max(x,na.rm=T)})
# Observed occurrence as starting values for z
inits <- function() list(z = zst, beta.psi=runif(5,-3,3), beta.p=runif(2,-3,3))
# Parameters monitored
params <- c("beta.psi", "beta.p")
# MCMC settings
nc <- 3   ;   ni <- 5000   ;   nb <- 2000   ;   nt <- 2
out.occ.cov <- jags(occ.data, inits, params, "occ_ibis.txt", n.chains=nc, n.iter=ni, n.burn = nb, n.thin=nt)

Processing function input....... 

Done. 
 
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 29614
   Unobserved stochastic nodes: 3227
   Total graph size: 88133

Initializing model

Adaptive phase..... 
Adaptive phase complete 
 

 Burn-in phase, 2000 iterations x 3 chains 
 

  |                                                        
  |                                                  |   0%
  |                                                        
  |*                                                 |   2%
  |                                                        
  |**                                                |   4%
  |                                                        
  |***                                               |   6%
  |                                                        
  |****                                              |   8%
  |                                                        
  |*****                                             |  10%
  |                                                        
  |******                                            |  12%
  |                                                        
  |*******                                           |  14%
  |                                                        
  |********                                          |  16%
  |                                                        
  |*********                                         |  18%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |***********                                       |  22%
  |                                                        
  |************                                      |  24%
  |                                                        
  |*************                                     |  26%
  |                                                        
  |**************                                    |  28%
  |                                                        
  |***************                                   |  30%
  |                                                        
  |****************                                  |  32%
  |                                                        
  |*****************                                 |  34%
  |                                                        
  |******************                                |  36%
  |                                                        
  |*******************                               |  38%
  |                                                        
  |********************                              |  40%
  |                                                        
  |*********************                             |  42%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |***********************                           |  46%
  |                                                        
  |************************                          |  48%
  |                                                        
  |*************************                         |  50%
  |                                                        
  |**************************                        |  52%
  |                                                        
  |***************************                       |  54%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |*****************************                     |  58%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |*******************************                   |  62%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |*********************************                 |  66%
  |                                                        
  |**********************************                |  68%
  |                                                        
  |***********************************               |  70%
  |                                                        
  |************************************              |  72%
  |                                                        
  |*************************************             |  74%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |***************************************           |  78%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |*****************************************         |  82%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |*******************************************       |  86%
  |                                                        
  |********************************************      |  88%
  |                                                        
  |*********************************************     |  90%
  |                                                        
  |**********************************************    |  92%
  |                                                        
  |***********************************************   |  94%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |************************************************* |  98%
  |                                                        
  |**************************************************| 100%

Sampling from joint posterior, 3000 iterations x 3 chains 
 

  |                                                        
  |                                                  |   0%
  |                                                        
  |*                                                 |   2%
  |                                                        
  |**                                                |   4%
  |                                                        
  |***                                               |   6%
  |                                                        
  |****                                              |   8%
  |                                                        
  |*****                                             |  10%
  |                                                        
  |******                                            |  12%
  |                                                        
  |*******                                           |  14%
  |                                                        
  |********                                          |  16%
  |                                                        
  |*********                                         |  18%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |***********                                       |  22%
  |                                                        
  |************                                      |  24%
  |                                                        
  |*************                                     |  26%
  |                                                        
  |**************                                    |  28%
  |                                                        
  |***************                                   |  30%
  |                                                        
  |****************                                  |  32%
  |                                                        
  |*****************                                 |  34%
  |                                                        
  |******************                                |  36%
  |                                                        
  |*******************                               |  38%
  |                                                        
  |********************                              |  40%
  |                                                        
  |*********************                             |  42%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |***********************                           |  46%
  |                                                        
  |************************                          |  48%
  |                                                        
  |*************************                         |  50%
  |                                                        
  |**************************                        |  52%
  |                                                        
  |***************************                       |  54%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |*****************************                     |  58%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |*******************************                   |  62%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |*********************************                 |  66%
  |                                                        
  |**********************************                |  68%
  |                                                        
  |***********************************               |  70%
  |                                                        
  |************************************              |  72%
  |                                                        
  |*************************************             |  74%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |***************************************           |  78%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |*****************************************         |  82%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |*******************************************       |  86%
  |                                                        
  |********************************************      |  88%
  |                                                        
  |*********************************************     |  90%
  |                                                        
  |**********************************************    |  92%
  |                                                        
  |***********************************************   |  94%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |************************************************* |  98%
  |                                                        
  |**************************************************| 100%

Calculating statistics....... 

Done. 
plot(out.occ.cov)

print(out.occ.cov$summary, digits=2)
               mean     sd     2.5%     25%     50%     75%   97.5% Rhat n.eff overlap0    f
beta.psi[1]   -2.63  0.158 -2.9e+00   -2.73   -2.62   -2.52   -2.32    1   255        0 1.00
beta.psi[2]   -1.34  0.308 -1.9e+00   -1.55   -1.33   -1.13   -0.75    1   781        0 1.00
beta.psi[3]   -0.82  0.322 -1.5e+00   -1.04   -0.82   -0.60   -0.21    1   397        0 1.00
beta.psi[4]    3.89  0.515  2.9e+00    3.55    3.89    4.23    4.91    1  1080        0 1.00
beta.psi[5]   -2.09  0.434 -2.9e+00   -2.39   -2.08   -1.80   -1.25    1  1863        0 1.00
beta.p[1]     -1.48  0.100 -1.7e+00   -1.55   -1.49   -1.42   -1.29    1   810        0 1.00
beta.p[2]      0.16  0.081  4.8e-03    0.11    0.16    0.22    0.32    1   969        0 0.98
deviance    4947.65 37.766  4.9e+03 4921.74 4946.99 4973.50 5022.08    1  2415        0 1.00
# fitted relationships
# --------------------
library(boot)
op<-par(mfrow=c(2,2), mar=c(4,1,1,2),oma=c(1,5,1,0))

MTCO.p <- inv.logit(out.occ.cov$summary[1,1] + out.occ.cov$summary[2,1] * MTCO)
od <- order(MTCO)
plot((MTCO[od]*attributes(MTCO)[[3]])+attributes(MTCO)[[2]],MTCO.p[od], type='l', axes=F, xlab='MTCO',ylim=c(0,1),ylab='',col.axis="orange",col="orange", fg="orange", col.lab="orange")
axis(1,seq(min((MTCO*attributes(MTCO)[[3]])+attributes(MTCO)[[2]]), max((MTCO*attributes(MTCO)[[3]])+attributes(MTCO)[[2]]), length=3),col.axis="orange",col="orange", fg="orange")
axis(2,c(0,0.5,1),las=1,col.axis="orange",col="orange", fg="orange")
MTWA.p <- inv.logit(out.occ.cov$summary[1,1] + out.occ.cov$summary[3,1] * MTWA)
od <- order(MTWA)
plot((MTWA[od]*attributes(MTWA)[[3]])+attributes(MTWA)[[2]],MTWA.p[od], type='l',axes=F, xlab='MTWA',ylim=c(0,1),ylab='',col.axis="orange",col="orange", fg="orange", col.lab="orange")
axis(1, seq(min((MTWA[od]*attributes(MTWA)[[3]])+attributes(MTWA)[[2]]), max((MTWA[od]*attributes(MTWA)[[3]])+attributes(MTWA)[[2]]), length=3),col.axis="orange",col="orange", fg="orange")
AET.PET.p <- inv.logit(out.occ.cov$summary[1,1] + out.occ.cov$summary[4,1] * AET.PET)
od <- order(AET.PET)
plot((AET.PET[od]*attributes(AET.PET)[[3]])+attributes(AET.PET)[[2]],AET.PET.p[od], type='l',axes=F, xlab='AET/PET',ylim=c(0,1),ylab='',col.axis="orange",col="orange", fg="orange", col.lab="orange")
axis(1,seq(min((AET.PET[od]*attributes(AET.PET)[[3]])+attributes(AET.PET)[[2]]), max((AET.PET[od]*attributes(AET.PET)[[3]])+attributes(AET.PET)[[2]]), length=3),col.axis="orange",col="orange", fg="orange")
axis(2,c(0,0.5,1),las=1,col.axis="orange",col="orange", fg="orange")
Wet.Intensity.p <- inv.logit(out.occ.cov$summary[1,1] + out.occ.cov$summary[5,1] * Wet.Intensity)
od <- order(Wet.Intensity)
plot((Wet.Intensity[od]*attributes(Wet.Intensity)[[3]])+attributes(Wet.Intensity)[[2]],Wet.Intensity.p[od], type='l',axes=F, xlab='Wet.Intensity',ylim=c(0,1),ylab='',col.axis="orange",col="orange", fg="orange", col.lab="orange")
axis(1,seq(min((Wet.Intensity[od]*attributes(Wet.Intensity)[[3]])+attributes(Wet.Intensity)[[2]]), max((Wet.Intensity[od]*attributes(Wet.Intensity)[[3]])+attributes(Wet.Intensity)[[2]]), length=3),col.axis="orange",col="orange", fg="orange")
mtext("Occupancy probability",side=2,line=3,outer=T,col.axis="orange",col="orange", fg="orange")
par(op)
# predictions for drawing the map
# -------------------------------
pred <- data.frame(Predicted = inv.logit(out.occ.cov$summary[1,1] + out.occ.cov$summary[2,1] * MTCO + out.occ.cov$summary[3,1] * MTWA +
        out.occ.cov$summary[4,1] * AET.PET + out.occ.cov$summary[5,1] * Wet.Intensity))
pred$Pentad <- aggregate(bi.m$GDD0,list(bi.m$Pentad),mean)$Group.1
# coordinates
pred$lat <- -(as.numeric(substr(pred$Pentad,1,2)) + (as.numeric(substr(pred$Pentad,3,4))+2.5)/60)
pred$long <- as.numeric(substr(pred$Pentad,6,7)) + (as.numeric(substr(pred$Pentad,8,9))+2.5)/60
# plot the bald ibis map
library(rgdal)
ZAmap <- readOGR(dsn = getwd(), layer ="RSA_Provinces_WGS84") # map of South Africa
OGR data source with driver: ESRI Shapefile 
Source: "/home/res/Dropbox/res/centre for statistical ecology/SEEC toolbox seminars/2018_05_31 Bayesian Occupancy", layer: "RSA_Provinces_WGS84"
with 11 features
It has 10 fields
Integer64 fields read as strings:  hectares 
proj4string(ZAmap)
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
coordinates(pred)= ~lat+long
proj4string(pred) <-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
nc <- as.numeric(cut(pred$Predicted,20))
par(mar=c(1,1,1,1))

plot(lat~long, data=pred, pch=15, cex=0.6, ann=F, axes=F, xlim=c(15,33), ylim = c(-35, -21), col=colorRampPalette(c("grey", "darkorange3"))(20)[nc])
plot(ZAmap, add=T)

…

LS0tCnRpdGxlOiAiU2luZ2xlLXNlYXNvbiBvY2N1cGFuY3kgbW9kZWxzIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpUaGlzIFtSIE1hcmtkb3duXShodHRwOi8vcm1hcmtkb3duLnJzdHVkaW8uY29tKSBOb3RlYm9vayBhY2NvbXBhbmllcyB0d28gW1NFRUNdKGh0dHA6Ly93d3cuc2VlYy51Y3QuYWMuemEvKSBbVG9vbGJveCBTZW1pbmFyc10oaHR0cDovL3d3dy5zZWVjLnVjdC5hYy56YS9zdGF0cy10b29sYm94LXNlbWluYXJzKSBnaXZlbiBieSBSZXMgQWx0d2VnZyBhdCB0aGUgVW5pdmVyc2l0eSBvZiBDYXBlIFRvd24sIG9uZSBvbiAyNSBBdWd1c3QgMjAxNiBhbmQgdGhlIG90aGVyIG9uZSBvbiAzMSBNYXkgMjAxOC4KCk9jY3VwYW5jeSBtb2RlbHMgYXJlIGJlY29taW5nIGluY3JlYXNpbmdseSBwb3B1bGFyIGFuZCB0aGVzZSB0d28gc2VtaW5hcnMgd2VyZSBpbnRlbmRlZCB0byBpbnRyb2R1Y2UgdGhlIGJhc2ljIGlkZWFzIG9mIHRoZSBzaW5nbGUtc2Vhc29uIG9jY3VwYW5jeSBtb2RlbC4gSSBzaG93IGhvdyB0byBmaXQgdGhlc2UgbW9kZWxzIGluIGEgZmVxdWVudGlzdCB3YXkgdXNpbmcgdGhlIFIgcGFja2FnZSBgdW5tYXJrZWQnLCBhbmQgaW4gYSBCYXllc2lhbiB3YXkgdXNpbmcgSkFHUy4KCk9jY3VwYW5jeSBjYW4gYmUgYSB1c2VmdWwgc3RhdGUgdmFyaWFibGUgaW4gbWFueSBzaXR1YXRpb25zLiBFLmcuIGlmIHdlIHdhbnQgdG8ga25vdyB3aGVyZSBhIHNwZWNpZXMgb2NjdXJzIGFuZCB3aGVyZSBpdCBkb2Vzbid0IG9jY3VyLCBidXQgd2UgYXJlIGxlc3Mgd29ycmllZCBhYm91dCBkZW5zaXR5IG9mIGluZGl2aWR1YWxzIGluIHBsYWNlcyB3aGVyZSBhIHNwZWNpZXMgZG9lcyBvY2N1ci4gRm9yIGV4YW1wbGUsIG1ldGFwb3B1bGF0aW9uIGVjb2xvZ3kgb2Z0ZW4gZGVhbHMgd2l0aCBwYXRjaCBvY2N1cGFuY3kgYXMgYSBzdGF0ZSB2YXJpYWJsZTsgc29tZSBvZiB0aGUgSVVDTiByZS1saXN0aW5nIGNyaXRlcmlhIGFyZSBiYXNlZCBvbiBvY2N1cGFuY3kgKHRoZSByYW5nZSBvZiBhIHNwZWNpZXMsIGFuZCBhcmVhIG9mIG9jY3VwYW5jeSk7IGFuZCBtYW5hZ2VycyBtYXkgbmVlZCB0byBrbm93IHdoaWNoIHNwZWNpZXMgb2NjdXIgaW4gd2hpY2ggb2YgdGhlIHNpdGVzIHRoZXkgYXJlIG1hbmFnaW5nLgoKT25lIHdheSB0byB0aGluayBvZiBvY2N1cGFuY3kgaXMgdG8gaW1hZ2luZSBhIHNldCBvZiAkbiQgZGlzY3JldGUgc2l0ZXMuIE9jY3VwYW5jeSAkXFBzaSQgaXMgdGhlbiB0aGUgcHJvcG9ydGlvbiBvZiBzaXRlcyB0aGF0IGFyZSBvY2N1cGllZCAkXGZyYWN7b2NjdXBpZWR9e259JC4gSWYgd2Uga25ldyBleGFjdGx5IHdoaWNoIHNpdGVzIGFyZSBvY2N1cGllZCBhbmQgbm90IG9jY3VwaWVkLCB3ZSBjb3VsZCB1c2UgJFxmcmFje29jY3VwaWVkfXtufSQgYXMgYW4gZXN0aW1hdG9yIG9mIG9jY3VwYW5jeS4gQW5kIHdlIGNvdWxkIHVzZSBtZXRob2RzIGxpa2UgbG9naXN0aWMgcmVncmVzc2lvbiwgJGxvZ2l0KFxQc2kpID0gZihjb3ZhcmlhdGVzKSQgdG8gZXhhbWluZSBwYXR0ZXJucyBpbiBvY2N1cGFuY3kuIFRoZSB0cm91YmxlIGlzIHRoYXQgd2UgdXN1YWxseSBvbmx5IGtub3cgd2hlcmUgYSBzcGVjaWVzIHdhcyBkZXRlY3RlZCBvciBub3QuIElmIHdlIHZpc2l0IGV2ZXJ5IHNpbmdsZSBzaXRlIHdlIGFyZSBpbnRlcmVzdGVkIGluIGFuZCByZWNvcmQgd2hldGhlciB3ZSBmaW5kIHRoZSBzcGVjaWVzLCB0aGVyZSBpcyBzdGlsbCBhIGdvb2QgY2hhbmNlIHRoYXQgb3VyIGluZm9ybWF0aW9uIGlzIGluY29tcGxldGUgYmVjYXVzZSB3ZSBtaXNzZWQgdGhlIHNwZWNpZXMgaW4gbWFueSBwbGFjZXMgd2hlcmUgaXQgYWN0dWFsbHkgb2NjdXJzLiBUaGF0IGlzIGJlY2F1c2UgdGhlIGRldGVjdGlvbiBwcm9iYWJpbGl0eSAkcCQgaXMgdXN1YWxseSBsZXNzIHRoYW4gMS4gVGhlIGZyYWN0aW9uIG9mIHNpdGVzIHRoYXQgYXJlIHtcaXRzaGFwZSBvYnNlcnZlZH0gdG8gYmUgb2NjdXBpZWQgaXMgdGh1cyBpbiByZWFsaXR5ICRcUHNpIFx0aW1lcyBwJCwgd2hpY2ggaXMgb25seSBjbG9zdCB0byAkXFBzaSQgaWYgJHAgXHNpbSAxJC4gQW5kIGEgbG9naXN0aWMgcmVncmVzc2lvbiBvbiB0aGVzZSBkZXRlY3Rpb25zIHdvdWxkIGFjdHVhbGx5IG1vZGVsICRsb2dpdChcUHNpIFx0aW1lcyBwKSA9IGYoY292YXJpYXRlcykkLiBBbnkgcGF0dGVybiBmb3VuZCBpbiB0aGlzIHdheSBjb3VsZCBiZSB0aGUgcmVzdWx0IG9mIHZhcmlhdGlvbiBpbiBvY2N1cGFuY3kgb3IgZGV0ZWN0aW9uLCBvciBib3RoLiBXZSBhcmUgaW50ZXJlc3RlZCBpbiBvY2N1cGFuY3kgYW5kIHNvIG5lZWQgdG8gZmluZCBhIHdheSB0byBzZXBhcmF0ZSB0aGUgZGV0ZWN0aW9uIHByb2Nlc3MgZnJvbSB0aGUgb2NjdXBhbmN5IHByb2Nlc3MuIEZvciB0aGlzLCB3ZSBuZWVkIHNvbWUga2luZCBvZiBpbmZvcm1hdGlvbiBvbiB0aGUgb2JzZXJ2YXRpb24gcHJvY2Vzcy4gT25lIHdheSB0byBnZXQgYXQgdGhlIG9ic2VydmF0aW9uIHByb2Nlc3MgaXMgYnkgdmlzaXRpbmcgc2l0ZXMgcmVwZWF0ZWRseSwgZWFjaCB0aW1lIHJlY29yZGluZyB3aGV0aGVyIHRoZSBzcGVjaWVzIHdhcyBzZWVuIG9yIG5vdC4KCiMjIyBTaW1wbGUgb2NjdXBhbmN5IG1vZGVsIHdpdGggY29uc3RhbnQgJFxQc2kkIGFuZCAkcCQKCkhlcmUgaXMgYSBzaW11bGF0aW9uIHRvIGRlbW9uc3RyYXRlIHRoZXNlIGlkZWFzIGFuZCBob3cgYSBzaW5nbGUtc2Vhc29uIG9jY3VwYW5jeSBtb2RlbCBjYW4gYmUgdXNlZCB0byBlc3RpbWF0ZSB0cnVlIG9jY3VwYW5jeS4gSSB1c2UgYW4gYXJ0aWZpY2FsIGxhbmRzY2FwZSBvZiAkMjAgXHRpbWVzIDIwJCBncmlkIGNlbGxzLCBlYWNoIG9mIHRoZW0gb2NjdXBpZWQgd2l0aCBwcm9iYWJpbGl0eSAkXFBzaSA9IDAuNSQuIFdpdGggZWFjaCB2aXNpdCwgd2UgaGF2ZSBhIGRldGVjdGlvbiBwcm9iYWJpbGl0eSAkcCA9IDAuMyQgb2YgZmluZGluZyB0aGUgc3BlY2llcyBpZiB0aGUgZ3JpZCBjZWxsIGlzIG9jY3VwaWVkLgoKYGBge3J9CgpzZXQuc2VlZCgyNDYpCgpvY2MucHJvYiA8LSAwLjUKZGV0LnByb2IgPC0gMC4zCgpncmlkIDwtIGV4cGFuZC5ncmlkKGxvbmc9YygxOjIwKSxsYXQ9YygxOjIwKSkKZ3JpZCRvY2MgPC0gcmJpbm9tKGRpbShncmlkKVsxXSwxLHByb2IgPSBvY2MucHJvYikgIyBvY2N1cGFuY3kKZ3JpZCRjb2wgPC0gaWZlbHNlKGdyaWQkb2NjPT0xLCdvcmFuZ2UnLCdncmV5JykKCmdyaWQkZGV0ZWN0ZWQxIDwtIHJiaW5vbShkaW0oZ3JpZClbMV0sMSxwcm9iID0gZ3JpZCRvY2MgKiBkZXQucHJvYikgIyBkZXRlY3Rpb25zIG9uIGZpcnN0IHN1cnZleQpncmlkJGRldGVjdGVkMiA8LSByYmlub20oZGltKGdyaWQpWzFdLDEscHJvYiA9IGdyaWQkb2NjICogZGV0LnByb2IpICMgZGV0ZWN0aW9ucyBvbiBzZWNvbmQgc3VydmV5CmdyaWQkZGV0ZWN0ZWQzIDwtIHJiaW5vbShkaW0oZ3JpZClbMV0sMSxwcm9iID0gZ3JpZCRvY2MgKiBkZXQucHJvYikgIyBkZXRlY3Rpb25zIG9uIHRoaXJkIHN1cnZleQoKYGBgCgpMZXQncyBwbG90IHRoZSBzaW11bGF0ZWQgbGFuZHNjYXBlLiBUcnVseSBvY2N1cGllZCBncmlkIGNlbGxzIGFyZSBvcmFuZ2U7IGRldGVjdGlvbnMgYmxhY2suIEEgZ3JleSBjZWxsIG1lYW5zIHRoYXQgdGhlIHNwZWNpZXMgd2FzIG5vdCBkZXRlY3RlZC4KCmBgYHtyfQpvcCA8LSBwYXIobWZyb3c9YygyLDIpLCBtYXI9YygxLDEsMSwxKSkKcGxvdChsYXR+bG9uZyxjb2w9Y29sLHBjaD0gMTUsIGRhdGE9Z3JpZCxhbm49RixheGVzPUYsY2V4PTEuNSkKCnBsb3QobGF0fmxvbmcsY29sPSdncmV5JyxwY2g9IDE1LCBkYXRhPWdyaWQsYW5uPUYsYXhlcz1GLGNleD0xLjUpCnBvaW50cyhsYXR+bG9uZyxjb2w9ZGV0ZWN0ZWQxLHBjaD0gMTUsIGRhdGE9Z3JpZFtncmlkJGRldGVjdGVkMT09MSxdLGNleD0xLjUpCgpwbG90KGxhdH5sb25nLGNvbD0nZ3JleScscGNoPSAxNSwgZGF0YT1ncmlkLGFubj1GLGF4ZXM9RixjZXg9MS41KQpwb2ludHMobGF0fmxvbmcsY29sPWRldGVjdGVkMixwY2g9IDE1LCBkYXRhPWdyaWRbZ3JpZCRkZXRlY3RlZDI9PTEsXSxjZXg9MS41KQoKcGxvdChsYXR+bG9uZyxjb2w9J2dyZXknLHBjaD0gMTUsIGRhdGE9Z3JpZCxhbm49RixheGVzPUYsY2V4PTEuNSkKcG9pbnRzKGxhdH5sb25nLGNvbD1kZXRlY3RlZDMscGNoPSAxNSwgZGF0YT1ncmlkW2dyaWQkZGV0ZWN0ZWQzPT0xLF0sY2V4PTEuNSkKcGFyKG9wKQpgYGAKCkZvciB0aGUgYW5hbHlzaXMsIHdlIG5lZWQgdGhlIGRhdGEgaW4gdGhlIGZvcm1hdCBvZiBzdXJ2ZXkgaGlzdG9yaWVzOgoKYGBge3J9CnkgPC0gZ3JpZFssNTo3XQoKaGVhZCh5KQpkaW0oeSkKYGBgCgpJbiBob3cgbWFueSBncmlkIGNlbGxzIGRpZCB3ZSBkZXRlY3QgdGhlIHNwZWNpZXMgYXQgbGVhc3Qgb25jZT8KCmBgYHtyfQpzdW0oYXBwbHkoeSwxLHN1bSk+MCkvNDAwIApgYGAKClRoaXMgYG5haXZlJyBlc3RpbWF0ZSBvZiBvY2N1cGFuY3kgaXMgbG93ZXIgdGhhbiB0aGUgdHJ1ZSBvY2N1cGFuY3kgcHJvYmFiaWxpdHksIHdoaWNoIC0tIGFzIHdlIGtub3cgaW4gdGhpcyBjYXNlIGJlY2F1c2Ugd2Ugc2ltdWxhdGVkIHRoZSBkYXRhIC0tIGlzIDAuNS4gTGV0J3Mgc2VlIHdoZXRoZXIgYW4gb2NjdXBhbmN5IG1vZGVsIHJldHJpZXZlcyBhIGJldHRlciBlc3RpbWF0ZS4KCmBgYHtyfQpyZXF1aXJlKHVubWFya2VkKQoKeSA8LSBncmlkWyw1OjddCmV4MSA8LSB1bm1hcmtlZEZyYW1lT2NjdSh5ID0geSkKc3VtbWFyeShleDEpCgptMSA8LSBvY2N1KH4xIH4xLCBkYXRhPWV4MSkKc3VtbWFyeShtMSkKCmJhY2tUcmFuc2Zvcm0obTEsIHR5cGU9InN0YXRlIikKYmFja1RyYW5zZm9ybShtMSwgdHlwZT0iZGV0IikKYGBgCgpZZXMsIGl0IGRvZXMhIE91ciBlc3RpbWF0ZWQgb2NjdXBhbmN5ICgwLjQ1MikgYW5kIGRldGVjdGlvbiBwcm9iYWJpbGl0aWVzICgwLjMyOCkgYXJlIGNsb3NlIChnaXZlbiB0aGUgc3RhbmRhcmQgZXJyb3JzKSB0byB0aGUgdmFsdWVzIHVzZWQgaW4gdGhlIHNpbXVsYXRpb24uIEJ1dCBkb2VzIG91ciBtb2RlbCBhZGVxdWF0ZWx5IHJlZmxlY3QgdGhlIHN0cnVjdHVyZSBpbiBvdXIgZGF0YT8gV2UgY2FuIGNoZWNrIHRoYXQgd2l0aCBhIEdvb2RuZXNzLW9mLWZpdCB0ZXN0IHVzaW5nIHRoZSBwYXJhbWV0cmljIGJvb3RzdHJhcCByb3V0aW5nIGluIGB1bm1hcmtlZCc6CgpgYGB7cn0KIyBnb29kbmVzcyBvZiBmaXQKCmNoaXNxIDwtIGZ1bmN0aW9uKGZtKSB7ICAjIGNhbGN1bGF0ZXMgY2hpLXNxdWFyZSB2YWx1ZQogIG9ic2VydmVkIDwtIGdldFkoZm1AZGF0YSkKICBleHBlY3RlZCA8LSBmaXR0ZWQoZm0pCiAgc3VtKChvYnNlcnZlZCAtIGV4cGVjdGVkKV4yL2V4cGVjdGVkKQogIH0KCnBiIDwtIHBhcmJvb3QobTEsIHN0YXRpc3RpYz1jaGlzcSwgbnNpbT0xMDApCgpwbG90KHBiLCBtYWluID0gIiIpCgpgYGAKCk5vdywgd2UgY2FuIGZpdCB0aGUgc2FtZSBtb2RlbCBpbiBhIEJheWVzaWFuIGFuYWx5c2lzIHVzaW5nIEpBR1M6CgpgYGB7cn0KIyBjb2RlIHRoZSBtb2RlbCBpbiBCVUdTCnNpbmsoIm9jYy50eHQiKQpjYXQoIgogICAgbW9kZWwgewogICAgIyBFY29sb2dpY2FsIG1vZGVsIGZvciB0cnVlIG9jY3VycmVuY2UKICAgIGZvciAoaSBpbiAxOm5zaXRlcykgewogICAgeltpXSB+IGRiZXJuKHBzaSkKCiAgICAjIE9ic2VydmF0aW9uIG1vZGVsIGZvciByZXBsaWNhdGVkIGRldGVjdGlvbi9ub25kZXRlY3Rpb24gb2JzZXJ2YXRpb25zCiAgICBwLmVmZltpXSA8LSB6W2ldICogcAogICAgZm9yIChqIGluIDE6bnZpc2l0cykgewogICAgeVtpLGpdIH4gZGJlcm4ocC5lZmZbaV0pCiAgICB9ICNqCiAgICB9ICNpCgogICAgIyBQcmlvcnMKICAgIHBzaSB+IGR1bmlmKDAsIDEpCiAgICBwIH4gZHVuaWYoMCwgMSkKCiAgICAjIERlcml2ZWQgcXVhbnRpdGllcwogICAgb2NjLmZzIDwtIHN1bSh6W10pCiAgICB9CiAgICAiLGZpbGwgPSBUUlVFKQpzaW5rKCkKCmxpYnJhcnkoamFnc1VJKQoKb2NjLmRhdGEgPC0gbGlzdCh5ID0geSwgbnNpdGVzID0gbnJvdyh5KSwgbnZpc2l0cyA9IG5jb2woeSkpCiMgSW5pdGlhbCB2YWx1ZXMKenN0IDwtIGFwcGx5KHksIDEsIG1heCkKIyBPYnNlcnZlZCBvY2N1cnJlbmNlIGFzIHN0YXJ0aW5nIHZhbHVlcyBmb3Igegppbml0cyA8LSBmdW5jdGlvbigpIGxpc3QoeiA9IHpzdCkKIyBQYXJhbWV0ZXJzIG1vbml0b3JlZApwYXJhbXMgPC0gYygicHNpIiwgInAiLCAib2NjLmZzIikKCiMgTUNNQyBzZXR0aW5ncwpuYyA8LSAzICAgOyAgIG5pIDwtIDUwMDAgICA7ICAgbmIgPC0gMjAwMCAgIDsgICBudCA8LSAyCgpvdXQub2NjIDwtIGphZ3Mob2NjLmRhdGEsIGluaXRzLCBwYXJhbXMsICJvY2MudHh0Iiwgbi5jaGFpbnM9bmMsIG4uaXRlcj1uaSwgbi5idXJuID0gbmIsIG4udGhpbj1udCkKcGxvdChvdXQub2NjKQpwcmludChvdXQub2NjJHN1bW1hcnksIGRpZ2l0cz0yKQpgYGAKClRoZSBCYXllc2lhbiBhbmFseXNpcyByZXR1cm5zIGVzdGltYXRlcyB0aGF0IGFyZSB2ZXJ5IHNpbWlsYXIgdG8gdGhlIGZyZXF1ZW50aXN0IGVzdGltYXRlcy4gVGhpcyBpcyBhcyBpdCBzaG91bGQgYmUgc2luY2Ugd2UgdXNlZCBub24taW5mb3JtYXRpdmUgcHJpb3JzLgoKIyMjIE9jY3VwYW5jeSBtb2RlbCB3aXRoIHNpdGUtbGV2ZWwgY292YXJpYXRlcwoKT2ssIHRoaXMgd2FzIGEgdmVyeSBzaW1wbGUgb2NjdXBhbmN5IG1vZGVsLCB3aXRoIGNvbnN0YW50IG9jY3VwYW5jeSBhbmQgZGV0ZWN0aW9uIHByb2JhYmlsaXRpZXMuIFVzdWFsbHksIHdlIGFyZSBpbnRlcmVzdGVkIGluIGhvdyBvY2N1cGFuY3kgdmFyaWVzIGFtb25nIHNpdGVzLiBJIG5leHQgc2ltdWxhdGUgYSBkYXRhIHNldCB3aGVyZSBvY2N1cGFuY3kgcHJvYmFiaWxpdGllcyBpbmNyZWFzZSBmcm9tIGxlZnQgdG8gcmlnaHQgLS0gb3IgZnJvbSB3ZXN0IHRvIGVhc3QsIGlmIHlvdSB3aXNoLiBBdCB0aGUgc2FtZSB0aW1lLCBkZXRlY3Rpb24gcHJvYmFiaWxpdGllcyBkZWNyZWFzZSBmcm9tIHdlc3QgdG8gZWFzdC4gU28gd2hlbiB3ZSBsb29rIGF0IHRoZSByYXcgZGV0ZWN0aW9uIGRhdGEsIGl0IGlzIG5vdCBvYnZpb3VzIHRoYXQgdGhlcmUgYXJlIGFueSB0cmVuZHMgaW4gdGhlc2UgZGF0YSBhdCBhbGwuIFRoZSB0cmVuZHMgaW4gdGhlIGVjb2xvZ2ljYWwgYW5kIG9ic2VydmF0aW9uIHByb2Nlc3MgY2FuY2VsIGVhY2ggb3RoZXIgb3V0LiBMZXQncyBzZWUgd2hldGhlciB3ZSBjYW4gZGlzZW50YW5nbGUgdGhlIHR3byBlZmZlY3RzIHVzaW5nIGFuIG9jY3VwYW5jeSBtb2RlbC4KCmBgYHtyfQpsaWJyYXJ5KGJvb3QpCgpzZXQuc2VlZCgyNDYpCgpncmlkIDwtIGV4cGFuZC5ncmlkKGxvbmc9YygxOjIwKSxsYXQ9YygxOjIwKSkKCm9jYy5wcm9iIDwtIGludi5sb2dpdCgtMiswLjIqZ3JpZCRsb25nKSAjIGdyYWRpZW50IGluIG9jY3VwYW5jeSBwcm9iYWJpbGl0eQpkZXQucHJvYiA8LSBpbnYubG9naXQoMi0wLjIqZ3JpZCRsb25nKSAjIG9wcG9zaXRlIGdyYWRpZW50IGluIGRldGVjdGlvbiBwcm9iYWJpbGl0eQoKZ3JpZCRvY2MgPC0gcmJpbm9tKGRpbShncmlkKVsxXSwxLHByb2IgPSBvY2MucHJvYikgIyBvY2N1cGFuY3kKZ3JpZCRjb2wgPC0gaWZlbHNlKGdyaWQkb2NjPT0xLCdvcmFuZ2UnLCdncmV5JykKCmdyaWQkZGV0ZWN0ZWQxIDwtIHJiaW5vbShkaW0oZ3JpZClbMV0sMSxwcm9iID0gZ3JpZCRvY2MgKiBkZXQucHJvYikgIyBkZXRlY3Rpb25zIG9uIGZpcnN0IHN1cnZleQpncmlkJGRldGVjdGVkMiA8LSByYmlub20oZGltKGdyaWQpWzFdLDEscHJvYiA9IGdyaWQkb2NjICogZGV0LnByb2IpICMgZGV0ZWN0aW9ucyBvbiBzZWNvbmQgc3VydmV5CmdyaWQkZGV0ZWN0ZWQzIDwtIHJiaW5vbShkaW0oZ3JpZClbMV0sMSxwcm9iID0gZ3JpZCRvY2MgKiBkZXQucHJvYikgIyBkZXRlY3Rpb25zIG9uIHRoaXJkIHN1cnZleQoKIyBwbG90IHRoZSB0cnVlIG9jY3VwYW5jeSBhbmQgZGV0ZWN0aW9ucyBkdXJpbmcgZWFjaCBvZiB0aGUgdGhyZWUgc3VydmV5cwpvcCA8LSBwYXIobWZyb3c9YygyLDIpLCBtYXI9YygxLDEsMSwxKSkKcGxvdChsYXR+bG9uZyxjb2w9Y29sLHBjaD0gMTUsIGRhdGE9Z3JpZCxhbm49RixheGVzPUYsY2V4PTEuNSkKCnBsb3QobGF0fmxvbmcsY29sPSdncmV5JyxwY2g9IDE1LCBkYXRhPWdyaWQsYW5uPUYsYXhlcz1GLGNleD0xLjUpCnBvaW50cyhsYXR+bG9uZyxjb2w9ZGV0ZWN0ZWQxLHBjaD0gMTUsIGRhdGE9Z3JpZFtncmlkJGRldGVjdGVkMT09MSxdLGNleD0xLjUpCgpwbG90KGxhdH5sb25nLGNvbD0nZ3JleScscGNoPSAxNSwgZGF0YT1ncmlkLGFubj1GLGF4ZXM9RixjZXg9MS41KQpwb2ludHMobGF0fmxvbmcsY29sPWRldGVjdGVkMixwY2g9IDE1LCBkYXRhPWdyaWRbZ3JpZCRkZXRlY3RlZDI9PTEsXSxjZXg9MS41KQoKcGxvdChsYXR+bG9uZyxjb2w9J2dyZXknLHBjaD0gMTUsIGRhdGE9Z3JpZCxhbm49RixheGVzPUYsY2V4PTEuNSkKcG9pbnRzKGxhdH5sb25nLGNvbD1kZXRlY3RlZDMscGNoPSAxNSwgZGF0YT1ncmlkW2dyaWQkZGV0ZWN0ZWQzPT0xLF0sY2V4PTEuNSkKcGFyKG9wKQoKIyB0aGUgdHJ1ZSBncmFkaWVudHMgbG9vayBsaWtlIHRoaXMKb3AgPC0gcGFyKG1mcm93PWMoMiwxKSwgbWFyPWMoMSw0LDEsMSksIG9tYT1jKDMsMCwwLDApKQpwbG90KDE6MjAsIGludi5sb2dpdCgtMiswLjIqMToyMCksIHlsaW09YygwLDEpLCB0eXBlPSdsJywgYXhlcz1GLCB5bGFiPSd0cnVlIG9jY3VwYW5jeScsIGx3ZD0yLCB4bGltPWMoMCwyMCkpICAjIHRydWUgb2NjdXBhbmN5CmF4aXMoMiwgYXQ9YygwLDAuNSwxKSwgbGFzPTEpCnBsb3QoMToyMCwgaW52LmxvZ2l0KDItMC4yKjE6MjApLCB5bGltPWMoMCwxKSwgdHlwZT0nbCcsIGF4ZXM9RiwgeWxhYj0ndHJ1ZSBkZXRlY3Rpb24nLCBsd2Q9Mix4bGltPWMoMCwyMCkpICAjIHRydWUgZGV0ZWN0aW9uIHByb2JhYmlsaXR5CmF4aXMoMiwgYXQ9YygwLDAuNSwxKSwgbGFzPTEpCmF4aXMoMSwgYXQ9YygwLDEwLDIwKSkKbXRleHQoIkxvbmdpdHVkZSIsIDEsIDEuNSwgb3V0ZXI9VCkKcGFyKG9wKQpgYGAKClRoZSBmaWd1cmVzIGFib3ZlIHNob3cgdGhlIGdyYWRpZW50cyBpbiBvY2N1cGFuY3kgYW5kIGRldGVjdGlvbi4gVGhleSBydW4gaW4gb3Bwb3NpdGUgZGlyZWN0aW9ucyBhbmQgdGhlcmVmb3JlIHdlIGNhbid0IHNlZSBhIGNsZWFyIHBhdHRlcm4gaW4gdGhlIHJhdyBkYXRhLgoKV2UgYW5hbHlzZSB0aGlzIGRhdGEgc2V0IHVzaW5nIGB1bm1hcmtlZCc6CgpgYGB7cn0KcmVxdWlyZSh1bm1hcmtlZCkKCnkgPC0gZ3JpZFssNTo3XQpzaXRlQ292cyA8LSBkYXRhLmZyYW1lKGxvbmc9Z3JpZFssImxvbmciXSkKZXgxIDwtIHVubWFya2VkRnJhbWVPY2N1KHkgPSB5LCBzaXRlQ292cyA9IHNpdGVDb3ZzKQpzdW1tYXJ5KGV4MSkKCm0xIDwtIG9jY3UofmxvbmcgfmxvbmcsIGRhdGE9ZXgxKQpzdW1tYXJ5KG0xKQpgYGAKClRoaXMgbW9kZWwgY29ycmVjdGx5IHJldHVybnMgdGhlIGdyYWRpZW50czogYSBwb3NpdGl2ZSBlZmZlY3Qgb2YgbG9uZ2l0dWRlIG9uIG9jY3VwYW5jeSBhbmQgYSBuZWdhdGl2ZSBvbmUgb24gZGV0ZWN0aW9uLiBJbiBmYWN0LCB0aGUgZXN0aW1hdGVkIHZhbHVlcyBhcmUgY2xvc2UgdG8gdGhlIG9uZXMgd2UgdXNlZCBpbiB0aGUgc2ltdWxhdGlvbnMgKCRsb2dpdChcUHNpKSA9IC0yICsgMC4yIFx0aW1lcyBsb25naXR1ZGUkIGFuZCAkbG9naXQocCkgPSAyIC0gMC4yIFx0aW1lcyBsb25naXR1ZGUkKS4KCkJ1dCB0aGlzIG1vZGVsIHdhcyBvZiBjb3Vyc2UgdGhlIHRydWUgZ2VuZXJhdGluZyBtb2RlbCwgd2hpY2ggd2Uga25ldyBpbiB0aGlzIGNhc2UuIExldCdzIGZpdCBzb21lIHNpbXBsZXIgbW9kZWxzIGFuZCBzZWUgd2hldGhlciBtb2RlbCBzZWxlY3Rpb24gcGlja3MgdGhlIGNvcnJlY3QgbW9kZWw6CgpgYGB7cn0KbTIgPC0gb2NjdSh+MSB+bG9uZywgZGF0YT1leDEpCnN1bW1hcnkobTIpCgptMyA8LSBvY2N1KH4xIH4xLCBkYXRhPWV4MSkKc3VtbWFyeShtMykKCm1sIDwtIGZpdExpc3QoInBzaShsb25nKXAobG9uZykiID0gbTEsICJwc2kobG9uZylwKC4pIj1tMiwgInBzaSguKXAoLikiID0gbTMpCm1vZFNlbChtbCkKYGBgCkl0IGRvZXMuIFdlIGhhdmUgY2xlYXIgZXZpZGVuY2UgZm9yIG1vZGVsIG0xIChpbmNsdWRpbmcgdGhlIGVmZmVjdCBvZiBsb25naXR1ZGUgb24gYm90aCAkXFBzaSQgYW5kICRwJCkgdG8gYmUgdGhlIGJlc3QgbW9kZWwgaW4gdGhlIHNldC4KCldoYXQgZG8gdGhlIGZpdHRlZCByZWxhdGlvbnNoaXBzIGxvb2sgbGlrZT8gV2UgdXNlIHRoZSBtb2RlbCB0byBwcmVkaWN0IG9jY3VwYW5jeSBhbmQgZGV0ZWN0aW9uIGFzIGEgZnVuY3Rpb24gb2YgbG9uZ2l0dWRlLgoKYGBge3J9CmJhY2tUcmFuc2Zvcm0obGluZWFyQ29tYihtMSwgY29lZmZpY2llbnRzID0gYygxLDApLCB0eXBlID0gInN0YXRlIikpICAjIG9jY3VwYW5jeSBhdCBsb25nID0gMApiYWNrVHJhbnNmb3JtKGxpbmVhckNvbWIobTEsIGNvZWZmaWNpZW50cyA9IGMoMSwyMCksIHR5cGUgPSAic3RhdGUiKSkgICMgb2NjdXBhbmN5IGF0IGxvbmcgPSAyMAoKbmV3RGF0YSA8LSBkYXRhLmZyYW1lKGxvbmc9MToyMCkKcHJlZGljdChtMSwgdHlwZSA9ICJzdGF0ZSIsIG5ld2RhdGEgPSBuZXdEYXRhLCBhcHBlbmREYXRhPVRSVUUpCnByLnMgPC0gcHJlZGljdChtMSwgdHlwZSA9ICJzdGF0ZSIsIG5ld2RhdGEgPSBuZXdEYXRhLCBhcHBlbmREYXRhPVRSVUUpCnByLmQgPC0gcHJlZGljdChtMSwgdHlwZSA9ICJkZXQiLCBuZXdkYXRhID0gbmV3RGF0YSwgYXBwZW5kRGF0YT1UUlVFKQoKb3AgPC0gcGFyKG1mcm93PWMoMiwxKSwgbWFyPWMoMSw0LDEsMSksIG9tYT1jKDMsMCwwLDApKQpwbG90KDE6MjAsIGludi5sb2dpdCgtMiswLjIqMToyMCksIHlsaW09YygwLDEpLCB0eXBlPSdsJywgYXhlcz1GLCB5bGFiPSd0cnVlIG9jY3VwYW5jeScsIGx3ZD0yLCB4bGltPWMoMCwyMCkpICAjIHRydWUgb2NjdXBhbmN5CnBvaW50cyhQcmVkaWN0ZWQgfiBsb25nLCBkYXRhPXByLnMpCmF4aXMoMiwgYXQ9YygwLDAuNSwxKSwgbGFzPTEpCnBsb3QoMToyMCwgaW52LmxvZ2l0KDItMC4yKjE6MjApLCB5bGltPWMoMCwxKSwgdHlwZT0nbCcsIGF4ZXM9RiwgeWxhYj0ndHJ1ZSBkZXRlY3Rpb24nLCBsd2Q9MiwgeGxpbT1jKDAsMjApKSAgIyB0cnVlIGRldGVjdGlvbiBwcm9iYWJpbGl0eQpwb2ludHMoUHJlZGljdGVkIH4gbG9uZywgZGF0YT1wci5kKQpheGlzKDIsIGF0PWMoMCwwLjUsMSksIGxhcz0xKQpheGlzKDEsIGF0PWMoMCwxMCwyMCkpCm10ZXh0KCJMb25naXR1ZGUiLCAxLCAxLjUsIG91dGVyPVQpCnBhcihvcCkKYGBgCgpXZSBzZWUgdGhhdCB0aGUgZml0dGVkIHJlbGF0aW9uc2hpcHMgKHNtYWxsIGNpcmNsZXMpIGNsb3NlbHkgZm9sbG93IHRoZSB0cnVlIHJlbGF0aW9uc2hpcHMgKHNvbGlkIGxpbmVzKS4gVGhlIHByZWRpY3RlZCBvY2N1cGFuY3kgcHJvYmFiaWxpdGllcyBmb3Igb3VyIHNpbXVsYXRlZCBsYW5kc2NhcGUgY29ycmVjdGx5IGRlc2NyaWJlIHRoZSBncmFkaWVudCBmcm9tIHdlc3QgdG8gZWFzdCBhbmQgbG9vayBsaWtlIHRoaXM6CgpgYGB7cn0KbmV3RGF0YSA8LSBkYXRhLmZyYW1lKGxvbmcgPSBncmlkJGxvbmcpCmdyaWQkUHJlZGljdGVkbTEgPC0gcHJlZGljdChtMSwgdHlwZSA9ICJzdGF0ZSIsIG5ld2RhdGEgPSBuZXdEYXRhKSRQcmVkaWN0ZWQKCm9wIDwtIHBhcihtZmNvbD1jKDEsMiksIG1hcj1jKDEsMSwxLDEpKQpwbG90KGxhdH5sb25nLGNvbD1jb2wscGNoPSAxNSwgZGF0YT1ncmlkLGFubj1GLGF4ZXM9RixjZXg9MS41KQpuYyA8LSBhcy5udW1lcmljKGN1dChncmlkJFByZWRpY3RlZG0xLDIwKSkKcGxvdChsYXR+bG9uZyxwY2g9IDE1LCBkYXRhPWdyaWQsYW5uPUYsYXhlcz1GLGNleD0xLjUsY29sPWNvbG9yUmFtcFBhbGV0dGUoYygid2hpdGUiLCAiZGFya29yYW5nZTMiKSkoMjApW25jXSkKCmBgYAoKQW5kIG5vdyBhbmFseXNlIHRoZSBzYW1lIG1vZGVsIGluIEpBR1M6CgpgYGB7cn0Kc2luaygib2NjX2Nvdi50eHQiKQpjYXQoIgogICAgbW9kZWwgewogICAgIyBFY29sb2dpY2FsIG1vZGVsIGZvciB0cnVlIG9jY3VycmVuY2UKICAgIGZvciAoaSBpbiAxOm5zaXRlcykgewogICAgeltpXSB+IGRiZXJuKHBzaVtpXSkKICAgIAogICAgIyBPYnNlcnZhdGlvbiBtb2RlbCBmb3IgcmVwbGljYXRlZCBkZXRlY3Rpb24vbm9uZGV0ZWN0aW9uIG9ic2VydmF0aW9ucwogICAgcC5lZmZbaV0gPC0geltpXSAqIHBbaV0KICAgIGZvciAoaiBpbiAxOm52aXNpdHMpIHsKICAgIHlbaSxqXSB+IGRiZXJuKHAuZWZmW2ldKQogICAgfSAjagogICAgfSAjaQogICAgCiAgICAjIGNvdmFyaWF0ZXMKICAgIGZvciAoaSBpbiAxOm5zaXRlcyl7CiAgICAgIGxvZ2l0KHBzaVtpXSkgPSBiZXRhLnBzaVsxXSArIGJldGEucHNpWzJdICogbG9uZ1tpXQogICAgICBsb2dpdChwW2ldKSA9IGJldGEucFsxXSArIGJldGEucFsyXSAqIGxvbmdbaV0KICAgIH0KCiAgICAjIFByaW9ycwogICAgZm9yIChiIGluIDE6Mil7CiAgICAgIGJldGEucHNpW2JdIH4gZG5vcm0oMCwgMC4wMSkKICAgICAgYmV0YS5wW2JdIH4gZG5vcm0oMCwgMC4wMSkKICAgIH0KICAgIAogICAgIyBEZXJpdmVkIHF1YW50aXRpZXMKICAgIG9jYy5mcyA8LSBzdW0oeltdKQogICAgfQogICAgIixmaWxsID0gVFJVRSkKc2luaygpCgpsaWJyYXJ5KGphZ3NVSSkKCm9jYy5kYXRhIDwtIGxpc3QoeSA9IGdyaWRbLDU6N10sIG5zaXRlcyA9IG5yb3coeSksIG52aXNpdHMgPSBuY29sKHkpLCBsb25nPWdyaWRbLCJsb25nIl0pCiMgSW5pdGlhbCB2YWx1ZXMKenN0IDwtIGFwcGx5KHksIDEsIG1heCkKIyBPYnNlcnZlZCBvY2N1cnJlbmNlIGFzIHN0YXJ0aW5nIHZhbHVlcyBmb3Igegppbml0cyA8LSBmdW5jdGlvbigpIGxpc3QoeiA9IHpzdCwgYmV0YS5wc2k9cnVuaWYoMiwtMywzKSwgYmV0YS5wPXJ1bmlmKDIsLTMsMykpCiMgUGFyYW1ldGVycyBtb25pdG9yZWQKcGFyYW1zIDwtIGMoImJldGEucHNpIiwgImJldGEucCIsICJvY2MuZnMiKQoKIyBNQ01DIHNldHRpbmdzCm5jIDwtIDMgICA7ICAgbmkgPC0gNTAwMCAgIDsgICBuYiA8LSAyMDAwICAgOyAgIG50IDwtIDIKCm91dC5vY2MuY292IDwtIGphZ3Mob2NjLmRhdGEsIGluaXRzLCBwYXJhbXMsICJvY2NfY292LnR4dCIsIG4uY2hhaW5zPW5jLCBuLml0ZXI9bmksIG4uYnVybiA9IG5iLCBuLnRoaW49bnQpCnBsb3Qob3V0Lm9jYy5jb3YpCnByaW50KG91dC5vY2MuY292JHN1bW1hcnksIGRpZ2l0cz0yKQpgYGAKCkFuZCBwbG90dGluZyB0aGUgZml0dGVkIHJlbGF0aW9uc2hpcHM6CgpgYGB7cn0KbmV3LmxvbmcgPC0gMToyMAoKcHIucyA8LSBpbnYubG9naXQob3V0Lm9jYy5jb3Ykc3VtbWFyeVsxLDFdICsgb3V0Lm9jYy5jb3Ykc3VtbWFyeVsyLDFdICogbmV3LmxvbmcpCnByLmQgPC0gaW52LmxvZ2l0KG91dC5vY2MuY292JHN1bW1hcnlbMywxXSArIG91dC5vY2MuY292JHN1bW1hcnlbNCwxXSAqIG5ldy5sb25nKQoKb3AgPC0gcGFyKG1mcm93PWMoMiwxKSwgbWFyPWMoMSw0LDEsMSksIG9tYT1jKDMsMCwwLDApKQpwbG90KDE6MjAsIGludi5sb2dpdCgtMiswLjIqMToyMCksIHlsaW09YygwLDEpLCB0eXBlPSdsJywgYXhlcz1GLCB5bGFiPSd0cnVlIG9jY3VwYW5jeScsIGx3ZD0yLCB4bGltPWMoMCwyMCkpICAjIHRydWUgb2NjdXBhbmN5CnBvaW50cyhuZXcubG9uZywgcHIucykKYXhpcygyLCBhdD1jKDAsMC41LDEpLCBsYXM9MSkKcGxvdCgxOjIwLCBpbnYubG9naXQoMi0wLjIqMToyMCksIHlsaW09YygwLDEpLCB0eXBlPSdsJywgYXhlcz1GLCB5bGFiPSd0cnVlIGRldGVjdGlvbicsIGx3ZD0yLCB4bGltPWMoMCwyMCkpICAjIHRydWUgZGV0ZWN0aW9uIHByb2JhYmlsaXR5CnBvaW50cyhuZXcubG9uZywgcHIuZCkKYXhpcygyLCBhdD1jKDAsMC41LDEpLCBsYXM9MSkKYXhpcygxLCBhdD1jKDAsMTAsMjApKQptdGV4dCgiTG9uZ2l0dWRlIiwgMSwgMS41LCBvdXRlcj1UKQpwYXIob3ApCgpuZXdEYXRhIDwtIGRhdGEuZnJhbWUobG9uZyA9IGdyaWQkbG9uZykKZ3JpZCRQcmVkaWN0ZWRtMSA8LSBpbnYubG9naXQob3V0Lm9jYy5jb3Ykc3VtbWFyeVsxLDFdICsgb3V0Lm9jYy5jb3Ykc3VtbWFyeVsyLDFdICogbmV3LmxvbmcpCgpvcCA8LSBwYXIobWZjb2w9YygxLDIpLCBtYXI9YygxLDEsMSwxKSkKcGxvdChsYXR+bG9uZyxjb2w9Y29sLHBjaD0gMTUsIGRhdGE9Z3JpZCxhbm49RixheGVzPUYsY2V4PTEuNSkKbmMgPC0gYXMubnVtZXJpYyhjdXQoZ3JpZCRQcmVkaWN0ZWRtMSwyMCkpCnBsb3QobGF0fmxvbmcscGNoPSAxNSwgZGF0YT1ncmlkLGFubj1GLGF4ZXM9RixjZXg9MS41LGNvbD1jb2xvclJhbXBQYWxldHRlKGMoIndoaXRlIiwgImRhcmtvcmFuZ2UzIikpKDIwKVtuY10pCmBgYAoKVGhhdCdzIGl0LgoKIyMjIEFuIGV4YW1wbGUgd2l0aCBhIHJlYWwgZGF0YSBzZXQgdGhlIFNBQkFQMiBkYXRhIGZvciBiYWxkIGliaXMKCkFib3ZlLCB3ZSBsb29rZWQgYXQgc2ltdWxhdGVkIGRhdGEsIHdoaWNoIGhhcyB0aGUgYWR2YW50YWdlIHRoYXQgd2Uga25vdyB0aGUgdHJ1dGggYW5kIHNvIHdlIGNhbiBlYXNpbHkgY2hlY2sgdGhhdCBvdXIgbW9kZWwgZG9lcyB0aGUgY29ycmVjdCB0aGluZy4gVXN1YWxseSwgb2YgY291cnNlLCB3ZSBkb24ndCBrbm93IHRydXRoIGJ1dCBhcmUgaW50ZXJlc3RlZCBpbiBleHRyYWN0aW5nIG9jY3VwYW5jeSBwYXR0ZXJucyBmcm9tIHJlYWwgZGF0YS4gSGVyZSBpcyBhbiBleGFtcGxlIHVzaW5nIGRhdGEgZnJvbSB0aGUgW1NlY29uZCBTb3V0aGVybiBBZnJpY2FuIEJpcmQgQXRsYXMgUHJvamVjdCBdKGh0dHA6Ly9zYWJhcDIuYWR1Lm9yZy56YS8pLiBJIGNob3NlIHRoZSBiYWxkIGliaXMgYXMgYW4gZXhhbXBsZS4KCmBgYHtyfQpiaS5tIDwtIHJlYWQuY3N2KCJiYWxkaWJpc2RhdGEuY3N2IikKCiMgc3RhcnQgb2NjdXBhbmN5IGFuYWx5c2lzCnJlcXVpcmUodW5tYXJrZWQpCgp3cml0ZS5jc3YoYmkubVssYygiUGVudGFkIiwiamRheSIsIlNwcCIsImxob3VycyIpXSwiYmkuY3N2Iixyb3cubmFtZXMgPSBGKQpiaS51bWYgPC0gY3N2VG9VTUYoImJpLmNzdiIsbG9uZz1UUlVFLCB0eXBlID0gInVubWFya2VkRnJhbWVPY2N1IikKc3VtbWFyeShiaS51bWYpCgojIGhhdmUgYSBsb29rIGF0IHRoZSBvYnNlcnZhdGlvbi1sZXZlbCBjb3ZhcmlhdGUKaGVhZChvYnNDb3ZzKGJpLnVtZixtYXRyaWNlcz1UKSRsaG91cnMpCgojIGFkZCBzaXRlLWxldmVsIGNvdmFyaWF0ZXMgKGZvciBzb21lIHJlYXNvbiwgdGhpcyBkaWRuJ3Qgd29yayBhdXRvbWF0aWNhbGx5IHdpdGggJ2NzdlRvVU1GJykKc2l0ZUNvdnMoYmkudW1mKSA8LSBkYXRhLmZyYW1lKEdERDA9YWdncmVnYXRlKGJpLm0kR0REMCxsaXN0KGJpLm0kUGVudGFkKSxtZWFuKSR4LAogICAgICAgICAgICAgICAgICAgICAgICAgICBHREQ1PWFnZ3JlZ2F0ZShiaS5tJEdERDUsbGlzdChiaS5tJFBlbnRhZCksbWVhbikkeCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgTVRDTz1hZ2dyZWdhdGUoYmkubSRNVENPLGxpc3QoYmkubSRQZW50YWQpLG1lYW4pJHgsCiAgICAgICAgICAgICAgICAgICAgICAgICAgIE1UV0E9YWdncmVnYXRlKGJpLm0kTVRXQSxsaXN0KGJpLm0kUGVudGFkKSxtZWFuKSR4LAogICAgICAgICAgICAgICAgICAgICAgICAgICBBRVQuUEVUPWFnZ3JlZ2F0ZShiaS5tJEFFVC5QRVQsbGlzdChiaS5tJFBlbnRhZCksbWVhbikkeCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgV2V0LkludGVuc2l0eT1hZ2dyZWdhdGUoYmkubSRXZXQuSW50ZW5zaXR5LGxpc3QoYmkubSRQZW50YWQpLG1lYW4pJHgsCiAgICAgICAgICAgICAgICAgICAgICAgICAgIERyeS5JbnRlbnNpdHk9YWdncmVnYXRlKGJpLm0kRHJ5LkludGVuc2l0eSxsaXN0KGJpLm0kUGVudGFkKSxtZWFuKSR4KQpzdW1tYXJ5KGJpLnVtZikKCm0xIDwtIG9jY3UofiBsaG91cnMgfiAxLCBkYXRhID0gYmkudW1mKQpzdW1tYXJ5KG0xKQoKbTIgPC0gb2NjdSh+IGxob3VycyB+IE1UQ08gKyBNVFdBICsgQUVULlBFVCArIFdldC5JbnRlbnNpdHksIGRhdGEgPSBiaS51bWYpCnN1bW1hcnkobTIpCgojIGZpdHRlZCByZWxhdGlvbnNoaXBzCiMgLS0tLS0tLS0tLS0tLS0tLS0tLS0KCm9wPC1wYXIobWZyb3c9YygyLDIpLCBtYXI9Yyg0LDEsMSwyKSxvbWE9YygxLDUsMSwwKSkKCk1UQ08ucCA8LSBkYXRhLmZyYW1lKE1UQ08gPSBzZXEobWluKGJpLm0kTVRDTyksIG1heChiaS5tJE1UQ08pLCBsZW5ndGg9MTAwKSwgTVRXQSA9IG1lYW4oYmkubSRNVFdBKSwgQUVULlBFVCA9IG1lYW4oYmkubSRBRVQuUEVUKSwgV2V0LkludGVuc2l0eSA9IG1lYW4oYmkubSRXZXQuSW50ZW5zaXR5KSkKcGxvdChNVENPLnAkTVRDTyxwcmVkaWN0KG0yLHR5cGU9J3N0YXRlJyxuZXdkYXRhPU1UQ08ucCkkUHJlZGljdGVkLCB0eXBlPSdsJywgYXhlcz1GLCB4bGFiPSdNVENPJyx5bGltPWMoMCwxKSx5bGFiPScnKQpheGlzKDEsc2VxKG1pbihiaS5tJE1UQ08pLCBtYXgoYmkubSRNVENPKSwgbGVuZ3RoPTMpKQpheGlzKDIsYygwLDAuNSwxKSxsYXM9MSkKCk1UV0EucCA8LSBkYXRhLmZyYW1lKE1UQ08gPSBtZWFuKGJpLm0kTVRDTyksIE1UV0EgPSBzZXEobWluKGJpLm0kTVRXQSksIG1heChiaS5tJE1UV0EpLCBsZW5ndGg9MTAwKSwgQUVULlBFVCA9IG1lYW4oYmkubSRBRVQuUEVUKSwgV2V0LkludGVuc2l0eSA9IG1lYW4oYmkubSRXZXQuSW50ZW5zaXR5KSkKcGxvdChNVFdBLnAkTVRXQSxwcmVkaWN0KG0yLHR5cGU9J3N0YXRlJyxuZXdkYXRhPU1UV0EucCkkUHJlZGljdGVkLCB0eXBlPSdsJyxheGVzPUYsIHhsYWI9J01UV0EnLHlsaW09YygwLDEpLHlsYWI9JycpCmF4aXMoMSwgc2VxKG1pbihiaS5tJE1UV0EpLCBtYXgoYmkubSRNVFdBKSwgbGVuZ3RoPTMpKQoKQUVULlBFVC5wIDwtIGRhdGEuZnJhbWUoTVRDTyA9IG1lYW4oYmkubSRNVENPKSwgTVRXQSA9IG1lYW4oYmkubSRNVFdBKSwgQUVULlBFVCA9IHNlcShtaW4oYmkubSRBRVQuUEVUKSwgbWF4KGJpLm0kQUVULlBFVCksIGxlbmd0aD0xMDApLCBXZXQuSW50ZW5zaXR5ID0gbWVhbihiaS5tJFdldC5JbnRlbnNpdHkpKQpwbG90KEFFVC5QRVQucCRBRVQuUEVULHByZWRpY3QobTIsdHlwZT0nc3RhdGUnLG5ld2RhdGE9QUVULlBFVC5wKSRQcmVkaWN0ZWQsIHR5cGU9J2wnLGF4ZXM9RiwgeGxhYj0nQUVUL1BFVCcseWxpbT1jKDAsMSkseWxhYj0nJykKYXhpcygxLHNlcShtaW4oYmkubSRBRVQuUEVUKSwgbWF4KGJpLm0kQUVULlBFVCksIGxlbmd0aD0zKSkKYXhpcygyLGMoMCwwLjUsMSksbGFzPTEpCgpXZXQuSW50ZW5zaXR5LnAgPC0gZGF0YS5mcmFtZShNVENPID0gbWVhbihiaS5tJE1UQ08pLCBNVFdBID0gbWVhbihiaS5tJE1UV0EpLCBBRVQuUEVUID0gbWVhbihiaS5tJEFFVC5QRVQpLCBXZXQuSW50ZW5zaXR5ID0gc2VxKG1pbihiaS5tJFdldC5JbnRlbnNpdHkpLCBtYXgoYmkubSRXZXQuSW50ZW5zaXR5KSwgbGVuZ3RoPTEwMCkpCnBsb3QoV2V0LkludGVuc2l0eS5wJFdldC5JbnRlbnNpdHkscHJlZGljdChtMix0eXBlPSdzdGF0ZScsbmV3ZGF0YT1XZXQuSW50ZW5zaXR5LnApJFByZWRpY3RlZCwgdHlwZT0nbCcsYXhlcz1GLCB4bGFiPSdXZXQuSW50ZW5zaXR5Jyx5bGltPWMoMCwxKSx5bGFiPScnKQpheGlzKDEsc2VxKG1pbihiaS5tJFdldC5JbnRlbnNpdHkpLCBtYXgoYmkubSRXZXQuSW50ZW5zaXR5KSwgbGVuZ3RoPTMpKQptdGV4dCgiT2NjdXBhbmN5IHByb2JhYmlsaXR5IixzaWRlPTIsbGluZT0zLG91dGVyPVQpCgpwYXIob3ApCgojIHByZWRpY3Rpb25zIGZvciBkcmF3aW5nIHRoZSBtYXAKIyAtLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCnByZWQgPC0gcHJlZGljdChtMiwgdHlwZT0nc3RhdGUnLCBuZXdEYXRhPXNpdGVDb3ZzKGJpLnVtZilbLDM6Nl0sIGFwcGVuZERhdGE9RikKcHJlZCRQZW50YWQgPC0gYWdncmVnYXRlKGJpLm0kR0REMCxsaXN0KGJpLm0kUGVudGFkKSxtZWFuKSRHcm91cC4xCgojIGNvb3JkaW5hdGVzCnByZWQkbGF0IDwtIC0oYXMubnVtZXJpYyhzdWJzdHIocHJlZCRQZW50YWQsMSwyKSkgKyAoYXMubnVtZXJpYyhzdWJzdHIocHJlZCRQZW50YWQsMyw0KSkrMi41KS82MCkKcHJlZCRsb25nIDwtIGFzLm51bWVyaWMoc3Vic3RyKHByZWQkUGVudGFkLDYsNykpICsgKGFzLm51bWVyaWMoc3Vic3RyKHByZWQkUGVudGFkLDgsOSkpKzIuNSkvNjAKCiMgcGxvdCB0aGUgYmFsZCBpYmlzIG1hcApsaWJyYXJ5KHJnZGFsKQpaQW1hcCA8LSByZWFkT0dSKGRzbiA9IGdldHdkKCksIGxheWVyID0iUlNBX1Byb3ZpbmNlc19XR1M4NCIpICMgbWFwIG9mIFNvdXRoIEFmcmljYQpwcm9qNHN0cmluZyhaQW1hcCkKCmNvb3JkaW5hdGVzKHByZWQpPSB+bGF0K2xvbmcKcHJvajRzdHJpbmcocHJlZCkgPC0iK3Byb2o9bG9uZ2xhdCArZGF0dW09V0dTODQgK25vX2RlZnMgK2VsbHBzPVdHUzg0ICt0b3dnczg0PTAsMCwwIgoKbmMgPC0gYXMubnVtZXJpYyhjdXQocHJlZCRQcmVkaWN0ZWQsMjApKQpwYXIobWFyPWMoMSwxLDEsMSkpCnBsb3QobGF0fmxvbmcsIGRhdGE9cHJlZCwgcGNoPTE1LCBjZXg9MC42LCBhbm49RiwgYXhlcz1GLCB4bGltPWMoMTUsMzMpLCB5bGltID0gYygtMzUsIC0yMSksIGNvbD1jb2xvclJhbXBQYWxldHRlKGMoImdyZXkiLCAiZGFya29yYW5nZTMiKSkoMjApW25jXSkKcGxvdChaQW1hcCwgYWRkPVQpCgpwYXIobWFyPWMoMCwwLDAsMCkpCnBsb3QobGF0fmxvbmcsIGRhdGE9YmkubSwgcGNoPTE1LCBjZXg9MC42LCBhbm49RiwgYXhlcz1GLCB4bGltPWMoMTUsMzMpLCB5bGltID0gYygtMzUsIC0yMSksIGNvbD0iZ3JleSIpCnBvaW50cyhsYXR+bG9uZywgZGF0YT1iaS5tW2JpLm0kU3BwPT0xLF0sIHBjaD0xNSwgY2V4PTAuNixjb2w9J2JsYWNrJykKcGxvdChaQW1hcCwgYWRkPVQpCgoKIyBjb2RlIHRoZSBtb2RlbCBpbiBCVUdTCnNpbmsoIm9jY19pYmlzLnR4dCIpCmNhdCgiCiAgICBtb2RlbCB7CiAgICAjIEVjb2xvZ2ljYWwgbW9kZWwgZm9yIHRydWUgb2NjdXJyZW5jZQogICAgZm9yIChpIGluIDE6bnNpdGVzKSB7CiAgICAgIHpbaV0gfiBkYmVybihwc2lbaV0pCiAgICAKICAgICAgIyBPYnNlcnZhdGlvbiBtb2RlbCBmb3IgcmVwbGljYXRlZCBkZXRlY3Rpb24vbm9uZGV0ZWN0aW9uIG9ic2VydmF0aW9ucwogICAgICBmb3IgKGogaW4gMTpuY2FyZHNbaV0pIHsKICAgICAgICBwLmVmZltpLGpdIDwtIHpbaV0gKiBwW2ksal0KICAgICAgICB5W2ksal0gfiBkYmVybihwLmVmZltpLGpdKQogICAgICB9ICNqCiAgICB9ICNpCiAgICAKICAgICMgY292YXJpYXRlcwogICAgZm9yIChpIGluIDE6bnNpdGVzKXsKICAgIGxvZ2l0KHBzaVtpXSkgPSBiZXRhLnBzaVsxXSArIGJldGEucHNpWzJdICogTVRDT1tpXSArIGJldGEucHNpWzNdICogTVRXQVtpXSArIGJldGEucHNpWzRdICogQUVULlBFVFtpXSArIGJldGEucHNpWzVdICogV2V0LkludGVuc2l0eVtpXQogICAgICBmb3IgKGogaW4gMTpuY2FyZHNbaV0pewogICAgICAgIGxvZ2l0KHBbaSxqXSkgPSBiZXRhLnBbMV0gKyBiZXRhLnBbMl0gKiBsaG91cnNbaSxqXQogICAgICB9ICNqCiAgICB9ICNpCgogIAogICAgIyBQcmlvcnMKICAgIGZvciAoYiBpbiAxOjUpewogICAgICBiZXRhLnBzaVtiXSB+IGRub3JtKDAsIDAuMDEpCiAgICB9CgogICAgZm9yIChiIGluIDE6Mil7CiAgICAgIGJldGEucFtiXSB+IGRub3JtKDAsIDAuMDEpCiAgICB9CiAgICB9CiAgICAiLGZpbGwgPSBUUlVFKQpzaW5rKCkKCmxpYnJhcnkoamFnc1VJKQoKeSA8LSBiaS51bWZAeQpkaW0oeSkgICMgMzIyMCBzaXRlcywgbWF4IDcxOSBjYXJkcwoKbmNhcmRzIDwtIGFzLm51bWVyaWModGFibGUoYmkubSRQZW50YWQpKQoKTVRDTz1zY2FsZShhZ2dyZWdhdGUoYmkubSRNVENPLGxpc3QoYmkubSRQZW50YWQpLG1lYW4pJHgpCk1UV0E9c2NhbGUoYWdncmVnYXRlKGJpLm0kTVRXQSxsaXN0KGJpLm0kUGVudGFkKSxtZWFuKSR4KQpBRVQuUEVUPXNjYWxlKGFnZ3JlZ2F0ZShiaS5tJEFFVC5QRVQsbGlzdChiaS5tJFBlbnRhZCksbWVhbikkeCkKV2V0LkludGVuc2l0eT1zY2FsZShhZ2dyZWdhdGUoYmkubSRXZXQuSW50ZW5zaXR5LGxpc3QoYmkubSRQZW50YWQpLG1lYW4pJHgpCgpvY2MuZGF0YSA8LSBsaXN0KHkgPSB5LCBuc2l0ZXMgPSBucm93KHkpLCBuY2FyZHMgPSBuY2FyZHMsCiAgICAgICAgICAgICAgICAgbGhvdXJzPW1hdHJpeChiaS51bWZAb2JzQ292cyRsaG91cnMsIG5yb3c9bnJvdyh5KSwgYnlyb3c9VCksCiAgICAgICAgICAgICAgICAgTVRDTz1hcy5udW1lcmljKE1UQ08pLAogICAgICAgICAgICAgICAgIE1UV0E9YXMubnVtZXJpYyhNVFdBKSwKICAgICAgICAgICAgICAgICBBRVQuUEVUPWFzLm51bWVyaWMoQUVULlBFVCksCiAgICAgICAgICAgICAgICAgV2V0LkludGVuc2l0eT1hcy5udW1lcmljKFdldC5JbnRlbnNpdHkpKQojIEluaXRpYWwgdmFsdWVzCnpzdCA8LSBhcHBseSh5LCAxLCBmdW5jdGlvbih4KXttYXgoeCxuYS5ybT1UKX0pCiMgT2JzZXJ2ZWQgb2NjdXJyZW5jZSBhcyBzdGFydGluZyB2YWx1ZXMgZm9yIHoKaW5pdHMgPC0gZnVuY3Rpb24oKSBsaXN0KHogPSB6c3QsIGJldGEucHNpPXJ1bmlmKDUsLTMsMyksIGJldGEucD1ydW5pZigyLC0zLDMpKQojIFBhcmFtZXRlcnMgbW9uaXRvcmVkCnBhcmFtcyA8LSBjKCJiZXRhLnBzaSIsICJiZXRhLnAiKQoKIyBNQ01DIHNldHRpbmdzCm5jIDwtIDMgICA7ICAgbmkgPC0gNTAwMCAgIDsgICBuYiA8LSAyMDAwICAgOyAgIG50IDwtIDIKCm91dC5vY2MuY292IDwtIGphZ3Mob2NjLmRhdGEsIGluaXRzLCBwYXJhbXMsICJvY2NfaWJpcy50eHQiLCBuLmNoYWlucz1uYywgbi5pdGVyPW5pLCBuLmJ1cm4gPSBuYiwgbi50aGluPW50KQpwbG90KG91dC5vY2MuY292KQpwcmludChvdXQub2NjLmNvdiRzdW1tYXJ5LCBkaWdpdHM9MikKCgojIGZpdHRlZCByZWxhdGlvbnNoaXBzCiMgLS0tLS0tLS0tLS0tLS0tLS0tLS0KbGlicmFyeShib290KQoKb3A8LXBhcihtZnJvdz1jKDIsMiksIG1hcj1jKDQsMSwxLDIpLG9tYT1jKDEsNSwxLDApKQoKTVRDTy5wIDwtIGludi5sb2dpdChvdXQub2NjLmNvdiRzdW1tYXJ5WzEsMV0gKyBvdXQub2NjLmNvdiRzdW1tYXJ5WzIsMV0gKiBNVENPKQpvZCA8LSBvcmRlcihNVENPKQpwbG90KChNVENPW29kXSphdHRyaWJ1dGVzKE1UQ08pW1szXV0pK2F0dHJpYnV0ZXMoTVRDTylbWzJdXSxNVENPLnBbb2RdLCB0eXBlPSdsJywgYXhlcz1GLCB4bGFiPSdNVENPJyx5bGltPWMoMCwxKSx5bGFiPScnLGNvbC5heGlzPSJvcmFuZ2UiLGNvbD0ib3JhbmdlIiwgZmc9Im9yYW5nZSIsIGNvbC5sYWI9Im9yYW5nZSIpCmF4aXMoMSxzZXEobWluKChNVENPKmF0dHJpYnV0ZXMoTVRDTylbWzNdXSkrYXR0cmlidXRlcyhNVENPKVtbMl1dKSwgbWF4KChNVENPKmF0dHJpYnV0ZXMoTVRDTylbWzNdXSkrYXR0cmlidXRlcyhNVENPKVtbMl1dKSwgbGVuZ3RoPTMpLGNvbC5heGlzPSJvcmFuZ2UiLGNvbD0ib3JhbmdlIiwgZmc9Im9yYW5nZSIpCmF4aXMoMixjKDAsMC41LDEpLGxhcz0xLGNvbC5heGlzPSJvcmFuZ2UiLGNvbD0ib3JhbmdlIiwgZmc9Im9yYW5nZSIpCgpNVFdBLnAgPC0gaW52LmxvZ2l0KG91dC5vY2MuY292JHN1bW1hcnlbMSwxXSArIG91dC5vY2MuY292JHN1bW1hcnlbMywxXSAqIE1UV0EpCm9kIDwtIG9yZGVyKE1UV0EpCnBsb3QoKE1UV0Fbb2RdKmF0dHJpYnV0ZXMoTVRXQSlbWzNdXSkrYXR0cmlidXRlcyhNVFdBKVtbMl1dLE1UV0EucFtvZF0sIHR5cGU9J2wnLGF4ZXM9RiwgeGxhYj0nTVRXQScseWxpbT1jKDAsMSkseWxhYj0nJyxjb2wuYXhpcz0ib3JhbmdlIixjb2w9Im9yYW5nZSIsIGZnPSJvcmFuZ2UiLCBjb2wubGFiPSJvcmFuZ2UiKQpheGlzKDEsIHNlcShtaW4oKE1UV0Fbb2RdKmF0dHJpYnV0ZXMoTVRXQSlbWzNdXSkrYXR0cmlidXRlcyhNVFdBKVtbMl1dKSwgbWF4KChNVFdBW29kXSphdHRyaWJ1dGVzKE1UV0EpW1szXV0pK2F0dHJpYnV0ZXMoTVRXQSlbWzJdXSksIGxlbmd0aD0zKSxjb2wuYXhpcz0ib3JhbmdlIixjb2w9Im9yYW5nZSIsIGZnPSJvcmFuZ2UiKQoKQUVULlBFVC5wIDwtIGludi5sb2dpdChvdXQub2NjLmNvdiRzdW1tYXJ5WzEsMV0gKyBvdXQub2NjLmNvdiRzdW1tYXJ5WzQsMV0gKiBBRVQuUEVUKQpvZCA8LSBvcmRlcihBRVQuUEVUKQpwbG90KChBRVQuUEVUW29kXSphdHRyaWJ1dGVzKEFFVC5QRVQpW1szXV0pK2F0dHJpYnV0ZXMoQUVULlBFVClbWzJdXSxBRVQuUEVULnBbb2RdLCB0eXBlPSdsJyxheGVzPUYsIHhsYWI9J0FFVC9QRVQnLHlsaW09YygwLDEpLHlsYWI9JycsY29sLmF4aXM9Im9yYW5nZSIsY29sPSJvcmFuZ2UiLCBmZz0ib3JhbmdlIiwgY29sLmxhYj0ib3JhbmdlIikKYXhpcygxLHNlcShtaW4oKEFFVC5QRVRbb2RdKmF0dHJpYnV0ZXMoQUVULlBFVClbWzNdXSkrYXR0cmlidXRlcyhBRVQuUEVUKVtbMl1dKSwgbWF4KChBRVQuUEVUW29kXSphdHRyaWJ1dGVzKEFFVC5QRVQpW1szXV0pK2F0dHJpYnV0ZXMoQUVULlBFVClbWzJdXSksIGxlbmd0aD0zKSxjb2wuYXhpcz0ib3JhbmdlIixjb2w9Im9yYW5nZSIsIGZnPSJvcmFuZ2UiKQpheGlzKDIsYygwLDAuNSwxKSxsYXM9MSxjb2wuYXhpcz0ib3JhbmdlIixjb2w9Im9yYW5nZSIsIGZnPSJvcmFuZ2UiKQoKV2V0LkludGVuc2l0eS5wIDwtIGludi5sb2dpdChvdXQub2NjLmNvdiRzdW1tYXJ5WzEsMV0gKyBvdXQub2NjLmNvdiRzdW1tYXJ5WzUsMV0gKiBXZXQuSW50ZW5zaXR5KQpvZCA8LSBvcmRlcihXZXQuSW50ZW5zaXR5KQpwbG90KChXZXQuSW50ZW5zaXR5W29kXSphdHRyaWJ1dGVzKFdldC5JbnRlbnNpdHkpW1szXV0pK2F0dHJpYnV0ZXMoV2V0LkludGVuc2l0eSlbWzJdXSxXZXQuSW50ZW5zaXR5LnBbb2RdLCB0eXBlPSdsJyxheGVzPUYsIHhsYWI9J1dldC5JbnRlbnNpdHknLHlsaW09YygwLDEpLHlsYWI9JycsY29sLmF4aXM9Im9yYW5nZSIsY29sPSJvcmFuZ2UiLCBmZz0ib3JhbmdlIiwgY29sLmxhYj0ib3JhbmdlIikKYXhpcygxLHNlcShtaW4oKFdldC5JbnRlbnNpdHlbb2RdKmF0dHJpYnV0ZXMoV2V0LkludGVuc2l0eSlbWzNdXSkrYXR0cmlidXRlcyhXZXQuSW50ZW5zaXR5KVtbMl1dKSwgbWF4KChXZXQuSW50ZW5zaXR5W29kXSphdHRyaWJ1dGVzKFdldC5JbnRlbnNpdHkpW1szXV0pK2F0dHJpYnV0ZXMoV2V0LkludGVuc2l0eSlbWzJdXSksIGxlbmd0aD0zKSxjb2wuYXhpcz0ib3JhbmdlIixjb2w9Im9yYW5nZSIsIGZnPSJvcmFuZ2UiKQptdGV4dCgiT2NjdXBhbmN5IHByb2JhYmlsaXR5IixzaWRlPTIsbGluZT0zLG91dGVyPVQsY29sLmF4aXM9Im9yYW5nZSIsY29sPSJvcmFuZ2UiLCBmZz0ib3JhbmdlIikKCnBhcihvcCkKCiMgcHJlZGljdGlvbnMgZm9yIGRyYXdpbmcgdGhlIG1hcAojIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0KcHJlZCA8LSBkYXRhLmZyYW1lKFByZWRpY3RlZCA9IGludi5sb2dpdChvdXQub2NjLmNvdiRzdW1tYXJ5WzEsMV0gKyBvdXQub2NjLmNvdiRzdW1tYXJ5WzIsMV0gKiBNVENPICsgb3V0Lm9jYy5jb3Ykc3VtbWFyeVszLDFdICogTVRXQSArCiAgICAgICAgb3V0Lm9jYy5jb3Ykc3VtbWFyeVs0LDFdICogQUVULlBFVCArIG91dC5vY2MuY292JHN1bW1hcnlbNSwxXSAqIFdldC5JbnRlbnNpdHkpKQpwcmVkJFBlbnRhZCA8LSBhZ2dyZWdhdGUoYmkubSRHREQwLGxpc3QoYmkubSRQZW50YWQpLG1lYW4pJEdyb3VwLjEKCiMgY29vcmRpbmF0ZXMKcHJlZCRsYXQgPC0gLShhcy5udW1lcmljKHN1YnN0cihwcmVkJFBlbnRhZCwxLDIpKSArIChhcy5udW1lcmljKHN1YnN0cihwcmVkJFBlbnRhZCwzLDQpKSsyLjUpLzYwKQpwcmVkJGxvbmcgPC0gYXMubnVtZXJpYyhzdWJzdHIocHJlZCRQZW50YWQsNiw3KSkgKyAoYXMubnVtZXJpYyhzdWJzdHIocHJlZCRQZW50YWQsOCw5KSkrMi41KS82MAoKIyBwbG90IHRoZSBiYWxkIGliaXMgbWFwCmxpYnJhcnkocmdkYWwpClpBbWFwIDwtIHJlYWRPR1IoZHNuID0gZ2V0d2QoKSwgbGF5ZXIgPSJSU0FfUHJvdmluY2VzX1dHUzg0IikgIyBtYXAgb2YgU291dGggQWZyaWNhCnByb2o0c3RyaW5nKFpBbWFwKQoKY29vcmRpbmF0ZXMocHJlZCk9IH5sYXQrbG9uZwpwcm9qNHN0cmluZyhwcmVkKSA8LSIrcHJvaj1sb25nbGF0ICtkYXR1bT1XR1M4NCArbm9fZGVmcyArZWxscHM9V0dTODQgK3Rvd2dzODQ9MCwwLDAiCgpuYyA8LSBhcy5udW1lcmljKGN1dChwcmVkJFByZWRpY3RlZCwyMCkpCnBhcihtYXI9YygxLDEsMSwxKSkKcGxvdChsYXR+bG9uZywgZGF0YT1wcmVkLCBwY2g9MTUsIGNleD0wLjYsIGFubj1GLCBheGVzPUYsIHhsaW09YygxNSwzMyksIHlsaW0gPSBjKC0zNSwgLTIxKSwgY29sPWNvbG9yUmFtcFBhbGV0dGUoYygiZ3JleSIsICJkYXJrb3JhbmdlMyIpKSgyMClbbmNdKQpwbG90KFpBbWFwLCBhZGQ9VCkKYGBgCgouLi4=