Mercurial > repos > xuebing > sharplabtool
comparison tools/mytools/cdf-old-not-used/cdf2.py @ 0:9071e359b9a3
Uploaded
| author | xuebing |
|---|---|
| date | Fri, 09 Mar 2012 19:37:19 -0500 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:9071e359b9a3 |
|---|---|
| 1 ''' | |
| 2 python wrapper for cdf plot | |
| 3 ''' | |
| 4 | |
| 5 import os,sys | |
| 6 | |
| 7 filename1 = sys.argv[1] | |
| 8 filename2 = sys.argv[2] | |
| 9 output = sys.argv[3] | |
| 10 c1 = sys.argv[4].strip() | |
| 11 c2 = sys.argv[5].strip() | |
| 12 label1,label2 = sys.argv[6].strip().split(',') | |
| 13 title = sys.argv[8].strip() | |
| 14 log = sys.argv[7] | |
| 15 | |
| 16 | |
| 17 | |
| 18 rf = open('tmp.r','w') | |
| 19 | |
| 20 rf.write("options(warn=-1)\n") | |
| 21 | |
| 22 rf.write("x1 <- read.table('"+filename1+"')\n") | |
| 23 rf.write("x2 <- read.table('"+filename2+"')\n") | |
| 24 | |
| 25 | |
| 26 rf.write("mycdf <- function(x1,x2) {\n") | |
| 27 rf.write(" sx1 <- sort(x1)\n") | |
| 28 rf.write(" sy1 <- c(1:length(sx1))/length(sx1)\n") | |
| 29 rf.write(" sx2 <- sort(x2)\n") | |
| 30 rf.write(" sy2 <- c(1:length(sx2))/length(sx2)\n") | |
| 31 rf.write(" mi <- min(c(x1,x2))\n") | |
| 32 rf.write(" ma <- max(c(x1,x2))\n") | |
| 33 rf.write(" sx1 <- c(mi,sx1,ma)\n") | |
| 34 rf.write(" sx2 <- c(mi,sx2,ma)\n") | |
| 35 rf.write(" sy1 <- c(0,sy1,1)\n") | |
| 36 rf.write(" sy2 <- c(0,sy2,1)\n") | |
| 37 rf.write(" pv <- ks.test(x1,x2)$p.value\n") | |
| 38 rf.write(" list(sx1,sy1,sx2,sy2,pv)\n") | |
| 39 rf.write("}\n") | |
| 40 | |
| 41 rf.write("list1 <- x1[,"+c1+"]\n") | |
| 42 rf.write("list2 <- x2[,"+c2+"]\n") | |
| 43 | |
| 44 rf.write("list1[list1 > 1e10] <- 1e10\n") | |
| 45 rf.write("list2[list2 > 1e10] <- 1e10\n") | |
| 46 rf.write("list1[list1 < -1e10] <- -1e10\n") | |
| 47 rf.write("list2[list2 < -1e10] <- -1e10\n") | |
| 48 | |
| 49 rf.write("max(list1)\n") | |
| 50 rf.write("max(list2)\n") | |
| 51 rf.write("min(list1)\n") | |
| 52 rf.write("min(list2)\n") | |
| 53 | |
| 54 | |
| 55 | |
| 56 if log == 'log': | |
| 57 rf.write("list1 <- log2(list1) \n") | |
| 58 rf.write("list2 <- log2(list2) \n") | |
| 59 rf.write("res <- mycdf(list1,list2)\n") | |
| 60 | |
| 61 rf.write("pdf('"+output+"')\n") | |
| 62 rf.write("plot(res[[1]],res[[2]],type='l',lty=1,lwd=2,col='red',main=paste('"+title+"',format(res[[5]], digit=2,scientific = T),sep=' , p='),xlab='value',ylab='cumulative frequency')\n") | |
| 63 rf.write("lines(res[[3]],res[[4]],type='l',lty=2,lwd=2,col='blue')\n") | |
| 64 rf.write("legend('bottomright',c('"+label1+"','"+label2+"'),col=c('red','blue'), lty=1:2, lwd=2, bty='n')\n") | |
| 65 rf.write("dev.off()\n") | |
| 66 | |
| 67 rf.close() | |
| 68 os.system("R --vanilla < tmp.r") | |
| 69 os.system('rm tmp.r') |