Repository: popsim-consortium/stdpopsim (language: python)

Issue:
make intervals and annotations code consistent and make all intervals integer-valued
There's a number of operations on intervals in `utils.py` that are public but mostly undocumented. These are:
- `check_intervals_array_shape`
- `check_intervals_validity`
- `build_intervals_array`
- `mask_intervals`

Of these, I think that `check_intervals_array_shape` and probably `check_intervals_validity` should be private, and we should document the others.

Also, right now it's easy to get very confusing errors trying to pass in interval arrays having floats rather than ints. I think we should require that intervals are *always* of integer type, and give informative errors if they are not. We should also make sure we give an informative error in the case where nested lists are passed in instead of a numpy array. 

I'm of two minds whether we want to convert, say `[[0, 2.0], [4.0, 10]]` to the integer-valued numpy array.



Discussion / hints:
Furthermore, right now `Annotations` returns an array of floats. Is this desired, or do we want to *always* have intervals be ints? Options are:
- intervals are floats: easier, as RateMaps have float coordinates and contig.length is a floa and indeed all other coordinates in msprime/tskit, Conversion to int happens in the SLiM code.
- intervals are ints: more consistent as we're always doing things with integer coordinates actually.
Also - I think that `contig.add_DFE( )` should take an Annotation object for the `intervals=` argument, in which case it Does the Right Thing.
I vote for integers. I think floats are a legacy that we should try to remove. In particular, contig.length is a float to support the `length_multiplier` params, which should be deprecated as discussed in #670.
> Also - I think that `contig.add_DFE( )` should take an Annotation object for the `intervals=` argument, in which case it Does the Right Thing.

+1 -- i think we should definitely do this

Fix patch:
diff --git a/stdpopsim/annotations.py b/stdpopsim/annotations.py
index 15dedb282..e7787df88 100644
--- a/stdpopsim/annotations.py
+++ b/stdpopsim/annotations.py
@@ -83,7 +83,7 @@ def get_chromosome_annotations(self, id):
         if not self.is_cached():
             self.download()
         file_path = os.path.join(self.cache_path, self.file_pattern.format(id=chrom.id))
-        ret = np.loadtxt(file_path)
+        ret = np.loadtxt(file_path, dtype="int32")
         if len(ret) == 0:
             raise ValueError(f"No annotations found for {id}")
         return ret
diff --git a/stdpopsim/genomes.py b/stdpopsim/genomes.py
index 1611fee5d..22c619b0b 100644
--- a/stdpopsim/genomes.py
+++ b/stdpopsim/genomes.py
@@ -380,7 +380,7 @@ def add_dfe(self, intervals, DFE):
         :param array intervals: A valid set of intervals.
         :param DFE dfe: A DFE object.
         """
-        stdpopsim.utils.check_intervals_validity(intervals)
+        stdpopsim.utils._check_intervals_validity(intervals)
         for j, ints in enumerate(self.interval_list):
             self.interval_list[j] = stdpopsim.utils.mask_intervals(ints, intervals)
         self.dfe_list.append(DFE)
diff --git a/stdpopsim/utils.py b/stdpopsim/utils.py
index 20ef9dab0..1d3641179 100644
--- a/stdpopsim/utils.py
+++ b/stdpopsim/utils.py
@@ -9,7 +9,6 @@
 import tarfile
 import contextlib
 import numpy as np
-import tskit
 
 
 def is_valid_demographic_model_id(model_id):
@@ -184,19 +183,23 @@ def add_if_unique(chrom, syn):
             add_if_unique(chrom, "chr" + "M")
 
 
-def check_intervals_array_shape(intervals):
-    if (
-        (not isinstance(intervals, np.ndarray))
-        or (len(intervals.shape) != 2)
-        or (intervals.shape[1] < 2)
-    ):
+def _check_intervals_array_shape(intervals):
+    if (len(intervals.shape) != 2) or (intervals.shape[1] < 2):
         raise ValueError(
-            "Intervals must be 2D objects with at least 2 columns " "[left, right)."
+            "Intervals must be 2D numpy arrays with 2 columns " "[left, right)."
         )
 
 
-def check_intervals_validity(intervals, start=0, end=np.inf):
-    check_intervals_array_shape(intervals)
+def _check_intervals_validity(intervals, start=0, end=np.inf):
+    # is numpy?
+    if not isinstance(intervals, np.ndarray):
+        raise ValueError("Intervals must be a numpy array, did you try to pass a list?")
+    # are the intervals ints?
+    if (intervals.dtype != "int32") and (intervals.dtype != "int64"):
+        raise ValueError(
+            "Intervals must be arrays of type integer" "(i.e., with dtype='int')."
+        )
+    _check_intervals_array_shape(intervals)
     if np.any(intervals[:, 0] >= intervals[:, 1]):
         raise ValueError(
             "Left positions cannot be greater than or equal to right positions."
@@ -207,20 +210,6 @@ def check_intervals_validity(intervals, start=0, end=np.inf):
         raise ValueError("Intervals must be non-overlapping.")
 
 
-def build_intervals_array(intervals, start=0, end=np.inf):
-    """
-    Converts a 2D list or numpy.array to an np.int32 array, which is sorted by
-    the first axis. It also checks for the validity of intervals and non-overlappingness.
-    """
-    # TODO: allow floats that can be converted to ints?
-    intervals = tskit.util.safe_np_int_cast(intervals, dtype=np.int64)
-    check_intervals_array_shape(intervals)
-    sorter = intervals[:, 0].argsort()
-    intervals = intervals[sorter]
-    check_intervals_validity(intervals, start, end)
-    return intervals
-
-
 def mask_intervals(intervals, mask):
     """
     Removes the intervals of ``mask`` from those in ``intervals``.
