beat.py

r"""
The aim of a beat detection algorithm is to report the times at which a typical
human listener might tap their foot to a piece of music. As a result, most
metrics for evaluating the performance of beat tracking systems involve
computing the error between the estimated beat times and some reference list of
beat locations. Many metrics additionally compare the beat sequences at
different metric levels in order to deal with the ambiguity of tempo.

Based on the methods described in:
    Matthew E. P. Davies,  Norberto Degara, and Mark D. Plumbley.
    "Evaluation Methods for Musical Audio Beat Tracking Algorithms",
    Queen Mary University of London Technical Report C4DM-TR-09-06
    London, United Kingdom, 8 October 2009.

See also the Beat Evaluation Toolbox:
    https://code.soundsoftware.ac.uk/projects/beat-evaluation/

Conventions
-----------

Beat times should be provided in the form of a 1-dimensional array of beat
times in seconds in increasing order.  Typically, any beats which occur before
5s are ignored; this can be accomplished using
:func:`mir_eval.beat.trim_beats()`.

Metrics
-------

* :func:`mir_eval.beat.f_measure`: The F-measure of the beat sequence, where an
  estimated beat is considered correct if it is sufficiently close to a
  reference beat
* :func:`mir_eval.beat.cemgil`: Cemgil's score, which computes the sum of
  Gaussian errors for each beat
* :func:`mir_eval.beat.goto`: Goto's score, a binary score which is 1 when at
  least 25\% of the estimated beat sequence closely matches the reference beat
  sequence
* :func:`mir_eval.beat.p_score`: McKinney's P-score, which computes the
  cross-correlation of the estimated and reference beat sequences represented
  as impulse trains
* :func:`mir_eval.beat.continuity`: Continuity-based scores which compute the
  proportion of the beat sequence which is continuously correct
* :func:`mir_eval.beat.information_gain`: The Information Gain of a normalized
  beat error histogram over a uniform distribution

"""

import numpy as np
import collections
from . import util
import warnings


# The maximum allowable beat time
MAX_TIME = 30000.0


def trim_beats(beats, min_beat_time=5.0):
    """Remove beats before min_beat_time.  A common preprocessing step.

    Parameters
    ----------
    beats : np.ndarray
        Array of beat times in seconds.
    min_beat_time : float
        Minimum beat time to allow
        (Default value = 5.)

    Returns
    -------
    beats_trimmed : np.ndarray
        Trimmed beat array.
    """
    # Remove beats before min_beat_time
    return beats[beats >= min_beat_time]


def validate(reference_beats, estimated_beats):
    """Check that the input annotations to a metric look like valid beat time
    arrays, and throws helpful errors if not.

    Parameters
    ----------
    reference_beats : np.ndarray
        reference beat times, in seconds
    estimated_beats : np.ndarray
        estimated beat times, in seconds
    """
    # If reference or estimated beats are empty,
    # warn because metric will be 0
    if reference_beats.size == 0:
        warnings.warn("Reference beats are empty.")
    if estimated_beats.size == 0:
        warnings.warn("Estimated beats are empty.")
    for beats in [reference_beats, estimated_beats]:
        util.validate_events(beats, MAX_TIME)


def _get_reference_beat_variations(reference_beats):
    """Return metric variations of the reference beats

    Parameters
    ----------
    reference_beats : np.ndarray
        beat locations in seconds

    Returns
    -------
    reference_beats : np.ndarray
        Original beat locations
    off_beat : np.ndarray
        180 degrees out of phase from the original beat locations
    double : np.ndarray
        Beats at 2x the original tempo
    half_odd : np.ndarray
        Half tempo, odd beats
    half_even : np.ndarray
        Half tempo, even beats
    """
    # Create annotations at twice the metric level
    interpolated_indices = np.arange(0, reference_beats.shape[0] - 0.5, 0.5)
    original_indices = np.arange(0, reference_beats.shape[0])
    double_reference_beats = np.interp(
        interpolated_indices, original_indices, reference_beats
    )
    # Return metric variations:
    # True, off-beat, double tempo, half tempo odd, and half tempo even
    return (
        reference_beats,
        double_reference_beats[1::2],
        double_reference_beats,
        reference_beats[::2],
        reference_beats[1::2],
    )


