annotate python/cluster.py @ 5:672621ab4db4

trim points we check to a circle per radius
author Henry S. Thompson <ht@inf.ed.ac.uk>
date Tue, 25 May 2021 17:19:00 -0400
parents 56508a6033a9
children 3ee329b129c6
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
4
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
1 #!/usr/bin/python3
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
2 '''Read a minutor block lockation tsv and resort it to show clusters'''
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
3 import sys
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
4
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
5 from util import *
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
6
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
7 if len(sys.argv)==1 or sys.argv[1]=='-h':
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
8 print("""Usage: cluster.py [-h] [-c n] infile.tsv [outfile.tsv]
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
9 n is cluster diameter, default is 5
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
10 default outfile is [infile]_c[n].tsv""")
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
11 exit(1)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
12 if sys.argv[1]=='-c':
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
13 sys.argv.pop(1)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
14 n=float(sys.argv.pop(1))
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
15 else:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
16 n=5.0
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
17
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
18 infile_name=sys.argv.pop(1)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
19 if len(sys.argv)>1:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
20 outfile_name=sys.argv.pop(1)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
21 else:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
22 outfile_name="%s_c%s.tsv"%(infile_name.split('.')[0],n)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
23
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
24 cc=[]
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
25
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
26 with open(infile_name,'r') as infile:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
27 with open(outfile_name,'w') as outfile:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
28 l=infile.readline().rstrip()
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
29 print(l,file=outfile)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
30 ff=PPAT.split(l)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
31 (nr,ox,oy,oz)=intsMaybe(ff)
5
672621ab4db4 trim points we check to a circle per radius
Henry S. Thompson <ht@inf.ed.ac.uk>
parents: 4
diff changeset
32 home=(float(ox),float(oy),float(oz))
4
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
33 et=ff[9]
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
34 l=infile.readline().rstrip()
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
35 print(l,file=outfile)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
36 (orad,ymin,ymax)=intsMaybe(PPAT.split(l))
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
37 _=infile.readline()
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
38 for l in infile:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
39 found=False
5
672621ab4db4 trim points we check to a circle per radius
Henry S. Thompson <ht@inf.ed.ac.uk>
parents: 4
diff changeset
40 (_,td,qd)=l.rstrip().split('\t')
672621ab4db4 trim points we check to a circle per radius
Henry S. Thompson <ht@inf.ed.ac.uk>
parents: 4
diff changeset
41 q=[float(i) for i in qd.split(',')]
672621ab4db4 trim points we check to a circle per radius
Henry S. Thompson <ht@inf.ed.ac.uk>
parents: 4
diff changeset
42 td=float(td)
672621ab4db4 trim points we check to a circle per radius
Henry S. Thompson <ht@inf.ed.ac.uk>
parents: 4
diff changeset
43 if td>orad:
672621ab4db4 trim points we check to a circle per radius
Henry S. Thompson <ht@inf.ed.ac.uk>
parents: 4
diff changeset
44 # Enforce a circular region, not square
672621ab4db4 trim points we check to a circle per radius
Henry S. Thompson <ht@inf.ed.ac.uk>
parents: 4
diff changeset
45 break
4
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
46 for c in cc:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
47 for p in c:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
48 if d(p,q)<=n:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
49 c.append(q)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
50 found=True
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
51 break
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
52 if found:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
53 break
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
54 if not found:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
55 cc.append([q])
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
56 oc=cc
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
57 cc=[] # lose
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
58 w=0
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
59 ow=-1
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
60 nc=[] # win
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
61 while True:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
62 for i,c in enumerate(oc):
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
63 win=False
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
64 for p in c:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
65 for g in oc[i+1:]:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
66 for q in g:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
67 if d(p,q)<=n:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
68 win=True
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
69 w+=1
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
70 nc.append(c+g)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
71 break
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
72 if win:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
73 break
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
74 if win:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
75 break
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
76 if not win:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
77 cc.append(c)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
78 print(len(cc),len(nc),ow,w,file=sys.stderr)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
79 if ow==w:
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
80 break
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
81 ow=w
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
82 oc=nc
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
83 nc=[]
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
84 for c in sorted(cc,reverse=True,key=lambda x:len(x)):
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
85 print(c,file=outfile)
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
86
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
87
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
88
56508a6033a9 minutor chunk hacking
Henry S. Thompson <ht@inf.ed.ac.uk>
parents:
diff changeset
89