@@ -229,8 +218,8 @@ def mask_intervals(intervals, mask):
     be the the intervals in ``intervals`` except with those in
     ``mask`` removed.
     """
-    check_intervals_validity(intervals)
-    check_intervals_validity(mask)
+    _check_intervals_validity(intervals)
+    _check_intervals_validity(mask)
     out = []
     last_mask_right = -1 * np.inf
     j = 0


Test patch:
diff --git a/tests/test_utils.py b/tests/test_utils.py
index e843f3269..6bce2e273 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -447,17 +447,7 @@ def test_intervals_array_shape(self):
         bad_arrays = (np.array([10, 20]), np.array([[[10, 20]]]))
         for bad_array in bad_arrays:
             with pytest.raises(ValueError):
-                utils.check_intervals_array_shape(intervals=bad_array)
-            with pytest.raises(ValueError):
-                utils.build_intervals_array(intervals=bad_array)
-        good_arrays = (
-            np.array([[10, 20], [30, 40]]),
-            np.array([[10, 20, 1], [30, 40, 2]]),
-        )
-        for good_array in good_arrays:
-            utils.check_intervals_array_shape(intervals=good_array)
-            built_array = utils.build_intervals_array(intervals=good_array)
-            assert (built_array == good_array).all()
+                utils._check_intervals_array_shape(intervals=bad_array)
 
     def test_invalid_intervals(self):
         # right <= left, overlapping, uncastable float->int
@@ -468,41 +458,17 @@ def test_invalid_intervals(self):
         )
         for invalid_array in invalid_arrays:
             with pytest.raises(ValueError):
-                utils.build_intervals_array(intervals=invalid_array)
-            with pytest.raises(ValueError):
-                utils.check_intervals_validity(intervals=invalid_array)
+                utils._check_intervals_validity(intervals=invalid_array)
         # interval outside [start, end)
         with pytest.raises(ValueError):
-            utils.build_intervals_array(intervals=[[10, 50]], start=20)
+            utils._check_intervals_validity(intervals=np.array([[10, 50]]), end=30)
         with pytest.raises(ValueError):
-            utils.check_intervals_validity(intervals=np.array([[10, 50]]), end=30)
-
-    def test_intervals_casting(self):
-        castable_intervals = (
-            [[10, 20], [20, 40]],
-            [[1, 2, 2], [2, 3, 2]],
-            np.array([[10, 20]]),
-        )
-        for intervals in castable_intervals:
-            casted = utils.build_intervals_array(intervals)
-            assert isinstance(casted, np.ndarray) and casted.dtype == np.int64
-            assert casted.shape[0] == len(intervals)
-            assert casted.shape[1] >= 2
-
-    def test_interval_sorting(self):
-        unsorted_intervals = (
-            np.array([[20, 30], [10, 20]]),
-            np.array([[10, 20], [20, 30], [2, 3]]),
-        )
-        for intervals in unsorted_intervals:
-            casted = utils.build_intervals_array(intervals)
-            assert not (np.all(np.diff(intervals[:, 0]) >= 0))
-            assert np.all(np.diff(casted[:, 0]) >= 0)
+            utils._check_intervals_validity(intervals=np.array([[10.2, 50.3]]))
 
     def intervals_to_set(self, intervals):
         # note: we may want to use intervaltree for more complicated operations
         # in the future
-        utils.check_intervals_validity(intervals)
+        utils._check_intervals_validity(intervals)
         return functools.reduce(
             lambda a, b: a | b, [set(range(a, b)) for a, b in intervals], set()
         )
@@ -537,8 +503,8 @@ def test_mask_intervals(self):
             ),
         ]:
             for u, v in [(a, b), (b, a), (a, empty_array), (empty_array, a)]:
-                u = np.array(u)
-                v = np.array(v)
+                u = np.array(u, dtype="int32")
+                v = np.array(v, dtype="int32")
                 umv = utils.mask_intervals(u, v)
                 assert u.dtype == umv.dtype
                 x = self.intervals_to_set(umv)