def stratify_by_scores(scores, goal_n_strata='auto', method='cum_sqrt_F',
                       n_bins = 'auto'):
    available_methods = ['equal_size', 'cum_sqrt_F']
    if method not in available_methods:
        raise ValueError("method argument is invalid")
    if (method == 'cum_sqrt_F') or (goal_n_strata == 'auto'):
        if n_bins == 'auto':
            width_score = _heuristic_bin_width(scores)
            n_bins = np.ceil(sp.ptp(scores)/width_score).astype(int)
            print("Automatically setting n_bins = {}.".format(n_bins))
        counts, score_bins = np.histogram(scores, bins=n_bins)
        sqrt_counts = np.sqrt(counts)
        csf = np.cumsum(sqrt_counts)
        if goal_n_strata == 'auto':
            width_csf = _heuristic_bin_width(csf)
            goal_n_strata = np.ceil(sp.ptp(csf)/width_csf).astype(int)
            print("Automatically setting goal_n_strata = {}.".format(goal_n_strata))
        elif method == 'cum_sqrt_F':
            width_csf = csf[-1]/goal_n_strata
    if method == 'equal_size':
        sorted_ids = scores.argsort()
        n_items = len(sorted_ids)
        quotient = n_items // goal_n_strata
        remainder = n_items % goal_n_strata
        allocations = np.empty(n_items, dtype='int')
        st_pops = [quotient for i in range(goal_n_strata)]
        for i in range(remainder):
            st_pops[i] += 1
        j = 0
        for k,nk in enumerate(st_pops):
            start = j
            end = j + nk
            allocations[sorted_ids[start:end]] = k
            j = end
    if method == 'cum_sqrt_F':
        if goal_n_strata > n_bins:
            warnings.warn("goal_n_strata > n_bins. "
                          "Consider increasing n_bins.")
        csf_bins = [x * width_csf for x in np.arange(goal_n_strata + 1)]
        j = 0
        new_bins = []
        for (idx,value) in enumerate(csf):
            if j == (len(csf_bins) - 1) or idx == (len(csf) - 1):
                new_bins.append(score_bins[-1])
                break
            if value >= csf_bins[j]:
                new_bins.append(score_bins[idx])
                j += 1
        new_bins[0] -= 0.01
        new_bins[-1] += 0.01
        allocations = np.digitize(scores, bins=new_bins, right=True) - 1
        nonempty_ids = np.unique(allocations)
        n_strata = len(nonempty_ids)
        indices = np.arange(n_strata)
        allocations = np.digitize(allocations, nonempty_ids, right=True)
        if n_strata < goal_n_strata:
            warnings.warn("Failed to create {} strata".format(goal_n_strata))
    return Strata(allocations)