def f_measure(reference_beats, estimated_beats, f_measure_threshold=0.07):
    """Compute the F-measure of correct vs incorrectly predicted beats.
    "Correctness" is determined over a small window.

    Examples
    --------
    >>> reference_beats = mir_eval.io.load_events('reference.txt')
    >>> reference_beats = mir_eval.beat.trim_beats(reference_beats)
    >>> estimated_beats = mir_eval.io.load_events('estimated.txt')
    >>> estimated_beats = mir_eval.beat.trim_beats(estimated_beats)
    >>> f_measure = mir_eval.beat.f_measure(reference_beats,
                                            estimated_beats)

    Parameters
    ----------
    reference_beats : np.ndarray
        reference beat times, in seconds
    estimated_beats : np.ndarray
        estimated beat times, in seconds
    f_measure_threshold : float
        Window size, in seconds
        (Default value = 0.07)

    Returns
    -------
    f_score : float
        The computed F-measure score

    """
    validate(reference_beats, estimated_beats)
    # When estimated beats are empty, no beats are correct; metric is 0
    if estimated_beats.size == 0 or reference_beats.size == 0:
        return 0.0
    # Compute the best-case matching between reference and estimated locations
    matching = util.match_events(reference_beats, estimated_beats, f_measure_threshold)

    precision = float(len(matching)) / len(estimated_beats)
    recall = float(len(matching)) / len(reference_beats)
    return util.f_measure(precision, recall)


def cemgil(reference_beats, estimated_beats, cemgil_sigma=0.04):
    """Cemgil's score, computes a gaussian error of each estimated beat.
    Compares against the original beat times and all metrical variations.

    Examples
    --------
    >>> reference_beats = mir_eval.io.load_events('reference.txt')
    >>> reference_beats = mir_eval.beat.trim_beats(reference_beats)
    >>> estimated_beats = mir_eval.io.load_events('estimated.txt')
    >>> estimated_beats = mir_eval.beat.trim_beats(estimated_beats)
    >>> cemgil_score, cemgil_max = mir_eval.beat.cemgil(reference_beats,
                                                        estimated_beats)

    Parameters
    ----------
    reference_beats : np.ndarray
        reference beat times, in seconds
    estimated_beats : np.ndarray
        query beat times, in seconds
    cemgil_sigma : float
        Sigma parameter of gaussian error windows
        (Default value = 0.04)

    Returns
    -------
    cemgil_score : float
        Cemgil's score for the original reference beats
    cemgil_max : float
        The best Cemgil score for all metrical variations
    """
    validate(reference_beats, estimated_beats)
    # When estimated beats are empty, no beats are correct; metric is 0
    if estimated_beats.size == 0 or reference_beats.size == 0:
        return 0.0, 0.0
    # We'll compute Cemgil's accuracy for each variation
    accuracies = []
    for reference_beats in _get_reference_beat_variations(reference_beats):
        accuracy = 0
        # Cycle through beats
        for beat in reference_beats:
            # Find the error for the closest beat to the reference beat
            beat_diff = np.min(np.abs(beat - estimated_beats))
            # Add gaussian error into the accuracy
            accuracy += np.exp(-(beat_diff**2) / (2.0 * cemgil_sigma**2))
        # Normalize the accuracy
        accuracy /= 0.5 * (estimated_beats.shape[0] + reference_beats.shape[0])
        # Add it to our list of accuracy scores
        accuracies.append(accuracy)
    # Return raw accuracy with non-varied annotations
    # and maximal accuracy across all variations
    return accuracies[0], np.max(accuracies)


