diff --git a/tests/test_metrics.py b/tests/test_metrics.py index c9981e89..2024a248 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -30,18 +30,10 @@ def test_dtw_accumulated_matrix(): s1 = cast([1, 2, 3], array_type) s2 = cast([1.0, 2.0, 2.0, 3.0], array_type) mask = tslearn.metrics.compute_mask(s1, s2, be=be) - matrix_1 = tslearn.metrics.dtw_accumulated_matrix(s1, s2, mask, be=be) - with pytest.deprecated_call(): - matrix_2 = tslearn.metrics.dtw_variants.accumulated_matrix( - to_time_series(s1), - to_time_series(s2), - tslearn.metrics.compute_mask(s1, s2, be=be), - be=be - ) + matrix = tslearn.metrics.dtw_accumulated_matrix(s1, s2, mask, be=be) expected = backend.array([[0.0, 1.0, 2.0, 6.0], [1.0, 0.0, 0.0, 1.0], [5.0, 1.0, 1.0, 0.0]]) - for matrix in [matrix_1, matrix_2]: - assert backend.belongs_to_backend(matrix) - np.testing.assert_array_equal(matrix, expected) + assert backend.belongs_to_backend(matrix) + np.testing.assert_array_equal(matrix, expected) def test_dtw(): @@ -69,24 +61,10 @@ def test_dtw(): if not backend.is_numpy: assert backend.belongs_to_backend(dist) - with pytest.deprecated_call(): - deprecated_path, deprecated_dist = tslearn.metrics.dtw_variants.dtw_path( - cast([1, 2, 3], array_type), - cast([1.0, 2.0, 2.0, 3.0, 4.0], array_type), - be=be - ) - assert deprecated_path == path - assert deprecated_dist == dist - with pytest.raises(ValueError): tslearn.metrics.dtw_path([], [], be=be) with pytest.raises(ValueError): tslearn.metrics.dtw_path([1], [[1, 2]], be=be) - with warnings.catch_warnings(): - with pytest.raises(ValueError): - tslearn.metrics.dtw_variants.dtw_path([], [], be=be) - with pytest.raises(ValueError): - tslearn.metrics.dtw_variants.dtw_path([1], [[1, 2]], be=be) # dtw n1, n2, d = 15, 10, 3 @@ -103,23 +81,10 @@ def test_dtw(): wrong_feature_size_y = cast(np.random.randn(n2, 2), array_type) tslearn.metrics.dtw(x, wrong_feature_size_y, be=be) - with pytest.deprecated_call(): - deprecated_dist = tslearn.metrics.dtw_variants.dtw( - cast([1, 2, 3], array_type), - cast([1.0, 2.0, 2.0, 3.0, 4.0], array_type), - be=be - ) - assert deprecated_dist == dist - with pytest.raises(ValueError): tslearn.metrics.dtw([], [], be=be) with pytest.raises(ValueError): tslearn.metrics.dtw([1], [[1, 2]], be=be) - with warnings.catch_warnings(): - with pytest.raises(ValueError): - tslearn.metrics.dtw_variants.dtw([], [], be=be) - with pytest.raises(ValueError): - tslearn.metrics.dtw_variants.dtw([1], [[1, 2]], be=be) # cdist_dtw dists = tslearn.metrics.cdist_dtw( @@ -146,15 +111,6 @@ def test_dtw(): np.testing.assert_array_equal(dists, parallel_dists) assert backend.belongs_to_backend(parallel_dists) - with pytest.deprecated_call(): - deprecated_dists = tslearn.metrics.dtw_variants.cdist_dtw( - cast([[1, 2, 2, 3], [1.0, 2.0, 3.0, 4.0]], array_type), - [[1, 2, 3], [2, 3, 4, 5]], # The second dataset can not be cast to array because of its shape - be=be - ) - np.testing.assert_allclose(deprecated_dists, dists) - assert backend.belongs_to_backend(deprecated_dists) - def test_ctw(): for be in backends: @@ -431,97 +387,90 @@ def test_masks(): for array_type in array_types: backend = instantiate_backend(be) - for sakoe_chiba_mask in [tslearn.metrics.sakoe_chiba_mask, - tslearn.metrics.dtw_variants.sakoe_chiba_mask]: + sk_mask = tslearn.metrics.sakoe_chiba_mask(4, 4, 1, be=be) + reference_mask = np.array( + [ + [True, True, False, False], + [True, True, True, False], + [False, True, True, True], + [False, False, True, True], + ] + ) + np.testing.assert_allclose(sk_mask, reference_mask) + assert backend.belongs_to_backend(sk_mask) + + sk_mask = tslearn.metrics.sakoe_chiba_mask(7, 3, 1, be=be) + reference_mask = np.array( + [ + [True, True, False], + [True, True, True], + [True, True, True], + [True, True, True], + [True, True, True], + [True, True, True], + [False, True, True], + ] + ) + np.testing.assert_allclose(sk_mask, reference_mask) + assert backend.belongs_to_backend(sk_mask) + + i_mask = tslearn.metrics.itakura_mask(6, 6, be=be) + reference_mask = np.array( + [ + [True, False, False, False, False, False], + [False, True, True, False, False, False], + [False, True, True, True, False, False], + [False, False, True, True, True, False], + [False, False, False, True, True, False], + [False, False, False, False, False, True], + ] + ) + np.testing.assert_allclose(i_mask, reference_mask) + assert backend.belongs_to_backend(i_mask) - sk_mask = sakoe_chiba_mask(4, 4, 1, be=be) - reference_mask = np.array( - [ - [True, True, False, False], - [True, True, True, False], - [False, True, True, True], - [False, False, True, True], - ] - ) - np.testing.assert_allclose(sk_mask, reference_mask) - assert backend.belongs_to_backend(sk_mask) + backend = instantiate_backend(be, array_type) - sk_mask = sakoe_chiba_mask(7, 3, 1, be=be) - reference_mask = np.array( - [ - [True, True, False], - [True, True, True], - [True, True, True], - [True, True, True], - [True, True, True], - [True, True, True], - [False, True, True], - ] - ) - np.testing.assert_allclose(sk_mask, reference_mask) - assert backend.belongs_to_backend(sk_mask) + # Test masks for different combinations of global_constraints / + # sakoe_chiba_radius / itakura_max_slope + sz = 10 + mask_no_constraint = tslearn.metrics.compute_mask(sz, sz, be=be) + np.testing.assert_array_equal(mask_no_constraint, np.full((sz, sz), True)) - for itakura_mask in [tslearn.metrics.itakura_mask, - tslearn.metrics.dtw_variants.itakura_mask]: - i_mask = itakura_mask(6, 6, be=be) - reference_mask = np.array( - [ - [True, False, False, False, False, False], - [False, True, True, False, False, False], - [False, True, True, True, False, False], - [False, False, True, True, True, False], - [False, False, False, True, True, False], - [False, False, False, False, False, True], - ] - ) - np.testing.assert_allclose(i_mask, reference_mask) - assert backend.belongs_to_backend(i_mask) - - for compute_mask in [tslearn.metrics.compute_mask, - tslearn.metrics.dtw_variants.compute_mask]: - backend = instantiate_backend(be, array_type) - - # Test masks for different combinations of global_constraints / - # sakoe_chiba_radius / itakura_max_slope - sz = 10 - mask_no_constraint = compute_mask(sz, sz, be=be) - np.testing.assert_array_equal(mask_no_constraint, np.full((sz, sz), True)) - - ts0 = cast(np.empty((sz, 1)), array_type) - ts1 = cast(np.empty((sz, 1)), array_type) - mask_no_constraint = compute_mask( - ts0, ts1, global_constraint=0, be=be - ) - np.testing.assert_array_equal(mask_no_constraint, np.full((sz, sz), True)) - assert backend.belongs_to_backend(mask_no_constraint) + ts0 = cast(np.empty((sz, 1)), array_type) + ts1 = cast(np.empty((sz, 1)), array_type) + mask_no_constraint = tslearn.metrics.compute_mask( + ts0, ts1, global_constraint=0, be=be + ) + np.testing.assert_array_equal(mask_no_constraint, np.full((sz, sz), True)) + assert backend.belongs_to_backend(mask_no_constraint) - mask_itakura = compute_mask( - ts0, ts1, global_constraint=1, be=be - ) - mask_itakura_bis = compute_mask( - ts0, ts1, itakura_max_slope=2.0, be=be - ) - np.testing.assert_allclose(mask_itakura, mask_itakura_bis) - assert backend.belongs_to_backend(mask_itakura) + mask_itakura = tslearn.metrics.compute_mask( + ts0, ts1, global_constraint=1, be=be + ) + mask_itakura_bis = tslearn.metrics.compute_mask( + ts0, ts1, itakura_max_slope=2.0, be=be + ) + np.testing.assert_allclose(mask_itakura, mask_itakura_bis) + assert backend.belongs_to_backend(mask_itakura) - mask_sakoe = compute_mask( - ts0, ts1, global_constraint=2, be=be - ) - mask_sakoe_bis = compute_mask( - ts0, ts1, sakoe_chiba_radius=1, be=be - ) - np.testing.assert_allclose(mask_sakoe, mask_sakoe_bis) - assert backend.belongs_to_backend(mask_sakoe) - - np.testing.assert_raises( - RuntimeWarning, - compute_mask, - ts0, - ts1, - sakoe_chiba_radius=1, - itakura_max_slope=2.0, - be=be, - ) + mask_sakoe = tslearn.metrics.compute_mask( + ts0, ts1, global_constraint=2, be=be + ) + mask_sakoe_bis = tslearn.metrics.compute_mask( + ts0, ts1, sakoe_chiba_radius=1, be=be + ) + np.testing.assert_allclose(mask_sakoe, mask_sakoe_bis) + assert backend.belongs_to_backend(mask_sakoe) + + np.testing.assert_raises( + RuntimeWarning, + tslearn.metrics.compute_mask, + ts0, + ts1, + sakoe_chiba_radius=1, + itakura_max_slope=2.0, + be=be, + ) # Tests for estimators that can set masks through metric_params n, sz, d = 15, 10, 3 @@ -580,34 +529,14 @@ def test_gak(): be=be ) np.testing.assert_allclose(s, 2.0, atol=1e-5) - with pytest.deprecated_call(): - backend.testing.assert_allclose( - s, - tslearn.metrics.softdtw_variants.sigma_gak( - dataset, - n_samples=200, - random_state=0, - be=be - ) - ) # GAK s1, s2 = cast([1, 2, 2, 3], array_type), cast([1.0, 2.0, 3.0, 4.0], array_type) g = tslearn.metrics.gak(s1, s2, sigma=2, be=be) np.testing.assert_allclose(g, 0.656297, atol=1e-5) - with pytest.deprecated_call(): - backend.testing.assert_allclose( - g, - tslearn.metrics.softdtw_variants.gak(s1, s2, sigma=2, be=be) - ) g = tslearn.metrics.unnormalized_gak(s1, s2, be=be) np.testing.assert_allclose(g, 3.745675, atol=1e-5) - with pytest.deprecated_call(): - backend.testing.assert_allclose( - g, - tslearn.metrics.softdtw_variants.unnormalized_gak(s1, s2, be=be) - ) g = tslearn.metrics.cdist_gak( dataset, sigma=2.0, be=be @@ -616,11 +545,7 @@ def test_gak(): g, np.array([[1.0, 0.656297], [0.656297, 1.0]]), atol=1e-5 ) assert backend.belongs_to_backend(g) - with pytest.deprecated_call(): - backend.testing.assert_allclose( - g, - tslearn.metrics.softdtw_variants.cdist_gak(dataset, sigma=2, be=be) - ) + g = tslearn.metrics.cdist_gak( [[1, 2, 2], [1.0, 2.0, 3.0, 4.0]], # Can not be cast to array because of its shape cast([[1, 2, 2, 3], [1.0, 2.0, 3.0, 4.0]], array_type), @@ -744,13 +669,6 @@ def test_dtw_path_from_metric(): # Use dtw_path as a reference path_ref, dist_ref = tslearn.metrics.dtw_path(s1, s2, be=be) - with pytest.deprecated_call(): - deprecated_path, deprecated_dist = tslearn.metrics.dtw_variants.dtw_path_from_metric( - s1, s2, metric="sqeuclidean", be=be - ) - np.testing.assert_equal(path_ref, deprecated_path) - np.testing.assert_allclose(backend.sqrt(deprecated_dist), dist_ref) - # Test of using a scipy distance function path, dist = tslearn.metrics.dtw_path_from_metric( s1, s2, metric="sqeuclidean", be=be diff --git a/tslearn/metrics/dtw_variants.py b/tslearn/metrics/dtw_variants.py index 32c1b8d2..b73c8b12 100644 --- a/tslearn/metrics/dtw_variants.py +++ b/tslearn/metrics/dtw_variants.py @@ -1,14 +1,11 @@ -import warnings - import numpy from numba import njit, prange from tslearn.backend import instantiate_backend -from tslearn.utils import to_time_series, to_time_series_dataset +from tslearn.utils import to_time_series -from .utils import _cdist_generic -from . _masks import compute_mask as compute_mask_ +from ._masks import compute_mask as compute_mask_ __author__ = "Romain Tavenard romain.tavenard[at]univ-rennes2.fr" @@ -55,779 +52,19 @@ def _local_squared_dist(x, y, be=None): If `be` is `None`, the backend is determined by the input arrays. See our :ref:`dedicated user-guide page ` for more information. - Returns - ------- - dist : float - Squared distance between x and y. - """ - be = instantiate_backend(be, x, y) - x = be.array(x) - y = be.array(y) - dist = 0.0 - for di in range(be.shape(x)[0]): - diff = x[di] - y[di] - dist += diff * diff - return dist - - -@njit() -def njit_accumulated_matrix(s1, s2, mask): - """Compute the accumulated cost matrix score between two time series. - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) - First time series. - s2 : array-like, shape=(sz2, d) - Second time series. - mask : array-like, shape=(sz1, sz2) - Mask. Unconsidered cells must have False values. - - Returns - ------- - mat : array-like, shape=(sz1, sz2) - Accumulated cost matrix. - """ - l1 = s1.shape[0] - l2 = s2.shape[0] - cum_sum = numpy.full((l1 + 1, l2 + 1), numpy.inf) - cum_sum[0, 0] = 0.0 - - for i in range(l1): - for j in range(l2): - if mask[i, j]: - cum_sum[i + 1, j + 1] = _njit_local_squared_dist(s1[i], s2[j]) - cum_sum[i + 1, j + 1] += min( - cum_sum[i, j + 1], cum_sum[i + 1, j], cum_sum[i, j] - ) - return cum_sum[1:, 1:] - - -def accumulated_matrix(s1, s2, mask, be=None): - """Compute the accumulated cost matrix score between two time series. - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) - First time series. - s2 : array-like, shape=(sz2, d) - Second time series. - mask : array-like, shape=(sz1, sz2) - Mask. Unconsidered cells must have False values. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - mat : array-like, shape=(sz1, sz2) - Accumulated cost matrix. - """ - warnings.warn( - "This method is deprecated, use tslearn.metrics.accumulated_matrix instead.", - DeprecationWarning - ) - - be = instantiate_backend(be, s1, s2) - s1 = be.array(s1) - s2 = be.array(s2) - l1 = be.shape(s1)[0] - l2 = be.shape(s2)[0] - cum_sum = be.full((l1 + 1, l2 + 1), be.inf) - cum_sum[0, 0] = 0.0 - - for i in range(l1): - for j in range(l2): - if mask[i, j]: - cum_sum[i + 1, j + 1] = _local_squared_dist(s1[i], s2[j], be=be) - cum_sum[i + 1, j + 1] += min( - cum_sum[i, j + 1], cum_sum[i + 1, j], cum_sum[i, j] - ) - return cum_sum[1:, 1:] - - -@njit(nogil=True) -def _njit_dtw(s1, s2, mask): - """Compute the dynamic time warping score between two time series. - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) - First time series. - s2 : array-like, shape=(sz2, d) - Second time series. - mask : array-like, shape=(sz1, sz2) - Mask. Unconsidered cells must have False values. - - Returns - ------- - dtw_score : float - Dynamic Time Warping score between both time series. - - """ - cum_sum = njit_accumulated_matrix(s1, s2, mask) - return numpy.sqrt(cum_sum[-1, -1]) - - -def _dtw(s1, s2, mask, be=None): - """Compute the dynamic time warping score between two time series. - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) - First time series. - s2 : array-like, shape=(sz2, d) - Second time series. - mask : array-like, shape=(sz1, sz2) - Mask. Unconsidered cells must have False values. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - dtw_score : float - Dynamic Time Warping score between both time series. - - """ - be = instantiate_backend(be, s1, s2) - s1 = be.array(s1) - s2 = be.array(s2) - cum_sum = accumulated_matrix(s1, s2, mask, be=be) - return be.sqrt(cum_sum[-1, -1]) - - -@njit() -def _njit_return_path(acc_cost_mat): - """Return path from accumulated cost matrix. - - Parameters - ---------- - acc_cost_mat : array-like, shape=(sz1, sz2) - Accumulated cost matrix. - - Returns - ------- - path : list of integer pairs - Matching path represented as a list of index pairs. In each pair, the - first index corresponds to a first time series s1 and the second one - corresponds to a second time series s2. - """ - sz1, sz2 = acc_cost_mat.shape - path = [(sz1 - 1, sz2 - 1)] - while path[-1] != (0, 0): - i, j = path[-1] - if i == 0: - path.append((0, j - 1)) - elif j == 0: - path.append((i - 1, 0)) - else: - arr = numpy.array( - [ - acc_cost_mat[i - 1][j - 1], - acc_cost_mat[i - 1][j], - acc_cost_mat[i][j - 1], - ] - ) - argmin = numpy.argmin(arr) - if argmin == 0: - path.append((i - 1, j - 1)) - elif argmin == 1: - path.append((i - 1, j)) - else: - path.append((i, j - 1)) - return path[::-1] - - -def _return_path(acc_cost_mat, be=None): - """Return path from accumulated cost matrix. - - Parameters - ---------- - acc_cost_mat : array-like, shape=(sz1, sz2) - Accumulated cost matrix. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - path : list of integer pairs - Matching path represented as a list of index pairs. In each pair, the - first index corresponds to a first time series s1 and the second one - corresponds to a second time series s2. - """ - be = instantiate_backend(be, acc_cost_mat) - sz1, sz2 = be.shape(acc_cost_mat) - path = [(sz1 - 1, sz2 - 1)] - while path[-1] != (0, 0): - i, j = path[-1] - if i == 0: - path.append((0, j - 1)) - elif j == 0: - path.append((i - 1, 0)) - else: - arr = be.array( - [ - acc_cost_mat[i - 1][j - 1], - acc_cost_mat[i - 1][j], - acc_cost_mat[i][j - 1], - ] - ) - argmin = be.argmin(arr) - if argmin == 0: - path.append((i - 1, j - 1)) - elif argmin == 1: - path.append((i - 1, j)) - else: - path.append((i, j - 1)) - return path[::-1] - - -def dtw_path( - s1, - s2, - global_constraint=None, - sakoe_chiba_radius=None, - itakura_max_slope=None, - be=None, -): - r"""Compute Dynamic Time Warping (DTW) similarity measure between - (possibly multidimensional) time series and return both the path and the - similarity. - - DTW is computed as the Euclidean distance between aligned time series, - i.e., if :math:`\pi` is the alignment path: - - .. math:: - - DTW(X, Y) = \sqrt{\sum_{(i, j) \in \pi} (X_{i} - Y_{j})^2} - - It is not required that both time series share the same size, but they must - be the same dimension. DTW was originally presented in [1]_ and is - discussed in more details in our :ref:`dedicated user-guide page `. - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) or (sz1,) - A time series. If shape is (sz1,), the time series is assumed to be univariate. - s2 : array-like, shape=(sz2, d) or (sz2,) - Another time series. If shape is (sz2,), the time series is assumed to be univariate. - global_constraint : {"itakura", "sakoe_chiba"} or None (default: None) - Global constraint to restrict admissible paths for DTW. - sakoe_chiba_radius : int or None (default: None) - Radius to be used for Sakoe-Chiba band global constraint. - The Sakoe-Chiba radius corresponds to the parameter :math:`\delta` mentioned in [1]_, - it controls how far in time we can go in order to match a given - point from one time series to a point in another time series. - If None and `global_constraint` is set to "sakoe_chiba", a radius of - 1 is used. - If both `sakoe_chiba_radius` and `itakura_max_slope` are set, - `global_constraint` is used to infer which constraint to use among the - two. In this case, if `global_constraint` corresponds to no global - constraint, a `RuntimeWarning` is raised and no global constraint is - used. - itakura_max_slope : float or None (default: None) - Maximum slope for the Itakura parallelogram constraint. - If None and `global_constraint` is set to "itakura", a maximum slope - of 2. is used. - If both `sakoe_chiba_radius` and `itakura_max_slope` are set, - `global_constraint` is used to infer which constraint to use among the - two. In this case, if `global_constraint` corresponds to no global - constraint, a `RuntimeWarning` is raised and no global constraint is - used. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - list of integer pairs - Matching path represented as a list of index pairs. In each pair, the - first index corresponds to s1 and the second one corresponds to s2. - - float - Similarity score - - Examples - -------- - >>> path, dist = dtw_path([1, 2, 3], [1., 2., 2., 3.]) - >>> path - [(0, 0), (1, 1), (1, 2), (2, 3)] - >>> float(dist) - 0.0 - >>> float(dtw_path([1, 2, 3], [1., 2., 2., 3., 4.])[1]) - 1.0 - - See Also - -------- - dtw : Get only the similarity score for DTW - cdist_dtw : Cross similarity matrix between time series datasets - dtw_path_from_metric : Compute a DTW using a user-defined distance metric - - References - ---------- - .. [1] H. Sakoe, S. Chiba, "Dynamic programming algorithm optimization for - spoken word recognition," IEEE Transactions on Acoustics, Speech and - Signal Processing, vol. 26(1), pp. 43--49, 1978. - - """ - warnings.warn( - "This method is deprecated, use tslearn.metrics.dtw_path instead.", - DeprecationWarning - ) - - be = instantiate_backend(be, s1, s2) - s1 = to_time_series(s1, remove_nans=True, be=be) - s2 = to_time_series(s2, remove_nans=True, be=be) - - if len(s1) == 0 or len(s2) == 0: - raise ValueError( - "One of the input time series contains only nans or has zero length." - ) - - if be.shape(s1)[1] != be.shape(s2)[1]: - raise ValueError("All input time series must have the same feature size.") - - mask = compute_mask_( - s1, - s2, - GLOBAL_CONSTRAINT_CODE[global_constraint], - sakoe_chiba_radius, - itakura_max_slope, - be=be, - ) - if be.is_numpy: - acc_cost_mat = njit_accumulated_matrix(s1, s2, mask=mask) - path = _njit_return_path(acc_cost_mat) - else: - acc_cost_mat = accumulated_matrix(s1, s2, mask=mask, be=be) - path = _return_path(acc_cost_mat, be=be) - return path, be.sqrt(acc_cost_mat[-1, -1]) - - -@njit() -def njit_accumulated_matrix_from_dist_matrix(dist_matrix, mask): - """Compute the accumulated cost matrix score between two time series using - a precomputed distance matrix. - - Parameters - ---------- - dist_matrix : array-like, shape=(sz1, sz2) - Array containing the pairwise distances. - mask : array-like, shape=(sz1, sz2) - Mask. Unconsidered cells must have False values. - - Returns - ------- - mat : array-like, shape=(sz1, sz2) - Accumulated cost matrix. - """ - l1, l2 = dist_matrix.shape - cum_sum = numpy.full((l1 + 1, l2 + 1), numpy.inf) - cum_sum[0, 0] = 0.0 - - for i in prange(l1): - for j in prange(l2): - if mask[i, j]: - cum_sum[i + 1, j + 1] = dist_matrix[i, j] - cum_sum[i + 1, j + 1] += min( - cum_sum[i, j + 1], cum_sum[i + 1, j], cum_sum[i, j] - ) - return cum_sum[1:, 1:] - - -def accumulated_matrix_from_dist_matrix(dist_matrix, mask, be=None): - """Compute the accumulated cost matrix score between two time series using - a precomputed distance matrix. - - Parameters - ---------- - dist_matrix : array-like, shape=(sz1, sz2) - Array containing the pairwise distances. - mask : array-like, shape=(sz1, sz2) - Mask. Unconsidered cells must have False values. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - mat : array-like, shape=(sz1, sz2) - Accumulated cost matrix. - """ - be = instantiate_backend(be, dist_matrix) - dist_matrix = be.array(dist_matrix) - l1, l2 = be.shape(dist_matrix) - cum_sum = be.full((l1 + 1, l2 + 1), be.inf) - cum_sum[0, 0] = 0.0 - - for i in range(l1): - for j in range(l2): - if mask[i, j]: - cum_sum[i + 1, j + 1] = dist_matrix[i, j] - cum_sum[i + 1, j + 1] += min( - cum_sum[i, j + 1], cum_sum[i + 1, j], cum_sum[i, j] - ) - return cum_sum[1:, 1:] - - -def dtw_path_from_metric( - s1, - s2=None, - metric="euclidean", - global_constraint=None, - sakoe_chiba_radius=None, - itakura_max_slope=None, - be=None, - **kwds -): - r"""Compute Dynamic Time Warping (DTW) similarity measure between - (possibly multidimensional) time series using a distance metric defined by - the user and return both the path and the similarity. - - Similarity is computed as the cumulative cost along the aligned time - series. - - It is not required that both time series share the same size, but they must - be the same dimension. DTW was originally presented in [1]_. - - Valid values for metric are the same as for scikit-learn - `pairwise_distances`_ function i.e. a string (e.g. "euclidean", - "sqeuclidean", "hamming") or a function that is used to compute the - pairwise distances. See `scikit`_ and `scipy`_ documentations for more - information about the available metrics. - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) or (sz1,) if metric!="precomputed", (sz1, sz2) otherwise - A time series or an array of pairwise distances between samples. - If shape is (sz1,), the time series is assumed to be univariate. - - s2 : array-like, shape=(sz2, d) or (sz2,), optional (default: None) - A second time series, only allowed if metric != "precomputed". - If shape is (sz2,), the time series is assumed to be univariate. - - metric : string or callable (default: "euclidean") - Function used to compute the pairwise distances between each points of - `s1` and `s2`. - - If metric is "precomputed", `s1` is assumed to be a distance matrix. - - If metric is an other string, it must be one of the options compatible - with sklearn.metrics.pairwise_distances. - - Alternatively, if metric is a callable function, it is called on pairs - of rows of `s1` and `s2`. The callable should take two 1 dimensional - arrays as input and return a value indicating the distance between - them. - - global_constraint : {"itakura", "sakoe_chiba"} or None (default: None) - Global constraint to restrict admissible paths for DTW. - - sakoe_chiba_radius : int or None (default: None) - Radius to be used for Sakoe-Chiba band global constraint. - The Sakoe-Chiba radius corresponds to the parameter :math:`\delta` mentioned in [1]_, - it controls how far in time we can go in order to match a given - point from one time series to a point in another time series. - If None and `global_constraint` is set to "sakoe_chiba", a radius of - 1 is used. - If both `sakoe_chiba_radius` and `itakura_max_slope` are set, - `global_constraint` is used to infer which constraint to use among the - two. In this case, if `global_constraint` corresponds to no global - constraint, a `RuntimeWarning` is raised and no global constraint is - used. - - itakura_max_slope : float or None (default: None) - Maximum slope for the Itakura parallelogram constraint. - If None and `global_constraint` is set to "itakura", a maximum slope - of 2. is used. - If both `sakoe_chiba_radius` and `itakura_max_slope` are set, - `global_constraint` is used to infer which constraint to use among the - two. In this case, if `global_constraint` corresponds to no global - constraint, a `RuntimeWarning` is raised and no global constraint is - used. - - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - **kwds - Additional arguments to pass to sklearn pairwise_distances to compute - the pairwise distances. - - Returns - ------- - list of integer pairs - Matching path represented as a list of index pairs. In each pair, the - first index corresponds to s1 and the second one corresponds to s2. - - float - Similarity score (sum of metric along the wrapped time series). - - Examples - -------- - Lets create 2 numpy arrays to wrap: - - >>> import numpy as np - >>> rng = np.random.RandomState(0) - >>> s1, s2 = rng.rand(5, 2), rng.rand(6, 2) - - The wrapping can be done by passing a string indicating the metric to pass - to scikit-learn pairwise_distances: - - >>> x, y = dtw_path_from_metric(s1, s2, - ... metric="sqeuclidean") # doctest: +ELLIPSIS - >>> x, float(y) - ([(0, 0), (0, 1), (1, 2), (2, 3), (3, 4), (4, 5)], 1.117...) - - Or by defining a custom distance function: - - >>> sqeuclidean = lambda x, y: np.sum((x-y)**2) - >>> x, y = dtw_path_from_metric(s1, s2, metric=sqeuclidean) # doctest: +ELLIPSIS - >>> x, float(y) - ([(0, 0), (0, 1), (1, 2), (2, 3), (3, 4), (4, 5)], 1.117...) - - Or by using a precomputed distance matrix as input: - - >>> from sklearn.metrics.pairwise import pairwise_distances - >>> dist_matrix = pairwise_distances(s1, s2, metric="sqeuclidean") - >>> x, y = dtw_path_from_metric(dist_matrix, - ... metric="precomputed") # doctest: +ELLIPSIS - >>> x, float(y) - ([(0, 0), (0, 1), (1, 2), (2, 3), (3, 4), (4, 5)], 1.117...) - - Notes - ----- - By using a squared euclidean distance metric as shown above, the output - path is the same as the one obtained by using dtw_path but the similarity - score is the sum of squared distances instead of the euclidean distance. - - See Also - -------- - dtw_path : Get both the matching path and the similarity score for DTW - - References - ---------- - .. [1] H. Sakoe, S. Chiba, "Dynamic programming algorithm optimization for - spoken word recognition," IEEE Transactions on Acoustics, Speech and - Signal Processing, vol. 26(1), pp. 43--49, 1978. - - .. _pairwise_distances: https://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise_distances.html - - .. _scikit: https://scikit-learn.org/stable/modules/metrics.html - - .. _scipy: https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.pdist.html - - """ # noqa: E501 - warnings.warn( - "This method is deprecated, use tslearn.metrics.dtw_path_from_metric instead.", - DeprecationWarning - ) - - be = instantiate_backend(be, s1, s2) - if metric == "precomputed": # Pairwise distance given as input - s1 = be.array(s1) - sz1, sz2 = be.shape(s1) - mask = compute_mask_( - sz1, - sz2, - GLOBAL_CONSTRAINT_CODE[global_constraint], - sakoe_chiba_radius, - itakura_max_slope, - be=be, - ) - dist_mat = s1 - else: - s1 = to_time_series(s1, remove_nans=True, be=be) - s2 = to_time_series(s2, remove_nans=True, be=be) - mask = compute_mask_( - s1, - s2, - GLOBAL_CONSTRAINT_CODE[global_constraint], - sakoe_chiba_radius, - itakura_max_slope, - be=be, - ) - dist_mat = be.pairwise_distances(s1, s2, metric=metric, **kwds) - - if be.is_numpy: - acc_cost_mat = njit_accumulated_matrix_from_dist_matrix(dist_mat, mask) - path = _njit_return_path(acc_cost_mat) - else: - dist_mat = be.array(dist_mat) - acc_cost_mat = accumulated_matrix_from_dist_matrix(dist_mat, mask, be=be) - path = _return_path(acc_cost_mat, be=be) - return path, acc_cost_mat[-1, -1] - - -def dtw( - s1, - s2, - global_constraint=None, - sakoe_chiba_radius=None, - itakura_max_slope=None, - be=None, -): - r"""Compute Dynamic Time Warping (DTW) similarity measure between - (possibly multidimensional) time series and return it. - - DTW is computed as the Euclidean distance between aligned time series, - i.e., if :math:`\pi` is the optimal alignment path: - - .. math:: - - DTW(X, Y) = \sqrt{\sum_{(i, j) \in \pi} \|X_{i} - Y_{j}\|^2} - - Note that this formula is still valid for the multivariate case. - - It is not required that both time series share the same size, but they must - be the same dimension. DTW was originally presented in [1]_ and is - discussed in more details in our :ref:`dedicated user-guide page `. - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) or (sz1,) - A time series. If shape is (sz1,), the time series is assumed to be univariate. - - s2 : array-like, shape=(sz2, d) or (sz2,) - Another time series. If shape is (sz2,), the time series is assumed to be univariate. - - global_constraint : {"itakura", "sakoe_chiba"} or None (default: None) - Global constraint to restrict admissible paths for DTW. - - sakoe_chiba_radius : int or None (default: None) - Radius to be used for Sakoe-Chiba band global constraint. - The Sakoe-Chiba radius corresponds to the parameter :math:`\delta` mentioned in [1]_, - it controls how far in time we can go in order to match a given - point from one time series to a point in another time series. - If None and `global_constraint` is set to "sakoe_chiba", a radius of - 1 is used. - If both `sakoe_chiba_radius` and `itakura_max_slope` are set, - `global_constraint` is used to infer which constraint to use among the - two. In this case, if `global_constraint` corresponds to no global - constraint, a `RuntimeWarning` is raised and no global constraint is - used. - - itakura_max_slope : float or None (default: None) - Maximum slope for the Itakura parallelogram constraint. - If None and `global_constraint` is set to "itakura", a maximum slope - of 2. is used. - If both `sakoe_chiba_radius` and `itakura_max_slope` are set, - `global_constraint` is used to infer which constraint to use among the - two. In this case, if `global_constraint` corresponds to no global - constraint, a `RuntimeWarning` is raised and no global constraint is - used. - - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - float - Similarity score - - Examples - -------- - >>> dtw([1, 2, 3], [1., 2., 2., 3.]) - 0.0 - >>> dtw([1, 2, 3], [1., 2., 2., 3., 4.]) - 1.0 - - The PyTorch backend can be used to compute gradients: - - >>> import torch - >>> s1 = torch.tensor([[1.0], [2.0], [3.0]], requires_grad=True) - >>> s2 = torch.tensor([[3.0], [4.0], [-3.0]]) - >>> sim = dtw(s1, s2, be="pytorch") - >>> print(sim) - tensor(6.4807, grad_fn=) - >>> sim.backward() - >>> print(s1.grad) - tensor([[-0.3086], - [-0.1543], - [ 0.7715]]) - - >>> s1_2d = torch.tensor([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0]], requires_grad=True) - >>> s2_2d = torch.tensor([[3.0, 3.0], [4.0, 4.0], [-3.0, -3.0]]) - >>> sim = dtw(s1_2d, s2_2d, be="pytorch") - >>> print(sim) - tensor(9.1652, grad_fn=) - >>> sim.backward() - >>> print(s1_2d.grad) - tensor([[-0.2182, -0.2182], - [-0.1091, -0.1091], - [ 0.5455, 0.5455]]) - - See Also - -------- - dtw_path : Get both the matching path and the similarity score for DTW - cdist_dtw : Cross similarity matrix between time series datasets - - References - ---------- - .. [1] H. Sakoe, S. Chiba, "Dynamic programming algorithm optimization for - spoken word recognition," IEEE Transactions on Acoustics, Speech and - Signal Processing, vol. 26(1), pp. 43--49, 1978. - - """ # noqa: E501 - warnings.warn( - "This method is deprecated, use tslearn.metrics.dtw instead.", - DeprecationWarning - ) - - be = instantiate_backend(be, s1, s2) - s1 = to_time_series(s1, remove_nans=True, be=be) - s2 = to_time_series(s2, remove_nans=True, be=be) - - if len(s1) == 0 or len(s2) == 0: - raise ValueError( - "One of the input time series contains only nans or has zero length." - ) - - if be.shape(s1)[1] != be.shape(s2)[1]: - raise ValueError("All input time series must have the same feature size.") - - mask = compute_mask_( - s1, - s2, - GLOBAL_CONSTRAINT_CODE[global_constraint], - sakoe_chiba_radius=sakoe_chiba_radius, - itakura_max_slope=itakura_max_slope, - be=be, - ) - if be.is_numpy: - return _njit_dtw(s1, s2, mask=mask) - return _dtw(s1, s2, mask=mask, be=be) + Returns + ------- + dist : float + Squared distance between x and y. + """ + be = instantiate_backend(be, x, y) + x = be.array(x) + y = be.array(y) + dist = 0.0 + for di in range(be.shape(x)[0]): + diff = x[di] - y[di] + dist += diff * diff + return dist def _max_steps(i, j, max_length, length_1, length_2): @@ -1483,528 +720,6 @@ def dtw_subsequence_path(subseq, longseq, be=None): return path, be.sqrt(acc_cost_mat[-1, :][global_optimal_match]) -@njit() -def njit_sakoe_chiba_mask(sz1, sz2, radius=1): - """Compute the Sakoe-Chiba mask. - - Parameters - ---------- - sz1 : int - The size of the first time series - sz2 : int - The size of the second time series. - radius : int - The radius of the band. - - Returns - ------- - mask : array-like, shape=(sz1, sz2) - Sakoe-Chiba mask. - - Examples - -------- - >>> njit_sakoe_chiba_mask(4, 4, 1) - array([[ True, True, False, False], - [ True, True, True, False], - [False, True, True, True], - [False, False, True, True]]) - >>> njit_sakoe_chiba_mask(7, 3, 1) - array([[ True, True, False], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [False, True, True]]) - """ - mask = numpy.full((sz1, sz2), False) - if sz1 > sz2: - width = sz1 - sz2 + radius - for i in prange(sz2): - lower = max(0, i - radius) - upper = min(sz1, i + width) + 1 - mask[lower:upper, i] = True - else: - width = sz2 - sz1 + radius - for i in prange(sz1): - lower = max(0, i - radius) - upper = min(sz2, i + width) + 1 - mask[i, lower:upper] = True - return mask - - -def sakoe_chiba_mask(sz1, sz2, radius=1, be=None): - """Compute the Sakoe-Chiba mask. - - Parameters - ---------- - sz1 : int - The size of the first time series - sz2 : int - The size of the second time series. - radius : int - The radius of the band. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - mask : array-like, shape=(sz1, sz2) - Sakoe-Chiba mask. - - Examples - -------- - >>> sakoe_chiba_mask(4, 4, 1) - array([[ True, True, False, False], - [ True, True, True, False], - [False, True, True, True], - [False, False, True, True]]) - >>> sakoe_chiba_mask(7, 3, 1) - array([[ True, True, False], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [False, True, True]]) - """ - warnings.warn( - "This method is deprecated, use tslearn.metrics.sakoe_chiba_mask instead.", - DeprecationWarning - ) - - be = instantiate_backend(be) - mask = be.full((sz1, sz2), False) - if sz1 > sz2: - width = sz1 - sz2 + radius - for i in range(sz2): - lower = max(0, i - radius) - upper = min(sz1, i + width) + 1 - mask[lower:upper, i] = True - else: - width = sz2 - sz1 + radius - for i in range(sz1): - lower = max(0, i - radius) - upper = min(sz2, i + width) + 1 - mask[i, lower:upper] = True - return mask - - -@njit() -def _njit_itakura_mask(sz1, sz2, max_slope=2.0): - """Compute the Itakura mask without checking that the constraints - are feasible. In most cases, you should use itakura_mask instead. - - Parameters - ---------- - sz1 : int - The size of the first time series - - sz2 : int - The size of the second time series. - - max_slope : float (default = 2) - The maximum slope of the parallelogram. - - Returns - ------- - mask : array-like, shape=(sz1, sz2) - Itakura mask. - """ - min_slope = 1 / float(max_slope) - max_slope *= float(sz1) / float(sz2) - min_slope *= float(sz1) / float(sz2) - - lower_bound = numpy.empty((2, sz2)) - lower_bound[0] = min_slope * numpy.arange(sz2) - lower_bound[1] = (sz1 - 1) - max_slope * (sz2 - 1) + max_slope * numpy.arange(sz2) - lower_bound_ = numpy.empty(sz2) - for i in prange(sz2): - lower_bound_[i] = max(round(lower_bound[0, i], 2), round(lower_bound[1, i], 2)) - lower_bound_ = numpy.ceil(lower_bound_) - - upper_bound = numpy.empty((2, sz2)) - upper_bound[0] = max_slope * numpy.arange(sz2) - upper_bound[1] = (sz1 - 1) - min_slope * (sz2 - 1) + min_slope * numpy.arange(sz2) - upper_bound_ = numpy.empty(sz2) - for i in prange(sz2): - upper_bound_[i] = min(round(upper_bound[0, i], 2), round(upper_bound[1, i], 2)) - upper_bound_ = numpy.floor(upper_bound_ + 1) - - mask = numpy.full((sz1, sz2), False) - for i in prange(sz2): - mask[int(lower_bound_[i]) : int(upper_bound_[i]), i] = True - return mask - - -def _itakura_mask(sz1, sz2, max_slope=2.0, be=None): - """Compute the Itakura mask without checking that the constraints - are feasible. In most cases, you should use itakura_mask instead. - - Parameters - ---------- - sz1 : int - The size of the first time series - sz2 : int - The size of the second time series. - max_slope : float (default = 2) - The maximum slope of the parallelogram. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - mask : array-like, shape=(sz1, sz2) - Itakura mask. - """ - min_slope = 1 / float(max_slope) - max_slope *= float(sz1) / float(sz2) - min_slope *= float(sz1) / float(sz2) - - lower_bound = be.empty((2, sz2)) - lower_bound[0] = min_slope * be.arange(sz2) - lower_bound[1] = (sz1 - 1) - max_slope * (sz2 - 1) + max_slope * be.arange(sz2) - lower_bound_ = be.empty(sz2) - for i in range(sz2): - lower_bound_[i] = max( - be.round(lower_bound[0, i], decimals=2), - be.round(lower_bound[1, i], decimals=2), - ) - lower_bound_ = be.ceil(lower_bound_) - - upper_bound = be.empty((2, sz2)) - upper_bound[0] = max_slope * be.arange(sz2) - upper_bound[1] = (sz1 - 1) - min_slope * (sz2 - 1) + min_slope * be.arange(sz2) - upper_bound_ = be.empty(sz2) - for i in range(sz2): - upper_bound_[i] = min( - be.round(upper_bound[0, i], decimals=2), - be.round(upper_bound[1, i], decimals=2), - ) - upper_bound_ = be.floor(upper_bound_ + 1) - - mask = be.full((sz1, sz2), False) - for i in range(sz2): - mask[int(lower_bound_[i]) : int(upper_bound_[i]), i] = True - return mask - - -def itakura_mask(sz1, sz2, max_slope=2.0, be=None): - """Compute the Itakura mask. - - Parameters - ---------- - sz1 : int - The size of the first time series - sz2 : int - The size of the second time series. - max_slope : float (default = 2) - The maximum slope of the parallelogram. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - mask : array-like, shape=(sz1, sz2) - Itakura mask. - - Examples - -------- - >>> itakura_mask(6, 6) - array([[ True, False, False, False, False, False], - [False, True, True, False, False, False], - [False, True, True, True, False, False], - [False, False, True, True, True, False], - [False, False, False, True, True, False], - [False, False, False, False, False, True]]) - """ - warnings.warn( - "This method is deprecated, use tslearn.metrics.itakura_mask instead.", - DeprecationWarning - ) - - be = instantiate_backend(be) - - if be.is_numpy: - mask = _njit_itakura_mask(sz1, sz2, max_slope=max_slope) - else: - mask = _itakura_mask(sz1, sz2, max_slope=max_slope, be=be) - - # Post-check - raise_warning = False - for i in range(sz1): - if not be.any(mask[i]): - raise_warning = True - break - if not raise_warning: - for j in range(sz2): - if not be.any(mask[:, j]): - raise_warning = True - break - if raise_warning: - warnings.warn( - "'itakura_max_slope' constraint is unfeasible " - "(ie. leads to no admissible path) for the " - "provided time series sizes", - RuntimeWarning, - ) - - return mask - - -def compute_mask( - s1, - s2, - global_constraint=0, - sakoe_chiba_radius=None, - itakura_max_slope=None, - be=None, -): - r"""Compute the mask (region constraint). - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) or (sz1,) - A time series or integer. - If shape is (sz1,), the time series is assumed to be univariate. - s2 : array-like, shape=(sz2, d) or (sz2,) - Another time series or integer. - If shape is (sz2,), the time series is assumed to be univariate. - global_constraint : {0, 1, 2} (default: 0) - Global constraint to restrict admissible paths for DTW: - - "itakura" if 1 - - "sakoe_chiba" if 2 - - no constraint otherwise - sakoe_chiba_radius : int or None (default: None) - Radius to be used for Sakoe-Chiba band global constraint. - The Sakoe-Chiba radius corresponds to the parameter :math:`\delta` mentioned in [1]_, - it controls how far in time we can go in order to match a given - point from one time series to a point in another time series. - If None and `global_constraint` is set to 2 (sakoe-chiba), a radius of - 1 is used. - If both `sakoe_chiba_radius` and `itakura_max_slope` are set, - `global_constraint` is used to infer which constraint to use among the - two. In this case, if `global_constraint` corresponds to no global - constraint, a `RuntimeWarning` is raised and no global constraint is - used. - itakura_max_slope : float or None (default: None) - Maximum slope for the Itakura parallelogram constraint. - If None and `global_constraint` is set to 1 (itakura), a maximum slope - of 2. is used. - If both `sakoe_chiba_radius` and `itakura_max_slope` are set, - `global_constraint` is used to infer which constraint to use among the - two. In this case, if `global_constraint` corresponds to no global - constraint, a `RuntimeWarning` is raised and no global constraint is - used. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - mask : array-like, shape=(sz1, sz2) - Constraint region. - """ - warnings.warn( - "This method is deprecated, use tslearn.metrics.compute_mask instead.", - DeprecationWarning - ) - - be = instantiate_backend(be, s1, s2) - # The output mask will be of shape (sz1, sz2) - if isinstance(s1, int) and isinstance(s2, int): - sz1, sz2 = s1, s2 - else: - s1 = be.array(s1) - s2 = be.array(s2) - sz1 = be.shape(s1)[0] - sz2 = be.shape(s2)[0] - if ( - global_constraint == 0 - and sakoe_chiba_radius is not None - and itakura_max_slope is not None - ): - raise RuntimeWarning( - "global_constraint is not set for DTW, but both " - "sakoe_chiba_radius and itakura_max_slope are " - "set, hence global_constraint cannot be inferred " - "and no global constraint will be used." - ) - if global_constraint == 2 or ( - global_constraint == 0 and sakoe_chiba_radius is not None - ): - if sakoe_chiba_radius is None: - sakoe_chiba_radius = 1 - if be.is_numpy: - mask = njit_sakoe_chiba_mask(sz1, sz2, radius=sakoe_chiba_radius) - else: - mask = sakoe_chiba_mask(sz1, sz2, radius=sakoe_chiba_radius, be=be) - elif global_constraint == 1 or ( - global_constraint == 0 and itakura_max_slope is not None - ): - if itakura_max_slope is None: - itakura_max_slope = 2.0 - mask = itakura_mask(sz1, sz2, max_slope=itakura_max_slope, be=be) - else: - mask = be.full((sz1, sz2), True) - return mask - - -def cdist_dtw( - dataset1, - dataset2=None, - global_constraint=None, - sakoe_chiba_radius=None, - itakura_max_slope=None, - n_jobs=None, - verbose=0, - be=None, -): - r"""Compute cross-similarity matrix using Dynamic Time Warping (DTW) - similarity measure. - - DTW is computed as the Euclidean distance between aligned time series, - i.e., if :math:`\pi` is the alignment path: - - .. math:: - - DTW(X, Y) = \sqrt{\sum_{(i, j) \in \pi} \|X_{i} - Y_{j}\|^2} - - Note that this formula is still valid for the multivariate case. - - It is not required that time series share the same size, but they - must be the same dimension. - DTW was originally presented in [1]_ and is - discussed in more details in our :ref:`dedicated user-guide page `. - - Parameters - ---------- - dataset1 : array-like, shape=(n_ts1, sz1, d) or (n_ts1, sz1) or (sz1,) - A dataset of time series. - If shape is (n_ts1, sz1), the dataset is composed of univariate time series. - If shape is (sz1,), the dataset is composed of a unique univariate time series. - - dataset2 : None or array-like, shape=(n_ts2, sz2, d) or (n_ts2, sz2) or (sz2,) (default: None) - Another dataset of time series. If `None`, self-similarity of - `dataset1` is returned. - If shape is (n_ts2, sz2), the dataset is composed of univariate time series. - If shape is (sz2,), the dataset is composed of a unique univariate time series. - - global_constraint : {"itakura", "sakoe_chiba"} or None (default: None) - Global constraint to restrict admissible paths for DTW. - - sakoe_chiba_radius : int or None (default: None) - Radius to be used for Sakoe-Chiba band global constraint. - The Sakoe-Chiba radius corresponds to the parameter :math:`\delta` mentioned in [1]_, - it controls how far in time we can go in order to match a given - point from one time series to a point in another time series. - If None and `global_constraint` is set to "sakoe_chiba", a radius of - 1 is used. - If both `sakoe_chiba_radius` and `itakura_max_slope` are set, - `global_constraint` is used to infer which constraint to use among the - two. In this case, if `global_constraint` corresponds to no global - constraint, a `RuntimeWarning` is raised and no global constraint is - used. - - itakura_max_slope : float or None (default: None) - Maximum slope for the Itakura parallelogram constraint. - If None and `global_constraint` is set to "itakura", a maximum slope - of 2. is used. - If both `sakoe_chiba_radius` and `itakura_max_slope` are set, - `global_constraint` is used to infer which constraint to use among the - two. In this case, if `global_constraint` corresponds to no global - constraint, a `RuntimeWarning` is raised and no global constraint is - used. - - n_jobs : int or None, optional (default=None) - The number of jobs to run in parallel. - ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. - ``-1`` means using all processors. See scikit-learns' - `Glossary `__ - for more details. - - verbose : int, optional (default=0) - The verbosity level: if non zero, progress messages are printed. - Above 50, the output is sent to stdout. - The frequency of the messages increases with the verbosity level. - If it more than 10, all iterations are reported. - `Glossary `__ - for more details. - - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - cdist : array-like, shape=(n_ts1, n_ts2) - Cross-similarity matrix. - - Examples - -------- - >>> cdist_dtw([[1, 2, 2, 3], [1., 2., 3., 4.]]) - array([[0., 1.], - [1., 0.]]) - >>> cdist_dtw([[1, 2, 2, 3], [1., 2., 3., 4.]], [[1, 2, 3], [2, 3, 4, 5]]) - array([[0. , 2.44948974], - [1. , 1.41421356]]) - - See Also - -------- - dtw : Get DTW similarity score - - References - ---------- - .. [1] H. Sakoe, S. Chiba, "Dynamic programming algorithm optimization for - spoken word recognition," IEEE Transactions on Acoustics, Speech and - Signal Processing, vol. 26(1), pp. 43--49, 1978. - """ # noqa: E501 - warnings.warn( - "This method is deprecated, use tslearn.metrics.cdist_dtw instead.", - DeprecationWarning - ) - - be = instantiate_backend(be, dataset1, dataset2) - dataset1 = to_time_series_dataset(dataset1, be=be) - if dataset2 is not None: - dataset2 = to_time_series_dataset(dataset2, be=be) - return _cdist_generic( - dist_fun=dtw, - dataset1=dataset1, - dataset2=dataset2, - n_jobs=n_jobs, - verbose=verbose, - compute_diagonal=False, - global_constraint=global_constraint, - sakoe_chiba_radius=sakoe_chiba_radius, - itakura_max_slope=itakura_max_slope, - be=be, - ) - - def lb_keogh(ts_query, ts_candidate=None, radius=1, envelope_candidate=None): r"""Compute LB_Keogh. diff --git a/tslearn/metrics/softdtw_variants.py b/tslearn/metrics/softdtw_variants.py index cc1986d1..0909b92b 100644 --- a/tslearn/metrics/softdtw_variants.py +++ b/tslearn/metrics/softdtw_variants.py @@ -1,11 +1,7 @@ """Soft-DTW toolbox""" import math -import warnings import numpy as np -from joblib import Parallel, delayed -from numba import njit -from sklearn.utils import check_random_state from tslearn.backend import instantiate_backend from tslearn.utils import ( @@ -18,6 +14,7 @@ _ts_size ) +from . import sigma_gak from ._dtw import dtw, dtw_path from .soft_dtw_fast import ( _jacobian_product_sq_euc, @@ -27,7 +24,6 @@ _soft_dtw, _soft_dtw_grad, ) -from .utils import _cdist_generic __author__ = "Romain Tavenard romain.tavenard[at]univ-rennes2.fr" @@ -36,441 +32,6 @@ VARIABLE_LENGTH_METRICS = ["dtw", "gak", "softdtw", "sax"] -def _gak(gram, be=None): - """Compute Global Alignment Kernel (GAK) between (possibly - multidimensional) time series and return it. - - Parameters - ---------- - gram : array-like, shape=(sz1, sz2) - Gram matrix. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - float - Kernel value - """ - be = instantiate_backend(be, gram) - sz1, sz2 = gram.shape - - cum_sum = be.zeros((sz1 + 1, sz2 + 1), dtype=gram.dtype) - cum_sum[0, 0] = 1.0 - - for i in range(sz1): - for j in range(sz2): - cum_sum[i + 1, j + 1] = ( - cum_sum[i, j + 1] + cum_sum[i + 1, j] + cum_sum[i, j] - ) * gram[i, j] - - return cum_sum[sz1, sz2] - - -@njit(nogil=True) -def _njit_gak(gram): - """Compute Global Alignment Kernel (GAK) between (possibly - multidimensional) time series and return it. - - Parameters - ---------- - gram : array-like, shape=(sz1, sz2) - Gram matrix. - - Returns - ------- - float - Kernel value - """ - sz1, sz2 = gram.shape - - cum_sum = np.zeros((sz1 + 1, sz2 + 1)) - cum_sum[0, 0] = 1.0 - - for i in range(sz1): - for j in range(sz2): - cum_sum[i + 1, j + 1] = ( - cum_sum[i, j + 1] + cum_sum[i + 1, j] + cum_sum[i, j] - ) * gram[i, j] - - return cum_sum[sz1, sz2] - - -def _gak_gram(s1, s2, sigma=1.0, be=None): - """Compute Global Alignment Kernel (GAK) Gram matrix between (possibly - multidimensional) time series and return it. - Raises ZeroDivisionError when sigma is close to zero. - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) - A time series. - s2 : array-like, shape=(sz2, d) - Another time series. - sigma : float (default 1.) - Bandwidth of the internal gaussian kernel used for GAK. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - gram : array-like, shape=(sz1, sz2) - Gram matrix. - """ - if math.isclose(sigma, 0.0): - raise ZeroDivisionError() - be = instantiate_backend(be, s1) - gram = -be.cdist(s1, s2, "sqeuclidean") / (2 * sigma**2) - gram -= be.log(2 - be.exp(gram)) - return be.exp(gram) - - -def unnormalized_gak(s1, s2, sigma=1.0, be=None): - r"""Compute Global Alignment Kernel (GAK) between (possibly - multidimensional) time series and return it. - - It is not required that both time series share the same size, but they must - be the same dimension. GAK was - originally presented in [1]_. - This is an unnormalized version. - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) or (sz1,) - A time series. - If shape is (sz1,), the time series is assumed to be univariate. - s2 : array-like, shape=(sz2, d) or (sz2,) - Another time series. - If shape is (sz2,), the time series is assumed to be univariate. - sigma : float (default 1.) - Bandwidth of the internal gaussian kernel used for GAK. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - float - Kernel value - - Examples - -------- - >>> unnormalized_gak([1, 2, 3], - ... [1., 2., 2., 3.], - ... sigma=2.) # doctest: +ELLIPSIS - 15.358... - >>> unnormalized_gak([1, 2, 3], - ... [1., 2., 2., 3., 4.]) # doctest: +ELLIPSIS - 3.166... - - See Also - -------- - gak : normalized version of GAK that ensures that k(x,x) = 1 - cdist_gak : Compute cross-similarity matrix using Global Alignment kernel - - References - ---------- - .. [1] M. Cuturi, "Fast global alignment kernels," - ICML 2011. - """ - warnings.warn( - "This method is deprecated, use tslearn.metrics.unnormalized_gak instead.", - DeprecationWarning - ) - - be = instantiate_backend(be, s1, s2) - s1 = to_time_series(s1, remove_nans=True, be=be) - s2 = to_time_series(s2, remove_nans=True, be=be) - - gram = _gak_gram(s1, s2, sigma=sigma, be=be) - - if be.is_numpy: - return _njit_gak(gram) - return _gak(gram, be=be) - - -def gak(s1, s2, sigma=1.0, be=None): # TODO: better doc (formula for the kernel) - r"""Compute Global Alignment Kernel (GAK) between (possibly - multidimensional) time series and return it. - - It is not required that both time series share the same size, but they must - be the same dimension. GAK was - originally presented in [1]_. - This is a normalized version that ensures that :math:`k(x,x)=1` for all - :math:`x` and :math:`k(x,y) \in [0, 1]` for all :math:`x, y`. - - Parameters - ---------- - s1 : array-like, shape=(sz1, d) or (sz1,) - A time series. - If shape is (sz1,), the time series is assumed to be univariate. - s2 : array-like, shape=(sz2, d) or (sz2,) - Another time series. - If shape is (sz2,), the time series is assumed to be univariate. - sigma : float (default 1.) - Bandwidth of the internal gaussian kernel used for GAK. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - float - Kernel value - - Examples - -------- - >>> float(gak([1, 2, 3], [1., 2., 2., 3.], sigma=2.)) # doctest: +ELLIPSIS - 0.839... - >>> float(gak([1, 2, 3], [1., 2., 2., 3., 4.])) # doctest: +ELLIPSIS - 0.273... - - See Also - -------- - cdist_gak : Compute cross-similarity matrix using Global Alignment kernel - - References - ---------- - .. [1] M. Cuturi, "Fast global alignment kernels," - ICML 2011. - """ - warnings.warn( - "This method is deprecated, use tslearn.metrics.gak instead.", - DeprecationWarning - ) - - be = instantiate_backend(be, s1, s2) - s1 = be.array(s1) - s2 = be.array(s2) - denom = be.sqrt( - unnormalized_gak(s1, s1, sigma=sigma, be=be) - * unnormalized_gak(s2, s2, sigma=sigma, be=be) - ) - return unnormalized_gak(s1, s2, sigma=sigma, be=be) / denom - - -def cdist_gak( - dataset1, - dataset2=None, - sigma=1.0, - n_jobs=None, - verbose=0, - be=None -): - r"""Compute cross-similarity matrix using Global Alignment kernel (GAK). - - GAK was originally presented in [1]_. - - Parameters - ---------- - dataset1 : array-like, shape=(n_ts1, sz1, d) or (n_ts1, sz1) or (sz1,) - A dataset of time series. - If shape is (n_ts1, sz1), the dataset is composed of univariate time series. - If shape is (sz1,), the dataset is composed of a unique univariate time series. - dataset2 : None or array-like, shape=(n_ts2, sz2, d) or (n_ts2, sz2) or (sz2,) (default: None) - Another dataset of time series. - If `None`, self-similarity of `dataset1` is returned. - If shape is (n_ts2, sz2), the dataset is composed of univariate time series. - If shape is (sz2,), the dataset is composed of a unique univariate time series. - sigma : float (default 1.) - Bandwidth of the internal gaussian kernel used for GAK - n_jobs : int or None, optional (default=None) - The number of jobs to run in parallel. - ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. - ``-1`` means using all processors. See scikit-learns' - `Glossary `__ - for more details. - verbose : int, optional (default=0) - The verbosity level: if non zero, progress messages are printed. - Above 50, the output is sent to stdout. - The frequency of the messages increases with the verbosity level. - If it more than 10, all iterations are reported. - `Glossary `__ - for more details. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - array-like, shape=(n_ts1, n_ts2) - Cross-similarity matrix. - - Examples - -------- - >>> cdist_gak([[1, 2, 2, 3], [1., 2., 3., 4.]], sigma=2.) - array([[1. , 0.65629661], - [0.65629661, 1. ]]) - >>> cdist_gak([[1, 2, 2], [1., 2., 3., 4.]], - ... [[1, 2, 2, 3], [1., 2., 3., 4.], [1, 2, 2, 3]], - ... sigma=2.) - array([[0.71059484, 0.29722877, 0.71059484], - [0.65629661, 1. , 0.65629661]]) - - See Also - -------- - gak : Compute Global Alignment kernel - - References - ---------- - .. [1] M. Cuturi, "Fast global alignment kernels," - ICML 2011. - """ # noqa: E501 - warnings.warn( - "This method is deprecated, use tslearn.metrics.unnormalized_gak instead.", - DeprecationWarning - ) - - be = instantiate_backend(be, dataset1, dataset2) - dataset1 = to_time_series_dataset(dataset1, be=be) - if dataset2 is not None: - dataset2 = to_time_series_dataset(dataset2, be=be) - return _cdist_gak( - dataset1=dataset1, - dataset2=dataset2, - sigma=sigma, - n_jobs=n_jobs, - verbose=verbose, - be=be, - ) - - -def _cdist_gak( - dataset1, - dataset2=None, - sigma=1.0, - n_jobs=None, - verbose=0, - be=None -): - if be is None: - be = instantiate_backend(dataset1, dataset2) - unnormalized_matrix = _cdist_generic( - dist_fun=unnormalized_gak, - dataset1=dataset1, - dataset2=dataset2, - n_jobs=n_jobs, - verbose=verbose, - sigma=sigma, - compute_diagonal=True, - be=be, - ) - if dataset2 is None: - diagonal = be.diag(be.sqrt(1.0 / be.diag(unnormalized_matrix))) - diagonal_left = diagonal_right = diagonal - else: - diagonal_left = Parallel(n_jobs=n_jobs, prefer="threads", verbose=verbose)( - delayed(unnormalized_gak)(dataset1[i], dataset1[i], sigma=sigma, be=be) - for i in range(len(dataset1)) - ) - diagonal_right = Parallel(n_jobs=n_jobs, prefer="threads", verbose=verbose)( - delayed(unnormalized_gak)(dataset2[j], dataset2[j], sigma=sigma, be=be) - for j in range(len(dataset2)) - ) - diagonal_left = be.diag(1.0 / be.sqrt(diagonal_left)) - diagonal_right = be.diag(1.0 / be.sqrt(diagonal_right)) - return diagonal_left @ unnormalized_matrix @ diagonal_right - - -def sigma_gak(dataset, n_samples=100, random_state=None, be=None): - r"""Compute sigma value to be used for GAK. - - This method was originally presented in [1]_. - - Parameters - ---------- - dataset : array-like, shape=(n_ts, sz, d) or (n_ts, sz1) or (sz,) - A dataset of time series. - If shape is (n_ts, sz), the dataset is composed of univariate time series. - If shape is (sz,), the dataset is composed of a unique univariate time series. - n_samples : int (default: 100) - Number of samples on which median distance should be estimated. - random_state : integer or numpy.RandomState or None (default: None) - The generator used to draw the samples. If an integer is given, it - fixes the seed. Defaults to the global numpy random number generator. - be : Backend object or string or None - Backend. If `be` is an instance of the class `NumPyBackend` or the string `"numpy"`, - the NumPy backend is used. - If `be` is an instance of the class `PyTorchBackend` or the string `"pytorch"`, - the PyTorch backend is used. - If `be` is `None`, the backend is determined by the input arrays. - See our :ref:`dedicated user-guide page ` for more information. - - Returns - ------- - float - Suggested bandwidth (:math:`\sigma`) for the Global Alignment kernel. - - Examples - -------- - >>> dataset = [[1, 2, 2, 3], [1., 2., 3., 4.]] - >>> float(sigma_gak(dataset=dataset, - ... n_samples=200, - ... random_state=0)) # doctest: +ELLIPSIS - 2.0... - - See Also - -------- - gak : Compute Global Alignment kernel - cdist_gak : Compute cross-similarity matrix using Global Alignment kernel - - References - ---------- - .. [1] M. Cuturi, "Fast global alignment kernels," - ICML 2011. - """ - warnings.warn( - "This method is deprecated, use tslearn.metrics.sigma_gak instead.", - DeprecationWarning - ) - - be = instantiate_backend(be, dataset) - - random_state = check_random_state(random_state) - dataset = to_time_series_dataset(dataset, be=be) - n_ts, sz, d = be.shape(dataset) - - # Remove points with nans from dataset - dataset = dataset.reshape((-1, d)) - mask = be.isfinite(be.sum(dataset, axis=-1)) - dataset = dataset[mask] - - nb_valid_samples = be.shape(dataset)[0] - replace = nb_valid_samples < n_samples - sample_indices = random_state.choice( - nb_valid_samples, - size=n_samples, - replace=replace - ) - dists = be.pdist( - dataset[sample_indices], - metric="euclidean", - ) - return be.median(dists) * be.sqrt(sz) - - def gamma_soft_dtw(dataset, n_samples=100, random_state=None, be=None): r"""Compute gamma value to be used for GAK/Soft-DTW. @@ -517,7 +78,6 @@ def gamma_soft_dtw(dataset, n_samples=100, random_state=None, be=None): .. [1] M. Cuturi, "Fast global alignment kernels," ICML 2011. """ - be = instantiate_backend(be, dataset) return ( 2.0 * sigma_gak(