#Upslope migration of snow avalanches in a warming climate 
#Complete data and model source files
##15 August 2021
#
#Florie Giacona1, Nicolas Eckert1, Christophe Corona2, Robin Mainieri1, Samuel Morin3, Markus Stoffel4,5,6, Brice Martin7, Mohamed Naaim8
#1Univ. Grenoble Alpes, INRAE, UR ETNA, Grenoble, France
#2GEOLAB, UMR6042 CNRS / Université Clermont-Auvergne, Maison des sciences de l’homme, Clermont-Ferrand, France
#3Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, CEN, 38000 Grenoble, France
#4Institute for Environmental Sciences - University of Geneva, Geneva, Switzerland
#5Department F.-A. Forel for Environmental and Aquatic Sciences - University of Geneva, Geneva, Switzerland
#6Department for Earth Sciences - University of Geneva, Geneva, Switzerland
#7CRESAT, Université de Haute Alsace, Mulhouse, France
#8Département AQUA, INRAE, Paris, France
#
#Model source file, to be run under Winbugs. The model is fully described in supplement S2
#
#Partial reference for the source paper:
#Giacona, F., Eckert, N., Corona, C., Mainieri, R., Morin, S., Stoffel, M., Martin, B., Naaim, M. (2021). 
#Upslope migration of snow avalanches in a warming climate. Proceedings of the National Academy of Sciences America, DOI: 10.1073/pnas.2107306



Model{  

#Observation level and latent variable specification
  for (i in 1 : N) {
  	for (t in 1:T) {
       y[i,t]  ~ dpois(lambda[i,t])
	   log(lambda[i,t]) <- log(e[i,t])+V[i]+beta[t]}
	   V[i] ~ dnorm(0, tauv)}

#Temporel term
for (t in 1:T) {
beta[t]<-G[t]+Z[t]  
Z[t]~ dnorm(0, tauz)}

#Smooth trend
G[1:T]~car.normal(adj[],weights[],num[],taug)

#Weigths and adjacency matrix for the tmooth trend
 for(t in 1:1) {
 weights[t] <- 2; adj[t] <- t+1
 weights[t+1] <- -1; adj[t+1] <- t+2; num[t] <- 2
 }
for(t in 2:2) {
weights[t+1] <- 2; adj[t+1] <- t-1
weights[t+2] <- 4; adj[t+2] <- t+1
weights[t+3] <- -1; adj[t+3] <- t+2; num[t] <- 3
}
for(t in 3:(T-2)) {
weights[6+(t-3)*4] <- -1; adj[6+(t-3)*4] <- t-2
weights[7+(t-3)*4] <- 4; adj[7+(t-3)*4] <- t-1
weights[8+(t-3)*4] <- 4; adj[8+(t-3)*4] <- t+1
weights[9+(t-3)*4] <- -1; adj[9+(t-3)*4] <- t+2; num[t] <- 4
}
for(t in (T-1):(T-1)) {
weights[(T-4)*4 + 6] <- 2; adj[(T-4)*4 + 6] <- t+1
weights[(T-4)*4 + 7] <- 4; adj[(T-4)*4 + 7] <- t-1
weights[(T-4)*4 + 8] <- -1; adj[(T-4)*4 + 8] <- t-2; num[t] <- 3
}
for(t in T:T) {
weights[(T-4)*4 + 9] <- 2; adj[(T-4)*4 + 9] <- t-1
weights[(T-4)*4 + 10] <- -1; adj[(T-4)*4 + 10] <- t-2; num[t] <- 2
}

#Priors
taug~dgamma(0.01,0.01)
tauv  ~ dgamma(0.01, 0.001)
tauz  ~ dgamma(0.01, 0.001)

#Standard deviations
sigmav <- sqrt(1 / tauv)
sigmaz <- sqrt(1 / tauz)
sigmag<- sqrt(1 / taug)

#Global means
vmoy<- (1/N)*sum(V[1:N])
gmoy<- (1/T)*sum(G[1:T])
zmoy<- (1/T)*sum(Z[1:T])
alpha<-gmoy+vmoy+zmoy

#Variance ratios
var.G <- sd(G[])*sd(G[])
f.temp <- var.G/(var.G+pow(sigmaz,2))
frac.expliquee.temps<-(var.G)/(var.G+pow(sigmaz,2)+pow(sigmav,2))
frac.bruit.espace<-(pow(sigmav,2))/(var.G+pow(sigmaz,2)+pow(sigmav,2))
frac.bruit.temps<-(pow(sigmaz,2))/(var.G+pow(sigmaz,2)+pow(sigmav,2))


#Estimation of annual frequencies and real counts

for (i in 1 : N) {
yireal[i]<- T*estat[i]*exp(gmoy+zmoy+V[i])
fireal[i]<-yireal[i]/(T*c[i])
yi[i]<-sum(y[i,1:T])}

 for (t in 1 : T) {
ytreal[t]<-sum(estat[1:N])*exp(vmoy+beta[t])
ftreal[t]<-ytreal[t]/sum(c[1:N])
ytrealliss[t]<-sum(estat[1:N])*exp(vmoy+G[t])
ftrealliss[t]<-ytrealliss[t]/sum(c[1:N])}

ytotreal1 <- sum(yireal[1:N])
ytotreal2 <- sum(ytreal[1:T])
ytotreal3 <- sum(ytrealliss[1:T])
ftotreal1 <- ytotreal1/(T*sum(c[1:N]))
ftotreal2 <- ytotreal2/(T*sum(c[1:N]))
ftotreal3 <- ytotreal3/(T*sum(c[1:N]))

#Estimation under stationarity assumption
ytotstat <- T*sum(estat[1:N])*exp(vmoy+gmoy+zmoy)
ftotstat<- ytotstat/(T*sum(c[1:N]))

#Estimation of total number y (just taking into account missing values)
ytot<-sum(yi[1:N])

}