#*******************See file functions*************************************
#**************************************************************************
#***If problems occur with trouser trawl cov functions then numerical******
#***accuracy enhancements used in cccovrich will need to be added**********

#Information matrix and diagnostic calculations for p=0.5 fit.
cov2par_function(x,type="logit",p=0.5) {
  ntotal_nfine+nwide
  cat("\n","cov2par: Fixed split p= ",p)
  fullfithood_sum ( nwide*log(ifelse(nwide>0.001,nwide/ntotal,1)) +
          nfine*log(ifelse(nfine>0.001,nfine/ntotal,1)) )
  cat("\n"," Selectivity curve is",type,".   Likelihood of full model is ",
      format(fullfithood))
  eta_x[1] + x[2]*lenclass
  if(type=="logit") {
    mu_p*exp(eta)/(1-p+exp(eta))
    dmudeta_p*(1-p)*exp(eta)/( (1-p+exp(eta))^2 ) }
  if(type=="cl") {
    dblexpo_exp(-exp(eta))
    mu_1-(1-p)/(1-p*dblexpo)
    dmudeta_(1-p)*p*exp(eta)*dblexpo/ ( (1-p*dblexpo)^2 ) }
  if(type=="logit" | type=="cl") {
    info_matrix(0,nrow=2,ncol=2)
    info[1,1]_sum( ntotal*dmudeta^2 / (mu*(1-mu)) )
    info[1,2]_sum( ntotal*lenclass*dmudeta^2 / (mu*(1-mu)) )
    info[2,1]_info[1,2]
    info[2,2]_sum( ntotal*(lenclass*dmudeta)^2 / (mu*(1-mu)) )
    covar_solve(info) 
    list(covar=covar)  }
  else return(NA)  }

#Information matrix and diagnostic calculations for estimated p fit.
cov3par_function(x,type="logit") {
  ntotal_nfine+nwide
  fullfithood_sum ( nwide*log(ifelse(nwide>0.001,nwide/ntotal,1)) +
          nfine*log(ifelse(nfine>0.001,nfine/ntotal,1)) )
  cat("\n"," Selectivity curve is",type,".   Likelihood of full model is ",
      format(fullfithood))
  ntotal_nfine+nwide
  eta_x[1] + x[2]*lenclass
  if(type=="logit") {
    mu_(x[3]*exp(eta))/(1-x[3]+exp(eta))
    dmudeta_ x[3]*(1-x[3])*exp(eta) / ( (1-x[3]+exp(eta))^2 )
    dmudp_ exp(eta)*(1+exp(eta)) / ( (1-x[3]+exp(eta))^2 ) }
  if(type=="cl") {
    dblexpo_exp(-exp(eta))
    mu_1-(1-x[3])/(1-x[3]*dblexpo)
    dmudeta_(1-x[3])*x[3]*exp(eta)*dblexpo/ ( (1-x[3]*dblexpo)^2 )
    dmudp_(1-dblexpo)/( (1-x[3]*dblexpo)^2 ) }
  if(type=="logit" | type=="cl") {
    info_matrix(0,nrow=3,ncol=3)
    info[1,1]_sum( ntotal*dmudeta^2 / (mu*(1-mu)) )
    info[1,2]_sum( ntotal*lenclass*dmudeta^2 / (mu*(1-mu)) )
    info[1,3]_sum( ntotal*dmudp*dmudeta / (mu*(1-mu)) )
    info[2,1]_info[1,2]
    info[2,2]_sum( ntotal*(lenclass*dmudeta)^2 / (mu*(1-mu)) )
    info[2,3]_sum( ntotal*lenclass*dmudeta*dmudp / (mu*(1-mu)) )
    info[3,1]_info[1,3]
    info[3,2]_info[2,3]
    info[3,3]_sum( ntotal*dmudp^2 / (mu*(1-mu)) )
    covar_solve(info) 
    list(covar=covar)  }
  else return(NA)  }


#Covariance matrix for logistic fit to covered codend data
cccov_function(x) {
  ntotal_nfine+nwide
# fullfithood_sum ( nwide*log(ifelse(nwide>0.001,nwide/ntotal,1)) +
#         nfine*log(ifelse(nfine>0.001,nfine/ntotal,1)) )
# cat("\n"," Selection curve is Richard's. Likelihood of full model is ",
#     format(fullfithood))
  eta_x[1] + x[2]*lenclass
  mu_exp(eta)/(1+exp(eta))
  dmudeta_mu/(1+exp(eta))
  info_matrix(0,nrow=2,ncol=2)
  info[1,1]_sum( ntotal*dmudeta^2 / (mu*(1-mu)) )
  info[1,2]_sum( ntotal*lenclass*dmudeta^2 / (mu*(1-mu)) )
  info[2,1]_info[1,2]
  info[2,2]_sum( ntotal*(lenclass*dmudeta)^2 / (mu*(1-mu)) )
  covar_solve(info)
  list(covar=covar) }

