###########################################################################
#The following GOF code is specific to 2 groups and 2 tags
###########################################################################

#Saturated Model: P_hist is a vector 
P_hist = nobs/nobs_total;
P_hist

#Number of unique tag histories in each group -- currently works only for two groups (ngroup=2)
#groupNum= as.vector(as.numeric(strsplit(data.inp$V4,split=";")))+1
#numUnique= vector(length=ngroup)
#for(g in 1:ngroup)
#{	numUnique[g]=0;
#	for(h in 1: length(groupNum))
#	{	if(groupNum[h]==g)
#		{	numUnique[g]=numUnique[g]+1;
#		}
#	}	
#}




######################################################################
# Estimate the probability of double-tagging at the 1st capture time #
######################################################################

#for each group, calculate the number tagged of each type (double or single) at each sample occassion
n_tagged<-array(data=0,dim=c(ngroup,nsample,ntag))
for(g in 1:ngroup)
{ for(i in 1:nobs[g])
  { if(sum(tag.hist[g,i,first[g,i],])==1)#single-tagged
  	{ n_tagged[g,first[g,i],1]=n_tagged[g,first[g,i],1]+1
  	}
  	else if(sum(tag.hist[g,i,first[g,i],])==2)#double-tagged
  	{ n_tagged[g,first[g,i],2]=n_tagged[g,first[g,i],2]+1
  	}
  }
}
n_tagged


# For each group, calculate the probability of being double-tagged at each sample occassion

#probability of double tagging at each sample time (& for each group)
p_double<-array(data=0,dim=c(ngroup,nsample))
for(g in 1:ngroup)
{ for(j in 1:nsample)
  { p_double[g,j]=n_tagged[g,j,2]/sum(n_tagged[g,j,])
  }
}
p_double



################################################################################
#The loglik function produces the log-likelihood for a specific observation i in group g
################################################################################

loglik_func<- function(group,observation_num,first,last,beta,capt.hist,tag.hist, numtag,ngroup,nobs_total,nobs,nsample,ntag,nstd_survival,nstd_capture, nstd_tagloss)
{ 
  std_parms= transform_beta(design_matrix,beta)
  #std_parms=c(0.6,0.24,0.8, 0.627, 1, 0.7, 1, 1, 0.55, 1, 0.7, 0.5, 0.5, 0.9, 0.7, 0.7, 0.16,0.375,1, 0.13, 0.375, 1) 
  std_survival=matrix(std_parms[1:nstd_survival],nrow=ngroup,ncol=(nsample-1),byrow=TRUE);
  std_capture=matrix(std_parms[(nstd_survival+1):(nstd_survival+nstd_capture)],nrow=ngroup,ncol=nsample,byrow=TRUE);
  std_tagloss=array(NA,dim=c(ngroup,(nsample-1),(nsample-1)));
  index=(nstd_survival+nstd_capture+1);
  for(g in 1:ngroup)
  { for(f in 1:(nsample-1))
  	{ for(j in 1:(nsample-1))
  	  { if(j>=f)
  	  	{ std_tagloss[g,f,j]=std_parms[index]
  	      index=(index+1);
  	    }
  	  }
  	}
  }
  std_bstar=matrix(std_parms[(nstd_survival+nstd_capture+nstd_tagloss+1):length(std_parms)],nrow=ngroup,ncol=nsample,byrow=TRUE);
  std_entry=matrix(nrow=ngroup,ncol=nsample);
  for(g in 1:ngroup)
  { for(j in 1:nsample)
  	{ std_entry[g,j]= bstar_to_b(g,j,std_bstar,nsample,ngroup)}
  }
  
  g=group;
  i=observation_num;
  
  lik=1 
  lik= lik*Psi(first[g,i],g,std_bstar,std_survival,std_capture, nsample,ngroup);
  lik= lik*main(g,i,first,last,std_survival,std_capture,std_tagloss, capt.hist,tag.hist,ngroup,nobs,nsample,ntag);
  if(numtag[g,i]==1)
  { lik= lik*chi1(first[g,i],last[g,i],g,std_survival,std_capture,std_tagloss, nsample)
  }
  if(numtag[g,i]==2)
  { lik= lik*chi2(first[g,i],last[g,i],g,std_survival,std_capture,std_tagloss, nsample)
  }
  P_0=P_0_func(std_survival,std_capture,std_tagloss,std_entry,ngroup,nsample);
  ## hat=Nhat_func(nobs,P_0,ngroup);
  loglik=log(lik)
  loglik
}





################################################################################
#Develop the conditional likelihood
################################################################################

loglik_cond=vector(length=ngroup,mode="numeric")
for(g in 1:ngroup)
{ for(i in 1:nobs[g])
  { if(sum(tag.hist[g,i,first[g,i],])==1) #single tag at first history
  	{P_0=P_0_func(std_survival,std_capture,std_tagloss,std_entry,ngroup,nsample)
       a=loglik_func(g,i,first,last,out$estimate,capt.hist,tag.hist,numtag,ngroup,nobs_total,nobs,nsample,ntag,nstd_survival,nstd_capture,nstd_tagloss)
       loglik_cond[g] = loglik_cond[g]+ a - log(max(1-P_0[g], 0.001)) + log(max(1 - p_double[g,first[g,i]],0.001))
  	}
  	
  	else #double tag at first history
  	{ P_0=P_0_func(std_survival,std_capture,std_tagloss,std_entry,ngroup,nsample)
        a=loglik_func(g,i,first,last,out$estimate,capt.hist,tag.hist,numtag,ngroup,nobs_total,nobs,nsample,ntag,nstd_survival,nstd_capture,nstd_tagloss)
        loglik_cond[g]=loglik_cond[g]+ a - log(max(1-P_0[g], 0.001)) +log(max(p_double[g,first[g,i]],0.001))
# print(loglik_cond[g])
	}
#a=loglik_cond
#  print(c(g, a[g]))
  }
}




################################################################################
## Saturated log-likelihood (ignoring factorial terms)
#################################################################################
freq= as.vector(as.numeric(strsplit(data.inp$V2,split=";")))
#groupNum= as.vector(as.numeric(strsplit(data.inp$V4,split=";")))+1

#numUnique= vector(length=ngroup)
#for(g in 1:ngroup)
#{	numUnique[g]=0;
#	for(h in 1: length(groupNum))
#	{	if(groupNum[h]==g)
#		{	numUnique[g]=numUnique[g]+1;
#		}
#	}	
#}


N_unique = vector (length=ngroup)
#for(g in 1:ngroup)
#{ N_unique[g]=numUnique[g+1]
# }
N_unique[ngroup]=length(data.inp$V3)


index=0
loglik_sat<- vector(length=ngroup,mode="numeric")
for(g in 1:ngroup)
{ for(i in 1:N_unique[g])
  {   
  index=index+1
  loglik_sat[g] = loglik_sat[g] + freq[index]*log(freq[index]/nobs[g]) 
  }
}

Tsat=sum(loglik_sat)
Tcum=sum(loglik_cond)

T = sum(loglik_sat) - sum(loglik_cond)
T

df<-sum(N_unique)-1-nbeta-ngroup*nsample
df

pv=1-pchisq(2*T,df)
pv