undetectable/voice-clone
1
1#!/usr/bin/env python2# -*- encoding: utf-8 -*-3"""4Sequential modeling5===================6 7Sequence alignment8------------------9.. autosummary::10 :toctree: generated/11 12 dtw13 rqa14 15Viterbi decoding16----------------17.. autosummary::18 :toctree: generated/19 20 viterbi21 viterbi_discriminative22 viterbi_binary23 24Transition matrices25-------------------26.. autosummary::27 :toctree: generated/28 29 transition_uniform30 transition_loop31 transition_cycle32 transition_local33"""34from __future__ import annotations35 36import numpy as np37from scipy.spatial.distance import cdist38from numba import jit39from .util import pad_center, fill_off_diagonal, is_positive_int, tiny, expand_to40from .util.exceptions import ParameterError41from .filters import get_window42from typing import Any, Iterable, List, Optional, Tuple, Union, overload43from typing_extensions import Literal44from ._typing import _WindowSpec, _IntLike_co45 46__all__ = [47 "dtw",48 "dtw_backtracking",49 "rqa",50 "viterbi",51 "viterbi_discriminative",52 "viterbi_binary",53 "transition_uniform",54 "transition_loop",55 "transition_cycle",56 "transition_local",57]58 59 60@overload61def dtw(62 X: np.ndarray,63 Y: np.ndarray,64 *,65 metric: str = ...,66 step_sizes_sigma: Optional[np.ndarray] = ...,67 weights_add: Optional[np.ndarray] = ...,68 weights_mul: Optional[np.ndarray] = ...,69 subseq: bool = ...,70 backtrack: Literal[False],71 global_constraints: bool = ...,72 band_rad: float = ...,73 return_steps: Literal[False] = ...,74) -> np.ndarray:75 ...76 77 78@overload79def dtw(80 *,81 C: np.ndarray,82 metric: str = ...,83 step_sizes_sigma: Optional[np.ndarray] = ...,84 weights_add: Optional[np.ndarray] = ...,85 weights_mul: Optional[np.ndarray] = ...,86 subseq: bool = ...,87 backtrack: Literal[False],88 global_constraints: bool = ...,89 band_rad: float = ...,90 return_steps: Literal[False] = ...,91) -> np.ndarray:92 ...93 94 95@overload96def dtw(97 X: np.ndarray,98 Y: np.ndarray,99 *,100 metric: str = ...,101 step_sizes_sigma: Optional[np.ndarray] = ...,102 weights_add: Optional[np.ndarray] = ...,103 weights_mul: Optional[np.ndarray] = ...,104 subseq: bool = ...,105 backtrack: Literal[False],106 global_constraints: bool = ...,107 band_rad: float = ...,108 return_steps: Literal[True],109) -> Tuple[np.ndarray, np.ndarray]:110 ...111 112 113@overload114def dtw(115 *,116 C: np.ndarray,117 metric: str = ...,118 step_sizes_sigma: Optional[np.ndarray] = ...,119 weights_add: Optional[np.ndarray] = ...,120 weights_mul: Optional[np.ndarray] = ...,121 subseq: bool = ...,122 backtrack: Literal[False],123 global_constraints: bool = ...,124 band_rad: float = ...,125 return_steps: Literal[True],126) -> Tuple[np.ndarray, np.ndarray]:127 ...128 129 130@overload131def dtw(132 X: np.ndarray,133 Y: np.ndarray,134 *,135 metric: str = ...,136 step_sizes_sigma: Optional[np.ndarray] = ...,137 weights_add: Optional[np.ndarray] = ...,138 weights_mul: Optional[np.ndarray] = ...,139 subseq: bool = ...,140 backtrack: Literal[True] = ...,141 global_constraints: bool = ...,142 band_rad: float = ...,143 return_steps: Literal[False] = ...,144) -> Tuple[np.ndarray, np.ndarray]:145 ...146 147 148@overload149def dtw(150 *,151 C: np.ndarray,152 metric: str = ...,153 step_sizes_sigma: Optional[np.ndarray] = ...,154 weights_add: Optional[np.ndarray] = ...,155 weights_mul: Optional[np.ndarray] = ...,156 subseq: bool = ...,157 backtrack: Literal[True] = ...,158 global_constraints: bool = ...,159 band_rad: float = ...,160 return_steps: Literal[False] = ...,161) -> Tuple[np.ndarray, np.ndarray]:162 ...163 164 165@overload166def dtw(167 X: np.ndarray,168 Y: np.ndarray,169 *,170 metric: str = ...,171 step_sizes_sigma: Optional[np.ndarray] = ...,172 weights_add: Optional[np.ndarray] = ...,173 weights_mul: Optional[np.ndarray] = ...,174 subseq: bool = ...,175 backtrack: Literal[True] = ...,176 global_constraints: bool = ...,177 band_rad: float = ...,178 return_steps: Literal[True],179) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:180 ...181 182 183@overload184def dtw(185 *,186 C: np.ndarray,187 metric: str = ...,188 step_sizes_sigma: Optional[np.ndarray] = ...,189 weights_add: Optional[np.ndarray] = ...,190 weights_mul: Optional[np.ndarray] = ...,191 subseq: bool = ...,192 backtrack: Literal[True] = ...,193 global_constraints: bool = ...,194 band_rad: float = ...,195 return_steps: Literal[True],196) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:197 ...198 199 200def dtw(201 X: Optional[np.ndarray] = None,202 Y: Optional[np.ndarray] = None,203 *,204 C: Optional[np.ndarray] = None,205 metric: str = "euclidean",206 step_sizes_sigma: Optional[np.ndarray] = None,207 weights_add: Optional[np.ndarray] = None,208 weights_mul: Optional[np.ndarray] = None,209 subseq: bool = False,210 backtrack: bool = True,211 global_constraints: bool = False,212 band_rad: float = 0.25,213 return_steps: bool = False,214) -> Union[215 np.ndarray, Tuple[np.ndarray, np.ndarray], Tuple[np.ndarray, np.ndarray, np.ndarray]216]:217 """Dynamic time warping (DTW).218 219 This function performs a DTW and path backtracking on two sequences.220 We follow the nomenclature and algorithmic approach as described in [#]_.221 222 .. [#] Meinard Mueller223 Fundamentals of Music Processing — Audio, Analysis, Algorithms, Applications224 Springer Verlag, ISBN: 978-3-319-21944-8, 2015.225 226 Parameters227 ----------228 X : np.ndarray [shape=(..., K, N)]229 audio feature matrix (e.g., chroma features)230 231 If ``X`` has more than two dimensions (e.g., for multi-channel inputs), all leading232 dimensions are used when computing distance to ``Y``.233 234 Y : np.ndarray [shape=(..., K, M)]235 audio feature matrix (e.g., chroma features)236 237 C : np.ndarray [shape=(N, M)]238 Precomputed distance matrix. If supplied, X and Y must not be supplied and239 ``metric`` will be ignored.240 241 metric : str242 Identifier for the cost-function as documented243 in `scipy.spatial.distance.cdist()`244 245 step_sizes_sigma : np.ndarray [shape=[n, 2]]246 Specifies allowed step sizes as used by the dtw.247 248 weights_add : np.ndarray [shape=[n, ]]249 Additive weights to penalize certain step sizes.250 251 weights_mul : np.ndarray [shape=[n, ]]252 Multiplicative weights to penalize certain step sizes.253 254 subseq : bool255 Enable subsequence DTW, e.g., for retrieval tasks.256 257 backtrack : bool258 Enable backtracking in accumulated cost matrix.259 260 global_constraints : bool261 Applies global constraints to the cost matrix ``C`` (Sakoe-Chiba band).262 263 band_rad : float264 The Sakoe-Chiba band radius (1/2 of the width) will be265 ``int(radius*min(C.shape))``.266 267 return_steps : bool268 If true, the function returns ``steps``, the step matrix, containing269 the indices of the used steps from the cost accumulation step.270 271 Returns272 -------273 D : np.ndarray [shape=(N, M)]274 accumulated cost matrix.275 D[N, M] is the total alignment cost.276 When doing subsequence DTW, D[N,:] indicates a matching function.277 wp : np.ndarray [shape=(N, 2)]278 Warping path with index pairs.279 Each row of the array contains an index pair (n, m).280 Only returned when ``backtrack`` is True.281 steps : np.ndarray [shape=(N, M)]282 Step matrix, containing the indices of the used steps from the cost283 accumulation step.284 Only returned when ``return_steps`` is True.285 286 Raises287 ------288 ParameterError289 If you are doing diagonal matching and Y is shorter than X or if an290 incompatible combination of X, Y, and C are supplied.291 292 If your input dimensions are incompatible.293 294 If the cost matrix has NaN values.295 296 Examples297 --------298 >>> import numpy as np299 >>> import matplotlib.pyplot as plt300 >>> y, sr = librosa.load(librosa.ex('brahms'), offset=10, duration=15)301 >>> X = librosa.feature.chroma_cens(y=y, sr=sr)302 >>> noise = np.random.rand(X.shape[0], 200)303 >>> Y = np.concatenate((noise, noise, X, noise), axis=1)304 >>> D, wp = librosa.sequence.dtw(X, Y, subseq=True)305 >>> fig, ax = plt.subplots(nrows=2, sharex=True)306 >>> img = librosa.display.specshow(D, x_axis='frames', y_axis='frames',307 ... ax=ax[0])308 >>> ax[0].set(title='DTW cost', xlabel='Noisy sequence', ylabel='Target')309 >>> ax[0].plot(wp[:, 1], wp[:, 0], label='Optimal path', color='y')310 >>> ax[0].legend()311 >>> fig.colorbar(img, ax=ax[0])312 >>> ax[1].plot(D[-1, :] / wp.shape[0])313 >>> ax[1].set(xlim=[0, Y.shape[1]], ylim=[0, 2],314 ... title='Matching cost function')315 """316 # Default Parameters317 default_steps = np.array([[1, 1], [0, 1], [1, 0]], dtype=np.uint32)318 default_weights_add = np.zeros(3, dtype=np.float64)319 default_weights_mul = np.ones(3, dtype=np.float64)320 321 if step_sizes_sigma is None:322 # Use the default steps323 step_sizes_sigma = default_steps324 325 # Use default weights if none are provided326 if weights_add is None:327 weights_add = default_weights_add328 329 if weights_mul is None:330 weights_mul = default_weights_mul331 else:332 # If we have custom steps but no weights, construct them here333 if weights_add is None:334 weights_add = np.zeros(len(step_sizes_sigma), dtype=np.float64)335 336 if weights_mul is None:337 weights_mul = np.ones(len(step_sizes_sigma), dtype=np.float64)338 339 # Make the default step weights infinite so that they are never340 # preferred over custom steps341 default_weights_add.fill(np.inf)342 default_weights_mul.fill(np.inf)343 344 # Append custom steps and weights to our defaults345 step_sizes_sigma = np.concatenate((default_steps, step_sizes_sigma))346 weights_add = np.concatenate((default_weights_add, weights_add))347 weights_mul = np.concatenate((default_weights_mul, weights_mul))348 349 # These asserts are bad, but mypy cannot trace the code paths properly350 assert step_sizes_sigma is not None351 assert weights_add is not None352 assert weights_mul is not None353 354 if np.any(step_sizes_sigma < 0):355 raise ParameterError("step_sizes_sigma cannot contain negative values")356 357 if len(step_sizes_sigma) != len(weights_add):358 raise ParameterError("len(weights_add) must be equal to len(step_sizes_sigma)")359 if len(step_sizes_sigma) != len(weights_mul):360 raise ParameterError("len(weights_mul) must be equal to len(step_sizes_sigma)")361 362 if C is None and (X is None or Y is None):363 raise ParameterError("If C is not supplied, both X and Y must be supplied")364 if C is not None and (X is not None or Y is not None):365 raise ParameterError("If C is supplied, both X and Y must not be supplied")366 367 c_is_transposed = False368 369 # calculate pair-wise distances, unless already supplied.370 # C_local will keep track of whether the distance matrix was supplied371 # by the user (False) or constructed locally (True)372 C_local = False373 if C is None:374 C_local = True375 # mypy can't figure out that this case does not happen376 assert X is not None and Y is not None377 # take care of dimensions378 X = np.atleast_2d(X)379 Y = np.atleast_2d(Y)380 381 # Perform some shape-squashing here382 # Put the time axes around front383 # Suppress types because mypy doesn't know these are ndarrays384 X = np.swapaxes(X, -1, 0) # type: ignore385 Y = np.swapaxes(Y, -1, 0) # type: ignore386 387 # Flatten the remaining dimensions388 # Use F-ordering to preserve columns389 X = X.reshape((X.shape[0], -1), order="F")390 Y = Y.reshape((Y.shape[0], -1), order="F")391 392 try:393 C = cdist(X, Y, metric=metric)394 except ValueError as exc:395 raise ParameterError(396 "scipy.spatial.distance.cdist returned an error.\n"397 "Please provide your input in the form X.shape=(K, N) "398 "and Y.shape=(K, M).\n 1-dimensional sequences should "399 "be reshaped to X.shape=(1, N) and Y.shape=(1, M)."400 ) from exc401 402 # for subsequence matching:403 # if N > M, Y can be a subsequence of X404 if subseq and (X.shape[0] > Y.shape[0]):405 C = C.T406 c_is_transposed = True407 408 C = np.atleast_2d(C)409 410 # if diagonal matching, Y has to be longer than X411 # (X simply cannot be contained in Y)412 if np.array_equal(step_sizes_sigma, np.array([[1, 1]])) and (413 C.shape[0] > C.shape[1]414 ):415 raise ParameterError(416 "For diagonal matching: Y.shape[-1] >= X.shape[-11] "417 "(C.shape[1] >= C.shape[0])"418 )419 420 max_0 = step_sizes_sigma[:, 0].max()421 max_1 = step_sizes_sigma[:, 1].max()422 423 # check C here for nans before building global constraints424 if np.any(np.isnan(C)):425 raise ParameterError("DTW cost matrix C has NaN values. ")426 427 if global_constraints:428 # Apply global constraints to the cost matrix429 if not C_local:430 # If C was provided as input, make a copy here431 C = np.copy(C)432 fill_off_diagonal(C, radius=band_rad, value=np.inf)433 434 # initialize whole matrix with infinity values435 D = np.ones(C.shape + np.array([max_0, max_1])) * np.inf436 437 # set starting point to C[0, 0]438 D[max_0, max_1] = C[0, 0]439 440 if subseq:441 D[max_0, max_1:] = C[0, :]442 443 # initialize step matrix with -1444 # will be filled in calc_accu_cost() with indices from step_sizes_sigma445 steps = np.zeros(D.shape, dtype=np.int32)446 447 # these steps correspond to left- (first row) and up-(first column) moves448 steps[0, :] = 1449 steps[:, 0] = 2450 451 # calculate accumulated cost matrix452 D: np.ndarray453 steps: np.ndarray454 D, steps = __dtw_calc_accu_cost(455 C, D, steps, step_sizes_sigma, weights_mul, weights_add, max_0, max_1456 )457 458 # delete infinity rows and columns459 D = D[max_0:, max_1:]460 steps = steps[max_0:, max_1:]461 462 return_values: List[np.ndarray]463 if backtrack:464 wp: np.ndarray465 if subseq:466 if np.all(np.isinf(D[-1])):467 raise ParameterError(468 "No valid sub-sequence warping path could "469 "be constructed with the given step sizes."470 )471 start = np.argmin(D[-1, :])472 _wp = __dtw_backtracking(steps, step_sizes_sigma, subseq, start)473 else:474 # perform warping path backtracking475 if np.isinf(D[-1, -1]):476 raise ParameterError(477 "No valid sub-sequence warping path could "478 "be constructed with the given step sizes."479 )480 481 _wp = __dtw_backtracking(steps, step_sizes_sigma, subseq)482 if _wp[-1] != (0, 0):483 raise ParameterError(484 "Unable to compute a full DTW warping path. "485 "You may want to try again with subseq=True."486 )487 488 wp = np.asarray(_wp, dtype=int)489 490 # since we transposed in the beginning, we have to adjust the index pairs back491 if subseq and (492 (X is not None and Y is not None and X.shape[0] > Y.shape[0])493 or c_is_transposed494 or C.shape[0] > C.shape[1]495 ):496 wp = np.fliplr(wp)497 return_values = [D, wp]498 else:499 return_values = [D]500 501 if return_steps:502 return_values.append(steps)503 504 if len(return_values) > 1:505 # Suppressing type check here because mypy can't506 # infer the exact length of the tuple507 return tuple(return_values) # type: ignore508 else:509 return return_values[0]510 511 512@jit(nopython=True, cache=False) # type: ignore513def __dtw_calc_accu_cost(514 C: np.ndarray,515 D: np.ndarray,516 steps: np.ndarray,517 step_sizes_sigma: np.ndarray,518 weights_mul: np.ndarray,519 weights_add: np.ndarray,520 max_0: int,521 max_1: int,522) -> Tuple[np.ndarray, np.ndarray]: # pragma: no cover523 """Calculate the accumulated cost matrix D.524 525 Use dynamic programming to calculate the accumulated costs.526 527 Parameters528 ----------529 C : np.ndarray [shape=(N, M)]530 pre-computed cost matrix531 D : np.ndarray [shape=(N, M)]532 accumulated cost matrix533 steps : np.ndarray [shape=(N, M)]534 Step matrix, containing the indices of the used steps from the cost535 accumulation step.536 step_sizes_sigma : np.ndarray [shape=[n, 2]]537 Specifies allowed step sizes as used by the dtw.538 weights_add : np.ndarray [shape=[n, ]]539 Additive weights to penalize certain step sizes.540 weights_mul : np.ndarray [shape=[n, ]]541 Multiplicative weights to penalize certain step sizes.542 max_0 : int543 maximum number of steps in step_sizes_sigma in dim 0.544 max_1 : int545 maximum number of steps in step_sizes_sigma in dim 1.546 547 Returns548 -------549 D : np.ndarray [shape=(N, M)]550 accumulated cost matrix.551 D[N, M] is the total alignment cost.552 When doing subsequence DTW, D[N,:] indicates a matching function.553 steps : np.ndarray [shape=(N, M)]554 Step matrix, containing the indices of the used steps from the cost555 accumulation step.556 557 See Also558 --------559 dtw560 """561 for cur_n in range(max_0, D.shape[0]):562 for cur_m in range(max_1, D.shape[1]):563 # accumulate costs564 for cur_step_idx, cur_w_add, cur_w_mul in zip(565 range(step_sizes_sigma.shape[0]), weights_add, weights_mul566 ):567 cur_D = D[568 cur_n - step_sizes_sigma[cur_step_idx, 0],569 cur_m - step_sizes_sigma[cur_step_idx, 1],570 ]571 cur_C = cur_w_mul * C[cur_n - max_0, cur_m - max_1]572 cur_C += cur_w_add573 cur_cost = cur_D + cur_C574 575 # check if cur_cost is smaller than the one stored in D576 if cur_cost < D[cur_n, cur_m]:577 D[cur_n, cur_m] = cur_cost578 579 # save step-index580 steps[cur_n, cur_m] = cur_step_idx581 582 return D, steps583 584 585@jit(nopython=True, cache=False) # type: ignore586def __dtw_backtracking(587 steps: np.ndarray,588 step_sizes_sigma: np.ndarray,589 subseq: bool,590 start: Optional[int] = None,591) -> List[Tuple[int, int]]: # pragma: no cover592 """Backtrack optimal warping path.593 594 Uses the saved step sizes from the cost accumulation595 step to backtrack the index pairs for an optimal596 warping path.597 598 Parameters599 ----------600 steps : np.ndarray [shape=(N, M)]601 Step matrix, containing the indices of the used steps from the cost602 accumulation step.603 step_sizes_sigma : np.ndarray [shape=[n, 2]]604 Specifies allowed step sizes as used by the dtw.605 subseq : bool606 Enable subsequence DTW, e.g., for retrieval tasks.607 start : int608 Start column index for backtraing (only allowed for ``subseq=True``)609 610 Returns611 -------612 wp : list [shape=(N,)]613 Warping path with index pairs.614 Each list entry contains an index pair615 (n, m) as a tuple616 617 See Also618 --------619 dtw620 """621 if start is None:622 cur_idx = (steps.shape[0] - 1, steps.shape[1] - 1)623 else:624 cur_idx = (steps.shape[0] - 1, start)625 626 wp = []627 # Set starting point D(N, M) and append it to the path628 wp.append((cur_idx[0], cur_idx[1]))629 630 # Loop backwards.631 # Stop criteria:632 # Setting it to (0, 0) does not work for the subsequence dtw,633 # so we only ask to reach the first row of the matrix.634 635 while (subseq and cur_idx[0] > 0) or (not subseq and cur_idx != (0, 0)):636 cur_step_idx = steps[(cur_idx[0], cur_idx[1])]637 638 # save tuple with minimal acc. cost in path639 cur_idx = (640 cur_idx[0] - step_sizes_sigma[cur_step_idx][0],641 cur_idx[1] - step_sizes_sigma[cur_step_idx][1],642 )643 644 # If we run off the side of the cost matrix, break here645 if min(cur_idx) < 0:646 break647 648 # append to warping path649 wp.append((cur_idx[0], cur_idx[1]))650 651 return wp652 653 654def dtw_backtracking(655 steps: np.ndarray,656 *,657 step_sizes_sigma: Optional[np.ndarray] = None,658 subseq: bool = False,659 start: Optional[Union[int, np.integer[Any]]] = None,660) -> np.ndarray:661 """Backtrack a warping path.662 663 Uses the saved step sizes from the cost accumulation664 step to backtrack the index pairs for a warping path.665 666 Parameters667 ----------668 steps : np.ndarray [shape=(N, M)]669 Step matrix, containing the indices of the used steps from the cost670 accumulation step.671 step_sizes_sigma : np.ndarray [shape=[n, 2]]672 Specifies allowed step sizes as used by the dtw.673 subseq : bool674 Enable subsequence DTW, e.g., for retrieval tasks.675 start : int676 Start column index for backtraing (only allowed for ``subseq=True``)677 678 Returns679 -------680 wp : list [shape=(N,)]681 Warping path with index pairs.682 Each list entry contains an index pair683 (n, m) as a tuple684 685 See Also686 --------687 dtw688 """689 if subseq is False and start is not None:690 raise ParameterError(691 f"start is only allowed to be set if subseq is True (start={start}, subseq={subseq})"692 )693 694 # Default Parameters695 default_steps = np.array([[1, 1], [0, 1], [1, 0]], dtype=np.uint32)696 697 if step_sizes_sigma is None:698 # Use the default steps699 step_sizes_sigma = default_steps700 else:701 # Append custom steps and weights to our defaults702 step_sizes_sigma = np.concatenate((default_steps, step_sizes_sigma))703 704 wp = __dtw_backtracking(steps, step_sizes_sigma, subseq, start)705 return np.asarray(wp, dtype=int)706 707 708@overload709def rqa(710 sim: np.ndarray,711 *,712 gap_onset: float = ...,713 gap_extend: float = ...,714 knight_moves: bool = ...,715 backtrack: Literal[False],716) -> np.ndarray:717 ...718 719 720@overload721def rqa(722 sim: np.ndarray,723 *,724 gap_onset: float = ...,725 gap_extend: float = ...,726 knight_moves: bool = ...,727 backtrack: Literal[True] = ...,728) -> Tuple[np.ndarray, np.ndarray]:729 ...730 731 732@overload733def rqa(734 sim: np.ndarray,735 *,736 gap_onset: float = ...,737 gap_extend: float = ...,738 knight_moves: bool = ...,739 backtrack: bool = ...,740) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:741 ...742 743 744def rqa(745 sim: np.ndarray,746 *,747 gap_onset: float = 1,748 gap_extend: float = 1,749 knight_moves: bool = True,750 backtrack: bool = True,751) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:752 """Recurrence quantification analysis (RQA)753 754 This function implements different forms of RQA as described by755 Serra, Serra, and Andrzejak (SSA). [#]_ These methods take as input756 a self- or cross-similarity matrix ``sim``, and calculate the value757 of path alignments by dynamic programming.758 759 Note that unlike dynamic time warping (`dtw`), alignment paths here are760 maximized, not minimized, so the input should measure similarity rather761 than distance.762 763 The simplest RQA method, denoted as `L` (SSA equation 3) and equivalent764 to the method described by Eckman, Kamphorst, and Ruelle [#]_, accumulates765 the length of diagonal paths with positive values in the input:766 767 - ``score[i, j] = score[i-1, j-1] + 1`` if ``sim[i, j] > 0``768 - ``score[i, j] = 0`` otherwise.769 770 The second method, denoted as `S` (SSA equation 4), is similar to the first,771 but allows for "knight moves" (as in the chess piece) in addition to strict772 diagonal moves:773 774 - ``score[i, j] = max(score[i-1, j-1], score[i-2, j-1], score[i-1, j-2]) + 1`` if ``sim[i, j] >775 0``776 - ``score[i, j] = 0`` otherwise.777 778 The third method, denoted as `Q` (SSA equations 5 and 6) extends this by779 allowing gaps in the alignment that incur some cost, rather than a hard780 reset to 0 whenever ``sim[i, j] == 0``.781 Gaps are penalized by two additional parameters, ``gap_onset`` and ``gap_extend``,782 which are subtracted from the value of the alignment path every time a gap783 is introduced or extended (respectively).784 785 Note that setting ``gap_onset`` and ``gap_extend`` to `np.inf` recovers the second786 method, and disabling knight moves recovers the first.787 788 .. [#] Serrà, Joan, Xavier Serra, and Ralph G. Andrzejak.789 "Cross recurrence quantification for cover song identification."790 New Journal of Physics 11, no. 9 (2009): 093017.791 792 .. [#] Eckmann, J. P., S. Oliffson Kamphorst, and D. Ruelle.793 "Recurrence plots of dynamical systems."794 World Scientific Series on Nonlinear Science Series A 16 (1995): 441-446.795 796 Parameters797 ----------798 sim : np.ndarray [shape=(N, M), non-negative]799 The similarity matrix to use as input.800 801 This can either be a recurrence matrix (self-similarity)802 or a cross-similarity matrix between two sequences.803 804 gap_onset : float > 0805 Penalty for introducing a gap to an alignment sequence806 807 gap_extend : float > 0808 Penalty for extending a gap in an alignment sequence809 810 knight_moves : bool811 If ``True`` (default), allow for "knight moves" in the alignment,812 e.g., ``(n, m) => (n + 1, m + 2)`` or ``(n + 2, m + 1)``.813 814 If ``False``, only allow for diagonal moves ``(n, m) => (n + 1, m + 1)``.815 816 backtrack : bool817 If ``True``, return the alignment path.818 819 If ``False``, only return the score matrix.820 821 Returns822 -------823 score : np.ndarray [shape=(N, M)]824 The alignment score matrix. ``score[n, m]`` is the cumulative value of825 the best alignment sequence ending in frames ``n`` and ``m``.826 path : np.ndarray [shape=(k, 2)] (optional)827 If ``backtrack=True``, ``path`` contains a list of pairs of aligned frames828 in the best alignment sequence.829 830 ``path[i] = [n, m]`` indicates that row ``n`` aligns to column ``m``.831 832 See Also833 --------834 librosa.segment.recurrence_matrix835 librosa.segment.cross_similarity836 dtw837 838 Examples839 --------840 Simple diagonal path enhancement (L-mode)841 842 >>> import numpy as np843 >>> import matplotlib.pyplot as plt844 >>> y, sr = librosa.load(librosa.ex('nutcracker'), duration=30)845 >>> chroma = librosa.feature.chroma_cqt(y=y, sr=sr)846 >>> # Use time-delay embedding to reduce noise847 >>> chroma_stack = librosa.feature.stack_memory(chroma, n_steps=10, delay=3)848 >>> # Build recurrence, suppress self-loops within 1 second849 >>> rec = librosa.segment.recurrence_matrix(chroma_stack, width=43,850 ... mode='affinity',851 ... metric='cosine')852 >>> # using infinite cost for gaps enforces strict path continuation853 >>> L_score, L_path = librosa.sequence.rqa(rec,854 ... gap_onset=np.inf,855 ... gap_extend=np.inf,856 ... knight_moves=False)857 >>> fig, ax = plt.subplots(ncols=2)858 >>> librosa.display.specshow(rec, x_axis='frames', y_axis='frames', ax=ax[0])859 >>> ax[0].set(title='Recurrence matrix')860 >>> librosa.display.specshow(L_score, x_axis='frames', y_axis='frames', ax=ax[1])861 >>> ax[1].set(title='Alignment score matrix')862 >>> ax[1].plot(L_path[:, 1], L_path[:, 0], label='Optimal path', color='c')863 >>> ax[1].legend()864 >>> ax[1].label_outer()865 866 Full alignment using gaps and knight moves867 868 >>> # New gaps cost 5, extending old gaps cost 10 for each step869 >>> score, path = librosa.sequence.rqa(rec, gap_onset=5, gap_extend=10)870 >>> fig, ax = plt.subplots(ncols=2, sharex=True, sharey=True)871 >>> librosa.display.specshow(rec, x_axis='frames', y_axis='frames', ax=ax[0])872 >>> ax[0].set(title='Recurrence matrix')873 >>> librosa.display.specshow(score, x_axis='frames', y_axis='frames', ax=ax[1])874 >>> ax[1].set(title='Alignment score matrix')875 >>> ax[1].plot(path[:, 1], path[:, 0], label='Optimal path', color='c')876 >>> ax[1].legend()877 >>> ax[1].label_outer()878 """879 880 if gap_onset < 0:881 raise ParameterError("gap_onset={} must be strictly positive")882 if gap_extend < 0:883 raise ParameterError("gap_extend={} must be strictly positive")884 885 score: np.ndarray886 pointers: np.ndarray887 score, pointers = __rqa_dp(sim, gap_onset, gap_extend, knight_moves)888 if backtrack:889 path = __rqa_backtrack(score, pointers)890 return score, path891 892 return score893 894 895@jit(nopython=True, cache=False) # type: ignore896def __rqa_dp(897 sim: np.ndarray, gap_onset: float, gap_extend: float, knight: bool898) -> Tuple[np.ndarray, np.ndarray]: # pragma: no cover899 """RQA dynamic programming implementation"""900 901 # The output array902 score = np.zeros(sim.shape, dtype=sim.dtype)903 904 # The backtracking array905 backtrack = np.zeros(sim.shape, dtype=np.int8)906 907 # These are place-holder arrays to limit the points being considered908 # at each step of the DP909 #910 # If knight moves are enabled, values are indexed according to911 # [(-1,-1), (-1, -2), (-2, -1)]912 #913 # If knight moves are disabled, then only the first entry is used.914 #915 # Using dummy vectors here makes the code a bit cleaner down below.916 sim_values = np.zeros(3)917 score_values = np.zeros(3)918 vec = np.zeros(3)919 920 if knight:921 # Initial limit is for the base case: diagonal + one knight922 init_limit = 2923 924 # Otherwise, we have 3 positions925 limit = 3926 else:927 init_limit = 1928 limit = 1929 930 # backtracking rubric:931 # 0 ==> diagonal move932 # 1 ==> knight move up933 # 2 ==> knight move left934 # -1 ==> reset without inclusion935 # -2 ==> reset with inclusion (ie positive value at init)936 937 # Initialize the first row and column with the data938 score[0, :] = sim[0, :]939 score[:, 0] = sim[:, 0]940 941 # backtracking initialization: the first row and column are all resets942 # if there's a positive link here, it's an inclusive reset943 for i in range(sim.shape[0]):944 if sim[i, 0]:945 backtrack[i, 0] = -2946 else:947 backtrack[i, 0] = -1948 949 for j in range(sim.shape[1]):950 if sim[0, j]:951 backtrack[0, j] = -2952 else:953 backtrack[0, j] = -1954 955 # Initialize the 1-1 case using only the diagonal956 if sim[1, 1] > 0:957 score[1, 1] = score[0, 0] + sim[1, 1]958 backtrack[1, 1] = 0959 else:960 link = sim[0, 0] > 0961 score[1, 1] = max(0, score[0, 0] - (link) * gap_onset - (~link) * gap_extend)962 if score[1, 1] > 0:963 backtrack[1, 1] = 0964 else:965 backtrack[1, 1] = -1966 967 # Initialize the second row with diagonal and left-knight moves968 i = 1969 for j in range(2, sim.shape[1]):970 score_values[:-1] = (score[i - 1, j - 1], score[i - 1, j - 2])971 sim_values[:-1] = (sim[i - 1, j - 1], sim[i - 1, j - 2])972 t_values = sim_values > 0973 if sim[i, j] > 0:974 backtrack[i, j] = np.argmax(score_values[:init_limit])975 score[i, j] = score_values[backtrack[i, j]] + sim[i, j] # or + 1 for binary976 else:977 vec[:init_limit] = (978 score_values[:init_limit]979 - t_values[:init_limit] * gap_onset980 - (~t_values[:init_limit]) * gap_extend981 )982 983 backtrack[i, j] = np.argmax(vec[:init_limit])984 score[i, j] = max(0, vec[backtrack[i, j]])985 # Is it a reset?986 if score[i, j] == 0:987 backtrack[i, j] = -1988 989 # Initialize the second column with diagonal and up-knight moves990 j = 1991 for i in range(2, sim.shape[0]):992 score_values[:-1] = (score[i - 1, j - 1], score[i - 2, j - 1])993 sim_values[:-1] = (sim[i - 1, j - 1], sim[i - 2, j - 1])994 t_values = sim_values > 0995 if sim[i, j] > 0:996 backtrack[i, j] = np.argmax(score_values[:init_limit])997 score[i, j] = score_values[backtrack[i, j]] + sim[i, j] # or + 1 for binary998 999 else:1000 vec[:init_limit] = (1001 score_values[:init_limit]1002 - t_values[:init_limit] * gap_onset1003 - (~t_values[:init_limit]) * gap_extend1004 )1005 1006 backtrack[i, j] = np.argmax(vec[:init_limit])1007 score[i, j] = max(0, vec[backtrack[i, j]])1008 # Is it a reset?1009 if score[i, j] == 0:1010 backtrack[i, j] = -11011 1012 # Now fill in the rest of the table1013 for i in range(2, sim.shape[0]):1014 for j in range(2, sim.shape[1]):1015 score_values[:] = (1016 score[i - 1, j - 1],1017 score[i - 1, j - 2],1018 score[i - 2, j - 1],1019 )1020 sim_values[:] = (sim[i - 1, j - 1], sim[i - 1, j - 2], sim[i - 2, j - 1])1021 t_values = sim_values > 01022 if sim[i, j] > 0:1023 # if knight is true, it's max of (-1,-1), (-1, -2), (-2, -1)1024 # otherwise, it's just the diagonal move (-1, -1)1025 # for backtracking purposes, if the max is 0 then it's the start of a new sequence1026 # if the max is non-zero, then we extend the existing sequence1027 backtrack[i, j] = np.argmax(score_values[:limit])1028 score[i, j] = (1029 score_values[backtrack[i, j]] + sim[i, j]1030 ) # or + 1 for binary1031 1032 else:1033 # if the max of our options is negative, then it's a hard reset1034 # otherwise, it's a skip move1035 vec[:limit] = (1036 score_values[:limit]1037 - t_values[:limit] * gap_onset1038 - (~t_values[:limit]) * gap_extend1039 )1040 1041 backtrack[i, j] = np.argmax(vec[:limit])1042 score[i, j] = max(0, vec[backtrack[i, j]])1043 # Is it a reset?1044 if score[i, j] == 0:1045 backtrack[i, j] = -11046 1047 return score, backtrack1048 1049 1050def __rqa_backtrack(score, pointers):1051 """RQA path backtracking1052 1053 Given the score matrix and backtracking index array,1054 reconstruct the optimal path.1055 """1056 1057 # backtracking rubric:1058 # 0 ==> diagonal move1059 # 1 ==> knight move up1060 # 2 ==> knight move left1061 # -1 ==> reset (sim = 0)1062 # -2 ==> start of sequence (sim > 0)1063 1064 # This array maps the backtracking values to the1065 # relative index offsets1066 offsets = [(-1, -1), (-1, -2), (-2, -1)]1067 1068 # Find the maximum to end the path1069 idx = list(np.unravel_index(np.argmax(score), score.shape))1070 1071 # Construct the path1072 path: List = []1073 while True:1074 bt_index = pointers[tuple(idx)]1075 1076 # A -1 indicates a non-inclusive reset1077 # this can only happen when sim[idx] == 0,1078 # and a reset with zero score should not be included1079 # in the path. In this case, we're done.1080 if bt_index == -1:1081 break1082 1083 # Other bt_index values are okay for inclusion1084 path.insert(0, idx)1085 1086 # -2 indicates beginning of sequence,1087 # so we can't backtrack any further1088 if bt_index == -2:1089 break1090 1091 # Otherwise, prepend this index and continue1092 idx = [idx[_] + offsets[bt_index][_] for _ in range(len(idx))]1093 1094 # If there's no alignment path at all, eg an empty cross-similarity1095 # matrix, return a properly shaped and typed array1096 if not path:1097 return np.empty((0, 2), dtype=np.uint)1098 1099 return np.asarray(path, dtype=np.uint)1100 1101 1102@jit(nopython=True, cache=False) # type: ignore1103def _viterbi(1104 log_prob: np.ndarray, log_trans: np.ndarray, log_p_init: np.ndarray1105) -> Tuple[np.ndarray, np.ndarray]: # pragma: no cover1106 """Core Viterbi algorithm.1107 1108 This is intended for internal use only.1109 1110 Parameters1111 ----------1112 log_prob : np.ndarray [shape=(T, m)]1113 ``log_prob[t, s]`` is the conditional log-likelihood1114 ``log P[X = X(t) | State(t) = s]``1115 log_trans : np.ndarray [shape=(m, m)]1116 The log transition matrix1117 ``log_trans[i, j] = log P[State(t+1) = j | State(t) = i]``1118 log_p_init : np.ndarray [shape=(m,)]1119 log of the initial state distribution1120 1121 Returns1122 -------1123 None1124 All computations are performed in-place on ``state, value, ptr``.1125 """1126 n_steps, n_states = log_prob.shape1127 1128 state = np.zeros(n_steps, dtype=np.uint16)1129 value = np.zeros((n_steps, n_states), dtype=np.float64)1130 ptr = np.zeros((n_steps, n_states), dtype=np.uint16)1131 1132 # factor in initial state distribution1133 value[0] = log_prob[0] + log_p_init1134 1135 for t in range(1, n_steps):1136 # Want V[t, j] <- p[t, j] * max_k V[t-1, k] * A[k, j]1137 # assume at time t-1 we were in state k1138 # transition k -> j1139 1140 # Broadcast over rows:1141 # Tout[k, j] = V[t-1, k] * A[k, j]1142 # then take the max over columns1143 # We'll do this in log-space for stability1144 1145 trans_out = value[t - 1] + log_trans.T1146 1147 # Unroll the max/argmax loop to enable numba support1148 for j in range(n_states):1149 ptr[t, j] = np.argmax(trans_out[j])1150 # value[t, j] = log_prob[t, j] + np.max(trans_out[j])1151 value[t, j] = log_prob[t, j] + trans_out[j, ptr[t][j]]1152 1153 # Now roll backward1154 1155 # Get the last state1156 state[-1] = np.argmax(value[-1])1157 1158 for t in range(n_steps - 2, -1, -1):1159 state[t] = ptr[t + 1, state[t + 1]]1160 1161 logp = value[-1:, state[-1]]1162 1163 return state, logp1164 1165 1166@overload1167def viterbi(1168 prob: np.ndarray,1169 transition: np.ndarray,1170 *,1171 p_init: Optional[np.ndarray] = ...,1172 return_logp: Literal[True],1173) -> Tuple[np.ndarray, np.ndarray]:1174 ...1175 1176 1177@overload1178def viterbi(1179 prob: np.ndarray,1180 transition: np.ndarray,1181 *,1182 p_init: Optional[np.ndarray] = ...,1183 return_logp: Literal[False] = ...,1184) -> np.ndarray:1185 ...1186 1187 1188def viterbi(1189 prob: np.ndarray,1190 transition: np.ndarray,1191 *,1192 p_init: Optional[np.ndarray] = None,1193 return_logp: bool = False,1194) -> Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]:1195 """Viterbi decoding from observation likelihoods.1196 1197 Given a sequence of observation likelihoods ``prob[s, t]``,1198 indicating the conditional likelihood of seeing the observation1199 at time ``t`` from state ``s``, and a transition matrix1200 ``transition[i, j]`` which encodes the conditional probability of