R/gn.R

In DNAseqtest: Generating and Testing DNA Sequences

gn <-
function(theta,merge2){
pix1<-theta[1:3]
piy1<-theta[4:6]
f00<-theta[7:9]
ss<-theta[10:15]
rho<-theta[16]
pix2<-c(pix1,1-sum(pix1))
Pix2<-diag(pix2)
piy2<-c(piy1,1-sum(piy1))
Piy2<-diag(piy2)
f0<-c(f00,1-sum(f00))
height<-theta[17:length(theta)]
if(sum(abs(pix2))>1){return(0)}
if(sum(abs(piy2))>1){return(0)}
if(sum(abs(f0))>1){return(0)}
if(rho<=0){return(0)}
if(any(theta<=0)){return(0)}
Sx1<-Smatrix(ss,pix2)
Sy1<-rho*Smatrix(ss,piy2)
k<-dim(merge2)[1]
if(k==1){F1<-Fmatrix(theta[17],theta[17],f0,Sx1,Sy1,Pix2,Piy2)
return(F1)
}
else{h<-c(0,height,1)
me<-merge2[k,]
f1<-f0
m<-me
t0<-h[k+1]
cc<-0
if(merge2[,1][k]>0){
  cc<-merge2[,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){
zzz<-0
iii<-k-i+1
if(iii==k){
xx<-0
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<-merge2[m[1],]
if(me[1]>0){
t1<-t0-h[me[1]+1]
}
else{
t1<-t0
}
if(me[2]>0){
t2<-t0-h[me[2]+1]
}
else{t2<-t0
}
if(t1<0||t1>1){
return(0)
}
if(t2<0||t2>1){
return(0)
}
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<-h[m[1]+1]
}
else{
return(f112)
}
me<-c(mm,m[-1])
}
  }
}