#*****************************************************************************************
# INTEGRATED POPULATION MODEL
#
#Model structure:
#  3 age classes (fledgling, yearling & adult)
#  Reproduction possible from second year of life
#  Post-breeding census
#  All rates are time-dependent, survival is age dependent and recapture is sex dependent  
#
#Data:
# 1) Population counts:  count.T = annual counts of fledglings (NAs) and yearling+adults 
# 2) Reproduction:	
#    nbf = annual number of females
#    nfem = annual number of females
#    J = annual number of fledglings per breeding female
# 3) Mark-recapture data = y (individual capture histories)
#*****************************************************************************************
    model{
#********************************
# 1: Priors and constraints
#*******************************
  # Recapture
    for (g in 1:2){
      for (t in 1:(n.occasions-1)){
        p[g,t] <- 1 / (1 + exp(-(mean.p[g] + epsilon.p[t])))
      } #t
      mean.p[g] ~ dunif(-5,5) # prior for mean, sex-specific recapture rate
    } #g 
    for (t in 1:(n.occasions-1)){epsilon.p[t] ~ dnorm(0,tau.p)} #t
    sigma.p ~ dunif(0,10)
    tau.p <- pow(sigma.p,-2)
    sigma2.p <- pow(sigma.p, 2) 
    
  # Survival 
    for (u in 1:3){
      for (t in 1:(n.occasions-1)){
        s[u,t] <- 1 / (1 + exp(-(mu.s[u] + epsilon[t,u])))
      } #t
    } #u
    
    for (u in 1:3){mu.s[u] ~ dunif(-5,5)} #u  # prior for mean age-specific survival rates
    
## Vital rate priors
    xi.1 ~ dunif(0,10)
    xi.2 ~ dunif(0,10)
    xi.3 ~ dunif(0,10)
    xi.4 ~ dunif(0,10)
    xi.5 ~ dunif(0,10)
    for (t in 1:(n.occasions-1)){
      eps.raw[t,1:5] ~ dmnorm(zeros[], tau.raw[,])
      epsilon[t,1] <- xi.1 * eps.raw[t,1] 
      epsilon[t,2] <- xi.2 * eps.raw[t,2]
      epsilon[t,3] <- xi.3 * eps.raw[t,3] 
      epsilon[t,4] <- xi.4 * eps.raw[t,4]
      epsilon[t,5] <- xi.5 * eps.raw[t,5]
    } #t
    
    tau.raw[1:5,1:5] ~ dwish(W[,], 6)
    sigma.raw[1:5,1:5] <- inverse(tau.raw[,])
    
  # Temporal variances
    sigma.s1 <- xi.1*sqrt(sigma.raw[1,1])
    sigma.s2 <- xi.2*sqrt(sigma.raw[2,2])
    sigma.s3 <- xi.3*sqrt(sigma.raw[3,3])
    sigma.s4 <- xi.4*sqrt(sigma.raw[4,4])
    sigma.s5 <- xi.5*sqrt(sigma.raw[5,5])

  # Temporal covariances 
    rho.s.12 <- sigma.raw[1,2]/sqrt(sigma.raw[1,1] * sigma.raw[2,2])
    rho.s.13 <- sigma.raw[1,3]/sqrt(sigma.raw[1,1] * sigma.raw[3,3])
    rho.s.23 <- sigma.raw[2,3]/sqrt(sigma.raw[2,2] * sigma.raw[3,3])
    rho.s.35 <- sigma.raw[3,5]/sqrt(sigma.raw[3,3] * sigma.raw[5,5])
    rho.s.45 <- sigma.raw[4,5]/sqrt(sigma.raw[4,4] * sigma.raw[5,5])
    
  # Unmarked rates
    for (t in 1:(n.occasions-1)){
      unm.A[t] ~ dunif(0,1) # prior for mean yearling and adults unmarked rate
      unm.F[t] ~ dunif(0,1) # prior for mean fledgling unmarked rate
    } #t
    
  # Observation error 
    # Peterson
    tauy.p ~ dgamma(0.001, 0.001)
    sigma2.y.p <- 1 / tauy.p
    
    # Fjord
    tauy.f ~ dgamma(0.001, 0.001)
    sigma2.y.f <- 1 / tauy.f
    
  # Initial population sizes
    nf ~ dnorm(100,0.0001)I(0,)  # fledglings
    Nf[1] <- round(nf)
    ny ~ dnorm(100,0.0001)I(0,)  # yearlings
    Ny[1] <- round(ny)
    na ~ dnorm(100,0.0001)I(0,)  # adults
    Na[1] <- round(na)
    UNMARKED.A[1] ~ dnorm(100,0.0001)I(0,) 	# Unmarked yearlings + unmarked adults
    UNMARKED.F[1] ~ dnorm(100,0.0001)I(0,) 	# Unmarked fledglings
    
#**************************************************
# 2: Likelihoods
#**************************************************
  # Reproduction data
    for (t in 1:(n.occasions-1)){
    # R
      nbf[t] ~ dbin(d[t],nfem[t])
      d[t] <-  1 / (1 + exp(-(mu.R + epsilon[t,5])))
    # Fecundity
      J[t] ~ dpois(R[t])
    R[t] <- exp(mu.f + epsilon[t,4] + log(nbf[t]))
    fec[t] <-  exp(mu.f + epsilon[t,4])
    } #t
    mu.R ~ dunif(-5,5) # prior for mean R 
    mu.f ~ dunif(-5,3) # prior for mean fec
    
  # Capture-recapture data
    for (i in 1:nind){
      z[i,f[i]] <- 1
      for (t in (f[i]+1):n.occasions){
    # State process
        z[i,t] ~ dbern(mu1[i,t])
        mu1[i,t] <- s[x1[i,t-1],t-1] * z[i,t-1]
    # Observation process
        y[i,t] ~ dbern(mu2[i,t])
        mu2[i,t] <- p[sex[i],t-1] * z[i,t]
      }#t
    }#i
    
  # Population count data
    # System process 
    for (t in 2:n.occasions){ 
      Nf.temp[t] <- (s[2,t-1]*fec[t-1]*0.5*d[t-1]*Ny[t-1]) + (s[3,t-1]*fec[t-1]*0.5*d[t-1]*Na[t-1]) 
      Nf[t] ~ dpois(Nf.temp[t])
      Ny[t] ~ dbin(s[1,t-1],Nf[t-1])
      mean2[t] ~ dbin(s[2,t-1],Ny[t-1]) 
      mean3[t] ~ dbin(s[3,t-1],Na[t-1])
      Na[t] <- mean2[t] + mean3[t]
      UNMARKED.A[t] <- Ntot.A[t-1]*unm.A[t-1]
      UNMARKED.F[t] <- Ntot.F[t-1]*unm.F[t-1]
    } #t
    
    for (t in 1:n.occasions){ 
      Ntot.A[t]<- Ny[t]+ Na[t] + UNMARKED.A[t]
      Ntot.F[t] <- Nf[t] + UNMARKED.F[t]
      Ntot[t] <- Ntot.A[t] + Ntot.F[t]
      log.Ntot.A[t] <- log(Ntot.A[t])
    }
    
    # Observation process
    # Petersons index
    for (t in 1:6){ 
    counts[t,1] ~ dlnorm(log.Ntot.A[t],tauy.p)
    counts[t,2] ~ dnorm(Ntot.F[t],tauy.p)
    }
    
    # Fjord counts
    for (t in 7:n.occasions){ 
    counts[t,1] ~ dlnorm(log.Ntot.A[t],tauy.f)
    counts[t,2] ~ dnorm(Ntot.F[t],tauy.f)
    }
    }