Compute the mean of correlation between bzt and vasp bandstructure on
    sym line, for all bands and locally (for each branches) the difference
    squared (%) if nb is specified.

def compare_sym_bands(bands_obj, bands_ref_obj, nb=None):
    """
    Compute the mean of correlation between bzt and vasp bandstructure on
    sym line, for all bands and locally (for each branches) the difference
    squared (%) if nb is specified.
    """

    nkpt = len(bands_obj.kpoints)
    if bands_ref_obj.is_spin_polarized:
        nbands = min(bands_obj.nb_bands, 2 * bands_ref_obj.nb_bands)
    else:
        # TODO: why is this needed? Shouldn't pmg take care of nb_bands?
        nbands = min(len(bands_obj.bands[Spin.up]),
                     len(bands_ref_obj.bands[Spin.up]))
    # print(nbands)
    arr_bands = np.array(bands_obj.bands[Spin.up][:nbands])
    # arr_bands_lavg = (arr_bands-np.mean(arr_bands,axis=1).reshape(nbands,1))

    if bands_ref_obj.is_spin_polarized:
        arr_bands_ref_up = np.array(bands_ref_obj.bands[Spin.up])
        arr_bands_ref_dw = np.array(bands_ref_obj.bands[Spin.down])
        # print(arr_bands_ref_up.shape)
        arr_bands_ref = np.vstack((arr_bands_ref_up, arr_bands_ref_dw))
        arr_bands_ref = np.sort(arr_bands_ref, axis=0)[:nbands]
        # print(arr_bands_ref.shape)
    else:
        arr_bands_ref = np.array(bands_ref_obj.bands[Spin.up][:nbands])

    # arr_bands_ref_lavg =
    # (arr_bands_ref-np.mean(arr_bands_ref,axis=1).reshape(nbands,1))

    # err = np.sum((arr_bands_lavg-arr_bands_ref_lavg)**2,axis=1)/nkpt
    corr = np.array(
        [distance.correlation(arr_bands[idx], arr_bands_ref[idx]) for idx in
         range(nbands)])

    if type(nb) == int: nb = [nb]

    bcheck = {}

    if max(nb) < nbands:
        branches = [[s['start_index'], s['end_index'], s['name']] for s in
                    bands_ref_obj.branches]

        if not bands_obj.is_metal() and not bands_ref_obj.is_metal():
            zero_ref = bands_ref_obj.get_vbm()['energy']
            zero = bands_obj.get_vbm()['energy']
            if not zero:
                vbm = bands_ref_obj.get_vbm()['band_index'][Spin.up][-1]
                zero = max(arr_bands[vbm])
        else:
            zero_ref = 0  # bands_ref_obj.efermi
            zero = 0  # bands_obj.efermi
            print(zero, zero_ref)

        for nbi in nb:
            bcheck[nbi] = {}

            bcheck[nbi]['Dist'] = np.mean(abs(arr_bands[nbi] - zero
                                              - arr_bands_ref[nbi] + zero_ref))
            bcheck[nbi]['Corr'] = corr[nbi]

            for start, end, name in branches:
                # werr.append((sum((arr_bands_corr[nb][start:end+1] -
                # arr_bands_ref_corr[nb][start:end+1])**2)/(end+1-start)*100,name))
                bcheck[nbi][name] = np.mean(abs(arr_bands[nbi][start:end + 1]
                                                - zero
                                                - arr_bands_ref[nbi][
                                                  start:end + 1] + zero_ref))
    else:
        bcheck = "No nb given"

    return bcheck