#Covariance matrix for Richard's fit to covered codend data
cccovrich_function(x) {
  ntotal_nfine+nwide
# fullfithood_sum ( nwide*log(ifelse(nwide>0.001,nwide/ntotal,1)) +
#         nfine*log(ifelse(nfine>0.001,nfine/ntotal,1)) )
# cat("\n"," Selection curve is Richard's. Likelihood of full model is ",
#     format(fullfithood))
  eta_x[1] + x[2]*lenclass
  nu_exp(eta)/(1+exp(eta)); mu_nu^(1/x[3])
  eta <- ifelse(mu == 1, Inf, eta) # To avoid numerical problems below
  dmudeta_(1/x[3])*mu/(1+exp(eta))
  dmudgamma_-x[3]^(-2)*mu*log(nu)
  info_matrix(0,nrow=3,ncol=3)
  info[1,1]_sum( ntotal*dmudeta^2 / (mu*(1-mu)) ,na.rm=T)
  info[1,2]_sum( ntotal*lenclass*dmudeta^2 / (mu*(1-mu)) ,na.rm=T)
  info[1,3]_sum( ntotal*dmudgamma*dmudeta / (mu*(1-mu)) ,na.rm=T)
  info[2,1]_info[1,2]
  info[2,2]_sum( ntotal*(lenclass*dmudeta)^2 / (mu*(1-mu)) ,na.rm=T)
  info[2,3]_sum( ntotal*lenclass*dmudeta*dmudgamma / (mu*(1-mu)) ,na.rm=T)
  info[3,1]_info[1,3]
  info[3,2]_info[2,3]
  info[3,3]_sum( ntotal*dmudgamma^2 / (mu*(1-mu)) ,na.rm=T)
  covar_solve(info)
  list(covar=covar) }

#Covariance matrix for Richard's fit to trouser trawl data
covrich_function(x,npars=3,p=0.5) {
  ntotal_nfine+nwide
  eta_x[1] + x[2]*lenclass 
  if(npars==2) cat("\n","covrich: Fixed split p= ",p) else p_x[4]
  nu_exp(eta)/(1+exp(eta)); r_nu^(1/x[3]); mu_p*r/(1-p+p*r)
  dmudr_p*(1-p)/((1-p+p*r)^2)
  drdeta_(1/x[3])*r/(1+exp(eta)); dmudeta_dmudr*drdeta
  drdgamma_-x[3]^(-2)*r*log(nu);  dmudgamma_dmudr*drdgamma
  dmudp_r/((1-p+p*r)^2)
  info_matrix(0,nrow=4,ncol=4)
  info[1,1]_sum( ntotal*dmudeta^2 / (mu*(1-mu)) )
  info[1,2]_sum( ntotal*lenclass*dmudeta^2 / (mu*(1-mu)) )
  info[1,3]_sum( ntotal*dmudgamma*dmudeta / (mu*(1-mu)) )
  info[1,4]_sum( ntotal*dmudeta*dmudp / (mu*(1-mu)) )
  info[2,1]_info[1,2]
  info[2,2]_sum( ntotal*(lenclass*dmudeta)^2 / (mu*(1-mu)) )
  info[2,3]_sum( ntotal*lenclass*dmudeta*dmudgamma / (mu*(1-mu)) )
  info[2,4]_sum( ntotal*lenclass*dmudeta*dmudp / (mu*(1-mu)) )
  info[3,1]_info[1,3]
  info[3,2]_info[2,3]
  info[3,3]_sum( ntotal*dmudgamma^2 / (mu*(1-mu)) )
  info[3,4]_sum( ntotal*dmudgamma*dmudp / (mu*(1-mu)) )
  info[4,1]_info[1,4]
  info[4,2]_info[2,4]
  info[4,3]_info[3,4]
  info[4,4]_sum( ntotal*dmudp*dmudp / (mu*(1-mu)) )
  if(npars==3) covar_solve(info) else covar_solve(info[1:3,1:3])
  list(covar=covar) }