Complete implementing convert_to_genetic_map_positions

python/tests/beagle_numba.py

@@ -47,27 +47,28 @@ def convert_to_genetic_map_positions(pos, genetic_map=None):
     :return: Genetic map positions (cM).
     :rtype: numpy.ndarray
     """
-    if genetic_map is None:
-        return pos / 1e6
     assert np.all(pos >= genetic_map.base_pos[0]) and np.all(
         pos <= genetic_map.base_pos[-1]
     ), "Physical positions outside of genetic map."
-    # Approximate genetic map distances between adjacent base positions.
+    if genetic_map is None:
+        return pos / 1e6
+    # Approximate genetic map distances by linear interpolation.
+    # Note np.searchsorted(a, v, side='right') returns i s.t. a[i-1] <= v < a[i].
     left_idx = np.searchsorted(genetic_map.base_pos, pos, side="right")
-    xMa = np.zeros(len(pos), dtype=np.float32)
-    bMa = np.zeros(len(pos), dtype=np.float32)
-    fbMfa = np.zeros(len(pos), dtype=np.float32)
+    est_gen_pos = np.zeros(len(pos), dtype=np.float32)
     for i in range(len(pos)):
-        xMa[i] = pos[i] - genetic_map.base_pos[left_idx[i] - 1]
-        bMa[i] = (
-            genetic_map.base_pos[left_idx[i]] - genetic_map.base_pos[left_idx[i] - 1]
-        )
-        fbMfa[i] = (
-            genetic_map.gen_pos[left_idx[i]] - genetic_map.gen_pos[left_idx[i] - 1]
-        )
-    est_cM_diff = xMa / bMa * fbMfa
-    # TODO: Calculate cumulative genetic map distances.
-    return est_cM_diff
+        a = genetic_map.base_pos[
+            left_idx[i] - 1
+        ]  # Physical position at left of marker i
+        b = genetic_map.base_pos[left_idx[i]]  # Physical position at right of marker i
+        fa = genetic_map.gen_pos[
+            left_idx[i] - 1
+        ]  # Genetic map position at left of marker i
+        fb = genetic_map.gen_pos[
+            left_idx[i]
+        ]  # Genetic map position at right of marker i
+        est_gen_pos[i] = fa + (pos[i] - a) / (b - a) * (fa - fb)
+    return est_gen_pos
 
 
 @njit