librosa/sequence.py
def viterbi_discriminative(prob, transition, p_state=None, p_init=None, return_logp=False):
    '''Viterbi decoding from discriminative state predictions.

    Given a sequence of conditional state predictions `prob[s, t]`,
    indicating the conditional likelihood of state `s` given the
    observation at time `t`, and a transition matrix `transition[i, j]`
    which encodes the conditional probability of moving from state `i`
    to state `j`, the Viterbi algorithm computes the most likely sequence
    of states from the observations.

    This implementation uses the standard Viterbi decoding algorithm
    for observation likelihood sequences, under the assumption that
    `P[Obs(t) | State(t) = s]` is proportional to
    `P[State(t) = s | Obs(t)] / P[State(t) = s]`, where the denominator
    is the marginal probability of state `s` occurring as given by `p_state`.

    Parameters
    ----------
    prob : np.ndarray [shape=(n_states, n_steps), non-negative]
        `prob[s, t]` is the probability of state `s` conditional on
        the observation at time `t`.
        Must be non-negative and sum to 1 along each column.

    transition : np.ndarray [shape=(n_states, n_states), non-negative]
        `transition[i, j]` is the probability of a transition from i->j.
        Each row must sum to 1.

    p_state : np.ndarray [shape=(n_states,)]
        Optional: marginal probability distribution over states,
        must be non-negative and sum to 1.
        If not provided, a uniform distribution is assumed.

    p_init : np.ndarray [shape=(n_states,)]
        Optional: initial state distribution.
        If not provided, it is assumed to be uniform.

    return_logp : bool
        If `True`, return the log-likelihood of the state sequence.

    Returns
    -------
    Either `states` or `(states, logp)`:

    states : np.ndarray [shape=(n_steps,)]
        The most likely state sequence.

    logp : scalar [float]
        If `return_logp=True`, the log probability of `states` given
        the observations.

    See Also
    --------
    viterbi : Viterbi decoding from observation likelihoods
    viterbi_binary: Viterbi decoding for multi-label, conditional state likelihoods

    Examples
    --------
    This example constructs a simple, template-based discriminative chord estimator,
    using CENS chroma as input features.

    .. note:: this chord model is not accurate enough to use in practice. It is only
            intended to demonstrate how to use discriminative Viterbi decoding.

    >>> # Create templates for major, minor, and no-chord qualities
    >>> maj_template = np.array([1,0,0, 0,1,0, 0,1,0, 0,0,0])
    >>> min_template = np.array([1,0,0, 1,0,0, 0,1,0, 0,0,0])
    >>> N_template   = np.array([1,1,1, 1,1,1, 1,1,1, 1,1,1.]) / 4.
    >>> # Generate the weighting matrix that maps chroma to labels
    >>> weights = np.zeros((25, 12), dtype=float)
    >>> labels = ['C:maj', 'C#:maj', 'D:maj', 'D#:maj', 'E:maj', 'F:maj',
    ...           'F#:maj', 'G:maj', 'G#:maj', 'A:maj', 'A#:maj', 'B:maj',
    ...           'C:min', 'C#:min', 'D:min', 'D#:min', 'E:min', 'F:min',
    ...           'F#:min', 'G:min', 'G#:min', 'A:min', 'A#:min', 'B:min',
    ...           'N']
    >>> for c in range(12):
    ...     weights[c, :] = np.roll(maj_template, c) # c:maj
    ...     weights[c + 12, :] = np.roll(min_template, c)  # c:min
    >>> weights[-1] = N_template  # the last row is the no-chord class
    >>> # Make a self-loop transition matrix over 25 states
    >>> trans = librosa.sequence.transition_loop(25, 0.9)

    >>> # Load in audio and make features
    >>> y, sr = librosa.load(librosa.util.example_audio_file())
    >>> chroma = librosa.feature.chroma_cens(y=y, sr=sr, bins_per_octave=36)
    >>> # Map chroma (observations) to class (state) likelihoods
    >>> probs = np.exp(weights.dot(chroma))  # P[class | chroma] proportional to exp(template' chroma)
    >>> probs /= probs.sum(axis=0, keepdims=True)  # probabilities must sum to 1 in each column
    >>> # Compute independent frame-wise estimates
    >>> chords_ind = np.argmax(probs, axis=0)
    >>> # And viterbi estimates
    >>> chords_vit = librosa.sequence.viterbi_discriminative(probs, trans)

    >>> # Plot the features and prediction map
    >>> import matplotlib.pyplot as plt
    >>> plt.figure(figsize=(10, 6))
    >>> plt.subplot(2,1,1)
    >>> librosa.display.specshow(chroma, x_axis='time', y_axis='chroma')
    >>> plt.colorbar()
    >>> plt.subplot(2,1,2)
    >>> librosa.display.specshow(weights, x_axis='chroma')
    >>> plt.yticks(np.arange(25) + 0.5, labels)
    >>> plt.ylabel('Chord')
    >>> plt.colorbar()
    >>> plt.tight_layout()

    >>> # And plot the results
    >>> plt.figure(figsize=(10, 4))
    >>> librosa.display.specshow(probs, x_axis='time', cmap='gray')
    >>> plt.colorbar()
    >>> times = librosa.frames_to_time(np.arange(len(chords_vit)))
    >>> plt.scatter(times, chords_ind + 0.75, color='lime', alpha=0.5, marker='+', s=15, label='Independent')
    >>> plt.scatter(times, chords_vit + 0.25, color='deeppink', alpha=0.5, marker='o', s=15, label='Viterbi')
    >>> plt.yticks(0.5 + np.unique(chords_vit), [labels[i] for i in np.unique(chords_vit)], va='center')
    >>> plt.legend(loc='best')
    >>> plt.tight_layout()

    '''

    n_states, n_steps = prob.shape

    if transition.shape != (n_states, n_states):
        raise ParameterError('transition.shape={}, must be '
                             '(n_states, n_states)={}'.format(transition.shape,
                                                              (n_states, n_states)))

    if np.any(transition < 0) or not np.allclose(transition.sum(axis=1), 1):
        raise ParameterError('Invalid transition matrix: must be non-negative '
                             'and sum to 1 on each row.')

    if np.any(prob < 0) or not np.allclose(prob.sum(axis=0), 1):
        raise ParameterError('Invalid probability values: each column must '
                             'sum to 1 and be non-negative')

    states = np.zeros(n_steps, dtype=int)
    values = np.zeros((n_steps, n_states), dtype=float)
    ptr = np.zeros((n_steps, n_states), dtype=int)

    # Compute log-likelihoods while avoiding log-underflow
    epsilon = np.finfo(prob.dtype).tiny

    # Compute marginal log probabilities while avoiding underflow
    if p_state is None:
        p_state = np.empty(n_states)
        p_state.fill(1./n_states)
    elif p_state.shape != (n_states,):
        raise ParameterError('Marginal distribution p_state must have shape (n_states,). '
                             'Got p_state.shape={}'.format(p_state.shape))
    elif np.any(p_state < 0) or not np.allclose(p_state.sum(axis=-1), 1):
        raise ParameterError('Invalid marginal state distribution: '
                             'p_state={}'.format(p_state))

    log_trans = np.log(transition + epsilon)
    log_marginal = np.log(p_state + epsilon)

    # By Bayes' rule, P[X | Y] * P[Y] = P[Y | X] * P[X]
    # P[X] is constant for the sake of maximum likelihood inference
    # and P[Y] is given by the marginal distribution p_state.
    #
    # So we have P[X | y] \propto P[Y | x] / P[Y]
    # if X = observation and Y = states, this can be done in log space as
    # log P[X | y] \propto \log P[Y | x] - \log P[Y]
    log_prob = np.log(prob.T + epsilon) - log_marginal

    if p_init is None:
        p_init = np.empty(n_states)
        p_init.fill(1./n_states)
    elif np.any(p_init < 0) or not np.allclose(p_init.sum(), 1):
        raise ParameterError('Invalid initial state distribution: '
                             'p_init={}'.format(p_init))

    log_p_init = np.log(p_init + epsilon)

    _viterbi(log_prob, log_trans, log_p_init, states, values, ptr)

    if return_logp:
        return states, values[-1, states[-1]]

    return states