def goto(
    reference_beats, estimated_beats, goto_threshold=0.35, goto_mu=0.2, goto_sigma=0.2
):
    """Calculate Goto's score, a binary 1 or 0 depending on some specific
    heuristic criteria

    Examples
    --------
    >>> reference_beats = mir_eval.io.load_events('reference.txt')
    >>> reference_beats = mir_eval.beat.trim_beats(reference_beats)
    >>> estimated_beats = mir_eval.io.load_events('estimated.txt')
    >>> estimated_beats = mir_eval.beat.trim_beats(estimated_beats)
    >>> goto_score = mir_eval.beat.goto(reference_beats, estimated_beats)

    Parameters
    ----------
    reference_beats : np.ndarray
        reference beat times, in seconds
    estimated_beats : np.ndarray
        query beat times, in seconds
    goto_threshold : float
        Threshold of beat error for a beat to be "correct"
        (Default value = 0.35)
    goto_mu : float
        The mean of the beat errors in the continuously correct
        track must be less than this
        (Default value = 0.2)
    goto_sigma : float
        The std of the beat errors in the continuously correct track must
        be less than this
        (Default value = 0.2)

    Returns
    -------
    goto_score : float
        Either 1.0 or 0.0 if some specific criteria are met
    """
    validate(reference_beats, estimated_beats)
    # When estimated beats are empty, no beats are correct; metric is 0
    if estimated_beats.size == 0 or reference_beats.size == 0:
        return 0.0
    # Error for each beat
    beat_error = np.ones(reference_beats.shape[0])
    # Flag for whether the reference and estimated beats are paired
    paired = np.zeros(reference_beats.shape[0])
    # Keep track of Goto's three criteria
    goto_criteria = 0
    for n in range(1, reference_beats.shape[0] - 1):
        # Get previous inner-reference-beat-interval
        previous_interval = 0.5 * (reference_beats[n] - reference_beats[n - 1])
        # Window start - in the middle of the current beat and the previous
        window_min = reference_beats[n] - previous_interval
        # Next inter-reference-beat-interval
        next_interval = 0.5 * (reference_beats[n + 1] - reference_beats[n])
        # Window end - in the middle of the current beat and the next
        window_max = reference_beats[n] + next_interval
        # Get estimated beats in the window
        beats_in_window = np.logical_and(
            (estimated_beats >= window_min), (estimated_beats < window_max)
        )
        # False negative/positive
        if beats_in_window.sum() == 0 or beats_in_window.sum() > 1:
            paired[n] = 0
            beat_error[n] = 1
        else:
            # Single beat is paired!
            paired[n] = 1
            # Get offset of the estimated beat and the reference beat
            offset = estimated_beats[beats_in_window] - reference_beats[n]
            # Scale by previous or next interval
            if offset < 0:
                beat_error[n] = offset[0] / previous_interval
            else:
                beat_error[n] = offset[0] / next_interval
    # Get indices of incorrect beats
    incorrect_beats = np.flatnonzero(np.abs(beat_error) > goto_threshold)
    # All beats are correct (first and last will be 0 so always correct)
    if incorrect_beats.shape[0] < 3:
        # Get the track of correct beats
        track = beat_error[incorrect_beats[0] + 1 : incorrect_beats[-1] - 1]
        goto_criteria = 1
    else:
        # Get the track of maximal length
        track_len = np.max(np.diff(incorrect_beats))
        track_start = np.flatnonzero(np.diff(incorrect_beats) == track_len)[0]
        # Is the track length at least 25% of the song?
        if track_len - 1 > 0.25 * (reference_beats.shape[0] - 2):
            goto_criteria = 1
            start_beat = incorrect_beats[track_start]
            end_beat = incorrect_beats[track_start + 1]
            track = beat_error[start_beat : end_beat + 1]
    # If we have a track
    if goto_criteria:
        # Are mean and std of the track less than the required thresholds?
        if np.mean(np.abs(track)) < goto_mu and np.std(track, ddof=1) < goto_sigma:
            goto_criteria = 3
    # If all criteria are met, score is 100%!
    return 1.0 * (goto_criteria == 3)


