# February 22 2005 This code will maximize the likelihood for the poisson
# migration model with tag loss.  The data is from the yellowtail
# flounder provided by DFO.  There are 2 strata, 5 years, 10 time periods.
# Working in directory: d:contract\Walsh\R Code

# Packages required to load for this program: boot





library(boot)

# Parameters in the model
# Lambda- tag retention rates dependent on time not strata
# S- survival (both natural and otherwise) and migration (between strata)
# lam- reporting rates depending on time, not strata and tag type
# p- exploitation rate, dependent on strata and time

# Statistics
# N- number of animals released by stratum, period and tag type
# R- number of animals recovered by stratum, period and tag type

# Functions


qaic<-function(observed, expected, over, nparm){
#Calculates the QAIC value for a given model.
   q<--2*loglikelihood(observed,expected)/over + 2*(nparm+1)
   q
} #End qaic


over <- function(obs, expected, cutoff, max.cutoff, nparm) {
   res <- (obs-expected)**2 / ifelse(expected>0,expected, 1)
   over<- sum(ifelse(expected>cutoff,ifelse(expected<max.cutoff,res,0),0))
   div <- sum(ifelse(expected>cutoff,ifelse(expected<max.cutoff,1,0),0))
   over <- over / (div - nparm)
   dev<-2*(loglikelihood(obs,obs)-loglikelihood(obs,expected))
#   dev<-2*(-(obs-expected)+obs*log(ifelse(obs>0, obs,1)/ifelse(expected>0,expected,1)))
#   sum.dev<-sum(ifelse(expected>cutoff,ifelse(expected<max.cutoff,dev,0),0))
   df<-div-nparm
    deviance<-dev/df
#   deviance<-sum.dev/df   
   return(over, df, div,nparm,  over*(df), deviance)
} #End over




create_residual<-function(beta.pack,nbeta, nsample,ntype,nstrata, R_observed, design,N)
# Create the residuals of a model from the expected and observed counts.
# (observed-expected) / sqrt(expected)
{  
   num_parm<-sum(nbeta$capture+nbeta$surv+nbeta$tag+nbeta$migrate+nbeta$report)
   #beta.pack<-beta_final
   Beta<-list()
   Beta$capture<-beta.pack[1:nbeta$capture]
   Beta$surv<-beta.pack[nbeta$capture+1:nbeta$surv]
   Beta$tag<-beta.pack[nbeta$capture+nbeta$surv+1:nbeta$tag]
   Beta$migrate<-beta.pack[nbeta$capture+nbeta$surv+nbeta$tag+1:nbeta$migrate]
   Beta$report<-beta.pack[nbeta$capture+nbeta$surv+nbeta$tag+nbeta$migrate+1:nbeta$report]

# Transform betas to standard parameters.  This gives me p*, surv*, tag*, migrate*  
# where the * indicates that there may not be as many standard parameters as for the
# saturated model.  For model 1, there will be only 1 standard parameter for each component.
   new<-list()
   new$capture<-transform_beta_to_std(Beta$capture)
   new$surv<-transform_beta_to_std(Beta$surv)
   new$migrate<-transform_beta_to_std(Beta$migrate)
   new$tag<-transform_beta_to_std(Beta$tag)
   new$report<-transform_beta_to_std(Beta$report)

# Make large matrices for multiplication
   big<-list()
   big$capture<-design$capture%*%new$capture
   big$surv<-design$surv%*%new$surv
   big$migrate<-design$migrate%*%new$migrate
   big$tag<-design$tag%*%new$tag
   big$report<-design$report%*%new$report

# Create matrices for calculating E(R)
   std<-list()
   std$surv<-create_parameter_matrix(nsample-1,nstrata, big$surv)
   std$capture<-create_parameter_matrix_capture(nsample, nstrata, big$capture)
   std$tag<-create_parameter_matrix(nsample-1,nstrata, big$tag)
   std$report <-create_report_matrix(nsample,nstrata,ntype,big$report)
   std$migrate <-create_migration_matrix(nsample,nstrata, big$migrate)  


# Calculate residual (o-e)/sqrt(e) and Overdisperison (o-e)**2/e  / df =number of cells-#estimated parameters

   R_expected<-Expected_Recovery(ntype,nsample,nstrata,N,std)
   residual<-array(rep(0,times=nsample*(nsample-1)*ntype*nstrata*nstrata))
   overdispersion<-array(rep(0,times=nsample*(nsample-1)*ntype*nstrata*nstrata))
   index<-1
    count<-0
   for (i in seq(1,nsample,by=2)){
      for (j in (i+1):nsample){
         for (k in 1:ntype){
            for (s in 1:nstrata){
               for (t in 1:nstrata){
                  #print(c(i,j,k,s,t,index,count))
                  # double tagged fish were not released in 2004 (period 9) Thus expected counts do not exist.
                  if (k==2 | k==3){
                     if(i!=9){residual[index]<-(R_observed[t,s,k,i,j]-R_expected[t,s,k,i,j])/sqrt(R_expected[t,s,k,i,j])                           
                              overdispersion[index]<-((R_observed[t,s,k,i,j]-R_expected[t,s,k,i,j])^2/R_expected[t,s,k,i,j])/num_parm
                              index<-index+1
                              count<-count+1}
                  }
                  else{ residual[index]<-(R_observed[t,s,k,i,j]-R_expected[t,s,k,i,j])/sqrt(R_expected[t,s,k,i,j])
                        overdispersion[index]<-((R_observed[t,s,k,i,j]-R_expected[t,s,k,i,j])^2/R_expected[t,s,k,i,j])/num_parm
                        index<-index+1
                        count<-count+1}     
                  if (abs(residual[index-1])>3){print(c(residual[index-1], R_observed[t,s,k,i,j], R_expected[t,s,k,i,j],t,s,k,i,j))}                                  
               }
            }
         }
      }
   }
   residual2<-residual[1:count]
   overdispersion2<-overdispersion[1:count]
   
   #pooled over all but tag type
   index<-1
   count<-0
   for (k in 1:ntype){
      residual[index]<-(sum(R_observed[1:nstrata,1:nstrata,k,1:nsample,1:nsample])-sum(R_expected[1:nstrata,1:nstrata,k,1:nsample,1:nsample]))/sqrt(sum(R_expected[1:nstrata,1:nstrata,k,1:nsample,1:nsample]))
      index<-index+1
      count<-count+1     
   }
  type<-residual[1:count] 
  
  #pooled over all but release period
  index<-1
   count<-0
   for (i in seq(1,nsample,by=2)){
      residual[index]<-(sum(R_observed[1:nstrata,1:nstrata,1:ntype,i,1:nsample])-sum(R_expected[1:nstrata,1:nstrata,1:ntype,i,1:nsample]))/sqrt(sum(R_expected[1:nstrata,1:nstrata,1:ntype,i,1:nsample]))
      index<-index+1
      count<-count+1     
   }
  mperiod<-residual[1:count]

   
   #pooled over all but recover period
  index<-1
   count<-0
   for (j in 2:nsample){
      residual[index]<-(sum(R_observed[1:nstrata,1:nstrata,1:ntype,1:nsample,j])-sum(R_expected[1:nstrata,1:nstrata,1:ntype,1:nsample,j]))/sqrt(sum(R_expected[1:nstrata,1:nstrata,1:ntype,1:nsample,j]))
      index<-index+1
      count<-count+1     
   }
  rperiod<-residual[1:count]
  
  #pooled over all but release strata
  index<-1
   count<-0
   for (t in 1:nstrata){
      residual[index]<-(sum(R_observed[t,1:nstrata,1:ntype,1:nsample,1:nsample])-sum(R_expected[t,1:nstrata,1:ntype,1:nsample,1:nsample]))/sqrt(sum(R_expected[t,1:nstrata,1:ntype,1:nsample,1:nsample]))
      index<-index+1
      count<-count+1     
   }
  mstrata<-residual[1:count]
  
  #pooled over all but recover strata
  index<-1
   count<-0
   for (s in 1:nstrata){
      residual[index]<-(sum(R_observed[1:nstrata,s,1:ntype,1:nsample,1:nsample])-sum(R_expected[1:nstrata,s,1:ntype,1:nsample,1:nsample]))/sqrt(sum(R_expected[1:nstrata,s,1:ntype,1:nsample,1:nsample]))
      index<-index+1
      count<-count+1     
   }
  rstrata<-residual[1:count]
  
   return(residual2,R_expected,overdispersion2,type,mperiod,rperiod,mstrata,rstrata)
} #End create_residual    



