<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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 -------- 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.
    # 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 /= .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)