# The Barbary lion model. The data is from pooling expert opinions across different factors, and different experts. See paper below to 
# see how the distribution was derived. 
# Lee, T. E., et al. (2015) "Assessing uncertainty in sighting records: an example of the Barbary lion"

# Input: random numbers drawn from the individual distributions for the quality of each sighting
# Output: The probability that the lion is extinct in each year after the last 'certain' sighting.

#==============================================================================================================================================
library(rjags)
library(R2jags)

Years = seq(1926,2016)
Tall <- 122
t0<-1895

N <- 91 # The model runs for every year after the last certain sighting (N years)
n <- Tall-N
A <- 2 # Fuzz factor

qDist = read.csv("Individual_q_Dist_SumW.csv", header=FALSE); # Random numbers drawn from the distribution for the quality of individual sightings
RR = dim(qDist)[1];
Randi = sample(1:RR, Tall,replace=T);

yyall <- c(1,  0,  0,  0, 
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	1,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0) # Certain sightings
yall <- matrix(0,1,Tall);
yall[1:length(yyall)] = yyall;
zzall <- c(0,  0,	0,	1,	
          0,	1,	1,	0,	
          0,	0,	1,	0,	         
          0,	0,	0,	1,	
          1,	1,	0,	0,	
          0,	0,	1,	0,	
          0,	1,	0,	1,	
          0,	0,	0,	0,	
          0,	0,	1,	1,	
          0,	0,	0,	1,	
          1,	0,	0,	0,	
          1,	0,	0,	1,	
          1,	0,	0,	0,	
          0,	0,	1,	0,	
          0,	0,	0,	0,	
          0,	1,	0,	0,	
          0,	0,	0,	0,	
          0,	0,	0) # Uncertain sightings


zall <- matrix(0,1,Tall);
zall[1:length(zzall)] = zzall;

# When an uncertain sighting occurs, the quality is the mean of the quality distribution. 
# Years following this sighting have a quality drawn randomly from the quality distribution of the most recent uncertain sighting. 

Thisqq <- c(1,1,1,1,
           1,2,3,3,
           3,3,4,4,
           4,4,4,5,
           6,7,7,7,
           7,7,8,8,
           8,9,9,10,
           10,10,10,10,
           10,10,11,12,
           12,12,12,13,
           14,14,14,14,
           15,15,15,16,
           17,17,17,17,
           17,17,18,18,
           18,18,18,18,
           18,19,19,19,
           19,19,19,19,
           19,19,19);

Thisq <- matrix(19,1,Tall);
Thisq[1:length(Thisqq)] = Thisqq;

qzall = matrix(nrow=Tall, ncol=1)
for (k in 1:Tall)
{
  qzall[k]<-qDist[Randi[k],Thisq[k]];
  if (zall[k]==1){qzall[k] = mean(qDist[,Thisq[k]])}
}

#==============================================================================================================================================


Extantmean = matrix(nrow=N, ncol=1)
Extantsd = matrix(nrow=N, ncol=1)
temean = matrix(nrow=N, ncol=1)
tesd = matrix(nrow=N, ncol=1)

for (k in 1:N)
{
  y = yall[0:(n+k)]
  z = zall[0:(n+k)]
  qz = qzall[0:(n+k)]
  T = n+k
  model1.string <-"
  model
  {
  
 # then program from Lee, T. E.,et al. (2013). Inferring extinctions from sighting records of variable reliability. Journal of Applied Ecology, 51 251-258
  te ~ dunif(0,T)  
  m1 ~ dunif(0,100)  
  f1 <-0  # no false sightings in certain sighting record
  m2 ~ dunif(0,100)  
  f2 ~ dunif(0,100) 
  Extant ~ dbern(0.5)  
  
  #Jeff~dbeta(0.5,0.5)
  #Extant ~ dbern(Jeff)  #  Jeff trick
  Extant ~ dbern(0.9) 
  
  for (i in 1:T)
  {
    mm1[i] <- Extant * m1+(1-Extant)*step(te-i) * m1 + f1  
    p1[i] <- 1- exp(-mm1[i])  
    y[i] ~ dbern(p1[i]) 
    
    #mm2[i] <- Extant * mp2 + (1-Extant)*step(te-i) * mp2 + f2
    #p2[i] <- 1 - exp(-mm2[i])
    #z[i] ~ dbern(p2[i])
  
    #added new sectionwarning
    h[i] <- 1/(-A*log(qz[i]))
    swell[i]~dgamma(h[i],h[i])
    m2[i] <- mp2*swell[i]
  
    mm2[i] <- Extant * m2[i] + (1-Extant)*step(te-i) * m2[i] + f2  
    p2[i] <- 1 - exp(-mm2[i])
    z[i] ~ dbern(p2[i])
  }  
}
 "
  
  model1.spec<-textConnection(model1.string)
  
  inits <- function() {list(te=10,m1=10,mp2=10,f2=10,Extant=1)}
    #list(te = 10, 
    #     m1 = runif(1, 0, 100),
    #     m2 = 10,
    #     f2 = 10,
    #     Extant = 1)}
  
  jags <- jags.model(model1.spec,
                     data = list('y' = y,'z' = z,'qz'=qz,'T'=T,'A'=A),
                     inits=inits,
                     n.chains=1)
  
  message(k, sprintf(" out of %s\n", N))
  
  res <- coda.samples(jags,
                      c('Extant', 'te'),
                      n.iter=60000,
                      n.burnin= 30000)
  #class(res)
  #summary(res)
  rmat <- as.matrix(res)
  Extantmean[k] = mean(rmat[,1])
  Extantsd[k] = sd(rmat[,1])
  temean[k] = mean(rmat[,2])
  tesd[k] = sd(rmat[,2])
}

Adf = data.frame(Extantmean,Extantsd, 
                 temean,tesd)
col_headings<-c('Extant_mean','Extant_sd',
                'te_mean','te_sd')
names(Adf)<-col_headings


write.csv(Adf, "Lion_QualityAsVar_9.csv")

#========== plotting================================
png('Lion_QualityAsVar2_9.png')

First0 = min(which(Extantmean <0.9));
Year1 = t0+round(temean[First0]);
Year2 = t0+round(temean[N]);

Extinctmean = 1- Extantmean;

plot(Years, Extinctmean,type="n",
     ylim=range(c(0,1)),
     pch=19, xlab="Years", ylab="Prob. extinct",cex.lab=1.4,
     xaxs="i",yaxs="i", # xaxs=i is to flush frame to axis limits
    # main="Scatter plot with std.dev error bars"
)
lines(Years, Extinctmean, type="l") # line for mean
# hack: draw arrows but with very special "arrowheads"
arrows(Years, Extinctmean-Extantsd, Years, Extinctmean+Extantsd, length=0.05, angle=90, code=3) # errorbars
for (k in 1:N){
  if (zall[n+k]==1) {points(Years[k],Extinctmean[k],pch = 19, bg = "black", cex = 1, lwd = 2)}
}

polygon(c(Years[1],Year1,Year1,Years[1]), c(1,1,0,0),col=rgb(0.5, 0.5, 0.5,0.25), border=NA)
polygon(c(Years[1],Year2,Year2,Years[1]), c(1,1,0,0),col=rgb(0.5, 0.5, 0.5,0.25), border=NA)

dev.off()