transform_to_std<-function(X,beta)
# Transform the beta estimates to standard estimates by first multiplying by the large design matrix.
{
   ans<-X%*%beta
   ans<-transform_beta_to_std(ans)
   ans
}

make.large.design<-function(design,nbeta,nstd)
# Make a large design matrix consiting of all design matrices for all parameters. 
# I need this to find the covariances of the standard parameters.
{
   sum.nbeta<-sum(nbeta$capture,nbeta$surv,nbeta$tag,nbeta$migrate,nbeta$report)
   sum.nstd<-sum(nstd$capture,nstd$surv,nstd$tag,nstd$migrate,nstd$report)
   big.design<-array(rep(0,times=sum.nbeta*sum.nstd),dim=c(sum.nstd,sum.nbeta))
   
   big.design[1:nstd$capture,1:nbeta$capture]<-design$capture
   big.design[(nstd$capture+1):(nstd$capture+nstd$surv), (nbeta$capture+1):(nbeta$capture+nbeta$surv)]<-design$surv
   big.design[(nstd$capture+nstd$surv+1):(nstd$capture+nstd$surv+nstd$tag), (nbeta$capture+nbeta$surv+1):(nbeta$capture+nbeta$surv+nbeta$tag)]<-design$tag      
   big.design[(nstd$capture+nstd$surv+nstd$tag+1):(nstd$capture+nstd$surv+nstd$tag+nstd$migrate), (nbeta$capture+nbeta$surv+nbeta$tag+1):(nbeta$capture+nbeta$surv+nbeta$tag+nbeta$migrate)]<-design$migrate
   big.design[(nstd$capture+nstd$surv+nstd$tag+nstd$migrate+1):sum.nstd, (nbeta$capture+nbeta$surv+nbeta$tag+nbeta$migrate+1):sum.nbeta]<-design$report
   big.design
 }        


