#####  ####### #       #######  #####  #######
#     # #       #       #       #     #    #
#       #       #       #       #          #
 #####  #####   #       #####   #          #
      # #       #       #       #          #
#     # #       #       #       #     #    #
 #####  ####### ####### #######  #####     #
#********************************************************************
#Copyright Russell Millar, ******************************************
#          Department of Statistics, ********************************
#          University of Auckland, **********************************
#          P.O. Box 92019, Auckland.  08 Jan 1998********************
#
#====================================================================
#THIS SOFTWARE IS FREE TO USE.
#
#THIS SOFTWARE MAY BE REDISTRIBUTED, BUT MUST BE REDISTRIBUTED FREE
#OF CHARGE.
#
#PLEASE ACKNOWLEDGE THE USE OF THIS SOFTWARE WHERE APPROPRIATE.
#
#NO WARRANTIES ARE GIVEN, AND NO RESPONSIBILITY IS ACCEPTED FOR 
#ANY FAULTS OR CONSEQUENCES ARISING FROM THE USE OF THIS SOFTWARE.
#====================================================================
#
#*********This is the file for loading functions*********************
#*********You also must read read in covfunctions***26 June 1996*****
#***********Sample program code in haddock.example**08 Jan  1998*****
#********************************************************************
#********************************************************************
#To use, create the vectors:
#lenclass   midpoint of lengthclass
#nfine      numbers in control codend (trouser trawl) or cover 
#nwide      numbers in exptl codend
#
#For a logistic fit to covered codend    ccfit()
#For a Richards fit to covered codend    ccfit(type="rich")
#For a logistic fit to trouser trawl     fit()
#For a Richards fit to trouser trawl     fit(type="rich")
#To fix split parameter to 0.5, say      psplit_0.5; fit(npars=2)
#
#Plots of fit and deviance resids are produced unless  plots=F  
#is specified.
#********************************************************************


#Logistic selection curve
#lselect_function(x) {
#  eta_x[1]+x[2]*lenclass
#  1-1/(1+exp(ifelse(eta>80,80,eta))) }

lselect_function(x) {
   expo_exp(pmin(500,x[1]+x[2]*lenclass))
   expo/(1+expo) }

cchood_function(x) {
  select_lselect(x)
  -sum(nwide*log(select) + nfine*log(1-select)) }

#ccrichhood_function(x) {
#  select_lselect(x[1:2])^(1/x[3])
#  Tmp_nwide*log(select)+nfine*log(1-select)
#  hood_-sum(Tmp) }

ccrichhood_function(x) {
  select_lselect(x[1:2])^(1/x[3])
  -sum( nwide*ifelse(select>0,log(select),-1e+06) +
      nfine*ifelse(select<1,log(1-select),-1e+06) ) }

hood2par_function(x) {
  expo_exp(x[1]+x[2]*lenclass)
  cselect_psplit*expo/( 1-psplit + expo)
  -sum( nwide*log(cselect) + nfine*log(1-cselect) ) }

hood3par_function(x) {
  expo_exp(x[1]+x[2]*lenclass)
  cselect_x[3]*expo/( (1-x[3]) + expo)
  -sum( nwide*log(cselect) + nfine*log(1-cselect) ) }

richhood_function(x) {
  select_lselect(x[1:2])^(1/x[3])
  cselect_x[4]*select/(1-x[4] + x[4]*select)
  -sum(nwide*log(cselect) + nfine*log(1-cselect)) }

richhood2par_function(x) {
  select_lselect(x[1:2])^(1/x[3])
  cselect_psplit*select/(1-psplit + psplit*select)
  -sum(nwide*log(cselect) + nfine*log(1-cselect)) }


#Returns the Pearson and deviance residuals
devres_function(y,yhat,n,suff.big=3) {
  if( any((n==0)&(y>0 | yhat>0)) ) stop("Wrong data in function devres")
  if( any((yhat==0 & y >0) | (yhat==n & y<n)) )
    stop("Impossibility in function devres")
  p_ifelse(n>0,y/n,0.5); phat_ifelse(n>0,yhat/n,0.5)
  sign_ifelse(y>=yhat,1,-1)
  Pearson_(y-yhat)/ifelse(n*phat*(1-phat)>0,sqrt(n*phat*(1-phat)),1)
  l_ifelse(y>0,y*(log(p)-log(phat)),0) +
     ifelse(y<n,(n-y)*(log(1-p)-log(1-phat)),0)
  cat("\n"," Pearson Chisquare=",round(sum(Pearson^2),4),
      ", Model deviance=",round(2*sum(l),4),"#lens used=",sum(n>0))
  suff.dat_(yhat>1 & (n-yhat)>suff.big)
  cat("\n"," Pearson Chisquare=",round(sum(Pearson[suff.dat]^2),4),
      ", Model deviance=",round(2*sum(l[suff.dat]),4),
      "#lens used=",sum(suff.dat),"\n")
  list(Pearson=Pearson,devres=sign*sqrt(2*l),suff.dat=suff.dat) }

