@@ -228,6 +228,25 @@ def compute_backward_probability_matrix(ref_h, query_h, rho, mu):
228
228
return bm
229
229
230
230
231
+ def _get_ref_hap_seg (a , b ):
232
+ """
233
+ Assuming all biallelic sites, get the index of a reference haplotype segment (i)
234
+ defined by the alleles a and b at two adjacent genotyped markers.
235
+
236
+ #i, a, b
237
+ 0, 0, 0
238
+ 1, 1, 0
239
+ 2, 0, 1
240
+ 3, 1, 1
241
+
242
+ See https://github.com/tskit-dev/tskit/issues/2802#issuecomment-1698767169
243
+
244
+ This is a helper function for `compute_state_probability_matrix`.
245
+ """
246
+ ref_hap_idx = a * 1 + b * 2
247
+ return ref_hap_idx
248
+
249
+
231
250
def compute_state_probability_matrix (fm , bm , ref_h , query_h , rho , mu ):
232
251
"""
233
252
Implement the HMM forward-backward algorithm to compute:
@@ -278,23 +297,6 @@ def compute_state_probability_matrix(fm, bm, ref_h, query_h, rho, mu):
278
297
assert np .all (np .isin (np .unique (ref_h ), [0 , 1 ]))
279
298
assert np .all (np .isin (np .unique (query_h ), [- 1 , 0 , 1 ]))
280
299
sm = np .zeros ((m , h ), dtype = np .float64 ) # HMM state probability matrix
281
-
282
- def _get_ref_hap_seg (a , b ):
283
- """
284
- Assuming all biallelic sites, get the index of a reference haplotype segment (i)
285
- defined by the alleles a and b at two adjacent genotyped markers.
286
-
287
- #i, a, b
288
- 0, 0, 0
289
- 1, 1, 0
290
- 2, 0, 1
291
- 3, 1, 1
292
-
293
- See https://github.com/tskit-dev/tskit/issues/2802#issuecomment-1698767169
294
- """
295
- ref_hap_idx = a * 1 + b * 2
296
- return ref_hap_idx
297
-
298
300
fwd_hap_probs = np .zeros ((m , 4 ), dtype = np .float64 )
299
301
bwd_hap_probs = np .zeros ((m , 4 ), dtype = np .float64 )
300
302
for i in np .arange (m - 1 , 0 , - 1 ):
0 commit comments