def get_spin_link_dict(peaklist):
    """
    Map each unique spin link to all of its corresponding peaks.

    NOESY peak lists represent spin links between Hydrogen atoms. Whether
    2D, 3D or 4D, each peak in a NOESY peak list has exactly two Hydrogen
    spins. Here, a spin link is represented by a frozenset containing the
    spin.assignment tuples for each Hydrogen atom. This function returns a
    dictionary mapping each unique spin link to a list of the Peaks in the
    PeakList that contain those two Hydrogen atoms.

    Examples
    --------
    >>> spin_link_dict = peaklist.spin_link_dict()
    >>> spin_link, peaks = spin_link_dict.popitem()
    >>> spin_link
    frozenset([Assignment(res_type='Q', res_num=21, atom='HN'), Assignment(
    res_type='G', res_num=17, atom='HN')])
    >>> print(peaks[0])
        Peak(spins=[
        Spin(res_type=G, res_num=17, atom=HN),
        Spin(res_type=G, res_num=17, atom=N),
        Spin(res_type=Q, res_num=21, atom=HN)])
    >>> print(peaks[1])
        Peak(spins=[
        Spin(res_type=Q, res_num=21, atom=HN),
        Spin(res_type=Q, res_num=21, atom=N),
        Spin(res_type=G, res_num=17, atom=HN)])

    Returns
    -------
    out : dictionary
        Each key is a frozenset made from two Hydrogen spin.assignment
        tuples. Each value is a list of peaks containing those two Hydrogen
        spins.
    """
    spin_link_dict = {}
    for peak in peaklist:
        spins = [spin for spin in peak
                 if spin.atom is not None and spin.atom[0] == 'H']
        if len(spins) != 2:
            err = ('expected 2 Hydrogens in each peak, '
                   'found %d' % len(spins))
            raise ValueError(err)
        link = frozenset(spin.assignment for spin in spins)
        spin_link_dict.setdefault(link, []).append(peak)
    return spin_link_dict