Mercurial > hg > lib > markup
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 |
rev | line source |
---|---|
4 | 1 #!/usr/bin/python3 |
2 '''Read a minutor block lockation tsv and resort it to show clusters''' | |
3 import sys | |
4 | |
5 from util import * | |
6 | |
7 if len(sys.argv)==1 or sys.argv[1]=='-h': | |
8 print("""Usage: cluster.py [-h] [-c n] infile.tsv [outfile.tsv] | |
9 n is cluster diameter, default is 5 | |
10 default outfile is [infile]_c[n].tsv""") | |
11 exit(1) | |
12 if sys.argv[1]=='-c': | |
13 sys.argv.pop(1) | |
14 n=float(sys.argv.pop(1)) | |
15 else: | |
16 n=5.0 | |
17 | |
18 infile_name=sys.argv.pop(1) | |
19 if len(sys.argv)>1: | |
20 outfile_name=sys.argv.pop(1) | |
21 else: | |
22 outfile_name="%s_c%s.tsv"%(infile_name.split('.')[0],n) | |
23 | |
24 cc=[] | |
25 | |
26 with open(infile_name,'r') as infile: | |
27 with open(outfile_name,'w') as outfile: | |
28 l=infile.readline().rstrip() | |
29 print(l,file=outfile) | |
30 ff=PPAT.split(l) | |
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 | 33 et=ff[9] |
34 l=infile.readline().rstrip() | |
35 print(l,file=outfile) | |
36 (orad,ymin,ymax)=intsMaybe(PPAT.split(l)) | |
37 _=infile.readline() | |
38 for l in infile: | |
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 | 46 for c in cc: |
47 for p in c: | |
48 if d(p,q)<=n: | |
49 c.append(q) | |
50 found=True | |
51 break | |
52 if found: | |
53 break | |
54 if not found: | |
55 cc.append([q]) | |
56 oc=cc | |
57 cc=[] # lose | |
58 w=0 | |
59 ow=-1 | |
60 nc=[] # win | |
61 while True: | |
62 for i,c in enumerate(oc): | |
63 win=False | |
64 for p in c: | |
65 for g in oc[i+1:]: | |
66 for q in g: | |
67 if d(p,q)<=n: | |
68 win=True | |
69 w+=1 | |
70 nc.append(c+g) | |
71 break | |
72 if win: | |
73 break | |
74 if win: | |
75 break | |
76 if not win: | |
77 cc.append(c) | |
78 print(len(cc),len(nc),ow,w,file=sys.stderr) | |
79 if ow==w: | |
80 break | |
81 ow=w | |
82 oc=nc | |
83 nc=[] | |
84 for c in sorted(cc,reverse=True,key=lambda x:len(x)): | |
85 print(c,file=outfile) | |
86 | |
87 | |
88 | |
89 |