def p_score(reference_beats, estimated_beats, p_score_threshold=0.2):
    """Get McKinney's P-score.
    Based on the autocorrelation of the reference and estimated beats

    Examples
    --------
    >>> reference_beats = mir_eval.io.load_events('reference.txt')
    >>> reference_beats = mir_eval.beat.trim_beats(reference_beats)
    >>> estimated_beats = mir_eval.io.load_events('estimated.txt')
    >>> estimated_beats = mir_eval.beat.trim_beats(estimated_beats)
    >>> p_score = mir_eval.beat.p_score(reference_beats, estimated_beats)

    Parameters
    ----------
    reference_beats : np.ndarray
        reference beat times, in seconds
    estimated_beats : np.ndarray
        query beat times, in seconds
    p_score_threshold : float
        Window size will be
        ``p_score_threshold*np.median(inter_annotation_intervals)``,
        (Default value = 0.2)

    Returns
    -------
    correlation : float
        McKinney's P-score

    """
    validate(reference_beats, estimated_beats)
    # Warn when only one beat is provided for either estimated or reference,
    # report a warning
    if reference_beats.size == 1:
        warnings.warn(
            "Only one reference beat was provided, so beat intervals"
            " cannot be computed."
        )
    if estimated_beats.size == 1:
        warnings.warn(
            "Only one estimated beat was provided, so beat intervals"
            " cannot be computed."
        )
    # When estimated or reference beats have <= 1 beats, can't compute the
    # metric, so return 0
    if estimated_beats.size <= 1 or reference_beats.size <= 1:
        return 0.0
    # Quantize beats to 10ms
    sampling_rate = int(1.0 / 0.010)
    # Shift beats so that the minimum in either sequence is zero
    offset = min(estimated_beats.min(), reference_beats.min())
    estimated_beats = np.array(estimated_beats - offset)
    reference_beats = np.array(reference_beats - offset)
    # Get the largest time index
    end_point = np.int64(
        np.ceil(np.max([np.max(estimated_beats), np.max(reference_beats)]))
    )
    # Make impulse trains with impulses at beat locations
    reference_train = np.zeros(end_point * sampling_rate + 1)
    beat_indices = np.ceil(reference_beats * sampling_rate).astype(np.int64)
    reference_train[beat_indices] = 1.0
    estimated_train = np.zeros(end_point * sampling_rate + 1)
    beat_indices = np.ceil(estimated_beats * sampling_rate).astype(np.int64)
    estimated_train[beat_indices] = 1.0
    # Window size to take the correlation over
    # defined as .2*median(inter-annotation-intervals)
    annotation_intervals = np.diff(np.flatnonzero(reference_train))
    win_size = int(np.round(p_score_threshold * np.median(annotation_intervals)))
    # Get full correlation
    train_correlation = np.correlate(reference_train, estimated_train, "full")
    # Get the middle element - note we are rounding down on purpose here
    middle_lag = train_correlation.shape[0] // 2
    # Truncate to only valid lags (those corresponding to the window)
    start = middle_lag - win_size
    end = middle_lag + win_size + 1
    train_correlation = train_correlation[start:end]
    # Compute and return the P-score
    n_beats = np.max([estimated_beats.shape[0], reference_beats.shape[0]])
    return np.sum(train_correlation) / n_beats


