def METHOD_NAME(bonds, n_atoms):
    """Returns an array mapping atoms to their bonded neighbors and an array
    that stores number of bonds made by each atom."""
    numbonds = np.bincount(bonds.reshape((bonds.shape[0] * 2)))
    bmap = np.zeros((n_atoms, numbonds.max()), int)
    bmap.fill(-1)
    index = np.zeros(n_atoms, int)
    for bond in bonds:
        a, b = bond
        bmap[a, index[a]] = b
        bmap[b, index[b]] = a
        index[bond] += 1
    return bmap, numbonds