make.large.design2<-function(design,nbeta,nstd)
# Make a large design matrix consiting of all design matrices for all parameters. 
# I need this to find the covariances of the standard parameters.
# THis adds one parameter for overdispersion
{
   sum.nbeta<-sum(nbeta$capture,nbeta$surv,nbeta$tag,nbeta$migrate,nbeta$report,nbeta$overdispersion)
   sum.nstd<-sum(nstd$capture,nstd$surv,nstd$tag,nstd$migrate,nstd$report,1)
   big.design<-array(rep(0,times=sum.nbeta*sum.nstd),dim=c(sum.nstd,sum.nbeta))
   
   big.design[1:nstd$capture,1:nbeta$capture]<-design$capture
   big.design[(nstd$capture+1):(nstd$capture+nstd$surv), (nbeta$capture+1):(nbeta$capture+nbeta$surv)]<-design$surv
   big.design[(nstd$capture+nstd$surv+1):(nstd$capture+nstd$surv+nstd$tag), (nbeta$capture+nbeta$surv+1):(nbeta$capture+nbeta$surv+nbeta$tag)]<-design$tag      
   big.design[(nstd$capture+nstd$surv+nstd$tag+1):(nstd$capture+nstd$surv+nstd$tag+nstd$migrate), (nbeta$capture+nbeta$surv+nbeta$tag+1):(nbeta$capture+nbeta$surv+nbeta$tag+nbeta$migrate)]<-design$migrate
   big.design[(nstd$capture+nstd$surv+nstd$tag+nstd$migrate+1):(sum.nstd-1), (nbeta$capture+nbeta$surv+nbeta$tag+nbeta$migrate+1):(sum.nbeta-1)]<-design$report
   big.design[sum.nstd,sum.nbeta]<-1
   big.design
 } # End make.large.design2


print_count <- function( counts, digits=1, width=5, format="f")
{
#     print out the count in a nice format
#     This function assumes s=t=2, 4 tag_types, and yrel=yrec=10
#
#     first convert counts to character strings to get ready for printing
#
       char_count <- formatC(counts, digits, width, format)
#
#     now to print out nicely by each tag type and then for each combination of s and t
       for(tag_type in 1:4){
          cat("\n\nTag Type", tag_type, "\n")
          cat("   ",formatC(c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10), digits=0, width, format="d"),"\n")
          cat("   ",formatC(c(1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2 ),  digits=0, width, format="d"),"\n")
          for(yrel in 1:10){
             for(s in 1:2) {
                cat(formatC(yrel,width=2,format="d"),
                    formatC(s,   width=1,format="d"),
                    char_count[s, ,tag_type, yrel, ], "\n")

             }
             cat("\n")  # blank line between each group of releases
          }
       }
}


loglikelihood<-function(observed,expected)
#Calculate the loglikilihood for the poisson model
{
   loglik<- -expected
   loglik<-loglik+observed*log(ifelse(expected>0, expected,1))
   sum.loglik<-sum(loglik)
   sum.loglik
} # End loglikelihood


quasi.likelihood<-function(observed,expected,overdispersion)
#Calculate the quasi-likelihood for the poisson model
{ 
   quasil<- -(observed-expected)
   quasil<-quasil+observed*log(ifelse(observed>0, observed,1)/ifelse(expected>0,expected,1))
   quasil<-2/overdispersion * quasil
   sum.quasi<-sum(quasil)
   sum.quasi
} # End quasi.likelihood


negative_loglikelihood<-function(beta.pack,ntype, nsample, nstrata, R_observed,N, design, nbeta)
# Calculate the negative log-likelihood function
{
# Unpack the beta's to make beta.p, beta.surv, beta.migrate, beta.tag, beta.report
#cat("betapack",beta.pack)
   Beta<-list()
   Beta$capture<-beta.pack[1:nbeta$capture]
   Beta$surv<-beta.pack[nbeta$capture+1:nbeta$surv]
   Beta$tag<-beta.pack[nbeta$capture+nbeta$surv+1:nbeta$tag]
   Beta$migrate<-beta.pack[nbeta$capture+nbeta$surv+nbeta$tag+1:nbeta$migrate]
   Beta$report<-beta.pack[nbeta$capture+nbeta$surv+nbeta$tag+nbeta$migrate+1:nbeta$report]
print(Beta)

# Transform betas to standard parameters.  This gives me p*, surv*, tag*, migrate*  
# where the * indicates that there may not be as many standard parameters as for the
# saturated model.  For model 1, there will be only 1 standard parameter for each component.
   new<-list()
   new$capture<-transform_beta_to_std(Beta$capture)
   new$surv<-transform_beta_to_std(Beta$surv)
   new$migrate<-transform_beta_to_std(Beta$migrate)
   new$tag<-transform_beta_to_std(Beta$tag)
   new$report<-transform_beta_to_std(Beta$report)

# Make large matrices for multiplication
   big<-list()
   big$capture<-design$capture%*%new$capture
   big$surv<-design$surv%*%new$surv
   big$migrate<-design$migrate%*%new$migrate
   big$tag<-design$tag%*%new$tag
   big$report<-design$report%*%new$report

# Create matrices for calculating E(R)
   std<-list()
   std$surv<-create_parameter_matrix(nsample-1,nstrata, big$surv)
   std$capture<-create_parameter_matrix_capture(nsample, nstrata, big$capture)
   std$tag<-create_parameter_matrix(nsample-1,nstrata, big$tag)
   std$report <-create_report_matrix(nsample,nstrata,ntype,big$report)
   std$migrate <-create_migration_matrix(nsample,nstrata, big$migrate)  

# Calculated the E(recovery)
   R_expected<-Expected_Recovery(ntype,nsample,nstrata,N,std)
#   print_count(R_expected)

# Calculate the negative log likelihood
# Calculate the loglikelihood
   loglik<-loglikelihood(R_observed, R_expected)
   neg_loglik<- -loglik


# Return the negative loglikelihood
#   cat("neg loglik",neg_loglik)
   neg_loglik

} # End negative_loglikelihood