def viterbi_discriminative(prob, transition, p_state=None, p_init=None, return_logp=False):
    '''Viterbi decoding from discriminative state predictions.

    Given a sequence of conditional state predictions `prob[s, t]`,
    indicating the conditional likelihood of state `s` given the
    observation at time `t`, and a transition matrix `transition[i, j]`
    which encodes the conditional probability of moving from state `i`
    to state `j`, the Viterbi algorithm computes the most likely sequence
    of states from the observations.

    This implementation uses the standard Viterbi decoding algorithm
    for observation likelihood sequences, under the assumption that
    `P[Obs(t) | State(t) = s]` is proportional to
    `P[State(t) = s | Obs(t)] / P[State(t) = s]`, where the denominator
    is the marginal probability of state `s` occurring as given by `p_state`.

    Parameters
    ----------
    prob : np.ndarray [shape=(n_states, n_steps), non-negative]
        `prob[s, t]` is the probability of state `s` conditional on
        the observation at time `t`.
        Must be non-negative and sum to 1 along each column.

    transition : np.ndarray [shape=(n_states, n_states), non-negative]
        `transition[i, j]` is the probability of a transition from i->j.
        Each row must sum to 1.

    p_state : np.ndarray [shape=(n_states,)]
        Optional: marginal probability distribution over states,
        must be non-negative and sum to 1.
        If not provided, a uniform distribution is assumed.

    p_init : np.ndarray [shape=(n_states,)]
        Optional: initial state distribution.
        If not provided, it is assumed to be uniform.

    return_logp : bool
        If `True`, return the log-likelihood of the state sequence.

    Returns
    -------
    Either `states` or `(states, logp)`:

    states : np.ndarray [shape=(n_steps,)]
        The most likely state sequence.

    logp : scalar [float]
        If `return_logp=True`, the log probability of `states` given
        the observations.

    See Also
    --------
    viterbi : Viterbi decoding from observation likelihoods
    viterbi_binary: Viterbi decoding for multi-label, conditional state likelihoods

    Examples
    --------
    This example constructs a simple, template-based discriminative chord estimator,
    using CENS chroma as input features.

    .. note:: this chord model is not accurate enough to use in practice. It is only
            intended to demonstrate how to use discriminative Viterbi decoding.

    >>> # Create templates for major, minor, and no-chord qualities
    >>> maj_template = np.array([1,0,0, 0,1,0, 0,1,0, 0,0,0])
    >>> min_template = np.array([1,0,0, 1,0,0, 0,1,0, 0,0,0])
    >>> N_template   = np.array([1,1,1, 1,1,1, 1,1,1, 1,1,1.]) / 4.
    >>> # Generate the weighting matrix that maps chroma to labels
    >>> weights = np.zeros((25, 12), dtype=float)
    >>> labels = ['C:maj', 'C#:maj', 'D:maj', 'D#:maj', 'E:maj', 'F:maj',
    ...           'F#:maj', 'G:maj', 'G#:maj', 'A:maj', 'A#:maj', 'B:maj',
    ...           'C:min', 'C#:min', 'D:min', 'D#:min', 'E:min', 'F:min',
    ...           'F#:min', 'G:min', 'G#:min', 'A:min', 'A#:min', 'B:min',
    ...           'N']
    >>> for c in range(12):
    ...     weights[c, :] = np.roll(maj_template, c) # c:maj
    ...     weights[c + 12, :] = np.roll(min_template, c)  # c:min
    >>> weights[-1] = N_template  # the last row is the no-chord class
    >>> # Make a self-loop transition matrix over 25 states
    >>> trans = librosa.sequence.transition_loop(25, 0.9)

    >>> # Load in audio and make features
    >>> y, sr = librosa.load(librosa.util.example_audio_file())
    >>> chroma = librosa.feature.chroma_cens(y=y, sr=sr, bins_per_octave=36)
    >>> # Map chroma (observations) to class (state) likelihoods
    >>> probs = np.exp(weights.dot(chroma))  # P[class | chroma] proportional to exp(template' chroma)
    >>> probs /= probs.sum(axis=0, keepdims=True)  # probabilities must sum to 1 in each column
    >>> # Compute independent frame-wise estimates
    >>> chords_ind = np.argmax(probs, axis=0)
    >>> # And viterbi estimates
    >>> chords_vit = librosa.sequence.viterbi_discriminative(probs, trans)

    >>> # Plot the features and prediction map
    >>> import matplotlib.pyplot as plt
    >>> plt.figure(figsize=(10, 6))
    >>> plt.subplot(2,1,1)
    >>> librosa.display.specshow(chroma, x_axis='time', y_axis='chroma')
    >>> plt.colorbar()
    >>> plt.subplot(2,1,2)
    >>> librosa.display.specshow(weights, x_axis='chroma')
    >>> plt.yticks(np.arange(25) + 0.5, labels)
    >>> plt.ylabel('Chord')
    >>> plt.colorbar()
    >>> plt.tight_layout()

    >>> # And plot the results
    >>> plt.figure(figsize=(10, 4))
    >>> librosa.display.specshow(probs, x_axis='time', cmap='gray')
    >>> plt.colorbar()
    >>> times = librosa.frames_to_time(np.arange(len(chords_vit)))
    >>> plt.scatter(times, chords_ind + 0.75, color='lime', alpha=0.5, marker='+', s=15, label='Independent')
    >>> plt.scatter(times, chords_vit + 0.25, color='deeppink', alpha=0.5, marker='o', s=15, label='Viterbi')
    >>> plt.yticks(0.5 + np.unique(chords_vit), [labels[i] for i in np.unique(chords_vit)], va='center')
    >>> plt.legend(loc='best')
    >>> plt.tight_layout()

    '''

    n_states, n_steps = prob.shape

    if transition.shape != (n_states, n_states):
        raise ParameterError('transition.shape={}, must be '
                             '(n_states, n_states)={}'.format(transition.shape,
                                                              (n_states, n_states)))

    if np.any(transition < 0) or not np.allclose(transition.sum(axis=1), 1):
        raise ParameterError('Invalid transition matrix: must be non-negative '
                             'and sum to 1 on each row.')

    if np.any(prob < 0) or not np.allclose(prob.sum(axis=0), 1):
        raise ParameterError('Invalid probability values: each column must '
                             'sum to 1 and be non-negative')

    states = np.zeros(n_steps, dtype=int)
    values = np.zeros((n_steps, n_states), dtype=float)
    ptr = np.zeros((n_steps, n_states), dtype=int)

    # Compute log-likelihoods while avoiding log-underflow
    epsilon = np.finfo(prob.dtype).tiny

    # Compute marginal log probabilities while avoiding underflow
    if p_state is None:
        p_state = np.empty(n_states)
        p_state.fill(1./n_states)
    elif p_state.shape != (n_states,):
        raise ParameterError('Marginal distribution p_state must have shape (n_states,). '
                             'Got p_state.shape={}'.format(p_state.shape))
    elif np.any(p_state < 0) or not np.allclose(p_state.sum(axis=-1), 1):
        raise ParameterError('Invalid marginal state distribution: '
                             'p_state={}'.format(p_state))

    log_trans = np.log(transition + epsilon)
    log_marginal = np.log(p_state + epsilon)

    # By Bayes' rule, P[X | Y] * P[Y] = P[Y | X] * P[X]
    # P[X] is constant for the sake of maximum likelihood inference
    # and P[Y] is given by the marginal distribution p_state.
    #
    # So we have P[X | y] \propto P[Y | x] / P[Y]
    # if X = observation and Y = states, this can be done in log space as
    # log P[X | y] \propto \log P[Y | x] - \log P[Y]
    log_prob = np.log(prob.T + epsilon) - log_marginal

    if p_init is None:
        p_init = np.empty(n_states)
        p_init.fill(1./n_states)
    elif np.any(p_init < 0) or not np.allclose(p_init.sum(), 1):
        raise ParameterError('Invalid initial state distribution: '
                             'p_init={}'.format(p_init))

    log_p_init = np.log(p_init + epsilon)

    _viterbi(log_prob, log_trans, log_p_init, states, values, ptr)

    if return_logp:
        return states, values[-1, states[-1]]

    return states
