@@ -138,26 +138,22 @@ def insert_branch_mutations(ts, mutations_per_branch=1):
138
138
139
139
def remove_mutation_times (ts ):
140
140
tables = ts .tables
141
- tables .mutations .clear ()
142
- for mut in ts .mutations ():
143
- tables .mutations .add_row (
144
- site = mut .site ,
145
- node = mut .node ,
146
- derived_state = mut .derived_state ,
147
- parent = mut .parent ,
148
- metadata = mut .metadata ,
149
- time = None ,
150
- )
141
+ tables .mutations .time = np .full_like (tables .mutations .time , tskit .UNKNOWN_TIME )
151
142
return tables .tree_sequence ()
152
143
153
144
154
- def insert_discrete_time_mutations (ts , num_times = 4 ):
145
+ def insert_discrete_time_mutations (ts , num_times = 4 , num_sites = 10 ):
146
+ """
147
+ Inserts mutations in the tree sequence at regularly-spaced num_sites
148
+ positions, at only a discrete set of times (the same for all trees): at
149
+ num_times times evenly spaced between 0 and the maximum tree height.
150
+ """
155
151
tables = ts .tables
156
152
tables .sites .clear ()
157
153
tables .mutations .clear ()
158
154
height = max ([t .time (t .roots [0 ]) for t in ts .trees ()])
159
- for pos in range ( int ( tables .sequence_length ) ):
160
- anc = "X" * pos
155
+ for j , pos in enumerate ( np . linspace ( 0 , tables .sequence_length , num_sites + 1 )[: - 1 ] ):
156
+ anc = "X" * j
161
157
tables .sites .add_row (position = pos , ancestral_state = anc )
162
158
t = ts .at (pos )
163
159
for k , s in enumerate (np .linspace (0 , height , num_times )):
@@ -166,11 +162,12 @@ def insert_discrete_time_mutations(ts, num_times=4):
166
162
(t .parent (n ) == tskit .NULL ) or (t .time (t .parent (n )) > s )
167
163
):
168
164
tables .mutations .add_row (
169
- site = pos , node = n , derived_state = anc + str (k ), time = s
165
+ site = j , node = n , derived_state = anc + str (k ), time = s
170
166
)
171
167
k += 1
172
168
tables .sort ()
173
- tables .compute_mutation_parents
169
+ tables .build_index ()
170
+ tables .compute_mutation_parents ()
174
171
return tables .tree_sequence ()
175
172
176
173
@@ -833,13 +830,14 @@ def shuffle_tables(
833
830
shuffle_individuals = True ,
834
831
shuffle_sites = True ,
835
832
shuffle_mutations = True ,
833
+ shuffle_migrations = True ,
836
834
keep_mutation_parent_order = False ,
837
835
):
838
836
"""
839
- Randomizes the order the Individual, Population, Edge, Site and Mutation Tables rows.
840
- Note that if mutations are completely shuffled, then TableCollection.sort() will not
841
- necessarily produce valid tables (unless all mutation times are present and
842
- distinct), since it only puts parent mutations before children if
837
+ Randomizes the order of rows in (possibly) all except the Node table. Note
838
+ that if mutations are completely shuffled, then TableCollection.sort() will
839
+ not necessarily produce valid tables (unless all mutation times are present
840
+ and distinct), since it only puts parent mutations before children if
843
841
canonical=True. However, setting keep_mutation_parent_order to True will
844
842
maintain the order of mutations within each site.
845
843
@@ -885,6 +883,20 @@ def shuffle_tables(
885
883
rng .shuffle (randomised_edges )
886
884
for e in randomised_edges :
887
885
tables .edges .add_row (e .left , e .right , e .parent , e .child , metadata = e .metadata )
886
+ # migrations
887
+ randomised_migrations = list (orig .migrations )
888
+ if shuffle_migrations :
889
+ rng .shuffle (randomised_migrations )
890
+ for m in randomised_migrations :
891
+ tables .migrations .add_row (
892
+ m .left ,
893
+ m .right ,
894
+ m .node ,
895
+ pop_id_map [m .source ],
896
+ pop_id_map [m .dest ],
897
+ m .time ,
898
+ m .metadata ,
899
+ )
888
900
# sites
889
901
randomised_sites = list (enumerate (orig .sites ))
890
902
if shuffle_sites :
@@ -934,7 +946,9 @@ def orig_cmp_site(i, j, tables):
934
946
935
947
936
948
def orig_cmp_mutation (i , j , tables , site_order , canonical = False , num_descendants = None ):
937
- ret = site_order [tables .mutations .site [i ]] - site_order [tables .mutations .site [j ]]
949
+ site_i = tables .mutations .site [i ]
950
+ site_j = tables .mutations .site [j ]
951
+ ret = site_order [site_i ] - site_order [site_j ]
938
952
if (
939
953
ret == 0
940
954
and (not tskit .is_unknown_time (tables .mutations .time [i ]))
@@ -1574,40 +1588,39 @@ def assert_table_collections_equal(t1, t2, ignore_provenance=False):
1574
1588
Checks for table collection equality, but step-by-step,
1575
1589
so it's easy to see what's different.
1576
1590
"""
1577
- if ignore_provenance :
1578
- t1 .provenances .clear ()
1579
- t2 .provenances .clear ()
1580
1591
assert_tables_equal (t1 .populations , t2 .populations , "populations" )
1581
1592
assert_tables_equal (t1 .individuals , t2 .individuals , "individuals" )
1582
1593
assert_tables_equal (t1 .nodes , t2 .nodes , "nodes" )
1583
1594
assert_tables_equal (t1 .edges , t2 .edges , "edges" )
1584
1595
assert_tables_equal (t1 .sites , t2 .sites , "sites" )
1585
1596
assert_tables_equal (t1 .mutations , t2 .mutations , "mutations" )
1586
1597
assert_tables_equal (t1 .migrations , t2 .migrations , "migrations" )
1587
- assert_tables_equal (t1 .provenances , t2 .provenances , "provenances" )
1598
+ if not ignore_provenance :
1599
+ assert_tables_equal (t1 .provenances , t2 .provenances , "provenances" )
1588
1600
assert t1 .metadata_schema == t2 .metadata_schema
1589
1601
assert t1 .metadata == t2 .metadata
1590
1602
assert t1 .metadata_bytes == t2 .metadata_bytes
1591
1603
assert t1 .sequence_length == t2 .sequence_length
1592
- assert t1 == t2
1604
+ assert t1 . equals ( t2 , ignore_provenance = ignore_provenance )
1593
1605
1594
1606
1595
1607
def assert_tables_equal (t1 , t2 , label = "" ):
1596
1608
if hasattr (t1 , "metadata_schema" ):
1597
1609
if t1 .metadata_schema != t2 .metadata_schema :
1598
- print (f"{ label } :::::::::: t1 ::::::::::::" )
1599
- print (t1 .metadata_schema )
1600
- print (f"{ label } :::::::::: t2 ::::::::::::" )
1601
- print (t2 .metadata_schema )
1602
- assert t1 .metadata_schema == t2 .metadata_schema
1603
- if t1 .num_rows != t2 .num_rows :
1604
- print (f"{ label } : t1.num_rows { t1 .num_rows } != { t2 .num_rows } t2.num_rows" )
1610
+ msg = (
1611
+ f"{ label } :::::::::: t1 ::::::::::::\n { t1 .metadata_schema } "
1612
+ f"{ label } :::::::::: t2 ::::::::::::\n { t1 .metadata_schema } "
1613
+ )
1614
+ raise AssertionError (msg )
1605
1615
for k , (e1 , e2 ) in enumerate (zip (t1 , t2 )):
1606
1616
if e1 != e2 :
1607
- print (f"{ label } :::::::::: t1 (row { k } ) ::::::::::::" )
1608
- print (e1 )
1609
- print (f"{ label } :::::::::: t2 (row { k } ) ::::::::::::" )
1610
- print (e2 )
1611
- assert e1 == e2
1612
- assert t1 .num_rows == t2 .num_rows
1617
+ msg = (
1618
+ f"{ label } :::::::::: t1 (row { k } ) ::::::::::::\n { e1 } "
1619
+ f"{ label } :::::::::: t2 (row { k } ) ::::::::::::\n { e2 } "
1620
+ )
1621
+ raise AssertionError (msg )
1622
+ if t1 .num_rows != t2 .num_rows :
1623
+ raise AssertionError (
1624
+ f"{ label } : t1.num_rows { t1 .num_rows } != { t2 .num_rows } t2.num_rows"
1625
+ )
1613
1626
assert t1 == t2
0 commit comments