#getcut.fun<-function(pp.org,nn=2*ngeno, plot=0){
# pp <-rev(sort(pp.org))
# ct <-round(pp*nn)
# dimen <-log((1:length(pp)), base=2 )/ct
# info <-cumsum(pp*log(1/pp))
# netinfo<-info-dimen
# cutpos <-netinfo==max(netinfo)
# if(plot==1){
# plot( netinfo, pp,type="b",cex=0.5)
# barplot(pp.org)
# abline(h=pp[cutpos])
# }
# return( (1:length(pp))[cutpos])
# }
final.BigMatB.fun<-function(mmle.cscn, hapreserv, loci, digit){
ini.cscn <-get.ini.possHap.hapReserv.fun(mmle.cscn, loci, digit, hapreserv)
mut.mat.cscn <-ini.cscn$prepsi
subPI.cscn <-get.subPI.fun(mmle.cscn, mut.mat=mut.mat.cscn)
## get BB matrix
jnk.cscn <-get.preBB.skip.fun(subPI.cscn, loci, digit, mmle.cscn)
preBB.cscn <-jnk.cscn$preBB.lst
indexBB.cscn <-jnk.cscn$indexBB
medBB.cscn <-get.medBB.skip.fun(subPI.lst=subPI.cscn, preBB.lst=preBB.cscn, indexBB.lst=indexBB.cscn,loci, digit, mmle.cscn)
fullBB.cscn <-get.fullBB.fun(medBB.lst=medBB.cscn, subPI.lst=subPI.cscn)
BB.cscn <-get.BB.fun(fullBB.cscn)
BigMatB.cscn <-getBigMatB.fun(BB.cscn, subPI.cscn)
return(list(fullBB.cscn=fullBB.cscn,
##BB.cscn=BB.cscn,
BigMatB.cscn =BigMatB.cscn ,
subPI.cscn = subPI.cscn
))
}
##--------------------------------------------------------------------------------------
## modified in "get.fullBB.fun": if(aa+bb==2 && is.null(medBB.lst[[rr]][[mm]])==0)
##
## on 03/10/05: modified in
## get.preBB.skip.fun: use mmle.cscn to construct matB
## get.medBB.skip.fun: use mmle.cscn to construct matB
## getPIstar.fun: remove the "if" statement
## on 03/11/05: add a new function:
## getBigMatB.fun: BB.lst=BB.cs, subPI.lst=subPI.cs
##
## on 04/06/05: add dimensionality check in function
## getBigMatB.fun, by using chkdim.Multiply.fun, and chkdim.Add.fun
## modified part is annotated by ******
##
##--------------------------------------------------------------------------------------
##----
get.step.mut.fun<-function(loci, digit, Ncol, Nrow,NAMEcol, NAMErow){
## Ncol & NAMEcol = col variable
## Nrow & NAMErow = row variable
UU<-matrix(0, ncol=loci, nrow=Ncol)
for(i in 1:loci){ UU[,i]<-as.numeric(substring(NAMEcol, (i-1)*digit+1,i*digit))}
VV<-matrix(0, ncol=loci, nrow=Nrow)
for(i in 1:loci){ VV[,i]<-as.numeric(substring(NAMErow,(i-1)*digit+1,i*digit))}
step.mut<-matrix(0, nrow=Nrow, ncol=Ncol)
for(v in 1:Nrow){
tmp<-rep(0,Ncol); for(i in 1:loci){tmp<-tmp+abs(VV[v,i]-UU[,i])}
step.mut[v,]<-tmp}
row.names(step.mut)<-NAMErow
dimnames(step.mut)[[2]]<-NAMEcol
return(step.mut)}
get.ini.possHap.hapReserv.fun<-function(em, loci, digit, prehapReserv){
possHap <-names(em)
RR <-length(possHap)
tmp <-prehapReserv[ match(possHap, prehapReserv)]
hapReserv<-tmp[is.na(tmp)==0]
Rstar <-length(hapReserv)
## cat("RR=", RR, " Rstar=", Rstar,"\n")
prepsi <-get.step.mut.fun(loci, digit, RR, Rstar, possHap, hapReserv )
return(list(possHap=possHap,
hapReserv=hapReserv,
RR=RR,
Rstar=Rstar,
prepsi=prepsi))
}
#################################################
## get matrix BB
#################################################
get.subPI.fun<-function(em, mut.mat){
subPI.lst<-list()
allhap <-names(em)
if(nrow(mut.mat)==ncol(mut.mat)){ subPI.lst[[1]]<-em }else{
mmin <-apply(mut.mat, 2, min)
tab <-table(mmin)
subPI.lst<-list()
subPI.lst[[1]]<-tmp<-em[mmin==0]
for(rr in 2:(1+max(mmin))){
del <-match(names(tmp), names(em))
subPI.lst[[rr]]<-tmp1<- (em[-del])[apply(matrix(mut.mat[,-del], ncol=(length(em)-length(del)))==(rr-1), 2, sum)>0]
tmp<-c(tmp1, tmp)
}
}
return(subPI.lst)
}
#WCC get.preBB.skip.fun<-function(subPI.lst=subPI.cs, loci, digit, mmle.cscn){
get.preBB.skip.fun<-function(subPI.lst, loci, digit, mmle.cscn){
preBB.lst<-list(); indexBB.lst<-list()
len <-length(subPI.lst)
preBB.lst[[1]] <-diag(rep(1, length(subPI.lst[[1]]))); dimnames(preBB.lst[[1]])<-list(names(subPI.lst[[1]]),names(subPI.lst[[1]]))
indexBB.lst[[1]]<-0
if(len>1)
{
sum.subPI.lst<-sapply(subPI.lst, sum)
for(rr in 2:len){
if(sum.subPI.lst[rr]==0){## no "rr"-step descendants
preBB.lst[[rr]]<-NA
} else
if(sum.subPI.lst[rr]> 0)
{
index.all<-0; ss<-0
while(index.all==0)
{## in order to find the right "ss"-step ancestors to obtain mut.step
ss<-ss+1
if(sum.subPI.lst[rr-ss]>0)
{
mut.step<-as.data.frame(get.step.mut.fun(loci, digit, length(subPI.lst[[rr-ss]]), length(subPI.lst[[rr]]),
names(subPI.lst[[rr-ss]]), names(subPI.lst[[rr]])))
mmin <-apply(mut.step, 1, min)
index.all<-max(mmin==ss)
} else {index.all<-0}
}
names(mmin)<-names(subPI.lst[[rr]])
indexBB.lst[[rr]]<-mmin
theta<-0
## changed on 03/10/2005, should use mmle.cscn instead of pcs or pcn to construct matB
poolPI<-mmle.cscn[match(names(subPI.lst[[rr-ss]]), names(mmle.cscn))]
if(length(poolPI)>1){
tmp.mut.step <- as.matrix(mut.step)
tmp <-(theta^(tmp.mut.step-ss))%*% diag(poolPI)
}else{
tmp <-(theta^(mut.step-ss)) * poolPI
}
preBB.lst[[rr]]<- tmp/ apply(tmp, 1,sum)
dimnames(preBB.lst[[rr]])<-dimnames(mut.step)
}
}## end of rr-loop
}## end of if(len>1)
return(list(preBB.lst=preBB.lst, indexBB.lst=indexBB.lst))
}
## subPI.lst<-subPI.cs[[j]][[ii]];preBB.lst<-preBB.cs[[j]][[ii]];indexBB.lst<-indexBB.cs[[j]][[ii]]
## rm(medBB.lst, subPI.lst, preBB.lst, indexBB.lst, len, rr, ss,tt,tmp, tmp1, tmp0, rowsum.tmp,pos,sstop, chk, prechk)
#WCC get.medBB.skip.fun<-function(subPI.lst=subPI.cn, preBB.lst=preBB.cn, indexBB.lst=indexBB.cn,loci, digit, mmle.cscn){
get.medBB.skip.fun<-function(subPI.lst, preBB.lst, indexBB.lst, loci, digit, mmle.cscn){
## if some hap couldn't find a direct 1-step ancestor, then need to use medBB to construct matB
medBB.lst<-list()
medBB.lst[[1]]<-list()
len <-length(subPI.lst)
medBB.lst[[1]][[1]]<-preBB.lst[[1]]
if(len>1)
{
for(rr in 2:len)
{
medBB.lst[[rr]]<-list()
if(prod(indexBB.lst[[rr]]==1))## i.e., each hap has its 1-step mutant ancestor
{
medBB.lst[[rr]][[1]]<-preBB.lst[[rr]]
} else
if(sum(preBB.lst[[rr]]=='NaN')==0)## i.e., from preBB we've found its closet mutant ancestor
{
medBB.lst[[rr]][[1]]<-preBB.lst[[rr]]
}else{## i.e., some hap dont have 1-step ancestors, need to prepare "BB-1step", "BB-2step", ...etc
##--- BB-1step ---
tmp0<-preBB.lst[[rr]];tmp0[tmp0=='NaN']<-0
medBB.lst[[rr]][[1]]<-tmp0
chk<-1; ss<-(min(indexBB.lst[[rr]])+1)
##for(ss in (min(indexBB.lst[[rr]])+1):(rr-1)){ ##--- BB-"ss"step ---
while(chk>0 & ss<rr)
{
pos<-(indexBB.lst[[rr]]>=ss )
if(sum(pos)==0) ## i.e., all hap have no "ss" step ancetors
{
medBB.lst[[rr]][[ss]]<-NULL
}else{## i.e., some hap have "ss" step ancestors
## need to get the BBmatrix for rr to rr-ss
## note rr=# of mutations needed for current group of hap to become the reserved hap
sstop<-0; tt<-ss
while(sstop==0)
{## in order to find the right "tt"-step ancestors to obtain mut.step
if(sum(subPI.lst[[rr-tt]])==0)
{
tt<-tt+1
}else{
mut.step<-as.data.frame(get.step.mut.fun(loci, digit, Nrow=length(pos),NAMErow=names(subPI.lst[[rr]]),
Ncol=length(subPI.lst[[rr-tt]]),NAMEcol = names(subPI.lst[[rr-tt]])))
mmin <-apply(mut.step, 1, min)
sstop <-prod(mmin[pos]==tt); if(sstop==0){tt<-tt+1}
## here I assume all hap will find their ancestor at the same tt-step, this need to be modified later
}
}## end of while
theta<-0
## changed on 03/10/2005, should use mmle.cscn instead of pcs or pcn to construct matB
poolPI<-mmle.cscn[match(names(subPI.lst[[rr-tt]]), names(mmle.cscn))]
if(length(poolPI)>1)
{
tmp.mut.step <- as.matrix(mut.step)
tmp<-(theta^(tmp.mut.step-tt))%*% diag(poolPI)
}else{
tmp<-(theta^(mut.step-tt)) * poolPI
}
rowsum.tmp<-apply(tmp, 1,sum); rowsum.tmp[rowsum.tmp==0]<-1
tmp1<- tmp/rowsum.tmp
medBB.lst[[rr]][[ss]]<-tmp1*0
medBB.lst[[rr]][[ss]][pos,]<-tmp1[pos,]
dimnames(medBB.lst[[rr]][[ss]])<-dimnames(mut.step)
}
## check if each hap finds its ancestor
if(length(medBB.lst[[rr]])>1){
tmp<-lapply(medBB.lst[[rr]],function(mat){if(1-is.null(mat)){apply(mat, 1, sum)}else{NA}})
prechk<-sapply(tmp, sum)
}else{ prechk<-apply(medBB.lst[[rr]][[1]],1, sum)}
chk<-sum(prechk==0, na.rm=T) ## chk=0 means all category can find its closest ancestor; chk>0 means some havent found yet.
if(chk>0){ss<-ss+1}
}## end of while
}
}## end of rr-loop
}
return(medBB.lst)
}
#WCC get.fullBB.fun<-function(medBB.lst=medBB.cs,subPI.lst=subPI.cs)
get.fullBB.fun<-function(medBB.lst,subPI.lst)
{
name.lst <-lapply(subPI.lst, names)
dim.vec <-sapply(subPI.lst, length)
len <-length(dim.vec)
fullBB.lst<-list()
fullBB.lst[[1]]<-medBB.lst[[1]]
if(len>1)
{
fullBB.lst[[2]]<-medBB.lst[[2]]
if(len>2)
{
for(rr in 3:length(medBB.lst))
{
fullBB.lst[[rr]]<-list()
for(kk in (1:(rr-1)))
{
mat<-matrix(0, nrow=dim.vec[rr], ncol=dim.vec[rr-kk])
row.names(mat) <-name.lst[[rr]]
dimnames(mat)[[2]]<-name.lst[[rr-kk]]
fullBB.lst[[rr]][[kk]]<-mat
}
for(kk in (1:(rr-1)))
{
for(mm in 1:length(medBB.lst[[rr]]))
{
aa<-sum(is.na(match(dimnames(fullBB.lst[[rr]][[kk]])[[2]], dimnames( medBB.lst[[rr]][[mm]])[[2]] )))==0
bb<-sum(is.na(match(dimnames( medBB.lst[[rr]][[mm]])[[2]], dimnames(fullBB.lst[[rr]][[kk]])[[2]] )))==0
if(aa+bb==2 && is.null(medBB.lst[[rr]][[mm]])==0)
## if(aa+bb==2)
{
fullBB.lst[[rr]][[kk]]<-medBB.lst[[rr]][[mm]]
}
}## end of mm-loop
}## end of kk-loop
}## end of rr-loop
}## end of if(len>2)
}## end of if (len>1)
return(fullBB.lst)
}
get.BB.fun<-function(fullBB.lst){
BB.lst <-list()
len <-c(1,(1:100))
cumlen <-cumsum(len)
BB.lst[[1]]<-fullBB.lst[[1]][[1]]
if(length(fullBB.lst)>1)
{
BB.lst[[ 2 ]]<-fullBB.lst[[2]][[1]]
if(length(fullBB.lst)>2)
{
for(rr in 3:length(fullBB.lst))
{
for(kk in 1:length(fullBB.lst[[rr]]))
{
BB.lst[[ cumlen[rr-1]+kk ]]<-fullBB.lst[[rr]][[kk]]
}
}
}
}
return(BB.lst)
}
get.allpossROUTE.fun<-function()
{
### when this file is modified, function "getPIstar.fun" need to be changed as well!!!!!!
m1 <-list(1)
m2 <- lapply(m1, 2, FUN = c)
m3 <-c(lapply(m2, 3, FUN = c),
lapply(m1, 4, FUN = c))
m4 <-c(lapply(m3, 5, FUN = c),
lapply(m2, 6, FUN = c),
lapply(m1, 7, FUN = c))
m5 <-c(lapply(m4, 8, FUN = c),
lapply(m3, 9, FUN = c),
lapply(m2, 10, FUN = c),
lapply(m1, 11, FUN = c))
m6 <-c(lapply(m5, 12, FUN = c),
lapply(m4, 13, FUN = c),
lapply(m3, 14, FUN = c),
lapply(m2, 15, FUN = c),
lapply(m1, 16, FUN = c))
m7 <-c(lapply(m6, 17, FUN = c),
lapply(m5, 18, FUN = c),
lapply(m4, 19, FUN = c),
lapply(m3, 20, FUN = c),
lapply(m2, 21, FUN = c),
lapply(m1, 22, FUN = c))
m8 <-c(lapply(m7, 23, FUN = c),
lapply(m6, 24, FUN = c),
lapply(m5, 25, FUN = c),
lapply(m4, 26, FUN = c),
lapply(m3, 27, FUN = c),
lapply(m2, 28, FUN = c),
lapply(m1, 29, FUN = c))
m9 <-c(lapply(m8, 30, FUN = c),
lapply(m7, 31, FUN = c),
lapply(m6, 32, FUN = c),
lapply(m5, 33, FUN = c),
lapply(m4, 34, FUN = c),
lapply(m3, 35, FUN = c),
lapply(m2, 36, FUN = c),
lapply(m1, 37, FUN = c))
m10<-c(lapply(m9, 38, FUN = c),
lapply(m8, 39, FUN = c),
lapply(m7, 30, FUN = c),
lapply(m6, 41, FUN = c),
lapply(m5, 42, FUN = c),
lapply(m4, 43, FUN = c),
lapply(m3, 44, FUN = c),
lapply(m2, 45, FUN = c),
lapply(m1, 46, FUN = c))
m = c(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10)
d = 10
r = c(1, rep(2, d-1))^(c(0,0:(d-2)))
di = d * (d-1) / 2 + 1
l = rep(1:d, r)
mat = matrix(0, sum(r), di)
for(i in 1:sum(r)) {
mat[i,m[[i]]] = m[[i]]
}
row.names(mat)<-paste("Step ",
c(1,2,rep(3,2), rep(4,4), rep(5,8), rep(6,16), rep(7,32), rep(8,64), rep(9,128), rep(10, 256)),
".",
c(1,1, (1:2), (1:4), (1:8), (1:16), (1:32), (1:64), (1:128), (1:256)), sep="")
dimnames(mat)[[2]]<-paste("mat",c(1,2,rep(3,2), rep(4,3), rep(5,4), rep(6,5),rep(7,6)),".(",(1:di),")",sep="")
return(mat)
}
chkdim.Multiply.fun<-function(mat){
chk1<-length(mat)>0
chk2<-all(dim(mat)>0)
chk3<-sum(is.na(mat))==0
return(list("chk1"=chk1,"chk2"=chk2, "chk3"=chk3))
}
#WCC chkdim.Add.fun<-function(mat, ncolShouldBe=RRstar){
chkdim.Add.fun<-function(mat, ncolShouldBe){
if(length(mat)==0)
{
chk1=chk2=0
}else{
chk1<-all(dim(mat)>0)
chk2<-(ncol(mat) == ncolShouldBe)
}
return(chk1*chk2)
}
#WCC getBigMatB.fun<-function(BB.lst=BB.cs, subPI.lst=subPI.cs){
getBigMatB.fun<-function(BB.lst, subPI.lst){
## modified on Apr 6, 2005; now check for the compatibility of dimensionality before adding or mulitplying
## see ****** part. functions used: chkdim.Multiply.fun, and chkdim.Add.fun
mmat <-get.allpossROUTE.fun()
r.ele <-c( 1,1,2,4,8,16,32,64,128,256)
cumsum.r <-cumsum(r.ele)
lenPI <-length(subPI.lst) # 4
lenBB <-length(BB.lst) # 7=1+1+2+3
lenROUTE <-cumsum.r[lenPI] # 8=1+1+2+4
RRstar <-length(subPI.lst[[1]]) ##*********
## changed on 03/11/05 (remove "if" condition)
mat.route<-as.matrix(mmat[(1:lenROUTE),(1:lenBB)])
pi.ele <-rep((1:10), r.ele)
j.index <-pi.ele[1:cumsum.r[lenPI]] ## indexing a mat B is for the level of H^{(j)}
preBigBB.lst<-list()
for(kk in 1:lenROUTE)## it's fine to have kk start from "1"; when lenROUTE=1, BB=bb[[1]]=Identity Matrix
{
bb<-BB.lst[mat.route[kk,]]
BB<-bb[[length(bb)]]
if(length(bb)>1)
{
##***** check the compability of dimenationality before multiplying
chk<-sapply(bb, chkdim.Multiply.fun)
if(all(chk==T))
{
for(jj in max(1,(length(bb)-1)):1){BB<-BB%*%bb[[jj]]}### range of jj is modified on 03/12/05 as well
}
else
{
BB<-matrix(0, nrow=length(subPI.lst[[ j.index[[kk]] ]]), ncol=RRstar)
colnames(BB) = names(subPI.lst[[1]])
rownames(BB) = names(subPI.lst[[ j.index[[kk]] ]])
}
}
preBigBB.lst[[kk]]<-BB
}
indx<-split((1:length(preBigBB.lst)), j.index)
BigMatB<-NULL
for(jj in 1:length(indx))
{
tmpB<-preBigBB.lst[ indx[[jj]] ]
if(length(tmpB)==1)
{
BigMatB<-rbind(BigMatB, tmpB[[1]])
}else{
##***** check the compability of dimenationality before adding
key<-sapply(tmpB,chkdim.Add.fun, ncolShouldBe=RRstar); pos<-(1:length(tmpB))[key==1]
if(all(pos==0))
{
sumB = NULL
} else {
Add.tmpB<-tmpB[pos]
sumB <-Add.tmpB[[1]]
for(ll in 2:length(Add.tmpB)){sumB<-sumB+Add.tmpB[[ll]]}
}
BigMatB<-rbind(BigMatB, sumB)
}
}
return(BigMatB)
}
## getBigMatB.fun<-function(BB.lst=BB.cs, subPI.lst=subPI.cs)
## {
## mmat <-get.allpossROUTE.fun()
## r.ele <-c( 1,1,2,4,8,16,32,64,128,256)
## cumsum.r <-cumsum(r.ele)
## lenPI <-length(subPI.lst) # 4
## lenBB <-length(BB.lst) # 7=1+1+2+3
## lenROUTE <-cumsum.r[lenPI] # 8=1+1+2+4
##
## ## changed on 03/11/05 (remove "if" condition)
## mat.route<-as.matrix(mmat[(1:lenROUTE),(1:lenBB)])
## pi.ele <-rep((1:10), r.ele)
## j.index <-pi.ele[1:cumsum.r[lenPI]] ## indexing a mat B is for the level of H^{(j)}
##
## preBigBB.lst<-list()
## for(kk in 1:lenROUTE)
## {
## bb<-BB.lst[mat.route[kk,]]
## ## in "getPIstar.fun" there is such if, but it seems to be useless so i removed
## ## on 03/10/05
## ## if(sum(is.na(bb))>0){
## ## ## prePIstar.lst[[kk]]<-rep(0, length(subPI.lst[[1]]))
## ## cat("dd=",dd); print(" is.na(bb)>0")
## ## }else{ below...}
## BB<-bb[[length(bb)]]
## for(jj in max(1,(length(bb)-1)):1){BB<-BB%*%bb[[jj]]}
## preBigBB.lst[[kk]]<-BB
## }
## indx<-split((1:length(preBigBB.lst)), j.index)
## BigMatB<-NULL
## for(jj in 1:length(indx))
## {
## tmpB<-preBigBB.lst[ indx[[jj]] ]
## if(length(tmpB)==1)
## {
## BigMatB<-rbind(BigMatB, tmpB[[1]])
## }else{
## sumB<-0
## for(ll in 1:length(tmpB)){sumB<-sumB+tmpB[[ll]]}
## BigMatB<-rbind(BigMatB, sumB)
## }
## }
## return(BigMatB)
## }
##
##
##
##rm(BB.lst, subPI.lst, prePIstar.lst, mmat, r.ele, cumsum.r, lenPI, lenBB, lenROUTE, mat.route, pi.ele, PI,bb, BB, kk, jj,PIstar)
#WCC getPIstar.fun<-function(BB.lst=BB.cs, subPI.lst=subPI.cs)
getPIstar.fun<-function(BB.lst, subPI.lst)
{
mmat <-get.allpossROUTE.fun()
r.ele <-c( 1,1,2,4,8,16,32,64,128,256)
cumsum.r <-cumsum(r.ele)
lenPI <-length(subPI.lst) # 4
lenBB <-length(BB.lst) # 7=1+1+2+3
lenROUTE <-cumsum.r[lenPI] # 8=1+1+2+4
## changed on 03/11/05 (remove "if" condition)
mat.route<-as.matrix(mmat[(1:lenROUTE),(1:lenBB)])
pi.ele <-rep((1:10), r.ele)
PI <-subPI.lst[pi.ele[1:cumsum.r[lenPI]]]
prePIstar.lst<-list()
for(kk in 1:lenROUTE)
{
bb<-BB.lst[mat.route[kk,]]
if(sum(is.na(bb))>0){
prePIstar.lst[[kk]]<-rep(0, length(subPI.lst[[1]]))
}else{
BB<-bb[[length(bb)]]
for(jj in max(1,(length(bb)-1)):1){BB<-BB%*%bb[[jj]]}
prePIstar.lst[[kk]]<-PI[[kk]] %*% BB
}
}
PIstar<-colSums(matrix(unlist(prePIstar.lst), nrow=length(prePIstar.lst), byrow=T))
names(PIstar)<-names(subPI.lst[[1]])
return(PIstar)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.