def METHOD_NAME(
    train_features: np.ndarray,
    pred_features: np.ndarray,
    target_isnan_positions: np.ndarray,
    kernel: Callable[[np.ndarray, np.ndarray], float],
    do_exp: bool = True,
) -> Iterator[np.ndarray]:
    """
    Given the (logarithm of) kernel as well as training and prediction
    range features, this method returns an iterator over sampling weights
    for each time step in the prediction range.
    Note that the number of sampling weights for each time step vary since
    the prediction for time step `pred_t` samples from all the training
    targets as well as predictions until `pred_t` - 1.
    :return: iterator over sampling weights
    Parameters
    ----------
    train_features
        shape: (num_features, train_length)
    pred_features
        shape: (num_features, prediction_length)
    target_isnan_positions:
        an array of indices where the target is a NaN
    kernel
        kernel function that maps pairs of arrays to real numbers
    do_exp:
        exponentiate the weights in case of exponential kernel
        (for numerical stability we do this here)
    Returns
    -------
    iterator over sampling weights
    """
    assert len(np.shape(train_features)) == 2, (
        "Train features should be 2D-array where the rows represent "
        "features and columns the time points."
    )
    assert len(np.shape(pred_features)) == 2, (
        "Prediction features should be 2D-array where the rows represent "
        "features and columns the time points."
    )
    train_length = train_features.shape[1]
    prediction_length = pred_features.shape[1]
    for pred_t in range(prediction_length):
        # Prediction for `pred_t` samples from all the training targets
        # as well as predictions until `pred_t` - 1
        sampling_weights = np.zeros(train_length + pred_t)
        for train_t in range(train_length):
            sampling_weights[train_t] = kernel(
                train_features[:, train_t], pred_features[:, pred_t]
            )
        for t in range(pred_t):
            sampling_weights[train_length + t] = kernel(
                pred_features[:, t], pred_features[:, pred_t]
            )
        if do_exp:
            # To avoid numerical issues with exponentiation.
            sampling_weights -= max(sampling_weights)
            sampling_weights = np.exp(sampling_weights)
        # reset kernel at positions where the target is NaN
        sampling_weights[target_isnan_positions] = 0.0
        # Sometimes (e.g. for a for seasonal climatological kernel ) all
        # positions with non-zero probability are NaNs, so after resetting
        # the weights at these positions sampling_weights has only zeroes.
        # In this case, we want to sample uniformly from the observed
        # positions.
        if np.sum(sampling_weights) == 0:
            sampling_weights[target_isnan_positions] = -1.0
            sampling_weights += 1.0
        yield sampling_weights