retentionlens_function(x,covar,type="logit",probs=c(0.25,0.5,0.75),sr=F) {
  np_length(probs)
  rlens_rep(NA,np)
  rlencovar_matrix(NA,nrow=np,ncol=np)
  if(type=="logit") {
    rlens_( log(probs/(1-probs)) -x[1] ) / x[2]
    if(!missing(covar)) {
      work_matrix(0,nrow=2,ncol=np)
      for(i in 1:np) work[,i]_c(-1/x[2],-rlens[i]/x[2])
      rlencovar_ t(work) %*% covar %*% work 
      lens_matrix(c(rlens,sqrt(diag(rlencovar))),nrow=np,byrow=F)
      if(sr) {
        srlen_rlens[np]-rlens[1]
        SR_c(srlen,sqrt(rlencovar[np,np]+rlencovar[1,1]-2*rlencovar[1,np]))
        return(list(lens=lens,sr=SR)) } 
      else return(list(lens=lens)) } }
  if(type=="rich") {
    tmp_(probs^x[3])/(1-probs^x[3]); rlens_(log(tmp)-x[1])/x[2]
    if(!missing(covar)) {
     work_matrix(0,nrow=3,ncol=np)
     for(i in 1:np)
      work[,i]_c(-1/x[2],-rlens[i]/x[2],log(probs[i])/(x[2]*(1-probs[i]^x[3])))
     rlencovar_ t(work) %*% covar %*% work
     lens_matrix(c(rlens,sqrt(diag(rlencovar))),nrow=np,byrow=F)
     if(sr) {
       srlen_rlens[np]-rlens[1]
       SR_c(srlen,sqrt(rlencovar[np,np]+rlencovar[1,1]-2*rlencovar[1,np]))
       return(list(lens=lens,sr=SR)) }
     else return(list(lens=lens)) } } }

fit_function(type="logit",npars=3,probs=c(0.25,0.5,0.75),x0=c(-10,0.3,0.5),
                delta=1.0,plots=T,cex=0.8,mkh=0.07,error.bars=F) {
  ntotal_nfine+nwide
  if(npars==2) cat("\n","fit: Fixed split, p= ",psplit)
  fullfithood_sum ( nwide*log(ifelse(nwide>0.001,nwide/ntotal,1)) +
          nfine*log(ifelse(nfine>0.001,nfine/ntotal,1)) )
  cat("\n"," Likelihood of full model is ",format(fullfithood))
  propn_ifelse((nwide+nfine)>0.001,nwide/(nwide+nfine),2)
  if(type=="logit") 
   {
   if(npars==2) {
    Tfit_nlmin(hood2par,x0[1:2],max.iter=100,max.fcal=200)
    Tcov_cov2par(Tfit$x,type="logit",p=psplit) 
    select_lselect(Tfit$x); cselect_psplit*select/(psplit*select + 1-psplit)
    Tlens_retentionlens(Tfit$x,cov=Tcov$covar,probs=probs,sr=T) }
   if(npars==3) {
    Tfit_nlmin(hood3par,x0,max.iter=150,max.fcal=300)
    Tcov_cov3par(Tfit$x,type="logit") 
    select_lselect(Tfit$x) 
    cselect_Tfit$x[3]*select/(Tfit$x[3]*select + 1-Tfit$x[3])
    Tlens_retentionlens(Tfit$x,cov=Tcov$covar[1:2,1:2],probs=probs,sr=T) }
   Tdevres_devres(nwide,ntotal*cselect,ntotal)
   }
  if(type=="rich") 
   {
   if(npars==2) {
    Tfit_nlmin(richhood2par,c(x0[1:2],delta),max.iter=200,max.fcal=300)
    cat("\n"," Likelihood of fitted model is ",
        format(-richhood2par(Tfit$x)),"\n")
    select_lselect(Tfit$x[1:2])^(1/Tfit$x[3])
    cselect_psplit*select/(psplit*select + 1-psplit)
    Tcov_covrich(Tfit$x,npars=2,p=psplit)
    Tlens_retentionlens(Tfit$x,cov=Tcov$covar[1:3,1:3],type="rich",
                        probs=probs,sr=T) }
   if(npars==3) {
    Tfit_nlmin(richhood,c(x0[1:2],1,x0[3]),max.iter=200,max.fcal=300)
    cat("\n"," Likelihood of fitted model is ",format(-richhood(Tfit$x)),"\n")
    select_lselect(Tfit$x[1:2])^(1/Tfit$x[3])
    cselect_Tfit$x[4]*select/(1-Tfit$x[4] + Tfit$x[4]*select) 
    Tcov_covrich(Tfit$x)
    Tlens_retentionlens(Tfit$x,cov=Tcov$covar[1:3,1:3],type="rich",
                        probs=probs,sr=T) }
   Tdevres_devres(nwide,ntotal*cselect,ntotal)
   }
  if(plots) {
   xyticks_c(length(lenclass)-1,10,7)
   if(!error.bars) 
    plot(lenclass[propn!=2], propn[propn!=2], type = "b", pch = 5, mkh=mkh,
           lab = xyticks, xlab = "", ylab = "", xlim = c(lenclass[1],
           lenclass[length(lenclass)]), ylim = c(0,1),cex=cex)
   else {
    lower.bnds_pmax(propn[propn!=2] + qnorm(0.1/2)*0.5/sqrt(ntotal)[propn!=2],0)
    upper.bnds_pmin(propn[propn!=2] - qnorm(0.1/2)*0.5/sqrt(ntotal)[propn!=2],1)
    error.bar(lenclass[propn!=2],propn[propn!=2],lower.bnds,upper.bnds,incr=F,
          pch = 5,lab = xyticks,
          xlab = "Length (cm)",ylab = "Propn retained in codend",mkh=0.07,
          xlim = c(lenclass[1], lenclass[length(lenclass)]), ylim = c(0,1)) }
   title(xlab = "Length (cm)", ylab = "Propn in large mesh codend",
                 main="Proportion of catch in large mesh codend",cex=cex)
   lines(lenclass,cselect,type="l",lty=2)
   plot(lenclass,Tdevres$devres,type="h",lab=xyticks,xlab="",ylab="",cex=cex)
   abline(h=0)
   title(xlab="Length (cm)",ylab="Deviance residual",
              main="Deviance residuals",cex=cex) }
  list(x=Tfit$x,converged=Tfit$converged,covar=Tcov$covar,lens=Tlens$lens,
       sr=Tlens$sr,devres=Tdevres$devres,suff.dat=Tdevres$suff.dat)  }