def continuity(
    reference_beats,
    estimated_beats,
    continuity_phase_threshold=0.175,
    continuity_period_threshold=0.175,
):
    """Get metrics based on how much of the estimated beat sequence is
    continually correct.

    Examples
    --------
    >>> reference_beats = mir_eval.io.load_events('reference.txt')
    >>> reference_beats = mir_eval.beat.trim_beats(reference_beats)
    >>> estimated_beats = mir_eval.io.load_events('estimated.txt')
    >>> estimated_beats = mir_eval.beat.trim_beats(estimated_beats)
    >>> CMLc, CMLt, AMLc, AMLt = mir_eval.beat.continuity(reference_beats,
                                                          estimated_beats)

    Parameters
    ----------
    reference_beats : np.ndarray
        reference beat times, in seconds
    estimated_beats : np.ndarray
        query beat times, in seconds
    continuity_phase_threshold : float
        Allowable ratio of how far is the estimated beat
        can be from the reference beat
        (Default value = 0.175)
    continuity_period_threshold : float
        Allowable distance between the inter-beat-interval
        and the inter-annotation-interval
        (Default value = 0.175)

    Returns
    -------
    CMLc : float
        Correct metric level, continuous accuracy
    CMLt : float
        Correct metric level, total accuracy (continuity not required)
    AMLc : float
        Any metric level, continuous accuracy
    AMLt : float
        Any metric level, total accuracy (continuity not required)
    """
    validate(reference_beats, estimated_beats)
    # Warn when only one beat is provided for either estimated or reference,
    # report a warning
    if reference_beats.size == 1:
        warnings.warn(
            "Only one reference beat was provided, so beat intervals"
            " cannot be computed."
        )
    if estimated_beats.size == 1:
        warnings.warn(
            "Only one estimated beat was provided, so beat intervals"
            " cannot be computed."
        )
    # When estimated or reference beats have <= 1 beats, can't compute the
    # metric, so return 0
    if estimated_beats.size <= 1 or reference_beats.size <= 1:
        return 0.0, 0.0, 0.0, 0.0
    # Accuracies for each variation
    continuous_accuracies = []
    total_accuracies = []
    # Get accuracy for each variation
    for reference_beats in _get_reference_beat_variations(reference_beats):
        # Annotations that have been used
        n_annotations = np.max([reference_beats.shape[0], estimated_beats.shape[0]])
        used_annotations = np.zeros(n_annotations)
        # Whether or not we are continuous at any given point
        beat_successes = np.zeros(n_annotations)
        for m in range(estimated_beats.shape[0]):
            # Is this beat correct?
            beat_success = 0
            # Get differences for this beat
            beat_differences = np.abs(estimated_beats[m] - reference_beats)
            # Get nearest annotation index
            nearest = np.argmin(beat_differences)
            min_difference = beat_differences[nearest]
            # Have we already used this annotation?
            if used_annotations[nearest] == 0:
                # Is this the first beat or first annotation?
                # If so, look forward.
                if m == 0 or nearest == 0:
                    # How far is the estimated beat from the reference beat,
                    # relative to the inter-annotation-interval?
                    if nearest + 1 < reference_beats.shape[0]:
                        reference_interval = (
                            reference_beats[nearest + 1] - reference_beats[nearest]
                        )
                    else:
                        # Special case when nearest + 1 is too large - use the
                        # previous interval instead
                        reference_interval = (
                            reference_beats[nearest] - reference_beats[nearest - 1]
                        )
                    # Handle this special case when beats are not unique
                    if reference_interval == 0:
                        if min_difference == 0:
                            phase = 1
                        else:
                            phase = np.inf
                    else:
                        phase = np.abs(min_difference / reference_interval)
                    # How close is the inter-beat-interval
                    # to the inter-annotation-interval?
                    if m + 1 < estimated_beats.shape[0]:
                        estimated_interval = estimated_beats[m + 1] - estimated_beats[m]
                    else:
                        # Special case when m + 1 is too large - use the
                        # previous interval
                        estimated_interval = estimated_beats[m] - estimated_beats[m - 1]
                    # Handle this special case when beats are not unique
                    if reference_interval == 0:
                        if estimated_interval == 0:
                            period = 0
                        else:
                            period = np.inf
                    else:
                        period = np.abs(1 - estimated_interval / reference_interval)
                    if (
                        phase < continuity_phase_threshold
                        and period < continuity_period_threshold
                    ):
                        # Set this annotation as used
                        used_annotations[nearest] = 1
                        # This beat is matched
                        beat_success = 1
                # This beat/annotation is not the first
                else:
                    # How far is the estimated beat from the reference beat,
                    # relative to the inter-annotation-interval?
                    reference_interval = (
                        reference_beats[nearest] - reference_beats[nearest - 1]
                    )
                    phase = np.abs(min_difference / reference_interval)
                    # How close is the inter-beat-interval
                    # to the inter-annotation-interval?
                    estimated_interval = estimated_beats[m] - estimated_beats[m - 1]
                    reference_interval = (
                        reference_beats[nearest] - reference_beats[nearest - 1]
                    )
                    period = np.abs(1 - estimated_interval / reference_interval)
                    if (
                        phase < continuity_phase_threshold
                        and period < continuity_period_threshold
                    ):
                        # Set this annotation as used
                        used_annotations[nearest] = 1
                        # This beat is matched
                        beat_success = 1
            # Set whether this beat is matched or not
            beat_successes[m] = beat_success
        # Add 0s at the begnning and end
        # so that we at least find the beginning/end of the estimated beats
        beat_successes = np.append(np.append(0, beat_successes), 0)
        # Where is the beat not a match?
        beat_failures = np.nonzero(beat_successes == 0)[0]
        # Take out those zeros we added
        beat_successes = beat_successes[1:-1]
        # Get the continuous accuracy as the longest track of successful beats
        longest_track = np.max(np.diff(beat_failures)) - 1
        continuous_accuracy = longest_track / (1.0 * beat_successes.shape[0])
        continuous_accuracies.append(continuous_accuracy)
        # Get the total accuracy - all sequences
        total_accuracy = np.sum(beat_successes) / (1.0 * beat_successes.shape[0])
        total_accuracies.append(total_accuracy)
    # Grab accuracy scores
    return (
        continuous_accuracies[0],
        total_accuracies[0],
        np.max(continuous_accuracies),
        np.max(total_accuracies),
    )