negative_quasi_likelihood<-function(beta.pack,ntype, nsample, nstrata, R_observed,N, design, nbeta)
# Calculate the negative log-likelihood function
{
# Unpack the beta's to make beta.p, beta.surv, beta.migrate, beta.tag, beta.report
cat("betapack",beta.pack)
   Beta<-list()
   Beta$capture<-beta.pack[1:nbeta$capture]
   Beta$surv<-beta.pack[nbeta$capture+1:nbeta$surv]
   Beta$tag<-beta.pack[nbeta$capture+nbeta$surv+1:nbeta$tag]
   Beta$migrate<-beta.pack[nbeta$capture+nbeta$surv+nbeta$tag+1:nbeta$migrate]
   Beta$report<-beta.pack[nbeta$capture+nbeta$surv+nbeta$tag+nbeta$migrate+1:nbeta$report]
   overdispersion<-beta.pack[length(beta.pack)]

# Transform betas to standard parameters.  This gives me p*, surv*, tag*, migrate*  
# where the * indicates that there may not be as many standard parameters as for the
# saturated model.  For model 1, there will be only 1 standard parameter for each component.
   new<-list()
   new$capture<-transform_beta_to_std(Beta$capture)
   new$surv<-transform_beta_to_std(Beta$surv)
   new$migrate<-transform_beta_to_std(Beta$migrate)
   new$tag<-transform_beta_to_std(Beta$tag)
   new$report<-transform_beta_to_std(Beta$report)


# Make large matrices for multiplication
   big<-list()
   big$capture<-design$capture%*%new$capture
   big$surv<-design$surv%*%new$surv
   big$migrate<-design$migrate%*%new$migrate
   big$tag<-design$tag%*%new$tag
   big$report<-design$report%*%new$report

# Create matrices for calculating E(R)
   std<-list()
   std$surv<-create_parameter_matrix(nsample-1,nstrata, big$surv)
   std$capture<-create_parameter_matrix_capture(nsample, nstrata, big$capture)
   std$tag<-create_parameter_matrix(nsample-1,nstrata, big$tag)
   std$report <-create_report_matrix(nsample,nstrata,ntype,big$report)
   std$migrate <-create_migration_matrix(nsample,nstrata, big$migrate)  

# Calculated the E(recovery)
   R_expected<-Expected_Recovery(ntype,nsample,nstrata,N,std)
#   print_count(R_expected)

# Calculate the negative log likelihood
# Calculate the loglikelihood
   quasilik<-quasi.likelihood(R_observed, R_expected,overdispersion)
   neg_loglik<- -quasilik


# Return the negative loglikelihood
#   cat("neg loglik",neg_loglik)
   neg_loglik

} # End negative_loglikelihood




