def METHOD_NAME(
    nbins,
    state_labels,
    weights,
    parent_ids,
    bin_assignments,
    label_assignments,
    state_map,
    labeled_pops,
    all_lags=False,
    labeled_fluxes=None,
    labeled_rates=None,
    unlabeled_rates=None,
):
    '''Estimate fluxes and rates over multiple iterations. The number of iterations is determined by how many
    vectors of weights, parent IDs, bin assignments, and label assignments are passed.
    If ``all_lags`` is true, then the average is over all possible lags within the length-N window given, otherwise
    simply the length N lag.
    Returns labeled flux matrix, labeled rate matrix, and unlabeled rate matrix.'''
    assert len(weights) == len(parent_ids) == len(bin_assignments) == len(label_assignments)
    nstates = len(state_labels)
    nbins = labeled_pops.shape[1] - 1
    # Prepare output arrays
    if labeled_fluxes is None:
        labeled_fluxes = np.zeros((nstates, nstates, nbins, nbins), weight_dtype)
    else:
        labeled_fluxes.fill(0.0)
    if labeled_rates is None:
        labeled_rates = np.zeros_like(labeled_fluxes)
    else:
        labeled_rates.fill(0.0)
    if unlabeled_rates is None:
        unlabeled_rates = np.zeros((nbins, nbins), weight_dtype)
    else:
        unlabeled_rates.fill(0.0)
    # Loop over all possible windows to accumulate flux matrix
    # flux matrix is [initial_label][final_label][initial_bin][final_bin]
    if all_lags:
        twindow = calculate_labeled_fluxes_alllags(nstates, weights, parent_ids, bin_assignments, label_assignments, labeled_fluxes)
    else:
        twindow = calculate_labeled_fluxes(nstates, weights, parent_ids, bin_assignments, label_assignments, labeled_fluxes)
    labeled_fluxes /= twindow
    # Calculate rate matrix for this window, using populations from the last iteration (which correspond
    # to the weights that contribute to the flux matrix)
    labeled_flux_to_rate(labeled_fluxes, labeled_pops, labeled_rates)
    # Calculate an unlabeled rate matrix
    unlabeled_fluxes = np.sum(labeled_fluxes, axis=(0, 1))
    unlabeled_pops = labeled_pops.sum(axis=0)
    unlabeled_rates[...] = labeled_flux_to_rate(unlabeled_fluxes[None, None, :, :], unlabeled_pops[None, :])[0, 0]
    return labeled_fluxes, labeled_rates, unlabeled_rates