Viterbi decoding from discriminative state predictions.

    Given a sequence of conditional state predictions `prob[s, t]`,
    indicating the conditional likelihood of state `s` given the
    observation at time `t`, and a transition matrix `transition[i, j]`
    which encodes the conditional probability of moving from state `i`
    to state `j`, the Viterbi algorithm computes the most likely sequence
    of states from the observations.

    This implementation uses the standard Viterbi decoding algorithm
    for observation likelihood sequences, under the assumption that
    `P[Obs(t) | State(t) = s]` is proportional to
    `P[State(t) = s | Obs(t)] / P[State(t) = s]`, where the denominator
    is the marginal probability of state `s` occurring as given by `p_state`.

    Parameters
    ----------
    prob : np.ndarray [shape=(n_states, n_steps), non-negative]
        `prob[s, t]` is the probability of state `s` conditional on
        the observation at time `t`.
        Must be non-negative and sum to 1 along each column.

    transition : np.ndarray [shape=(n_states, n_states), non-negative]
        `transition[i, j]` is the probability of a transition from i->j.
        Each row must sum to 1.

    p_state : np.ndarray [shape=(n_states,)]
        Optional: marginal probability distribution over states,
        must be non-negative and sum to 1.
        If not provided, a uniform distribution is assumed.

    p_init : np.ndarray [shape=(n_states,)]
        Optional: initial state distribution.
        If not provided, it is assumed to be uniform.

    return_logp : bool
        If `True`, return the log-likelihood of the state sequence.

    Returns
    -------
    Either `states` or `(states, logp)`:

    states : np.ndarray [shape=(n_steps,)]
        The most likely state sequence.

    logp : scalar [float]
        If `return_logp=True`, the log probability of `states` given
        the observations.

    See Also
    --------
    viterbi : Viterbi decoding from observation likelihoods
    viterbi_binary: Viterbi decoding for multi-label, conditional state likelihoods

    Examples
    --------
    This example constructs a simple, template-based discriminative chord estimator,
    using CENS chroma as input features.

    .. note:: this chord model is not accurate enough to use in practice. It is only
            intended to demonstrate how to use discriminative Viterbi decoding.

    >>> # Create templates for major, minor, and no-chord qualities
    >>> maj_template = np.array([1,0,0, 0,1,0, 0,1,0, 0,0,0])
    >>> min_template = np.array([1,0,0, 1,0,0, 0,1,0, 0,0,0])
    >>> N_template   = np.array([1,1,1, 1,1,1, 1,1,1, 1,1,1.]) / 4.
    >>> # Generate the weighting matrix that maps chroma to labels
    >>> weights = np.zeros((25, 12), dtype=float)
    >>> labels = ['C:maj', 'C#:maj', 'D:maj', 'D#:maj', 'E:maj', 'F:maj',
    ...           'F#:maj', 'G:maj', 'G#:maj', 'A:maj', 'A#:maj', 'B:maj',
    ...           'C:min', 'C#:min', 'D:min', 'D#:min', 'E:min', 'F:min',
    ...           'F#:min', 'G:min', 'G#:min', 'A:min', 'A#:min', 'B:min',
    ...           'N']
    >>> for c in range(12):
    ...     weights[c, :] = np.roll(maj_template, c) # c:maj
    ...     weights[c + 12, :] = np.roll(min_template, c)  # c:min
    >>> weights[-1] = N_template  # the last row is the no-chord class
    >>> # Make a self-loop transition matrix over 25 states
    >>> trans = librosa.sequence.transition_loop(25, 0.9)

    >>> # Load in audio and make features
    >>> y, sr = librosa.load(librosa.util.example_audio_file())
    >>> chroma = librosa.feature.chroma_cens(y=y, sr=sr, bins_per_octave=36)
    >>> # Map chroma (observations) to class (state) likelihoods
    >>> probs = np.exp(weights.dot(chroma))  # P[class | chroma] proportional to exp(template' chroma)
    >>> probs /= probs.sum(axis=0, keepdims=True)  # probabilities must sum to 1 in each column
    >>> # Compute independent frame-wise estimates
    >>> chords_ind = np.argmax(probs, axis=0)
    >>> # And viterbi estimates
    >>> chords_vit = librosa.sequence.viterbi_discriminative(probs, trans)

    >>> # Plot the features and prediction map
    >>> import matplotlib.pyplot as plt
    >>> plt.figure(figsize=(10, 6))
    >>> plt.subplot(2,1,1)
    >>> librosa.display.specshow(chroma, x_axis='time', y_axis='chroma')
    >>> plt.colorbar()
    >>> plt.subplot(2,1,2)
    >>> librosa.display.specshow(weights, x_axis='chroma')
    >>> plt.yticks(np.arange(25) + 0.5, labels)
    >>> plt.ylabel('Chord')
    >>> plt.colorbar()
    >>> plt.tight_layout()

    >>> # And plot the results
    >>> plt.figure(figsize=(10, 4))
    >>> librosa.display.specshow(probs, x_axis='time', cmap='gray')
    >>> plt.colorbar()
    >>> times = librosa.frames_to_time(np.arange(len(chords_vit)))
    >>> plt.scatter(times, chords_ind + 0.75, color='lime', alpha=0.5, marker='+', s=15, label='Independent')
    >>> plt.scatter(times, chords_vit + 0.25, color='deeppink', alpha=0.5, marker='o', s=15, label='Viterbi')
    >>> plt.yticks(0.5 + np.unique(chords_vit), [labels[i] for i in np.unique(chords_vit)], va='center')
    >>> plt.legend(loc='best')
    >>> plt.tight_layout()