######################################################### ## [1] Function cal_hier_score26 -- calculate the hierarchical scores for L=2,3,4,5,6 cal_hier_score26<-function(data, kmax=10000, ptim=100, anneal.coeff=1e-6, myoutf = "result.txt") { mytf = unique(c(data[,1], data[,2])) conOut = file(myoutf, "w") close(conOut) for(lev in 2:6) { cat("\n") res = matrix(0, length(mytf), lev) row.names(res) = mytf colnames(res) = 1:lev hs = rep(0, ptim) for(p in 1:ptim) { ## Initiation cat("\r\r\r", lev, "->", p) tfs = sample(1:lev, length(mytf), replace = T) names(tfs) = mytf tmp1 = tfs[data[,1]] tmp2 = tfs[data[,2]] xx = sum(tmp1==tmp2) emax = (sum(tmp1>tmp2)+xx)/(sum(tmp11) break } xx = 1:lev old.s = tfs[tmp.tf] xx = xx[xx!=old.s] tfs[tmp.tf] = sample(xx, 1) tmp1 = tfs[data[,1]] tmp2 = tfs[data[,2]] xx = sum(tmp1==tmp2) etmp = (sum(tmp1>tmp2)+xx)/(sum(tmp1emax) { emax = etmp }else { tt = anneal.coeff*k prob = exp((etmp-emax)/tt) if(prob<0.9) { tfs[tmp.tf] = old.s }else { emax = etmp } } } for(i in 1:length(mytf)) { res[i, tfs[i]]= res[i, tfs[i]] + 1 } hs[p] = emax } res = res/ptim max.lev = apply(res, 1, which.max) tmp1 = max.lev[data[,1]] tmp2 = max.lev[data[,2]] xx = sum(tmp1==tmp2) hs1 = (sum(tmp1>tmp2)+xx)/(sum(tmp1tmp2)+xx)/(sum(tmp1tmp2)+xx)/(sum(tmp11) break } xx = 1:lev old.s = tfs[tmp.tf] xx = xx[xx!=old.s] tfs[tmp.tf] = sample(xx, 1) tmp1 = tfs[data[,1]] tmp2 = tfs[data[,2]] xx = sum(tmp1==tmp2) etmp = (sum(tmp1>tmp2)+xx)/(sum(tmp1emax) { emax = etmp }else { tt = anneal.coeff[p]*k prob = exp((etmp-emax)/tt) if(prob<0.9) { tfs[tmp.tf] = old.s }else { emax = etmp } } } hs[p] = emax } res = cbind(anneal.coeff, hs) write.table(res, myoutf, sep="\t") plot(-log10(res[,1]), res[,2], ylim=c(0, max(hs)), xlab="-log10(anneal.coeff)", ylab="Hierarchical score (HS)", main="", type="b") } ######################################################### ## [3] Function Opt_Kmax3 -- Plot HS vs. kmax Opt_kmax3<-function(data, kmax=10000, anneal.coeff=1e-6, lev=3, myoutf = "result.txt") { mytf = unique(c(data[,1], data[,2])) hs = rep(0, kmax) ## Initiation tfs = sample(1:lev, length(mytf), replace = T) names(tfs) = mytf tmp1 = tfs[data[,1]] tmp2 = tfs[data[,2]] xx = sum(tmp1==tmp2) emax = (sum(tmp1>tmp2)+xx)/(sum(tmp11) break } xx = 1:lev old.s = tfs[tmp.tf] xx = xx[xx!=old.s] tfs[tmp.tf] = sample(xx, 1) tmp1 = tfs[data[,1]] tmp2 = tfs[data[,2]] xx = sum(tmp1==tmp2) etmp = (sum(tmp1>tmp2)+xx)/(sum(tmp1emax) { emax = etmp }else { tt = anneal.coeff*k prob = exp((etmp-emax)/tt) if(prob<0.9) { tfs[tmp.tf] = old.s }else { emax = etmp } } hs[k] = emax } res = cbind(1:kmax, hs) write.table(res, myoutf, sep="\t") plot(log10(res[,1]), res[,2], ylim=c(0, max(hs)), xlab="log10(kmax)", ylab="Hierarchical score (HS)", main="", type="b") }