# NB the following codes approach the computations in a very direct manner, 
# with all markers brought into a single Bayesian network, and not exploiting 
# the ideas in section 3.4.2. For the UAF2 results, in particular, this leads
# to a network that is very large indeed (so the corresponding line is 
# commented out below).
# But it is not actually necessary to run this analysis in a single net, 
# since the product rule applies in UAF scenario.

# 8 markers simple paternity analysis, for data in Table 3.

markers<-c('D13','D3','D5','D7','FGA','TH01','TPOX','VWA') 

# Evidence: Mgt,Cgt,Pfgt
evidence<-list(
	D13=c('10-13','13-13','11-13'),
	D3=c('16-17','16-17','17-18'),
	D5=c('11-12','11-11','11-11'),
	D7=c('10-12','10-11','11-12'),
	FGA=c('23-23','21-23','21-23'),
	TH01=c('6-6','6-7','7-7'),
	TPOX=c('8-11','8-11','11-11'),
	VWA=c('18-18','18-18','17-18')
)

alleles <-list( 
	D13 = c("10","11","12","13","x"),
	D3 = c("15","16","17","18","x"),
	D5 = c("11","12","x"),
	D7 = c("10","11","12","x"),
	FGA = c("20","21","23","24","25","x"),
	TH01 = c("6","7","9","9.3","x"),
	TPOX = c("8","10","11","x"),
	VWA = c("17","18","20","x") 
	) 
freqsC <-list( 
	D13 = c(0.051,0.339,0.248,0.124,0.238),
	D3 = c(0.262,0.253,0.215,0.152,0.118),
	D5 = c(0.361,0.384,0.255),
	D7 = c(0.243,0.207,0.166,0.384),
	FGA = c(0.127,0.185,0.134,0.136,0.071,0.347),
	TH01 = c(0.232,0.19,0.114,0.368,0.096),
	TPOX = c(0.535,0.056,0.243,0.166),
	VWA = c(0.281,0.2,0.005,0.514)
	)

IBD<-function(patt,m,spg,smg,aspg,asmg) 
{ 
# pattern|relation distribution 
zibd<-matrix(0,7,5) 
zibd[1,1]<-1
zibd[c(2,5),2]<-zibd[c(2,4),3]<-zibd[c(1,3),4]<-zibd[c(1,2),5]<-0.5
#
qr<-c('none','1pa2','2pa1','hsma','hspa')%in%vs.relation
zibd<-zibd[,qr,drop=FALSE] 
qp<-apply(zibd,1,sum)>0
zibd<-zibd[qp,,drop=FALSE] 
# 
pool1<-cs('pool1',m);
pool2<-cs('pool2',m); pool3<-cs('pool3',m) 
# IBD patterns
tab(c(patt,'relation'),c(sum(qp),sum(qr)),zibd,
        c('none','1p2p','1m2m','1p2m','1m2p','1p1m','2p2m')[qp])
# correlated founding genes - suspect and alt suspect 
founder(spg)
founder(pool1) 
founder(pool2) 
founder(pool3)
which(aspg,patt,lw=c(2,1,2,2,2,2,2),spg,pool1)
which(smg,patt,lw=c(3,3,3,3,2,1,3),spg,aspg,pool2)
which(asmg,patt,lw=c(4,4,3,1,4,4,2),spg,aspg,smg,pool3) 
}

UAF<-function(M,m,g,Nc='c',Nd='d',Npool='pool',Ntemp='temp')
{
# M is Dirichlet parameter, m is marker suffix, g is vector
# of founder gene names (excluding suffix)
# NG correlated founding genes, Polya urn
# c's and d's are binary selection variables
# pool's and temp's are allele-valued
# for N = 2,3,...,NG
# cN is choice between a new draw from gene pool (poolN)
# and random choice among preceding g's (tempN1)
# d's create tempN1 from successive biased binary choices,
# with tempNM holding intermediate stages
# tempNj is random choice among first j genes out of N

NG<-length(g)

founder(cs(g[1],m))
for(i in 2:NG) founder(cs(Npool,i,m))

for(i in 2:NG) query(cs(Nc,i,m),c(i-1,M)/(M+i-1))

which(cs(g[2],m),cs(Nc,2,m),cs(g[1],m),cs(Npool,2,m))

for(i in 3:NG) 
{
for(j in 2:(i-1)) query(cs(Nd,i,j,m),c(j-1,1)/j)
which(cs(Ntemp,i,2,m),cs(Nd,i,2,m),cs(g[1],m),cs(g[2],m))
if(i>3) for(j in 3:(i-1)) 
	which(cs(Ntemp,i,j,m),cs(Nd,i,j,m),cs(Ntemp,i,j-1,m),cs(g[j],m))
which(cs(g[i],m),cs(Nc,i,m),cs(Ntemp,i,i-1,m),cs(Npool,i,m))
}

}