def information_gain(reference_beats, estimated_beats, bins=41):
    """Get the information gain - K-L divergence of the beat error histogram
    to a uniform histogram

    Examples
    --------
    >>> reference_beats = mir_eval.io.load_events('reference.txt')
    >>> reference_beats = mir_eval.beat.trim_beats(reference_beats)
    >>> estimated_beats = mir_eval.io.load_events('estimated.txt')
    >>> estimated_beats = mir_eval.beat.trim_beats(estimated_beats)
    >>> information_gain = mir_eval.beat.information_gain(reference_beats,
                                                          estimated_beats)

    Parameters
    ----------
    reference_beats : np.ndarray
        reference beat times, in seconds
    estimated_beats : np.ndarray
        query beat times, in seconds
    bins : int
        Number of bins in the beat error histogram
        (Default value = 41)

    Returns
    -------
    information_gain_score : float
        Entropy of beat error histogram
    """
    validate(reference_beats, estimated_beats)
    # If an even number of bins is provided,
    # there will be no bin centered at zero, so warn the user.
    if not bins % 2:
        warnings.warn(
            "bins parameter is even, " "so there will not be a bin centered at zero."
        )
    # Warn when only one beat is provided for either estimated or reference,
    # report a warning
    if reference_beats.size == 1:
        warnings.warn(
            "Only one reference beat was provided, so beat intervals"
            " cannot be computed."
        )
    if estimated_beats.size == 1:
        warnings.warn(
            "Only one estimated beat was provided, so beat intervals"
            " cannot be computed."
        )
    # When estimated or reference beats have <= 1 beats, can't compute the
    # metric, so return 0
    if estimated_beats.size <= 1 or reference_beats.size <= 1:
        return 0.0
    # Get entropy for reference beats->estimated beats
    # and estimated beats->reference beats
    forward_entropy = _get_entropy(reference_beats, estimated_beats, bins)
    backward_entropy = _get_entropy(estimated_beats, reference_beats, bins)
    # Pick the larger of the entropies
    norm = np.log2(bins)
    if forward_entropy > backward_entropy:
        # Note that the beat evaluation toolbox does not normalize
        information_gain_score = (norm - forward_entropy) / norm
    else:
        information_gain_score = (norm - backward_entropy) / norm
    return information_gain_score