Expected_Recovery<-function(ntype,nsample,nstrata,N,std){
# Multiply matrices in order to find E(R).  
   # First I will need an identity matrix
   I<-diag(nstrata)
   R_expected<-array(rep(0, times=ntype*nstrata*nsample*nstrata*nsample),c(nstrata,nstrata,ntype,nsample,nsample))

   for (k in 1:ntype){
      for (s in 1:(nsample-1)){
         for (t in (s+1):nsample){
            R_expected[,,k,s,t]<-diag(N[,k,s])

         # 1. Type 1- Single tag or high reward single tag
            if(k==1 | k==4) {
               i<-s
               while(i<=(t-2) &(t-2)>0){
                  R_expected[,,k,s,t]<-R_expected[,,k,s,t]%*% std$tag[,,i] %*% std$surv[,,i] %*% std$migrate[,,i]%*%(I-std$capture[,,i+1])
                  i<-i+1
                  }
               # 100% reporting rate for type 4    
               if(k==1){R_expected[,,k,s,t]<- R_expected[,,k,s,t] %*% std$tag[,,t-1]%*% std$surv[,,t-1] %*% std$migrate[,,t-1] %*% std$capture[,,t]%*%std$report[,,k,t]}
               else {R_expected[,,k,s,t]<- R_expected[,,k,s,t] %*% std$tag[,,t-1]%*% std$surv[,,t-1] %*% std$migrate[,,t-1] %*% std$capture[,,t]}
#              print(c(k,s,t))
#              print(R_expected[,,k,s,t])
            }
         # 2. Double tagged fish
            if(k==3){
               i<-s
               while(i<=(t-2) &(t-2)>0){
                  R_expected[,,k,s,t]<-R_expected[,,k,s,t]%*% std$tag[,,i]%*%std$tag[,,i] %*% std$surv[,,i] %*% std$migrate[,,i]%*%(I-std$capture[,,i+1])
                  i<-i+1
                  }
               R_expected[,,k,s,t]<- R_expected[,,k,s,t] %*% std$tag[,,t-1]%*% std$tag[,,t-1] %*% std$surv[,,t-1] %*% std$migrate[,,t-1] %*% std$capture[,,t]%*%std$report[,,k,t]
 #             print(c(k,s,t))
 #             print(R_expected[,,k,s,t])
            }
         # 3. Double tagged fish recaptured as single tag
            if(k==2){
               Lam_temp<-array(rep(1,times=nstrata*nstrata), dim=c(nstrata,nstrata))
                             
               # Tag retention
               for(i in s:(t-1)){ #Possible times of tag last are s:(t-1)
                  temp<-array(rep(1,times=nstrata))
                  temp<-diag(temp)                
                  time_tag_lost<-i
                  j<-s
                  #Before time of tag loss
                  while (j<time_tag_lost){
                     temp<-temp%*%std$tag[,,j]%*%std$tag[,,j]%*%std$surv[,,j] %*% std$migrate[,,j]%*%(I-std$capture[,,j+1])
                     j<-j+1
                  }
                  
                  done.last<-FALSE #Whether or not multiplied by last time period.
                  #At time of tag loss                  
                  if (time_tag_lost<(t-1)){
                     temp<-temp%*%std$tag[,,time_tag_lost]%*%(I-std$tag[,,time_tag_lost])%*%std$surv[,,time_tag_lost] %*% std$migrate[,,time_tag_lost]%*%(I-std$capture[,,time_tag_lost+1])
                     } else {
                     temp<-temp%*%std$tag[,,time_tag_lost]%*%(I-std$tag[,,time_tag_lost])%*%std$surv[,,time_tag_lost] %*% std$migrate[,,time_tag_lost]%*%std$capture[,,time_tag_lost+1]%*%std$report[,,k,time_tag_lost+1]
                     done.last<-TRUE}
                     
                  j<-time_tag_lost+1
                  #After time of tag loss
                  while(j<(t-1)){
                     temp<-temp%*%std$tag[,,j]%*%std$surv[,,j] %*% std$migrate[,,j]%*%(I-std$capture[,,j+1])
                     j<-j+1
                  }
                  #At last time
                  if(done.last==FALSE & j==(t-1)){
                     temp<-temp%*%std$tag[,,j]%*%std$surv[,,j] %*% std$migrate[,,j]%*%std$capture[,,j+1]%*%std$report[,,k,j+1]
                     } 
                  if(i==s){Lam_temp<-temp} else {Lam_temp<-Lam_temp+temp}                  
               }
               R_expected[,,k,s,t]<-diag(N[,3,s]) #Start out released as double tag, there are no type 2 releases
               R_expected[,,k,s,t]<-R_expected[,,k,s,t]%*%Lam_temp
#               print(c(k,s,t))
#               print(R_expected[,,k,s,t])   
            }   
         }
      }
   }
   R_expected
#   print(R_expected)

} #End Expected_Recovery

create_parameter_matrix_capture<-function(nsample,nstrata,big)
# Creates the 2x2 standard capture matrix where n= the number of components in the 
# matrix (eg 2x2x10), m is a column vector containing the dimensions,v is the number 
# of periods for which the standard parameter is defined (ie 9 for survival, 10 for 
# capture), w is the number of strata, big is the matrix that holds the parameter
# components (see big.p, big.surv, etc).  Note that this makes capture start at 
#time 2 and end at time 10 which is inline with how the matrices are formally written 
#(ie there is no P1).
{
   index=1
   x<-array(rep(0,times=nsample*nstrata*nstrata),dim=c(nstrata,nstrata,nsample))
   for(i in 2:nsample) {
      for (j in 1:nstrata){
         x[j,j,i]<-big[index] 
         index=index+1
      }
   }
   x      
} # End create_parameter_matrix_capture

create_parameter_matrix<-function(nsample,nstrata,big)
# Creates the 2x2 standard parameter matrix where n= the number of components in the 
# matrix (eg 2x2x10), m is a column vector containing the dimensions,v is the number 
# of periods for which the standard parameter is defined (ie 9 for survival, 10 for 
# capture), w is the number of strata, big is the matrix that holds the parameter
# components (see big.p, big.surv, etc).
{
   index=1
   x<-array(rep(0,times=nsample*nstrata*nstrata),dim=c(nstrata,nstrata,nsample))
   for(i in 1:nsample) {
      for (j in 1:nstrata){
         x[j,j,i]<-big[index] 
         index=index+1
      }
   }
   x      
} # End create_parameter_matrix


