#Bayes 1 model with up to 4 joinpoints
model{
  for(i in 1:n){
    Obs[i]~dpois(mean[i])

    #Modelling of the expected cases for every period
    #the maximum and minimum bounds are imposed to avoid numerical overflow.
    #x1 and x2 are the initial and final years of study.
    log(mean[i])<-log(Pob[i])+alpha+beta0*(anyo[i]-(x1+x2)/2)+
    min(max(delta[1]*beta[1]*BaseJoin[i,1]+
    delta[2]*beta[2]*BaseJoin[i,2]+
    delta[3]*beta[3]*BaseJoin[i,3]+
    delta[4]*beta[4]*BaseJoin[i,4],-10),10)
    
    #estimated rate at every period
    rate[i]<-exp(alpha+beta0*(anyo[i]-(x1+x2)/2)+min(max(delta[1]*beta[1]*BaseJoin[i,1]+
    delta[2]*beta[2]*BaseJoin[i,2]+delta[3]*beta[3]*BaseJoin[i,3]+delta[4]*beta[4]*BaseJoin[i,4],-10),10))
  }

  #APC at every period
  for(i in 1:(n-1)){APC[i]<-100*(rate[i+1]-rate[i])/rate[i]}

  #Prior distributions for joinpoints' locations
  for(j in 1:4){tau[j]~dunif(tauinf,tausup)}
  tauinf<-x1+2 
  tausup<-x2-2
  #This condition imposes an order restriction and a minimimum distance between joinpoints
  #one has to be given as data and necesarily has to be equal to 1.
  one~dbern(condicion)
  condicion<-step(tau[2]-tau[1]-2)*step(tau[3]-tau[2]-2)*step(tau[4]-tau[3]-2)

  #flat prior distributions for alpha and beta0  
  alpha~dflat()
  beta0~dflat()
  
  #Prior distribution for parameter beta
  beta[1:4]~dmnorm(mean.beta[1:4],prec.u[1:4,1:4])
  for(i in 1:4){
    mean.beta[i]<-0
    for(k in 1:n){BaseJoin2[k,i]<-BaseJoin[k,i]*pow(w[k],1/2)}
    #prior precision matrix for beta.
    #The parameter g makes the prior distribution of beta to be a Cauchy distribution
    for(j in 1:4){
      prec.u[i,j]<-g*inprod(BaseJoin2[,i],BaseJoin2[,j])*max(step(i-j)*step(j-i),delta[i]*delta[j])
    }
  }
  #Prior distribution for g. 
  g~dgamma(0.5,n2)
  n2<-n/2
  
  #Parameter w for the prior precision of beta
  for(j in 1:n){
    w[j]<-Pob[j]*exp(alpha+beta0*(anyo[j]-(x1+x2)/2))
  }
  
  #########################
  #Break-point definitions#
  #########################
  for(j in 1:4){
    for(i in 1:n){
      BaseJoin[i,j]<-1+(anyo[i]-tau[j])*gamma[i,j]
      gamma[i,j]<-gamma1[j]+(gamma2[j]-gamma1[j])*step(anyo[i]-tau[j])
    }
    gamma1[j]<-(d0*c2[j]-c0*d2[j])/(d2[j]*c1[j]-d1[j]*c2[j])
    gamma2[j]<-(d1[j]*c0-d0*c1[j])/(d2[j]*c1[j]-d1[j]*c2[j])
  
    c1[j]<-inprod(anyocent[,j],primerperiodo[,j])
    c2[j]<-inprod(anyocent[,j],segundoperiodo[,j])
    d1[j]<-inprod(anyocent2[,j],primerperiodo[,j])
    d2[j]<-inprod(anyocent2[,j],segundoperiodo[,j])

    for(i in 1:n){
      anyocent[i,j]<-anyo[i]-tau[j]
      anyocent2[i,j]<-anyo[i]*anyocent[i,j]
      primerperiodo[i,j]<-step(-anyocent[i,j])
      segundoperiodo[i,j]<-1-primerperiodo[i,j]
    }

  }
  c0<-n;d0<-sum(anyo[])

  ################################  
  # Prior distributions for delta#  
  ################################
  #posible configurations for vector delta
  auxdelta[1,1]<-1;  auxdelta[2,1]<-1;  auxdelta[3,1]<-1;  auxdelta[4,1]<-1
  auxdelta[1,2]<-1;  auxdelta[2,2]<-1;  auxdelta[3,2]<-1;  auxdelta[4,2]<-0
  auxdelta[1,3]<-1;  auxdelta[2,3]<-1;  auxdelta[3,3]<-0;  auxdelta[4,3]<-1
  auxdelta[1,4]<-1;  auxdelta[2,4]<-0;  auxdelta[3,4]<-1;  auxdelta[4,4]<-1    
  auxdelta[1,5]<-0;  auxdelta[2,5]<-1;  auxdelta[3,5]<-1;  auxdelta[4,5]<-1
  auxdelta[1,6]<-1;  auxdelta[2,6]<-1;  auxdelta[3,6]<-0;  auxdelta[4,6]<-0
  auxdelta[1,7]<-1;  auxdelta[2,7]<-0;  auxdelta[3,7]<-1;  auxdelta[4,7]<-0
  auxdelta[1,8]<-0;  auxdelta[2,8]<-1;  auxdelta[3,8]<-1;  auxdelta[4,8]<-0
  auxdelta[1,9]<-1;  auxdelta[2,9]<-0;  auxdelta[3,9]<-0;  auxdelta[4,9]<-1
  auxdelta[1,10]<-0;  auxdelta[2,10]<-1;  auxdelta[3,10]<-0;  auxdelta[4,10]<-1
  auxdelta[1,11]<-0;  auxdelta[2,11]<-0;  auxdelta[3,11]<-1;  auxdelta[4,11]<-1
  auxdelta[1,12]<-1;  auxdelta[2,12]<-0;  auxdelta[3,12]<-0;  auxdelta[4,12]<-0  
  auxdelta[1,13]<-0;  auxdelta[2,13]<-1;  auxdelta[3,13]<-0;  auxdelta[4,13]<-0
  auxdelta[1,14]<-0;  auxdelta[2,14]<-0;  auxdelta[3,14]<-1;  auxdelta[4,14]<-0
  auxdelta[1,15]<-0;  auxdelta[2,15]<-0;  auxdelta[3,15]<-0;  auxdelta[4,15]<-1
  auxdelta[1,16]<-0;  auxdelta[2,16]<-0;  auxdelta[3,16]<-0;  auxdelta[4,16]<-0
  #delta will have as a value one of the former configurations
  for(i in 1:4){delta[i]<-auxdelta[i,comp]}
  #The probability of every configuration is given by p
  comp~dcat(p[])

  #every value of p is defined as the probability of the corresponding configuration of deltas
  for(i in 1:16){p[i]<-q[i]/sum(q[])} 
  q[1]<-p1*p2*p3*p4;
  q[2]<-p1*p2*p3*(1-p4)/4;q[3]<-p1*p2*(1-p3)*p4/4;q[4]<-p1*(1-p2)*p3*p4/4;q[5]<-(1-p1)*p2*p3*p4/4;
  q[6]<-p1*p2*(1-p3)*(1-p4)/6;q[7]<-p1*(1-p2)*p3*(1-p4)/6;q[8]<-(1-p1)*p2*p3*(1-p4)/6;q[9]<-p1*(1-p2)*(1-p3)*p4/6;q[10]<-(1-p1)*p2*(1-p3)*p4/6;
  q[11]<-(1-p1)*(1-p2)*p3*p4/6
  q[12]<-p1*(1-p2)*(1-p3)*(1-p4)/4;q[13]<-(1-p1)*p2*(1-p3)*(1-p4)/4;q[14]<-(1-p1)*(1-p2)*p3*(1-p4)/4;q[15]<-(1-p1)*(1-p2)*(1-p3)*p4/4
  q[16]<-(1-p1)*(1-p2)*(1-p3)*(1-p4)

  #prior probabilities of any joinpoint being included in the model
  p1~dbeta(paux1,paux2)
  p2~dbeta(paux1,paux2)
  p3~dbeta(paux1,paux2)
  p4~dbeta(paux1,paux2)

  #Indicator variable pointing which of the models M_0,M_1,M_2,M_3 or M_4 is sampled
  ModelSampled<-sum(delta[])
}