Add min cM distance value to avoid division by zero

python/tests/beagle.py

@@ -101,6 +101,7 @@ def get_weights(genotyped_pos, ungenotyped_pos, genetic_map=None):
     :return: Weights and associated marker interval start indices.
     :rtype: tuple(numpy.ndarray, numpy.ndarray)
     """
+    MIN_CM_DIST = 1e-7  # See `ImpData.java` in BEAGLE 4.1 source code
     assert len(genotyped_pos) > 1
     assert len(ungenotyped_pos) > 0
     # Check that the two sets are mutually exclusive.
@@ -131,7 +132,8 @@ def get_weights(genotyped_pos, ungenotyped_pos, genetic_map=None):
             k = np.max(np.where(ungenotyped_pos[i] > genotyped_pos)[0])  # Inefficient
             cm_mP1_x = genotyped_cm[k + 1] - ungenotyped_cm[i]
             cm_mP1_m = genotyped_cm[k + 1] - genotyped_cm[k]
-            assert cm_mP1_m > 0, "Division by zero in weight calculation."
+            if cm_mP1_m == 0:
+                cm_mP1_m = MIN_CM_DIST
             weights[i] = cm_mP1_x / cm_mP1_m
         marker_interval_start[i] = k
     assert 0 <= np.min(weights)

python/tests/beagle_numba.py

@@ -67,6 +67,7 @@ def get_weights(genotyped_pos, ungenotyped_pos, genotyped_cm, ungenotyped_cm):
     :return: Weights and marker interval start indices.
     :rtype: tuple(numpy.ndarray, numpy.ndarray)
     """
+    MIN_CM_DIST = 1e-7  # See `ImpData.java` in BEAGLE 4.1 source code
     m = len(genotyped_pos)
     x = len(ungenotyped_pos)
     # Calculate weights for ungenotyped markers.
@@ -88,6 +89,8 @@ def get_weights(genotyped_pos, ungenotyped_pos, genotyped_cm, ungenotyped_cm):
             # Between two genotyped markers.
             cm_mP1_x = genotyped_cm[idx_m + 1] - ungenotyped_cm[idx_x]
             cm_mP1_m = genotyped_cm[idx_m + 1] - genotyped_cm[idx_m]
+            if cm_mP1_m == 0:
+                cm_mP1_m = MIN_CM_DIST
             weights[idx_x] = cm_mP1_x / cm_mP1_m
         marker_interval_start[idx_x] = idx_m
     return (weights, marker_interval_start)