R/simapp.R

In DNAseqtest: Generating and Testing DNA Sequences

simapp <-
function(theta,seqLength,merge1){
genseq3 <-
function(t1,t2,seq,Sx2,Sy2,Pix,Piy){
sequ<-matrix(0,ncol=2,nrow=length(seq))
seqx<-seq
seqy<-seq
SS<-diag(c(1,1,1,1))
p1<-SS+Sx2%*%Pix
p2<-SS+Sy2%*%Piy
print(p1)
for(i in 1:(length(seq)*t1)){
km<-sample(1:length(seq),1)
seqx<-replace(seqx,km,(sample(1:4,1,prob=p1[seqx[km],],replace=T)))
}
for(i in 1:(length(seq)*t2)){
km<-sample(1:length(seq),1)
seqy<-replace(seqy,km,(sample(1:4,1,prob=p2[seqy[km],],replace=T)))
}
sequ<-cbind(seqx,seqy)
return(sequ)
}
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))
Sx1<-Smatrix(ss,pix2)
Sy1<-rho*Smatrix(ss,piy2)
merge2<-merge1
k<-dim(merge2)[1]
height2<-c(theta[17:length(theta)],1)
height<-height2
h<-c(0,height)
me<-merge2[k,]
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]
}
}
}
seqr<-matrix(0,seqLength,k-1)
fs<-matrix(0,seqLength,k+1)
seq1<-cbind(sample(1:4,seqLength,prob=f0,replace=T))
for(i in 1:k){
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))
if(i!=1){
seq1<-seqr[,(k-i+1)]
}
seq2<-genseq3(t1,t2,seq1,Sx2,Sy2,Pix,Piy)
if(me[1]<0){
fs[,(abs(me[1]))]<-seq2[,1]
}
else{
seqr[,(abs(me[1]))]<-seq2[,1]
}
if(me[2]<0){
fs[,(abs(me[2]))]<-seq2[,2]
}
else{
seqr[,(abs(me[2]))]<-seq2[,2]
}
if(m[1]>0)
t0<-h[m[1]+1]
else{
return(fs)
}
me<-c(mm,m[-1])
}
}