create_migration_matrix<-function(nsample,nstrata,big)
# Creates the 2x2 migration matrices.  First I put in the off diagonals, then I 
# put in the diagonals as they are 1-sum(off diagonals). 
{
   index=1
   x<-array(rep(0,times=(nsample-1)*nstrata*nstrata),dim=c(nstrata,nstrata,(nsample-1)))
   for(i in 1:(nsample-1)) {
      for (j in 1:nstrata){
         for (k in 1:nstrata){
            if(j != k ){x[j,k,i]<-big[index] #note order of indices, R is "colujmn major order",  
            index=index+1}                   #first subscript moves fastest, last subscript slowest
         }
      }
   x
   }      

#Now add in the diagonal of the migration matrices, note that rows sum to 1.
   for (i in 1:(nsample-1)) {
      for (j in 1:nstrata){
         x[j,j,i]<-1-sum(x[j,,i])
      }
   }
   x
} #End create_migration_matrix

create_report_matrix<-function(nsample,nstrata,ntype,big)
# Same as create_parameter_matrix_report but it relies on ntype as well, one extra dimension.
# Creates the 2x2 reporting rate matrix 
{
   index=1
   x<-array(rep(0,times=nsample*nstrata*ntype),dim=c(nstrata,nstrata,ntype,nsample))
   for (l in 1:ntype){
      for(i in 2:nsample) {
         for (j in 1:nstrata){
           x[j,j,l,i]<-big[index] 
           index=index+1
         }  
      }  
   }
   x
} # End create_report_matrix        


#initial_beta<-function(n, x, design)
# Find the initial values of the betas given n=number of standard parameters, 
# x=standard parameters, design=design matrix.  This is done using a least squares
# fit of a logistic link ie. logistic(phi)=XB
# Note that I need a for loop year to transform each element one at a time because
# I evaluate whether or not I am at an upper or lower limit in "transform_std_to_beta".
#{
#   y=0
#   for (i in 1:n){
#      y[i]<-transform_std_to_beta(x[i])
#   }
#   ans<-lsfit(design, y,intercept=FALSE)
#   beta<-ans$coefficients
#   beta
#} # End initial_beta

report.initial<-function(ntype,nsample,nstrata,R_observed,N)
# Find the initial reporting rates.  This is the ratio of R(k)/N(k) * N(4)/R(4) or the type k tags 
# reported and released to the type high reward tags reported and released.
{
   report<-c(rep(0,times=ntype*(nsample-1)*nstrata))
   index<-1
   for (k in 1:ntype){
      for (i in 2:nsample){
         for (s in 1:nstrata){
            if(k==4){report[index]<-1} #100% reporting rate for high reward tags
            else {report[index]<-min(((sum(R_observed[,,k,,])+1)/(sum(N[,k,])+1)) * ((sum(N[,4,])+1)/(sum(R_observed[,,4,,])+1)),1)}
            index<-index+1
            
         }
      }
   }
   report
} # End report.initial


   
initial_beta<-function(x, design)
# Find the initial values of the betas given 
# x=standard parameters, design=design matrix.  This is done using a least squares
# fit of a logistic link ie. logistic(phi)=XB
{
   y<-mylogit(x)
   ans<-lsfit(design, y,intercept=FALSE)
   beta<-ans$coefficients
   beta
} #End initial_beta


mylogit <- function( x)
# Looks at a vector and if the values are between 0 and 1 then the logit function is evaluated.  
# The logit() function is found in the "boot" package, this package needs to be loaded before use.  
# Alternatively, replace logit (x) with log(x/(1-x)).
{
  mylogit <- ifelse( x<=0, -10,
     ifelse(x>=1,10, logit(x)))
} # End mylogit


transform_beta_to_std<-function(beta_val)
#This transforms the beta parameters to the standard parameters.
   {
   std_val=exp(beta_val)/(1+exp(beta_val))
   std_val
   } # End transform_beta_to_std



###########--------------------Begin Program--------------------

#Design matrices: There are a number of design matrices to choose from.   
#The list is as follows and last updated February 14, 2005
#*_c constant parameters
#*_s strata variable, constant over time

