Repository: SCECcode/pycsep

Issue:
Spatial binning doesn't properly account for floating point precision
I tested forecast models with a catalog whose events sometimes occur exactly on a cell boundary.
I found that [`region.get_index_of`](https://github.com/SCECcode/pycsep/blob/20be367d88bb81f64c8f3d5d433ffe856c6a9996/csep/core/regions.py#L618) generally does not assign them the spatial bin according to what is expected from

> _"inclusive on the lower bound and exclusive on the upper bound."_ [[`utils.calc.bin1d_vec` docstring](https://github.com/SCECcode/pycsep/blob/20be367d88bb81f64c8f3d5d433ffe856c6a9996/csep/utils/calc.py#L54)]

but the opposite (i.e., the neighboring bin). Although this function considers a tolerance (`numpy.finfo(numpy.float64).eps`) to (supposedly) account for numerical precision, I found that it is not sufficient to get the desired binning behavior, even after the revision in [commit 949d342c](https://github.com/SCECcode/pycsep/commit/949d342c181a984071e0f8e83a7460e631bb39c2).

I believe the corresponding [unit tests](https://github.com/SCECcode/pycsep/blob/20be367d88bb81f64c8f3d5d433ffe856c6a9996/tests/test_calc.py#L29) don't catch these literal "corner cases", as those tests are too constructed/synthetic and not based on real data such as a loaded CSEP region. Specifically, they do not consider issues due to _double_ precision (float64).


## Demo & Reasoning

### 1. The problem

Let's do a quick self check:

```python
region = csep.core.regions.california_relm_region()

idx = 3
origin = region.origins()[idx]
print(idx, origin)

idx2 = region.get_index_of(*origin)
print(idx2, region.origins()[idx2])  # not the same as `idx` and `origin`
```

outputs:
```
3 [-125.4   40.4]
2 [-125.4   40.3]
```

Ooops. Not even a bin's origin is considered to be inside its own bin (but the neighboring one)!

> Btw: the above output may represent the coordinates as more precise than they numerically are:
> ```
> print(idx, [x for x in origin])
> ```
> gives: 
> ```
> 3 [-125.39999999999999, 40.400000000000006]
> ```
> which highlights that our unit tests are currently too optimistic/simplistic.


### 2. The correction

Only when adding a tiny amount to the Latitude and Longitude of the event coordinates gives the desired behavior:

```python
tinyamount = 1e-12

idx3 = region.get_index_of(
    origin[0] + tinyamount,
    origin[1] + tinyamount
)
print(idx3, region.origins()[idx3])  # same as `idx` and `origin`
```

outputs:
```
3 [-125.4   40.4]
```

Now we can do this self check for _all_ origins in the region:

```python
# tinyamount = 0  # uncomment for raising an Error
for i, origin in enumerate(region.origins()):

    assert i == region.get_index_of(
        origin[0] + tinyamount,
        origin[1] + tinyamount
    )
```

Only if (at least) `tinyamount = 1e-12` is added, the loop runs without Error. 

> Side note: Actually, if `tinyamount = 0`, it wouldn't raise an `AssertionError`, but this [`ValueError`](https://github.com/SCECcode/pycsep/blob/20be367d88bb81f64c8f3d5d433ffe856c6a9996/csep/core/regions.py#L633) in `get_index_of` for the first bin, because `get_index_of` tries to assign a bin that lies East of it (which doesn't exist and is outside the region). You can get an `AssertionError` like so:
> ```python
> tinyamount = 0
> lons = set()
> for i, origin in enumerate(region.origins()):
>
>     if origin[0] not in lons:
>         lons.add(origin[0])
>         continue
>     
>     assert i == region.get_index_of(
>         origin[0] + tinyamount,
>         origin[1] + tinyamount
>     )
> ```

### 3. The complexity

But this numerical issue gets slightly more complicated than that: 

#### a) Influence of the region

The region influences (i) the `tinyamount` that is just sufficient to bin correctly, and (ii) which coordinate component(s) must be corrected. For instance for Italy, this is sufficient:

```python
region = csep.core.regions.italy_csep_region()
tinyamount = 2e-13  # Note #1: smaller than for California region
for i, origin in enumerate(region.origins()):

    assert i == region.get_index_of(
        origin[0],  # Note #2: adding `tinyamount` to the Longitude NOT needed (!)
        origin[1] + tinyamount
    )
```

I don't know why, but it's certainly related to numerical representations of values.

#### b) Influence of the floating-point format/precision

The just-sufficient `tinyamount` also depends on the floating-point format (i.e., precision). For instance, assuming the events have single precision (float32), this is necessary:

```python
tinyamount = 1.6e-6  # for Italy region
tinyamount = 3.1e-6  # for California region
for i, origin in enumerate(region.origins()):

    assert i == region.get_index_of(
        origin[0].astype(np.float32) + tinyamount,  # (adding `tinyamount` not necessary for Italy region)
        origin[1].astype(np.float32) + tinyamount
    )
```

#### Summary

I tried to find out the just-sufficient threshold also as function of the corresponding `.eps` and `.resolution` property of `np.finfo` for the used precision (e.g., `np.finfo(np.float64).eps`):

| Region     | Precision | `tinyamount` | `.resolution *` | `.eps *` |
|------------|-----------|--------------|----------------:|---------:|
| California | float64   | 1.0e-12      |          1000.0 |  ~4400.0 |
|    Italy   | float64   | 1.4e-13      |           140.0 |   ~630.0 |
| California | float32   | 3.1e-6       |             3.1 |    ~26.0 |
|    Italy   | float32   | 1.6e-6       |             1.6 |    ~13.0 |


Apparently, `tinyamount` not only depends on the region, but it also doesn't scale with the precision.



## Suggested solution(s)

Some options (with increasing personal preference):
 1. add the required `tinyamount` to event coordinates;
 2. use `bin1d_vec`'s `tol` argument, which is not used by the `get_index_of` wrapper;
 3. revise the [tolerance section of `bin1d_vec`](https://github.com/SCECcode/pycsep/blob/20be367d88bb81f64c8f3d5d433ffe856c6a9996/csep/utils/calc.py#L82-L84);
 4. simplify the tolerance section of `bin1d_vec` by only adding a `tinyamount`.

**Option 1** is only a temporary (quick-)fix.

**For Option 2**, note that `tol = tinyamount` is not sufficient due to the [follow-up 'absolute tolerance' section](https://github.com/SCECcode/pycsep/blob/20be367d88bb81f64c8f3d5d433ffe856c6a9996/csep/utils/calc.py#L87-L90C). It needs to be slightly higher, e.g., `tol = tinyamount * 1.1` in our cases.

**For Option 3**, we _could_ do:

```diff
-    a0_tol = numpy.abs(a0) * numpy.finfo(numpy.float64).eps
-    h_tol = numpy.abs(h) * numpy.finfo(numpy.float64).eps
-    p_tol = numpy.abs(p) * numpy.finfo(numpy.float64).eps
+#    tol_fact = numpy.finfo(numpy.float64).eps  # (current)
+    tol_fact = # NEW THRESHOLD NEEDED
+    a0_tol = numpy.abs(a0) * tol_fact
+    h_tol = numpy.abs(h) * tol_fact
+    p_tol = numpy.abs(p) * tol_fact
```

Compared to `tol`, `tol_fact` is relative; I tried to find out the necessary thresholds again:

| Region     | Precision | `tol_fact` | `.resolution *` | `.eps *` |
|------------|-----------|------------|----------------:|---------:|
| California | float64   | 4.14e-15   |            4.14 |     18.7 |
|    Italy   | float64   | 1.74e-15   |            1.74 |      7.9 |
| California | float32   | 2.38e-8    |          0.0238 |      0.2 |
|    Italy   | float32   | 2.35e-8    |          0.0235 |      0.2 |

It also doesn't scale with the precision of the floating-point format, so I don't see the benefit of this option.

**For Option 4**, we could revert the [whole current tolerance section](https://github.com/SCECcode/pycsep/blob/20be367d88bb81f64c8f3d5d433ffe856c6a9996/csep/utils/calc.py#L82C1-L93C75) that was introduced in [commit 949d342c](https://github.com/SCECcode/pycsep/commit/949d342c181a984071e0f8e83a7460e631bb39c2) to a [prior state](https://github.com/SCECcode/pycsep/blob/0b1feb94a7aa078f1f21f0994208e99449d31d77/csep/utils/calc.py#L74C1-L77C42) and extend it (keeping the `tol` argument):

```diff
-    a0_tol = numpy.abs(a0) * numpy.finfo(numpy.float64).eps
-    h_tol = numpy.abs(h) * numpy.finfo(numpy.float64).eps
-    p_tol = numpy.abs(p) * numpy.finfo(numpy.float64).eps
-
-
-    # absolute tolerance
-    if tol is None:
-        idx = numpy.floor((p + (p_tol + a0_tol) - a0) / (h - h_tol))
-    else:
-        idx = numpy.floor((p + (tol + a0_tol) - a0) / (h - h_tol))
+    tol = tol or 0
+
+    dtype = np.asarray(p).dtype
+    if dtype.name in ('float64', 'float32'):
+        float_tol = numpy.finfo(dtype).resolution * {
+            'float64': 2000,
+            'float32': 10,
+        }[dtype.name]
+        tol = max(tol, float_tol)
+
+    idx = numpy.floor((p + tol - a0) / h)
```

I added a safety factor of 2-3 compared to Table 1 just to be more on the safe side.

Btw: We could even _always_ add `numpy.finfo(numpy.float32).resolution * 10 = 1e-5` independent of the floating-point format, which is on the geographical order of ~1 meter. So instead of the above change, this one:
```diff
+    float_min_tol = 1e-5  # conservative (minimum) tolerance
+    tol = max(tol or 0, float_min_tol)
+    idx = numpy.floor((p + tol - a0) / h)
```
This would be the same as (and replace!) setting/fixing `tol=0.00001` when binning the magnitude [across various places in pyCSEP code](https://github.com/search?q=repo%3ASCECcode%2Fpycsep+tol%3D0.0000&type=code).


### The implication

_If a catalog contains such "corner case" events_, fixing this issue _will_ change the test statistic – but only slightly.
Nevertheless, it would lead to **irreproducible previous results** (e.g., in my case for IGPE, ~at the 3 to 5th decimal digit).
Exactly such slight difference (compared to an independent binning implementation) was the reason for spotting this issue.

<br>

---

**A personal note**

Considering that event locations are sometimes very rough, it might have been a better idea to define grid cells so that their centers lie on rough coordinates (e.g., 13.5°), rather than their boundaries (13.45° / 13.55°); i.e., shift the grid by half a bin size in each direction.

---

**P.S.**

@khawajasim: Could [`QuadtreeGrid2D._find_location`](https://github.com/SCECcode/pycsep/blob/20be367d88bb81f64c8f3d5d433ffe856c6a9996/csep/core/regions.py#L1087) be affected by the same issue?

Discussion:
@wsavran: Hi @mherrmann3, thanks for flagging this issue. This is something that came up and the reason for writing the ```bin1d_vec``` function. Sadly, the tests that I wrote surrounding this function were not sufficient to spot the condition that you've identified in this comment.

I think the first thing we should do is add in some test functions that recreate this issue and prevent builds from passing until this has been corrected.

We might also want to consider the uncertainty in the event locations, because in that case there might be some average rate contribution from the adjacent cells. For these 'border' events, the error ellipse for the location would likely intersect both neighboring cells. I would imagine that this would translate into smoother models scoring higher, which might not reflect nature.@wsavran: Hi again @mherrmann3, I'm looking at the [Decimal](https://docs.python.org/3/library/decimal.html) package to see if there are any solutions in here. It's a python module that allows for the precise representation of floating point arithmetic. @mherrmann3: I identified the root of the issue 🎉. No need to add `tinyamount`.

**tl;dr:** just round `h` to the intended bin size with
```diff
-        h = bins[1] - bins[0]
+        h = np.round(bins[1] - bins[0], 6)
```
in [this line](https://github.com/SCECcode/pycsep/blob/2feb9148/csep/utils/calc.py#L80). This change passes the "self-check loop" over all origins (see first comment) without errors.

This change is sufficient because the numerical precision of `h` was the culprit: it is imprecise beyond floating point resolution from the intended `0.1`, which originates from the floating-point operation in this very line – the simple subtraction of two floats – even though the bins, i.e., `region.xs`/`.ys`, are indeed "perfectly precise" (within the resolution of floats) thanks to [creating them with `csep.utils.calc.cleaner_range`](https://github.com/SCECcode/pycsep/blob/main/csep/core/regions.py#L784-L785). So perhaps some sort of [Catastrophic cancellation](https://en.wikipedia.org/wiki/Catastrophic_cancellation).
Eventually, `h_tol` doesn't sufficiently account for this imprecision.

- [ ] **I will create a PR with a more elaborate fix (not hardcoding `6`) + further considerations/improvements/optimizations (currently working on it).**

@wsavran: I also tried out the `decimal` package, but the performance is abysmal (at least 50x slower due to mapping via strings, which is the way the package works to circumvent numerical precision).


### How I got there / Aha-moment

I investigated the first origin in the self-check closer (`lon: -125.39999999999999, lat: 40.1`), which – as mentioned in the original comment – raises _"this [`ValueError`](https://github.com/SCECcode/pycsep/blob/20be367d88bb81f64c8f3d5d433ffe856c6a9996/csep/core/regions.py#L633) in `get_index_of`"_.

Independently determining the bin indices and corresponding locations with

```python
def bin_closest(p, bins):
    ind = np.argmin(np.abs(bins - p))
    print(f"ind: {ind} [{bins[ind]}]")

lon, lat = region.origins()[0]
bin_closest(lon, region.xs)
bin_closest(lat, region.ys)
```
yields
```
ind: 0 [-125.4]
ind: 86 [40.1]
```
which is the actual location of this origin (as expected).

Logging `idx` and `idy` inside `get_index_of` as calculated with the current `bin1d_vec` gives:

```
idx: 0 [-125.4]
idy: 85 [40.0]
```
🤔

As a reminder, the core of `bin1d_vec` is:

```
a0 = numpy.min(bins)
h = bins[1] - bins[0]
idx = numpy.floor((p + (p_tol + a0_tol) - a0) / (h - h_tol))
```

Specifically, `(p + (p_tol + a0_tol) - a0) / (h - h_tol)` yields for the `lat` coordinate:

```
85.99999999999896
```
which `numpy.floor` rounds down to `85`. But why, considering all `tol`erances? Let's log all terms:
```
                        p: 40.1
                       a0: 31.5
                   p - a0: 8.600000000000001
p + (p_tol + a0_tol) - a0: 8.600000000000016

                        h: 0.10000000000000142
                    h_tol: 0.00000000000000002
      (float64 resolution: 0.00000000000000100 - just for comparison)
                h - h_tol: 0.1000000000000014
```
**Aha!** This means `h - h_tol` is still too large! (numerically, it should be exactly or just "below" the actual binning of `0.1` to get a "`p/h` ratio" that is numerically `>=86.0` (so that when rounded down yields the desired `86.0`).

Apparently, the `h_tol`erance consideration doesn't account for numerical precision/resolution issues, especially when they derive from an additional preceding floating-point _operation_. The deviation of `h` from `0.1` is even greater than the `float64` resolution (logged above for comparison), which is yet 50x greater than `h_tol`.

> Btw: `tinyamount` just sufficiently nudged `p` to avoid these numerical effects of a tiny bit smaller "`p/h` ratio".

Solution patch:
diff --git a/csep/core/catalogs.py b/csep/core/catalogs.py
index b4257bb..dc1d880 100644
--- a/csep/core/catalogs.py
+++ b/csep/core/catalogs.py
@@ -395,7 +395,7 @@ class AbstractBaseCatalog(LoggingMixin):
     def get_mag_idx(self):
         """ Return magnitude index from region magnitudes """
         try:
-            return bin1d_vec(self.get_magnitudes(), self.region.magnitudes, tol=0.00001, right_continuous=True)
+            return bin1d_vec(self.get_magnitudes(), self.region.magnitudes, right_continuous=True)
         except AttributeError:
             raise CSEPCatalogException("Cannot return magnitude index without self.region.magnitudes")
 
@@ -699,16 +699,21 @@ class AbstractBaseCatalog(LoggingMixin):
         event_flag[idx] = 1
         return event_flag
 
-    def magnitude_counts(self, mag_bins=None, tol=0.00001, retbins=False):
-        """ Computes the count of events within mag_bins
-
+    def magnitude_counts(self, mag_bins=None, tol=None, retbins=False):
+        """ Computes the count of events within magnitude bins
 
         Args:
-            mag_bins: uses csep.utils.constants.CSEP_MW_BINS as default magnitude bins
-            retbins (bool): if this is true, return the bins used
+            mag_bins (array-like): magnitude bin edges, by default tries to use magnitude bin edges
+                                   associated with region, otherwise :data:`csep.utils.constants.CSEP_MW_BINS`
+            tol (float): overwrite numerical tolerance, by default determined automatically from the
+                         magnitudes' dtype to account for the limited precision of floating-point values.
+                         Only necessary to specify if the magnitudes were subject to some
+                         floating-point operations after loading or generating them
+                         (increased roundoff error, see :func:`csep.utils.calc.bin1d_vec`).
+            retbins (bool): if true, return the used bins in a tuple together with the counts.
 
         Returns:
-            numpy.ndarray: showing the counts of hte events in each magnitude bin
+            numpy.ndarray: events counts in each magnitude bin
         """
         # todo: keep track of events that are ignored
         if mag_bins is None:
@@ -734,18 +739,24 @@ class AbstractBaseCatalog(LoggingMixin):
         else:
             return out
 
-    def spatial_magnitude_counts(self, mag_bins=None, tol=0.00001):
-        """ Return counts of events in space-magnitude region.
+    def spatial_magnitude_counts(self, mag_bins=None, tol=None):
+        """ Return counts of events in space-magnitude bins.
 
-        We figure out the index of the polygons and create a map that relates the spatial coordinate in the
-        Cartesian grid with the polygon in region.
+        It figures out the index of the polygons and maps the spatial coordinates in the Cartesian
+        grid with the polygon in region.
 
         Args:
-            mag_bins (list, numpy.array): magnitude bins (optional), if empty tries to use magnitude bins associated with region
-            tol (float): tolerance for comparisons within magnitude bins
+            mag_bins (array-like): magnitude bin edges (optional), by default uses magnitude bin edges
+                                   associated with region.
+            tol (float): overwrite numerical tolerance, by default determined automatically from the
+                         magnitudes' dtype to account for the limited precision of floating-point values.
+                         Only necessary to specify if the magnitudes were subject to some
+                         floating-point operations after loading or generating them
+                         (increased roundoff error, see :func:`csep.utils.calc.bin1d_vec`).
 
         Returns:
-            output: unnormalized event count in each bin, 1d ndarray where index corresponds to midpoints
+            numpy.ndarray: unnormalized event count in each space-magnitude bin (2d, with indices
+                           corresponding to spatial midpoints and magnitude bin, respectively)
 
         """
         # make sure region is specified with catalog
diff --git a/csep/core/forecasts.py b/csep/core/forecasts.py
index 4b2f863..37860e9 100644
--- a/csep/core/forecasts.py
+++ b/csep/core/forecasts.py
@@ -213,16 +213,21 @@ class MarkedGriddedDataSet(GriddedDataSet):
         """ Returns counts of events in magnitude bins """
         return numpy.sum(self.data, axis=0)
 
-    def get_magnitude_index(self, mags, tol=0.00001):
+    def get_magnitude_index(self, mags, tol=None):
         """ Returns the indices into the magnitude bins of selected magnitudes
 
         Note: the right-most bin is treated as extending to infinity.
 
         Args:
-            mags (array-like): list of magnitudes
+            mags (array-like): magnitudes bin edges.
+            tol (float): overwrite numerical tolerance, by default determined automatically from the
+                         magnitudes' dtype to account for the limited precision of floating-point values.
+                         Only necessary to specify if the magnitudes were subject to some
+                         floating-point operations after loading or generating them
+                         (increased roundoff error, see :func:`csep.utils.calc.bin1d_vec`).
 
         Returns:
-            idm (array-like): indices corresponding to mags
+            numpy.ndarray: indices corresponding to the magnitude bins
 
         Raises:
             ValueError
diff --git a/csep/core/regions.py b/csep/core/regions.py
index 04795c0..21f7d8e 100644
--- a/csep/core/regions.py
+++ b/csep/core/regions.py
@@ -313,12 +313,7 @@ def magnitude_bins(start_magnitude, end_magnitude, dmw):
     Returns:
         bin_edges (numpy.ndarray)
     """
-    # convert to integers to prevent accumulating floating point errors
-    const = 10000
-    start = numpy.floor(const * start_magnitude)
-    end = numpy.floor(const * end_magnitude)
-    d = const * dmw
-    return numpy.arange(start, end + d / 2, d) / const
+    return cleaner_range(start_magnitude, end_magnitude, dmw)
 
 def create_space_magnitude_region(region, magnitudes):
     """Simple wrapper to create space-magnitude region """
@@ -495,7 +490,7 @@ def compute_vertices(origin_points, dh, tol=numpy.finfo(float).eps):
     """
     return list(map(lambda x: compute_vertex(x, dh, tol=tol), origin_points))
 
-def _bin_catalog_spatio_magnitude_counts(lons, lats, mags, n_poly, mask, idx_map, binx, biny, mag_bins, tol=0.00001):
+def _bin_catalog_spatio_magnitude_counts(lons, lats, mags, n_poly, mask, idx_map, binx, biny, mag_bins, tol=None):
     """
     Returns a list of event counts as ndarray with shape (n_poly, n_cat) where each value
     represents the event counts within the polygon.
@@ -625,8 +620,8 @@ class CartesianGrid2D:
         Returns:
             idx: ndarray-like
         """
-        idx = bin1d_vec(numpy.array(lons), self.xs)
-        idy = bin1d_vec(numpy.array(lats), self.ys)
+        idx = bin1d_vec(lons, self.xs)
+        idy = bin1d_vec(lats, self.ys)
         if numpy.any(idx == -1) or numpy.any(idy == -1):
             raise ValueError("at least one lon and lat pair contain values that are outside of the valid region.")
         if numpy.any(self.bbox_mask[idy, idx] == 1):
@@ -1061,7 +1056,7 @@ class QuadtreeGrid2D:
         self.cell_area = cell_area
         return self.cell_area
 
-    def get_index_of(self, lons, lats): 
+    def get_index_of(self, lons, lats):
         """ Returns the index of lons, lats in self.polygons
 
         Args:
@@ -1182,7 +1177,7 @@ class QuadtreeGrid2D:
         out = numpy.zeros([len(self.quadkeys), len(mag_bins)])
 
         idx_loc = self.get_index_of(lon, lat)
-        idx_mag = bin1d_vec(mag, mag_bins, tol=0.00001, right_continuous=True)
+        idx_mag = bin1d_vec(mag, mag_bins, right_continuous=True)
 
         numpy.add.at(out, (idx_loc, idx_mag), 1)
 
diff --git a/csep/utils/calc.py b/csep/utils/calc.py
index dc208cd..ad46b59 100644
--- a/csep/utils/calc.py
+++ b/csep/utils/calc.py
@@ -51,67 +51,75 @@ def discretize(data, bin_edges, right_continuous=False):
     x_new = bin_edges[idx]
     return x_new
 
+def _get_tolerance(v):
+    """Determine numerical tolerance due to limited precision of floating-point values.
+
+    ... to account for roundoff error.
+
+    In other words, returns a maximum possible difference that can be considered negligible.
+    Only relevant for floating-point values.
+    """
+    if issubclass(v.dtype.type, numpy.floating):
+        return numpy.abs(v) * numpy.finfo(v.dtype).eps
+    return 0  # assuming it's an int
+
 def bin1d_vec(p, bins, tol=None, right_continuous=False):
     """Efficient implementation of binning routine on 1D Cartesian Grid.
 
-    Returns the indices of the points into bins. Bins are inclusive on the lower bound
+    Bins are inclusive on the lower bound
     and exclusive on the upper bound. In the case where a point does not fall within the bins a -1
     will be returned. The last bin extends to infinity when right_continuous is set as true.
 
+    To correctly bin points that are practically on a bin edge, this function accounts for the
+    limited precision of floating-point numbers (the roundoff error) with a numerical tolerance.
+    If the provided points were subject to some floating-point operations after loading or
+    generating them, the roundoff error increases (which is not accounted for) and requires
+    overwriting the `tol` argument.
+
     Args:
-        p (array-like): Point(s) to be placed into b
-        bins (array-like): bins to considering for binning, must be monotonically increasing
-        right_continuous (bool): if true, consider last bin extending to infinity
+        p (array-like): Point(s) to be placed into bins.
+        bins (array-like): bin edges; must be sorted (monotonically increasing)
+        tol (float): overwrite numerical tolerance, by default determined automatically from the
+                     points' dtype to account for the roundoff error.
+        right_continuous (bool): if true, consider last bin extending to infinity.
 
     Returns:
-        idx (array-like): indexes hashed into grid
+        numpy.ndarray: indices for the points corresponding to the bins.
 
     Raises:
         ValueError:
     """
-    bins = numpy.array(bins)
-    p = numpy.array(p)
-    a0 = numpy.min(bins)
-    # if user supplies only a single bin, do 2 things: 1) fix right continuous to true, and use of h is arbitrary
+    bins = numpy.asarray(bins)
+    p = numpy.asarray(p)
+
+    # if not np.all(bins[:-1] <= bins[1:]):  # check for sorted bins, which is a requirement
+    #     raise ValueError("Bin edges are not sorted.")  # (pyCSEP always passes sorted bins)
+    a0 = bins[0]
     if bins.size == 1:
+        # for a single bin, set `right_continuous` to true; h is now arbitrary
         right_continuous = True
-        h = 1
+        h = 1.
     else:
         h = bins[1] - bins[0]
 
-    a0_tol = numpy.abs(a0) * numpy.finfo(numpy.float64).eps
-    h_tol = numpy.abs(h) * numpy.finfo(numpy.float64).eps
-    p_tol = numpy.abs(p) * numpy.finfo(numpy.float64).eps
-
-    # absolute tolerance
-    if tol is None:
-        idx = numpy.floor((p + (p_tol + a0_tol) - a0) / (h - h_tol))
-    else:
-        idx = numpy.floor((p + (tol + a0_tol) - a0) / (h - h_tol))
     if h < 0:
         raise ValueError("grid spacing must be positive and monotonically increasing.")
-    # account for floating point uncertainties by considering extreme case
+
+    # account for floating point precision
+    a0_tol = _get_tolerance(a0)
+    h_tol = a0_tol  # must be based on *involved* numbers
+    p_tol = tol or _get_tolerance(p)
+
+    idx = numpy.floor((p - a0 + p_tol + a0_tol) / (h - h_tol))
+    idx = numpy.asarray(idx)  # assure idx is an array
 
     if right_continuous:
         # set upper bin index to last
-        try:
-            idx[(idx < 0)] = -1
-            idx[(idx >= len(bins) - 1)] = len(bins) - 1
-        except TypeError:
-            if idx >= len(bins) - 1:
-                idx = len(bins) - 1
-            if idx < 0:
-                idx = -1
+        idx[idx < 0] = -1
+        idx[idx >= len(bins) - 1] = len(bins) - 1
     else:
-        try:
-            idx[((idx < 0) | (idx >= len(bins)))] = -1
-        except TypeError:
-            if idx < 0 or idx >= len(bins):
-                idx = -1
-    try:
-        idx = idx.astype(numpy.int64)
-    except AttributeError:
-        idx = int(idx)
+        idx[(idx < 0) | (idx >= len(bins))] = -1
+    idx = idx.astype(numpy.int64)
     return idx
 
 def _compute_likelihood(gridded_data, apprx_rate_density, expected_cond_count, n_obs):
@@ -215,12 +223,13 @@ def cleaner_range(start, end, h):
     Returns:
         bin_edges (numpy.ndarray)
     """
-    # convert to integers to prevent accumulating floating point errors
-    const = 100000
-    start = numpy.floor(const * start)
-    end = numpy.floor(const * end)
-    d = const * h
-    return numpy.arange(start, end + d / 2, d) / const
+    # determine required scaling to account for decimal places of bin edges and stepping
+    num_decimals_bins = len(str(float(start)).split('.')[1])
+    scale = max(10**num_decimals_bins, 1 / h)
+    start = numpy.round(scale * start)
+    end = numpy.round(scale * end)
+    d = scale * h
+    return numpy.arange(start, end + d / 2, d) / scale
 
 def first_nonnan(arr, axis=0, invalid_val=-1):
     mask = arr==arr


Test patch:
diff --git a/tests/test_calc.py b/tests/test_calc.py
index 618443d..000d14c 100644
--- a/tests/test_calc.py
+++ b/tests/test_calc.py
@@ -106,6 +106,13 @@ class TestBin1d(unittest.TestCase):
         expected = [-1, -1, 0, 0, 0, 0, 0, -1]
         self.assertListEqual(test.tolist(), expected)
 
+    def test_bin1d_single_bin2(self):
+        data = [4.0]
+        bin_edges = [3.0]
+        test = bin1d_vec(data, bin_edges)
+        expected = [0]
+        self.assertListEqual(test.tolist(), expected)
+
     def test_upper_limit_right_continuous(self):
         data = [40, 40, 40]
         bin_edges = [0, 10, 20, 30]
@@ -134,18 +141,27 @@ class TestBin1d(unittest.TestCase):
         expected = [-1, 3, -1]
         self.assertListEqual(test.tolist(), expected)
 
-    def test_scalar_outside(self):
-        from csep.utils.calc import bin1d_vec
-        mbins = numpy.arange(5.95, 9, 0.1)  # This gives bins from 5.95 to 8.95
-        idx = bin1d_vec(5.95, mbins, tol=0.00001, right_continuous=True)
-        self.assertEqual(idx, 0)
+    def test_scalar_inside(self):
+        mbins = numpy.arange(5.95, 9, 0.1)  # (magnitude) bins from 5.95 to 8.95
 
-        idx = bin1d_vec(6, mbins, tol=0.00001, right_continuous=True)  # This would give 0: Which is fine.
-        self.assertEqual(idx, 0)
+        for i, m in enumerate(mbins):
+            idx = bin1d_vec(m, mbins, right_continuous=True)  # corner cases
+            self.assertEqual(idx, i)
 
-        idx = bin1d_vec(5, mbins, tol=0.00001, right_continuous=True)
-        self.assertEqual(idx, -1)
+            idx = bin1d_vec(m + 0.05, mbins, right_continuous=True)  # center bins
+            self.assertEqual(idx, i)
 
-        idx = bin1d_vec(4, mbins, tol=0.00001, right_continuous=True)
+            idx = bin1d_vec(m + 0.099999999, mbins, right_continuous=True)  # end of bins
+            self.assertEqual(idx, i)
+
+        idx = bin1d_vec(10, mbins, right_continuous=True)  # larger than last bin edge
+        self.assertEqual(idx, mbins.size - 1)
+
+    def test_scalar_outside(self):
+        mbins = numpy.arange(5.95, 9, 0.1)  # (magnitude) bins from 5.95 to 8.95
+
+        idx = bin1d_vec(5, mbins, right_continuous=True)
         self.assertEqual(idx, -1)
 
+        idx = bin1d_vec(4, mbins, right_continuous=True)
+        self.assertEqual(idx, -1)
diff --git a/tests/test_regions.py b/tests/test_regions.py
index a1e2cce..ccf4fd3 100644
--- a/tests/test_regions.py
+++ b/tests/test_regions.py
@@ -3,7 +3,11 @@ import unittest
 
 import numpy
 
-from csep.core.regions import italy_csep_region, california_relm_region, nz_csep_region
+from csep.core.regions import (
+    italy_csep_region, italy_csep_collection_region,
+    california_relm_region, california_relm_collection_region,
+    nz_csep_region, nz_csep_collection_region,
+)
 
 def get_italy_region_fname():
     root_dir = os.path.dirname(os.path.abspath(__file__))
@@ -86,3 +90,76 @@ class TestNZRegion(unittest.TestCase):
         r = nz_csep_region()
         # they dont have to be in the same order, but they need
         numpy.testing.assert_array_equal(r.midpoints().sort(), self.from_dat.sort())
+
+class TestBin2d(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        super(TestBin2d, cls).setUpClass()
+
+        # (loading those is the bottleneck of this test case)
+        cls.regions = [
+            italy_csep_region(),
+            italy_csep_collection_region(),
+            california_relm_region(),
+            california_relm_collection_region(),
+            nz_csep_region(),
+            nz_csep_collection_region(),
+            # global_region()  # extreme slow-down (~2min loading + ~5min per loop + ~4s per vect)
+        ]
+
+    def test_bin2d_regions_origins(self):
+        """every origin must be inside its own bin
+        """
+
+        for region in self.regions:
+            origins = region.origins()
+            self._test_bin2d_region_loop(region, origins)
+            self._test_bin2d_region_vect(region, origins)
+            self._test_bin2d_region_vect(region, origins.astype(numpy.float32))
+
+    def test_bin2d_regions_midpoints(self):
+        """every midpoint must be inside its own bin
+        """
+
+        for region in self.regions:
+            midpoints = region.midpoints()
+            self._test_bin2d_region_loop(region, midpoints)
+            self._test_bin2d_region_vect(region, midpoints)
+            self._test_bin2d_region_vect(region, midpoints.astype(numpy.float32))
+
+    def test_bin2d_regions_endcorner(self):
+        """every corner point (~opposite end of the origin) must be inside its own bin
+        """
+
+        for region in self.regions:
+            frac = 0.9999999999
+            endcorners = region.origins() + frac*region.dh
+            self._test_bin2d_region_loop(region, endcorners)
+            self._test_bin2d_region_vect(region, endcorners)
+            frac = 0.999  # decrease frac for float32 due to its lower resolution
+            endcorners = region.origins() + frac*region.dh
+            self._test_bin2d_region_vect(region, endcorners.astype(numpy.float32))
+
+    def _test_bin2d_region_loop(self, region, coords):
+        """(slow) loop over origins; each time, calls bin1d_vec for lat & lon scalars
+        """
+
+        for i, origin in enumerate(coords):
+            idx = region.get_index_of(
+                origin[0],
+                origin[1],
+            )
+            self.assertEqual(i, idx)
+
+    def _test_bin2d_region_vect(self, region, coords):
+        """call bin1d_vec once for all lat origins & all lon origins
+
+        Besides, also tests if vectors with ndim=2 are consumed properly
+        (returns 2nd-order/nested list ([[...]]).
+        """
+
+        lons, lats = numpy.split(coords.T, 2)  # both have ndim=2!
+        test = region.get_index_of(lons, lats).tolist()  # nested list ([[...]])
+        expected = [numpy.arange(len(region.origins())).tolist()]  # embed in another list
+        self.assertListEqual(test, expected)