1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
#!/usr/bin/R CMD BATCH
Rreverse = function(x) {
tmp =c(1:length(x))
for(i in 1:length(x)) {
tmp[i] =x[length(x)-i+1]
}
tmp
}
argv = function(x) {
args = commandArgs()
offset = 0
for(i in 1:length(args)) {
if (args[i] == "--args")
offset=i;
}
args[offset+x]
}
#par(mfrow=c(2,3))
ncolors =5
mycolors =heat.colors(ncolors)
# mycolors =topo.colors(ncolors)
#mycolors =terrain.colors(ncolors)
#mycolors =cm.colors(ncolors)
mycolors =Rreverse(mycolors)
levelnames = c("0.10", "0.25", "0.50", "0.75", "0.95", "0.99")
levels = c(0.10 , 0.25 , 0.50 , 0.75 , 0.95 , 0.99)
mycolors =gray((ncolors-1):0 / ncolors)
# minmax =c( c(0.001,0.25), c(0.25,0.50), c(0.50,0.75), c(0.75,0.95), c(0.95,1) )
mins =c(0.05, 0.25, 0.50, 0.75, 0.95)
maxs =c(0.25, 0.50, 0.75, 0.95, 1.00)
mymax = 120
tsize =1.25
tlab =1.25
mymax=0
plotcolor1 = function(x)
{
x = x
x = 1-x
gray(1-x)
}
LOD = function(x)
{
log10(x/(1-x))
}
transformcolor = function(x)
{
x1 = sqrt(x+0.001)-sqrt(0.0005)
x2 = (x1 + x*2)/3
x3 = max(0,min(1,LOD((x+1)/2)/3))
(x2 + 2*x3)/3
}
plotcolor2 = function(x)
{
x = LOD(x)-LOD(0.25)/(LOD(100)-LOD(0.25))
}
plotcolor = function(x)
{
# x[x<0.97] = x[x<0.97]/3+0.05
# x = x # -x*x + x*x*x #- x*x*x*x + x*x*x*x*x
x = 1-x
# x = x*x*x
gray(x)
}
plotwidth = function(x)
{
x = log10(x/(1-x))
x[x < 0] = 0
x[x > 2] = 2
x = x/2 + 2
x;
}
plotsegs = function(normal,xa0, ya0, xa1, ya1, weight, lab1, lab2)
{
xmax=max(xa1)
ymax=max(ya1)
weight = weight / normal
plot(xa0,ya0,type="n",ylab=lab1,xlab=lab2,font.lab=3,cex.lab=tlab,cex.axis=tsize)
tmp = data.frame(xa0,ya0,xa1,ya1,weight)
tmp2 = tmp
tmp2 = tmp[tmp$weight>0.001, ]
o = sort.list(tmp2$weight);
tmp3 = tmp2[o,]
tmp2 = tmp3
for(i in 1:length(tmp2$weight)) {
tmp2$weight[i] = transformcolor(tmp2$weight[i])
}
colors = gray(1-tmp2$weight)
widths = plotwidth(tmp2$weight)
segments(tmp2$xa0,tmp2$ya0,tmp2$xa1,tmp2$ya1,col=colors,lwd=widths)
legend(xmax*0.82,ymax*0.27,legend=levelnames,
lwd=plotwidth(levels),
col=plotcolor(levels),
cex=tsize,bty="n",
y.intersp=1)
}
filename = argv(1)
filename2 = paste(filename,".pdf",sep="")
pdf(file=filename2,title=filename,height=8,width=8)
par(mar=c(4,4,0.05,0.05))
input = file(filename);
open(input,"r");
total = scan(input,n=1,sep=" ")
name1 = scan(input,what="",n=1,sep=" ")
name1 = scan(input,what="",n=1,sep=" ")
name2 = scan(input,what="",n=1,sep=" ")
close(input)
d = read.table(filename,skip=1)
plotsegs(total,d$V1,d$V2,d$V3,d$V4,d$V5,name1,name2);
# d =read.table("m02")
# plotsegs(d$V1,d$V2,d$V3,d$V4,d$V5,tnames[0+1],tnames[2+1])
# d =read.table("m03")
# plotsegs(d$V1,d$V2,d$V3,d$V4,d$V5,tnames[0+1],tnames[3+1])
# d =read.table("m12")
# plotsegs(d$V1,d$V2,d$V3,d$V4,d$V5,tnames[1+1],tnames[2+1])
# d =read.table("m13")
# plotsegs(d$V1,d$V2,d$V3,d$V4,d$V5,tnames[1+1],tnames[3+1])
# d =read.table("m23")
# plotsegs(d$V1,d$V2,d$V3,d$V4,d$V5,tnames[2+1],tnames[3+1])
#dev.print(device=postscript,file="test.ps",height=8,width=16,horizontal=T)
|