model<-14
# Create Model 1 (all parameters constant)- Read in Design matrices
if (model==1) {
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_c"),ncol=1,byrow=T)
   report<-matrix(scan("design_report_c"),ncol=1,byrow=T)
   migrate<-matrix(scan("design_migrate_c"),ncol=1,byrow=T)
   surv<-matrix(scan("design_surv_c"),ncol=1,byrow=T)
}
if (model==2){#variable over strata
   tag <- matrix(scan("design_tag_s"),ncol=2,byrow=T)
   capture<-matrix(scan("design_p_s"),ncol=2,byrow=T)
   report<-matrix(scan("design_report_s"),ncol=2,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_s"),ncol=2,byrow=T)
}
if (model==3){#variable over strata and year
   tag <- matrix(scan("design_tag_st"),ncol=10,byrow=T)
   capture<-matrix(scan("design_p_st"),ncol=10,byrow=T)
   report<-matrix(scan("design_report_st"),ncol=10,byrow=T)
   migrate<-matrix(scan("design_migrate_st"),ncol=10,byrow=T)
   surv<-matrix(scan("design_surv_st"),ncol=10,byrow=T)
}
if (model==4){
#reporting dependent on tag type, tag constant, capture variable over strata constant over year, 
#survival variable over strata constant over year, migratin variable over strata and year.
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_s"),ncol=2,byrow=T) 
   surv<-matrix(scan("design_surv_s"),ncol=2,byrow=T)
   migrate<-matrix(scan("design_migrate_st"),ncol=10,byrow=T)
   report<-matrix(scan("design_report_k"), ncol=4,byrow=T)
}
if (model==5){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_s"),ncol=2,byrow=T)
   report<-matrix(scan("design_report_k"),ncol=4,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_s"),ncol=2,byrow=T)
}
if (model==6){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_st"),ncol=10,byrow=T)
   report<-matrix(scan("design_report_c"),ncol=1,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_s"),ncol=2,byrow=T)
}
if (model==7){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_s"),ncol=2,byrow=T)
   report<-matrix(scan("design_report_k"),ncol=4,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_c"),ncol=1,byrow=T)
}
if (model==8){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_s"),ncol=2,byrow=T)
   report<-matrix(scan("design_report_s"),ncol=2,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_st"),ncol=10,byrow=T)
}
if (model==9){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_s"),ncol=2,byrow=T)
   report<-matrix(scan("design_report_c"),ncol=1,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_c"),ncol=1,byrow=T)
}
if (model==10){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_s"),ncol=2,byrow=T)
   report<-matrix(scan("design_report_k12"),ncol=3,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_st"),ncol=10,byrow=T)
}
if (model==11){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_s"),ncol=2,byrow=T)
   report<-matrix(scan("design_report_c"),ncol=1,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_st"),ncol=10,byrow=T)
}
if (model==12){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_c"),ncol=1,byrow=T)
   report<-matrix(scan("design_report_c"),ncol=1,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_st"),ncol=10,byrow=T)
}
if (model==13){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_t"),ncol=5,byrow=T)
   report<-matrix(scan("design_report_c"),ncol=1,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_st"),ncol=10,byrow=T)
}
if (model==14){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_st"),ncol=10,byrow=T)
   report<-matrix(scan("design_report_c"),ncol=1,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_st"),ncol=10,byrow=T)
}
if (model==15){
   tag <- matrix(scan("design_tag_c"),ncol=1,byrow=T)
   capture<-matrix(scan("design_p_t"),ncol=10,byrow=T)
   report<-matrix(scan("design_report_c"),ncol=1,byrow=T)
   migrate<-matrix(scan("design_migrate_s"),ncol=2,byrow=T)
   surv<-matrix(scan("design_surv_c"),ncol=10,byrow=T)
}
      
design<-list()
design$tag<-tag
design$capture<-capture
design$report<-report
design$migrate<-migrate
design$surv<-surv

# Global Variables
#1. Number of standard parameters
nstd<-list()
nstd$tag<-length(design$tag[,1])
nstd$capture<-length(design$capture[,1])
nstd$report<-length(design$report[,1])
nstd$migrate<-length(design$migrate[,1])
nstd$surv<-length(design$surv[,1])


#2. Number of beta parameters
nbeta<-list()
nbeta$capture<-length(design$capture[1,])
nbeta$tag<-length(design$tag[1,])
nbeta$report<-length(design$report[1,])
nbeta$migrate<-length(design$migrate[1,])
nbeta$surv<-length(design$surv[1,])

#3. Other variables
nstrata<-2  #number of strata
nsample<-10 #number of sample times
ntype<-4    #number of tag types

# Read in the number of fish recovered and the number of fish released.  
# Columns for released are: mperiod, mstratum, tagtype, N. 
# Columns for recovered are: mperiod, rperiod, mstratum, rstratum, tagtype, R.
num.recovered<-matrix(scan("NumberRecovered.dat"), ncol=6,byrow=T)
num.released<-matrix(scan("NumberReleased.dat"), ncol=4,byrow=T)

# Create 2*1*nsample*ntype matrices for the number released N[mstratum,tagtype,mperiod]
N<-array(rep(0,times=nsample*ntype*nstrata),dim=c(nstrata,ntype,nsample-1))
for(index in 1:length(num.released[,1]))
   {
   N[num.released[index,2],num.released[index,3],num.released[index,1]]<-num.released[index,4]
   }

# Create 4*4*nsample*ntype matrices for the number recovered R_observed[mstratum,rstratum,tagtype,mperiod,rperiod]
R_observed<-array(rep(0,times=nsample*nsample*ntype*nstrata*nstrata), dim=c(nstrata,nstrata,ntype,nsample,nsample))
for (index in 1:length(num.recovered[,1]))
   {
   R_observed[num.recovered[index,3],num.recovered[index,4], num.recovered[index,5], num.recovered[index,1], num.recovered[index,2]]<-num.recovered[index,6]
   }

# Initialize standard parameters (hard coded in)
std<-list()
std$tag<-matrix(scan("tag_initial3.txt"),ncol=1,byrow=T)
std$capture<-matrix(scan("p_initial3.txt"),ncol=1,byrow=T)
std$migrate<-matrix(scan("migrate_initial3.txt"),ncol=1,byrow=T)
std$surv<-matrix(scan("surv_initial3.txt"),ncol=1,byrow=T)
std$report<-report.initial(ntype,nsample,nstrata,R_observed,N)


# Find initial values for the betas using least squares
# If initial values are 0 or 1, then lower=-10 or upper=10 values are given for
# beta values.  Note that I am using a logistic link function.
beta.init<-list()
beta.init$tag<-initial_beta(std$tag,design$tag)
beta.init$capture<-initial_beta(std$capture,design$capture)
beta.init$report<-initial_beta(std$report,design$report)
beta.init$migrate<-initial_beta(std$migrate,design$migrate)
beta.init$surv<-initial_beta(std$surv,design$surv)

# Pack Betas into one long vector
beta.pack<-c(beta.init$capture, beta.init$surv, beta.init$tag,beta.init$migrate,beta.init$report)

# Read in Capture Histories
#This does not work as my tag histories are 111111111 rather than 1 1 1 1 1 1 1 so I need to be able to distinguish columns here.
#tag.history <- matrix(scan("TagHistory.dat"),what=integer, ncol=22,byrow=T)

# Maximize the likelihood

new.beta<-nlm(negative_loglikelihood, beta.pack, ntype=ntype,nsample=nsample,nstrata=nstrata, R_observed=R_observed,N=N,design=design, nbeta=nbeta, hessian=TRUE, iterlim=1000)

# Make a large design matrix
large_X<-make.large.design(design,nbeta,nstd)
               
# Calculate the Covariance(XB)=X'Cov(Beta)X  where Cov(beta) is the hessian from new.beta               
library(MASS)
inv_hessian<-ginv(new.beta$hessian)
cov_XB<-large_X%*%inv_hessian%*%t(large_X)

# Make long parameter vector for the standard parameters
parm<-transform_to_std(large_X,new.beta$estimate)

# Make dg/dXB vector which is as long as parm vector  for logit this is simply p(1-p)
dg_dXB<-parm*(1-parm)

# Induce the delta method to get Covariance of standard parameters

cov_std<-diag(c(dg_dXB)) %*% cov_XB %*% diag(c(dg_dXB))

stderr<-diag(cov_std)
stderr<-sqrt(stderr)



# Create residual plot
beta_final<-new.beta$estimate
residual<-create_residual(beta_final, nbeta, nsample,ntype,nstrata, R_observed, design,N)
predicted<-residual$residual2
R_expected<-residual$R_expected
index<-1
   for (i in seq(1,nsample,by=2)){
      for (j in (i+1):nsample){
         for (k in 1:ntype){
            for (s in 1:nstrata){
               for (t in 1:nstrata){
                  # double tagged fish were not released in 2004 (period 9) Thus expected counts do not exist.
                  if (k==2 | k==3){
                     if(i!=9){predicted[index]<-R_expected[t,s,k,i,j]                           
                           index<-index+1}
                     else{cat("release 9")}}
                  else{ predicted[index]<-R_expected[t,s,k,i,j]
                        index<-index+1}
                if (predicted[index-1]>5){print(c(t,s,k,i,j,R_expected[t,s,k,i,j]))} 
                                                      
               }
            }
         }
      }
   }

y<-seq(-5,110,1)
plot(predicted,residual$residual2, ylab="Residual", xlab="Expected Recovery", xlim=c(0,101), ylim=c(-4,8),pch=20)
lines(y,rep(0, times=length(y)))

plot(R_observed,R_expected)
plot(residual$overdispersion)

#plotting residuals by factor
strata<-seq(1,nstrata)
period<-seq(2,nsample)
mperiod<-seq(1,nsample,by=2)
type<-seq(1:ntype)

plot(type,residual$type,ylab="Residual", xlab="Tag Type",xlim=c(0,ntype),ylim=c(-4,4),pch=20)
abline(h=0)
plot(strata,residual$mstrata,ylab="Residual", xlab="Release Strata", xlim=c(1,nstrata),ylim=c(-4,4),pch=20)
abline(h=0)
plot(strata,residual$rstrata,ylab="Residual", xlab="Recovery Strata", xlim=c(1,nstrata),ylim=c(-4,4), pch=20)
abline(h=0)
plot(mperiod,residual$mperiod,ylab="Residual", xlab="Release Period", xlim=c(0,nsample),ylim=c(-4,4),pch=20)
abline(h=0)
plot(period,residual$rperiod,ylab="Residual", xlab="Recovery Period", xlim=c(0,nsample),ylim=c(-4,4),pch=20)
abline(h=0)


# Look at over-dispersion calculated as {(o-e)^2/e}/(count-nparm) and the deviance
overd<-over(R_observed, R_expected, 0, 100, length(beta.pack)) 
overd
# adjust standard errors by overdispersion.

stderr<-sqrt(overd$over)*stderr
stderr

# Calculate the goodness of fit statistic
# difference in parameters is the number of recovery counts (Rijstg) - number of parameters -1 for overdispersion parameter.
D<-2*(loglikelihood(R_observed,R_observed)-loglikelihood(R_observed,R_expected))
D
1-pchisq(D,overd$df) #p-value is chi-square on difference in parameters=df



# Calculated the QAIC value
q.aic<-qaic(R_observed, R_expected, overd$over, overd$nparm)
q.aic

q.sat<-qaic(R_observed,R_observed,1,overd$div)
q.sat

# Examples of how to use nlm
# 2parameters x1,x2
#f<-function(x){x[1]^2*(1-x[1])^2*x[2]^2*(1-x[2])^4}
#nlm(f,c(0.2,0.2))

#a=0.33
# 1 parameter x, and nuisance value y
#f<-function(x,y){x^2*(1-x)^2*y^2*(1-y)^4}
#nlm(f,0.5,y=a)