ccfit_function(type="logit",probs=c(0.25,0.5,0.75),x0=c(-10,0.3),plots=T,
                                                   error.bars=F) {
  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"," Likelihood of full model is ",format(fullfithood))
  if(type=="logit")
   {
    Tfit_nlmin(cchood,x0,max.iter=100,max.fcal=200)
    cat("\n"," Likelihood of logistic model is ",format(-cchood(Tfit$x)))
    Tcov_cccov(Tfit$x)
    Tlens_retentionlens(Tfit$x,Tcov$covar,sr=T)
    select_lselect(Tfit$x)
    Tdevres_devres(nwide,ntotal*select,ntotal)
   }
  if(type=="rich")
   {
    Tfit_nlmin(ccrichhood,c(x0,1),max.iter=200,max.fcal=300)
    cat("\n"," Likelihood of Richard's model is ",format(-ccrichhood(Tfit$x)))
    Tcov_cccovrich(Tfit$x)
    Tlens_retentionlens(Tfit$x,Tcov$covar,type="rich",sr=T)
    select_lselect(Tfit$x[1:2])^(1/Tfit$x[3])
    Tdevres_devres(nwide,ntotal*select,ntotal) 
   }
  if(plots) {
   propn_ifelse(ntotal>0.001,nwide/ntotal,2) 
   xyticks_c(length(lenclass)-1,10,7)
   if(!error.bars) {
    plot(lenclass[propn!=2], propn[propn!=2],type = "b",pch = 5,lab = xyticks, 
          xlab = "Length (cm)", ylab = "Propn retained in codend",mkh=0.07, 
          xlim = c(lenclass[1], lenclass[length(lenclass)]), ylim = c(0,1)) }
   else {
    lower.bnds_pmax(propn[propn!=2] + qnorm(0.1/2)*0.5/sqrt(ntotal)[propn!=2],0)
    upper.bnds_pmin(propn[propn!=2] - qnorm(0.1/2)*0.5/sqrt(ntotal)[propn!=2],1)
    error.bar(lenclass[propn!=2],propn[propn!=2],lower.bnds,upper.bnds,incr=F,
          pch = 5,lab = xyticks,
          xlab = "Length (cm)",ylab = "Propn retained in codend",mkh=0.07,
          xlim = c(lenclass[1], lenclass[length(lenclass)]), ylim = c(0,1)) }
    title(main="Proportion in codend")
    lines(lenclass,select,type="l",lty=2) 
    abline(h=c(0.25,0.5,0.75),lty=3)
    plot(lenclass,Tdevres$devres,type="h",xlab="Length (cm)",lab=xyticks,
                  ylab="Deviance residual")
    title("Residual plot") 
    abline(h=0) }
  list(x=Tfit$x,covar=Tcov$covar,converged=Tfit$converged,lens=Tlens$lens,
       sr=Tlens$sr,devres=Tdevres$devres,suff.dat=Tdevres$suff.dat) }