def _get_entropy(reference_beats, estimated_beats, bins):
    """Compute the entropy of the beat error histogram.

    This is a helper function for the information gain
    metric, and needs to be run twice: once backwards, once
    forwards.

    Parameters
    ----------
    reference_beats : np.ndarray
        reference beat times, in seconds
    estimated_beats : np.ndarray
        query beat times, in seconds
    bins : int
        Number of bins in the beat error histogram

    Returns
    -------
    entropy : float
        Entropy of beat error histogram
    """
    beat_error = np.zeros(estimated_beats.shape[0])
    for n in range(estimated_beats.shape[0]):
        # Get index of closest annotation to this beat
        beat_distances = estimated_beats[n] - reference_beats
        closest_beat = np.argmin(np.abs(beat_distances))
        absolute_error = beat_distances[closest_beat]
        # If the first annotation is closest...
        if closest_beat == 0:
            # Inter-annotation interval - space between first two beats
            interval = 0.5 * (reference_beats[1] - reference_beats[0])
        # If last annotation is closest...
        if closest_beat == (reference_beats.shape[0] - 1):
            interval = 0.5 * (reference_beats[-1] - reference_beats[-2])
        else:
            if absolute_error < 0:
                # Closest annotation is the one before the current beat
                # so look at previous inner-annotation-interval
                start = reference_beats[closest_beat]
                end = reference_beats[closest_beat - 1]
                interval = 0.5 * (start - end)
            else:
                # Closest annotation is the one after the current beat
                # so look at next inner-annotation-interval
                start = reference_beats[closest_beat + 1]
                end = reference_beats[closest_beat]
                interval = 0.5 * (start - end)
        # The actual error of this beat
        beat_error[n] = 0.5 * absolute_error / interval
    # Put beat errors in range (-.5, .5)
    beat_error = np.mod(beat_error + 0.5, -1) + 0.5
    # Note these are slightly different the beat evaluation toolbox
    # (they are uniform)
    histogram_bin_edges = np.linspace(-0.5, 0.5, bins + 1)
    # Get the histogram
    raw_bin_values = np.histogram(beat_error, histogram_bin_edges)[0]
    # Turn into a proper probability distribution
    raw_bin_values = raw_bin_values / (1.0 * np.sum(raw_bin_values))
    # Set zero-valued bins to 1 to make the entropy calculation well-behaved
    raw_bin_values[raw_bin_values == 0] = 1
    # Calculate entropy
    return -np.sum(raw_bin_values * np.log2(raw_bin_values))


def evaluate(reference_beats, estimated_beats, **kwargs):
    """Compute all metrics for the given reference and estimated annotations.

    Examples
    --------
    >>> reference_beats = mir_eval.io.load_events('reference.txt')
    >>> estimated_beats = mir_eval.io.load_events('estimated.txt')
    >>> scores = mir_eval.beat.evaluate(reference_beats, estimated_beats)

    Parameters
    ----------
    reference_beats : np.ndarray
        Reference beat times, in seconds
    estimated_beats : np.ndarray
        Query beat times, in seconds
    **kwargs
        Additional keyword arguments which will be passed to the
        appropriate metric or preprocessing functions.

    Returns
    -------
    scores : dict
        Dictionary of scores, where the key is the metric name (str) and
        the value is the (float) score achieved.

    """
    # Trim beat times at the beginning of the annotations
    reference_beats = util.filter_kwargs(trim_beats, reference_beats, **kwargs)
    estimated_beats = util.filter_kwargs(trim_beats, estimated_beats, **kwargs)

    # Now compute all the metrics

    scores = collections.OrderedDict()

    # F-Measure
    scores["F-measure"] = util.filter_kwargs(
        f_measure, reference_beats, estimated_beats, **kwargs
    )

    # Cemgil
    scores["Cemgil"], scores["Cemgil Best Metric Level"] = util.filter_kwargs(
        cemgil, reference_beats, estimated_beats, **kwargs
    )

    # Goto
    scores["Goto"] = util.filter_kwargs(
        goto, reference_beats, estimated_beats, **kwargs
    )

    # P-Score
    scores["P-score"] = util.filter_kwargs(
        p_score, reference_beats, estimated_beats, **kwargs
    )

    # Continuity metrics
    (
        scores["Correct Metric Level Continuous"],
        scores["Correct Metric Level Total"],
        scores["Any Metric Level Continuous"],
        scores["Any Metric Level Total"],
    ) = util.filter_kwargs(continuity, reference_beats, estimated_beats, **kwargs)

    # Information gain
    scores["Information gain"] = util.filter_kwargs(
        information_gain, reference_beats, estimated_beats, **kwargs
    )

    return scores