run<-function(model,use,pars) 
{ 
clean(T)

if(model=='IBD')
{
	alpha<-pars[1]; gamma<-pars[2]
	tab('relation',5,c(1-alpha-gamma,alpha/2,alpha/2,gamma/2,gamma/2),
        c('none','1pa2','2pa1','hsma','hspa'))
} else if(model=='UAF1' || model=='UAF2')
{
	M<-pars
} else stop('not a valid model')

query('TfeqPf')

for(m in markers[use]) 
{ 
vs('alleles',alleles[[m]])
gene.freq<<-freqsC[[m]]

Tfmg<-cs('Tfmg',m); Tfpg<-cs('Tfpg',m); Afmg<-cs('Afmg',m); Afpg<-cs('Afpg',m)
Pfmg<-cs('Pfmg',m); Pfpg<-cs('Pfpg',m); Pfgt<-cs('Pfgt',m)
Mmg<-cs('Mmg',m); Mpg<-cs('Mpg',m); Mgt<-cs('Mgt',m)
Cmg<-cs('Cmg',m); Cpg<-cs('Cpg',m); Cgt<-cs('Cgt',m)
patt<-cs('patt',m)

if(model=='IBD')
{
	IBD(patt,m,Pfpg,Pfmg,Afpg,Afmg)
	founder(Mpg); founder(Mmg)
} else if(model=='UAF1')
{
	UAF(M,m,c('Pfpg','Pfmg','Afpg','Afmg'))
	founder(Mpg); founder(Mmg)
} else if(model=='UAF2')
{
 	UAF(M,m,c('Pfpg','Pfmg','Afpg','Afmg','Mpg','Mmg'))
}

se(Tfmg,Pfmg,Afmg,'TfeqPf')
se(Tfpg,Pfpg,Afpg,'TfeqPf')
genotype(Pfgt,Pfpg,Pfmg)

genotype(Mgt,Mpg,Mmg)
mendel(Cpg,Tfpg,Tfmg); mendel(Cmg,Mpg,Mmg); genotype(Cgt,Cpg,Cmg)
}


for(m in markers[use]) 
{
Pfgt<-cs('Pfgt',m); Mgt<-cs('Mgt',m); Cgt<-cs('Cgt',m)
prop.evid(Mgt,evidence[[m]][1])
prop.evid(Cgt,evidence[[m]][2])
prop.evid(Pfgt,evidence[[m]][3])
}

x<-marg(tcq[[acliq('TfeqPf')]],'TfeqPf')
cat(x[1]/x[2],':',markers[use],'\n')
invisible(x)

}

########################

# Computations for 'Baseline' column of Table 9 ###################

# Overall analysis using all 8 markers
run('IBD',1:8,c(0,0)) 

# Marker-by-marker analysis
prod<-1 
for(m in 1:8) 
	{x<-run('IBD',m,c(0,0)); prod<-prod*x[1]/x[2]} 

# Print 'product rule' likelihood ratio
prod

# Computations for 'UAF1' column of Table 9 #######################

# Overall analysis using all 8 markers
run('UAF1',1:8,100) 

# Marker-by-marker analysis
prod<-1 
for(m in 1:8) 
	{x<-run('UAF1',m,100); prod<-prod*x[1]/x[2]} 

# Print 'product rule' likelihood ratio
prod

# Computations for 'UAF2' column of Table 9 #######################

# Overall analysis using all 8 markers (see comments at head of file)
# run('UAF2',1:8,100) 

# Marker-by-marker analysis
prod<-1 
for(m in 1:8) 
	{x<-run('UAF2',m,100); prod<-prod*x[1]/x[2]} 

# Print 'product rule' likelihood ratio
prod

# Computations for 'IBD' column of Table 9 ########################

# Overall analysis using all 8 markers
run('IBD',1:8,c(.05,.05)) 

# Marker-by-marker analysis
prod<-1 
for(m in 1:8) 
	{x<-run('IBD',m,c(.05,.05)); prod<-prod*x[1]/x[2]} 

# Print 'product rule' likelihood ratio
prod