R/gn2.R

In DNAseqtest: Generating and Testing DNA Sequences

gn2 <-
function(theta,merge2)
{
pix1<-theta[1:3]
piy1<-theta[4:6]
f00<-theta[7:9]
ss<-theta[10:15]
rho1<-matrix(theta[16:(15+length(merge2))],,2)
pix2<-c(pix1,1-sum(pix1))
Pix2<-diag(pix2)
piy2<-c(piy1,1-sum(piy1))
Piy2<-diag(piy2)
f0<-c(f00,1-sum(f00))
if(sum(abs(pix2))>1){
return(0)
}
if(sum(abs(piy2))>1){
return(0)
}
if(sum(abs(f0))>1){
return(0)
}
Sx1<-Smatrix(ss,pix2)
Sy1<-Smatrix(ss,piy2)
k<-dim(merge2)[1]
if(k==1){
F1<-Fmatrix(theta[16],theta[17],f0,(rho1[1,1]*Sx1),(rho1[1,2]*Sy1),Pix2,Piy2)
return(F1)
}
else{
merge<-matrix(0,k,2)
for(i in 1:k){
merge[i,]<-rev(sort(merge2[i,]))
}
me<-merge[k,]
f1<-f0
m<-me
cc<-0
if(merge2[,1][k]>0){
cc<-merge[,1][k]
v<-matrix(0,k,1)
v[1]<-merge2[k,1]
cur<-0
max<-1
while(cur<max){
cur<-cur+1
for(i in c(1,2)){
if(merge2[(v[cur]),i]>0){
max<-max+1
v[max]<-merge2[(v[cur]),i]
}
}
cc<-matrix(0,max,1)
for(i in 1:max){
cc[i]<-v[i]
}
}
}
for(i in 1:k){
iii<-k-i+1
if(iii==k){
if(any(merge2[k,1]==cc)){
pix<-pix2
Pix<-Pix2
piy<-piy2
Piy<-Piy2
Sx2<-Sx1
Sy2<-Sy1
}
else{
piy<-pix2
Piy<-Pix2
pix<-piy2
Pix<-Piy2
Sy2<-Sx1
Sx2<-Sy1
}
}
else{
if(any(iii==cc)){
pix<-pix2
Pix<-Pix2
piy<-pix2
Piy<-Pix2
Sx2<-Sx1
Sy2<-Sx1
}
else{
piy<-piy2
Piy<-Piy2
pix<-piy2
Pix<-Piy2
Sx2<-Sy1
Sy2<-Sy1
}
}
m<-rev(sort(me))
mm<-merge[m[1],]
t1<-rho1[iii,][1]
t2<-rho1[iii,][2]
hp<-rev(order(me))
for(j in 1:(4^(i-1))){
if(j==1){
F1<-Fmatrix(t1,t2,f1[1:4],Sx2,Sy2,Pix,Piy)
}
else{
F1<-Fmatrix(t1,t2,f1[(4*(j-1)+1):(4*(j-1)+4)],Sx2,Sy2,Pix,Piy)
}
if(j==1){
f<-c(as.vector(F1),f1[-1:-4])
}
else{
f<-c(f[1:(4^2*(j-1))],as.vector(F1),f1[-1:-(4*j)])
}
}
f12<-array(f,c(rep(4,i+1)))
f112<-aperm(f12,c(hp))
f1<-as.vector(f112)
if(m[1]>0)
t0<-9999
else{
return(f112)
}
me<-c(mm,m[-1])
}
  }
}