Repository: mie-lab/trackintel
Pull request #594 (2024-01-04 16:48:04)

BUG: `calculate_distance_matrix` computes all pairs
fixes  #592
One of the added tests has a very generous relative tolerance, this is due to #593.
I made this function probably quite a bit slower. I thought that correctness is more important and if nobody experienced this bug in 3 years maybe it isn't called that often either.
Tried to use pairwise distance for everything. That wasn't possible for the LineStrings as sklearn and scipy both enforce float as data type. Thus I copied the code out of sklearn and adapted it a bit. 

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/mie-lab/trackintel/pull/594?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mie-lab) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Comparison is base [(`e98354d`)](https://app.codecov.io/gh/mie-lab/trackintel/commit/e98354d662d6668316dba8d5f14c9448c9d8dc8b?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mie-lab) 93.15% compared to head [(`6da32a9`)](https://app.codecov.io/gh/mie-lab/trackintel/pull/594?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mie-lab) 93.29%.


<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master     #594      +/-   ##
==========================================
+ Coverage   93.15%   93.29%   +0.13%     
==========================================
  Files          33       33              
  Lines        2104     2072      -32     
  Branches      374      370       -4     
==========================================
- Hits         1960     1933      -27     
+ Misses        130      127       -3     
+ Partials       14       12       -2     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/mie-lab/trackintel/pull/594?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mie-lab).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=mie-lab).


Patch:
diff --git a/trackintel/geogr/distances.py b/trackintel/geogr/distances.py
index 38a4ce1..ab4beb6 100644
--- a/trackintel/geogr/distances.py
+++ b/trackintel/geogr/distances.py
@@ -1,3 +1,4 @@
+import itertools
 import math
 import multiprocessing
 import warnings
