-
Notifications
You must be signed in to change notification settings - Fork 0
/
max_quartet_dist.ml
58 lines (54 loc) · 1.83 KB
/
max_quartet_dist.ml
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
(* given a tree T, we want to see if an S maximizing d_Q(T,S) across all S can
* be taken to be a comb tree.
* Thus for every S we ask if it has set or tied the record. If it has, and it
* is a comb tree, then we record it.
*
* the output of this program is a collection of d[T,S]'s and the corresponding
* distances, where S here is the record as described above.
* *)
exception Not_comb
let min_n_tips = 4
and max_n_tips = 8
let () =
if not (!Sys.interactive) then
begin
for n_tips = min_n_tips to max_n_tips do
Format.fprintf Format.std_formatter "@\n%d tips@\n" n_tips;
let all = Array.of_list (Btree.make_all_unrooted n_tips) in
let n = Array.length all
and qsets = Array.map Quartet_btree.qset_of_t all
in
for i=0 to n-1 do
let max_qdist = ref 0
and maximally_distant = ref (-1)
in
for j=0 to n-1 do
if i <> j then begin
let d = Quartet.qset_half_symdiff qsets.(i) qsets.(j) in
if !max_qdist < d || (* store true record always *)
(!max_qdist = d && Btree.is_comb all.(j))
(* store tie if it's a comb *)
then begin
max_qdist := d;
maximally_distant := j;
end
end
done;
Format.fprintf Format.std_formatter "d[%a,@,%a] = %d@\n"
Btree.ppr_btree all.(i)
Btree.ppr_btree all.(!maximally_distant)
!max_qdist;
if not (Btree.is_comb all.(!maximally_distant)) then begin
Format.fprintf Format.std_formatter "@\nfound not comb for %a@\n"
Btree.ppr_btree all.(i);
for j=0 to n-1 do
Format.fprintf Format.std_formatter "%d : %a@\n"
(Quartet.qset_half_symdiff qsets.(i) qsets.(j))
Btree.ppr_btree all.(j)
done;
raise Not_comb
end
done
done
end
let x = Btree_io.of_newick_str