@@ -82,44 +83,43 @@ def point_haversine_dist(lon_1, lat_1, lon_2, lat_2, r=6371000, float_flag=False
     return r * np.arccos(cos_lat_d - cos_lat1 * cos_lat2 * (1 - cos_lon_d))
 
 
-def calculate_distance_matrix(X, Y=None, dist_metric="haversine", n_jobs=0, **kwds):
+def calculate_distance_matrix(X, Y=None, dist_metric="haversine", n_jobs=None, **kwds):
     """
-    Calculate a distance matrix based on a specific distance metric.
+    Compute the distance matrix from a vector array X and optional Y.
 
+    X and Y can either be Point geometries or LineString geometries.
     If only X is given, the pair-wise distances between all elements in X are calculated.
     If X and Y are given, the distances between all combinations of X and Y are calculated.
-    Distances between elements of X and X, and distances between elements of Y and Y are not calculated.
 
     Parameters
     ----------
-    X : Trackintel class
+    X : GeoDataFrame
 
-    Y : Trackintel class, optional
-        Should be of the same type as X
+    Y : GeoDataFrame, optional
+        Should have same geometry type as X
 
     dist_metric: {{'haversine', 'euclidean', 'dtw', 'frechet'}}, optional
         The distance metric to be used for calculating the matrix. By default 'haversine.
 
-        For staypoints or positionfixes, a common choice is 'haversine' or 'euclidean'. This function wraps around
-        the ``pairwise_distance`` function from scikit-learn if only `X` is given and wraps around the
-        ``scipy.spatial.distance.cdist`` function if X and Y are given.
+        For Point geometries we provide the 'haversine' metric.
+        This function then wraps `sklearn.metrics.pairwise_distances`.
         Therefore the following metrics are also accepted:
 
-        via ``scikit-learn``: `['cityblock', 'cosine', 'euclidean', 'l1', 'l2', 'manhattan']`
+        - via ``scikit-learn``: `['cityblock', 'cosine', 'euclidean', 'l1', 'l2', 'manhattan']`
 
-        via ``scipy.spatial.distance``: `['braycurtis', 'canberra', 'chebyshev', 'correlation', 'dice', 'hamming', 'jaccard',
+        - via ``scipy.spatial.distance``: `['braycurtis', 'canberra', 'chebyshev', 'correlation', 'dice', 'hamming', 'jaccard',
         'kulsinski', 'mahalanobis', 'minkowski', 'rogerstanimoto', 'russellrao', 'seuclidean', 'sokalmichener',
         'sokalsneath', 'sqeuclidean', 'yule']`
 
-        For triplegs, common choice is 'dtw' or 'frechet'. This function uses the implementation
-        from similaritymeasures.
+        For LineStrings, we provide the metrics {'dtw', 'frechet'} via the implementation from similaritymeasures.
 
     n_jobs: int, optional
-        Number of cores to use: 'dtw', 'frechet' and all distance metrics from `pairwise_distance` (only available
-        if only X is given) are parallelized. By default 1.
+        The number of jobs to use for the computation. Ignored for LineStrings.
+        None means 1 unless in a joblib.parallel_backend context. -1 means using all processors.
+        See `sklearn.metrics.pairwise_distances` for more informations.
 
     **kwds:
-        optional keywords passed to the distance functions.
+        Optional keywords passed to the distance functions.
 
     Returns
     -------
@@ -133,108 +133,53 @@ def calculate_distance_matrix(X, Y=None, dist_metric="haversine", n_jobs=0, **kw
     >>> pfs.calculate_distance_matrix(dist_metric="haversine")
     """
     geom_type = X.geometry.iat[0].geom_type
-    if Y is None:
-        Y = X
-    assert (
-        Y.geometry.iat[0].geom_type == Y.geometry.iat[0].geom_type
-    ), "x and y need same geometry type (only first column checked)"
+    if Y is not None and Y.geometry.iloc[0].geom_type != geom_type:
+        raise ValueError("X and Y need to have same geometry type.")
+    if geom_type not in ["Point", "LineString"]:
+        raise ValueError(f"We only support 'Point' and 'LineString'. Your geometry is {geom_type}")
 
     if geom_type == "Point":
-        x1 = X.geometry.x.values
-        y1 = X.geometry.y.values
-        x2 = Y.geometry.x.values
-        y2 = Y.geometry.y.values
-
         if dist_metric == "haversine":
-            # create point pairs for distance calculation
-            nx = len(X)
-            ny = len(Y)
-
-            # if y != x they could have different dimensions
-            if ny >= nx:
-                ix_1, ix_2 = np.triu_indices(nx, k=1, m=ny)
-                trilix = np.tril_indices(nx, k=-1, m=ny)
-            else:
-                ix_1, ix_2 = np.tril_indices(nx, k=-1, m=ny)
-                trilix = np.triu_indices(nx, k=1, m=ny)
+            # curry our haversine distance
+            def haversine_curry(a, b, **_):
+                return point_haversine_dist(*a, *b, float_flag=True)
 
-            x1 = x1[ix_1]
-            y1 = y1[ix_1]
-            x2 = x2[ix_2]
-            y2 = y2[ix_2]
+            dist_metric = haversine_curry
+        X = shapely.get_coordinates(X.geometry)
+        Y = shapely.get_coordinates(Y.geometry) if Y is not None else X
+        return pairwise_distances(X, Y, metric=dist_metric, n_jobs=n_jobs, **kwds)
 
-            d = point_haversine_dist(x1, y1, x2, y2)
+    # geom_type == "LineString"
+    # for LineStrings we cannot use pairwise_distance because it enforces float in its array
+    if dist_metric == "dtw":
 
-            D = np.zeros((nx, ny))
-            D[(ix_1, ix_2)] = d
+        def dtw(a, b, **kwds):
+            return similaritymeasures.dtw(a, b, **kwds)[0]
 
-            # mirror triangle matrix to be conform with scikit-learn format and to
-            # allow for non-symmetric distances in the future
-            D[trilix] = D.T[trilix]
-
-        else:
-            xy1 = np.concatenate((x1.reshape(-1, 1), y1.reshape(-1, 1)), axis=1)
-
-            if Y is not None:
-                xy2 = np.concatenate((x2.reshape(-1, 1), y2.reshape(-1, 1)), axis=1)
-                D = cdist(xy1, xy2, metric=dist_metric, **kwds)
-            else:
-                D = pairwise_distances(xy1, metric=dist_metric, n_jobs=n_jobs)
-
-        return D
-
-    elif geom_type == "LineString":
-        if dist_metric in ["dtw", "frechet"]:
-            # these are the preparation steps for all distance functions based only on coordinates
-
-            if dist_metric == "dtw":
-                d_fun = partial(similaritymeasures.dtw, **kwds)
-            else:
-                d_fun = partial(similaritymeasures.frechet_dist, **kwds)
-
-            # get combinations of distances that have to be calculated
-            nx = len(X)
-            ny = len(Y)
-
-            if ny >= nx:
-                ix_1, ix_2 = np.triu_indices(nx, k=1, m=ny)
-                trilix = np.tril_indices(nx, k=-1, m=ny)
-            else:
-                ix_1, ix_2 = np.tril_indices(nx, k=-1, m=ny)
-                trilix = np.triu_indices(nx, k=1, m=ny)
-
-            # get the coordinates as list of each LineString
-            left = list(X.iloc[ix_1].geometry.apply(lambda x: x.coords))
-            right = list(Y.iloc[ix_2].geometry.apply(lambda x: x.coords))
-
-            # map the combinations to the distance function
-            if n_jobs == -1 or n_jobs > 1:
-                if n_jobs == -1:
-                    n_jobs = multiprocessing.cpu_count()
-                with multiprocessing.Pool(processes=n_jobs) as pool:
-                    left_right = list(zip(left, right))
-                    res = list(pool.starmap(d_fun, left_right))
-            else:
-                res = list(map(d_fun, left, right))
-
-            if dist_metric == "dtw":
-                # the first return is the dtw distance, see docs of similaritymeasures.dtw
-                d = [dist[0] for dist in res]
-            else:
-                d = res
-
-            # write results to (symmetric) distance matrix
-            D = np.zeros((nx, ny))
-            D[(ix_1, ix_2)] = d
-            D[trilix] = D.T[trilix]
-            return D
-
-        else:
-            raise AttributeError(
-                "Metric unknown. We only support ['dtw', 'frechet'] for LineStrings. " f"You passed {dist_metric}"
-            )
+        dist_metric = dtw
+    elif dist_metric == "frechet":
+        dist_metric = similaritymeasures.frechet_dist
+    else:
+        raise ValueError(f"Metric '{dist_metric}' unknown. We only support ['dtw', 'frechet'] for LineStrings")
+    X = X.geometry.values
+    Y = Y.geometry.values if Y is not None else X
+
+    # the following code is adapted from scikit-learn pairwise_distance
+    # https://github.com/scikit-learn/scikit-learn/blob/3f89022fa04d293152f1d32fbc2a5bdaaf2df364/sklearn/metrics/pairwise.py#L1784
+    out = np.zeros((len(X), len(Y)), dtype="float")
+    if X is Y:
+        # Only calculate metric for upper triangle
+        iterator = itertools.combinations(range(len(X)), 2)
+        for i, j in iterator:
+            out[i, j] = dist_metric(X[i].coords, Y[j].coords, **kwds)
+        # Make symmetric
+        out = out + out.T
     else:
-        raise AttributeError(f"We only support 'Point' and 'LineString'. Your geometry is {geom_type}")
+        # Calculate all cells
+        iterator = itertools.product(range(len(X)), range(len(Y)))
+        for i, j in iterator:
+            out[i, j] = dist_metric(X[i].coords, Y[j].coords, **kwds)
+    return out
 
 
 def meters_to_decimal_degrees(meters, latitude):


Test patch:
diff --git a/tests/geogr/test_distances.py b/tests/geogr/test_distances.py
index c71d55b..2cf3f77 100644
--- a/tests/geogr/test_distances.py
+++ b/tests/geogr/test_distances.py
@@ -6,6 +6,7 @@ import numpy as np
 import pandas as pd
 import pytest
 from geopandas.testing import assert_geodataframe_equal
+import shapely
 from shapely import wkt
 from shapely.geometry import LineString, MultiLineString, Point
 from sklearn.metrics import pairwise_distances
@@ -82,6 +83,42 @@ def example_triplegs():
     return pfs, tpls
 
 
+@pytest.fixture
+def geolife_sp():
+    """Load stored geolife staypoints"""
+    sp_file = os.path.join("tests", "data", "geolife", "geolife_staypoints.csv")
+    return ti.read_staypoints_csv(sp_file, tz="utc", index_col="id")
+
+
+@pytest.fixture
+def two_pfs():
+    """Generate two Positionfixes with their euclidean distance"""
+    p1 = Point(0, 0)
+    p2 = Point(0, 4)
+    p3 = Point(3, 4)  # euclidean dist p1 to p3 = 5
+    p4 = Point(3, 0)  # euclidean dist p2 to p4 = 5
+
+    t = pd.Timestamp.now(tz="utc")  # time doesn't matter
+
+    data0 = [
+        {"user_id": 0, "tracked_at": t, "geom": p1},
+        {"user_id": 0, "tracked_at": t, "geom": p2},
+        {"user_id": 0, "tracked_at": t, "geom": p3},
+    ]
+    pfs0 = ti.Positionfixes(data=data0, geometry="geom", crs="EPSG:4326")
+    pfs0.index.name = "id"
+
+    data1 = [
+        {"user_id": 0, "tracked_at": t, "geom": p1},
+        {"user_id": 0, "tracked_at": t, "geom": p4},
+    ]
+    pfs1 = ti.Positionfixes(data=data1, geometry="geom", crs="EPSG:4326")
+    pfs1.index.name = "id"
+    euc00 = np.array([[0.0, 4.0, 5.0], [4.0, 0.0, 3.0], [5.0, 3.0, 0.0]])
+    euc01 = np.array([[0.0, 3.0], [4.0, 5.0], [5.0, 4.0]])
+    return pfs0, euc00, pfs1, euc01
+
+
 class TestHaversineDist:
     def test_haversine_dist(self):
         """
@@ -108,11 +145,9 @@ class TestHaversineDist:
             haversine_output = point_haversine_dist(latlng[0], latlng[1], latlng[2], latlng[3])
             assert np.isclose(haversine_output, haversine, atol=0.1)
 
-    def test_haversine_vectorized(self):
-        sp_file = os.path.join("tests", "data", "geolife", "geolife_staypoints.csv")
-        sp = ti.read_staypoints_csv(sp_file, tz="utc", index_col="id")
-        x = sp.geometry.x.values
-        y = sp.geometry.y.values
+    def test_haversine_vectorized(self, geolife_sp):
+        x = geolife_sp.geometry.x.values
+        y = geolife_sp.geometry.y.values
 
         n = len(x)
         # our distance
@@ -134,7 +169,7 @@ class TestHaversineDist:
         d_theirs = D_theirs[ix_1, ix_2]
         assert np.sum(np.abs(d_ours - d_theirs)) < 0.01  # 1cm for 58 should be good enough
 
-    def test_example_from_sklean(self):
+    def test_example_from_sklearn(self):
         bsas = [-34.83333, -58.5166646]
         paris = [49.0083899664, 2.53844117956]
         bsas_in_radians = [radians(_) for _ in bsas]
@@ -149,12 +184,9 @@ class TestHaversineDist:
 class TestCalculate_distance_matrix:
     """Tests for the calculate_distance_matrix() function."""
 
-    def test_shape_for_different_array_length(self):
-        sp_file = os.path.join("tests", "data", "geolife", "geolife_staypoints.csv")
-        sp = ti.read_staypoints_csv(sp_file, tz="utc", index_col="id")
-
-        x = sp.iloc[0:5]
-        y = sp.iloc[5:15]
+    def test_shape_for_different_array_length(self, geolife_sp):
+        x = geolife_sp.iloc[0:5]
+        y = geolife_sp.iloc[5:15]
 
         d_euc1 = calculate_distance_matrix(X=x, Y=y, dist_metric="euclidean")
         d_euc2 = calculate_distance_matrix(X=y, Y=x, dist_metric="euclidean")
@@ -166,12 +198,9 @@ class TestCalculate_distance_matrix:
         assert np.isclose(0, np.sum(np.abs(d_euc1 - d_euc2.T)))
         assert np.isclose(0, np.sum(np.abs(d_hav1 - d_hav2.T)))
 
-    def test_keyword_combinations(self):
-        sp_file = os.path.join("tests", "data", "geolife", "geolife_staypoints.csv")
-        sp = ti.read_staypoints_csv(sp_file, tz="utc", index_col="id")
-
-        x = sp.iloc[0:5]
-        y = sp.iloc[5:15]
+    def test_keyword_combinations(self, geolife_sp):
+        x = geolife_sp.iloc[0:5]
+        y = geolife_sp.iloc[5:15]
 
         _ = calculate_distance_matrix(X=x, Y=y, dist_metric="euclidean", n_jobs=-1)
         _ = calculate_distance_matrix(X=y, Y=x, dist_metric="haversine", n_jobs=-1)
@@ -180,25 +209,42 @@ class TestCalculate_distance_matrix:
         d_euc = calculate_distance_matrix(X=x, Y=x, dist_metric="euclidean")
 
         assert not np.array_equal(d_mink1, d_mink2)
-        assert np.array_equal(d_euc, d_mink2)
+        assert np.allclose(d_euc, d_mink2, rtol=1e-8)
 
-    def test_compare_haversine_to_scikit_xy(self):
+    def test_compare_haversine_to_scikit_xy(self, geolife_sp):
         """Test the results using our haversine function and scikit function."""
-        sp_file = os.path.join("tests", "data", "geolife", "geolife_staypoints.csv")
-        sp = ti.read_staypoints_csv(sp_file, tz="utc", index_col="id")
-        our_d_matrix = calculate_distance_matrix(X=sp, Y=sp, dist_metric="haversine")
+        our_d_matrix = calculate_distance_matrix(X=geolife_sp, Y=geolife_sp, dist_metric="haversine")
 
-        x = sp.geometry.x.values
-        y = sp.geometry.y.values
+        x = geolife_sp.geometry.x.values
+        y = geolife_sp.geometry.y.values
 
-        x_rad = np.asarray([radians(_) for _ in x])
-        y_rad = np.asarray([radians(_) for _ in y])
-        yx = np.concatenate((y_rad.reshape(-1, 1), x_rad.reshape(-1, 1)), axis=1)
+        x_rad = np.radians(x)
+        y_rad = np.radians(y)
+        yx = np.stack((y_rad, x_rad), axis=1)
 
         their_d_matrix = pairwise_distances(yx, metric="haversine") * 6371000
         # atol = 10mm
         assert np.allclose(our_d_matrix, their_d_matrix, atol=0.01)
 
+    def test_compare_haversine_to_scikit_xy_with_Y(self, two_pfs):
+        pfs0, _, pfs1, _ = two_pfs
+        res01 = calculate_distance_matrix(X=pfs0, Y=pfs1, dist_metric="haversine")
+        rad0 = np.radians(shapely.get_coordinates(pfs0.geometry))
+        rad1 = np.radians(shapely.get_coordinates(pfs1.geometry))
+        sol01 = pairwise_distances(rad0, rad1, metric="haversine") * 6371000
+        # TODO: increase precision of haversine dist such we can lower the tolerance
+        # see issue #593 for more information
+        assert np.allclose(res01, sol01, rtol=1e-2)
+
+    def test_known_euclidean_distance(self, two_pfs):
+        """Test the result comparing to known euclidean distances"""
+        pfs0, euc00, pfs1, euc01 = two_pfs
+        res00 = calculate_distance_matrix(X=pfs0, dist_metric="euclidean")
+        res01 = calculate_distance_matrix(X=pfs0, Y=pfs1, dist_metric="euclidean")
+        print(res00)
+        assert np.all(euc00 == res00)
+        assert np.all(euc01 == res01)
+
     def test_trajectory_distance_dtw(self, geolife_tpls):
         """Calculate Linestring length using dtw, single and multi core."""
         tpls = geolife_tpls
@@ -247,9 +293,9 @@ class TestCalculate_distance_matrix:
         """Test if the an error is raised when passing unknown 'dist_metric'."""
         tpls = geolife_tpls
 
-        with pytest.raises(AttributeError):
+        with pytest.raises(ValueError, match=f"Metric '{12345}' unknown."):
             tpls.iloc[0:4].as_triplegs.calculate_distance_matrix(dist_metric=12345, n_jobs=1)
-        with pytest.raises(AttributeError):
+        with pytest.raises(ValueError, match=f"Metric '{'random'}' unknown."):
             tpls.iloc[0:4].as_triplegs.calculate_distance_matrix(dist_metric="random", n_jobs=1)
 
     def test_distance_error(self, single_linestring):
@@ -260,9 +306,17 @@ class TestCalculate_distance_matrix:
         gdf = gpd.GeoDataFrame(a_list, columns=["id", "geometry"]).set_geometry("geometry")
         gdf = gdf.set_crs("wgs84")
 
-        with pytest.raises(AttributeError):
+        with pytest.raises(ValueError, match="We only support 'Point' and 'LineString'."):
             calculate_distance_matrix(X=gdf, dist_metric="dtw", n_jobs=1)
 
+    def test_same_geometry_type(self, single_linestring, geolife_tpls):
+        """Test if error is raised if X and Y have different geometries"""
+        multi = MultiLineString([single_linestring, single_linestring])
+        t = [(0, multi), (1, multi)]
+        gdf = gpd.GeoDataFrame(t, columns=["id", "geometry"], geometry="geometry", crs="wgs84")
+        with pytest.raises(ValueError, match="X and Y need to have same geometry type."):
+            calculate_distance_matrix(gdf, geolife_tpls)
+
 
 class TestCheck_gdf_planar:
     """Tests for check_gdf_planar() method."""