How to use the matscipy.neighbours.first_neighbours function in matscipy

def test_first_neighbours(self):
        i = [1,1,1,1,3,3,3]
        self.assertArrayAlmostEqual(first_neighbours(5, i), [-1,0,-1,4,-1,7])
        i = [0,1,2,3,4,5]
        self.assertArrayAlmostEqual(first_neighbours(6, i), [0,1,2,3,4,5,6])

def test_first_neighbours(self):
        i = [1,1,1,1,3,3,3]
        self.assertArrayAlmostEqual(first_neighbours(5, i), [-1,0,-1,4,-1,7])
        i = [0,1,2,3,4,5]
        self.assertArrayAlmostEqual(first_neighbours(6, i), [0,1,2,3,4,5,6])

)
            tmp_3_summed[:, x, y] = np.bincount(
                    bincount_bins, weights=weights, 
                    minlength=unique_pairs_i1_i2.shape[0]
            ) 
        nl_index_i1_j1 = None
        nl_index_i2_j1 = None
        weights = None
        tmp_3 = None
        bincount_bins = None
        if divide_by_masses:
            geom_mean_mass_i1_i2 = np.sqrt(
                np.take(masses_i, unique_pairs_i1_i2[:, 0]) * np.take(masses_i, unique_pairs_i1_i2[:, 1])
            )
            tmp_3_summed /= geom_mean_mass_i1_i2[:, np.newaxis, np.newaxis]
        index_ptr = first_neighbours(nat, unique_pairs_i1_i2[:, 0])
        term_8 = bsr_matrix((tmp_3_summed, unique_pairs_i1_i2[:, 1], index_ptr), shape=(3*nat, 3*nat))
        if symmetry_check:
            print("check term 8", np.linalg.norm(term_8.todense() - term_8.todense().T))
        return term_8

-------
    a : ase.Atoms
        Atomic configuration of the hydrogenated slab.
    """
    if b is None:
        b = a.copy()
        b.set_pbc(np.logical_not(mask))

    if exclude is None:
        exclude = np.zeros(len(a), dtype=bool)

    i_a, j_a, D_a, d_a = neighbour_list('ijDd', a, cutoff)
    i_b, j_b = neighbour_list('ij', b, cutoff)

    firstneigh_a = first_neighbours(len(a), i_a)
    firstneigh_b = first_neighbours(len(b), i_b)

    coord_a = np.bincount(i_a, minlength=len(a))
    coord_b = np.bincount(i_b, minlength=len(b))

    hydrogens = []
    # Surface atoms have coord_a != coord_b. Those need hydrogenation
    for k in np.arange(len(a))[np.logical_and(coord_a!=coord_b,
                                              np.logical_not(exclude))]:
        l1_a = firstneigh_a[k]
        l2_a = firstneigh_a[k+1]
        l1_b = firstneigh_b[k]
        l2_b = firstneigh_b[k+1]
        n_H = 0
        for l_a in range(l1_a, l2_a):
            assert i_a[l_a] == k
            bond_exists = False

# to memory consumption. If we would divide by masses afterwards,
        # we would have to create a sparse matrix with the same size as the
        # Hessian matrix, i.e. we would momentarily need twice the given memory.
        if divide_by_masses: 
            masses_i = atoms.get_masses().reshape(-1, 1, 1)
            geom_mean_mass_n = np.sqrt(
                np.take(masses_i, i_n) * np.take(masses_i, j_n)
            ).reshape(-1, 1, 1)
        else:
            masses_i = None
            geom_mean_mass_n = None

        #------------------------------------------------------------------------ 
        # Calculate pair contribution to the Hessian matrix
        #------------------------------------------------------------------------ 
        first_i = first_neighbours(nat, i_n)
        e_nc = (dr_nc.T / abs_dr_n).T # normalized distance vectors r_i^{\mu\nu}
        outer_e_ncc = e_nc.reshape(-1, 3, 1) * e_nc.reshape(-1, 1, 3)
        eye_minus_outer_e_ncc = np.eye(3, dtype=e_nc.dtype) - outer_e_ncc
        D = self._calculate_hessian_pair_term(
            nat, i_n, j_n, abs_dr_n, first_i, 
            drep_n, ddrep_n, outer_e_ncc, eye_minus_outer_e_ncc, 
            divide_by_masses, geom_mean_mass_n, masses_i
        )
        drep_n = None
        ddrep_n = None

        #------------------------------------------------------------------------ 
        # Calculate contribution of embedding term
        #------------------------------------------------------------------------ 
        # For each pair in the neighborlist, create arrays which store the 
        # derivatives of the embedding energy of the corresponding atoms.

try:
                from scipy.sparse import bsr_matrix
            except ImportError:
                raise ImportError(
                    "Import error: Can not output the hessian matrix since scipy.sparse could not be loaded!")

        f = self.f
        dict = self.dict
        df = self.df
        df2 = self.df2

        nat = len(atoms)
        atnums = atoms.numbers

        i_n, j_n, dr_nc, abs_dr_n = neighbour_list('ijDd', atoms, dict)
        first_i = first_neighbours(nat, i_n)

        e_n = np.zeros_like(abs_dr_n)
        de_n = np.zeros_like(abs_dr_n)
        dde_n = np.zeros_like(abs_dr_n)
        for params, pair in enumerate(dict):
            if pair[0] == pair[1]:
                mask1 = atnums[i_n] == pair[0]
                mask2 = atnums[j_n] == pair[0]
                mask = np.logical_and(mask1, mask2)

                e_n[mask] = f[pair](abs_dr_n[mask])
                de_n[mask] = df[pair](abs_dr_n[mask])
                dde_n[mask] = df2[pair](abs_dr_n[mask])

            if pair[0] != pair[1]:
                mask1 = np.logical_and(

Returns
    -------
    a : ase.Atoms
        Atomic configuration of the hydrogenated slab.
    """
    if b is None:
        b = a.copy()
        b.set_pbc(np.logical_not(mask))

    if exclude is None:
        exclude = np.zeros(len(a), dtype=bool)

    i_a, j_a, D_a, d_a = neighbour_list('ijDd', a, cutoff)
    i_b, j_b = neighbour_list('ij', b, cutoff)

    firstneigh_a = first_neighbours(len(a), i_a)
    firstneigh_b = first_neighbours(len(b), i_b)

    coord_a = np.bincount(i_a, minlength=len(a))
    coord_b = np.bincount(i_b, minlength=len(b))

    hydrogens = []
    # Surface atoms have coord_a != coord_b. Those need hydrogenation
    for k in np.arange(len(a))[np.logical_and(coord_a!=coord_b,
                                              np.logical_not(exclude))]:
        l1_a = firstneigh_a[k]
        l2_a = firstneigh_a[k+1]
        l1_b = firstneigh_b[k]
        l2_b = firstneigh_b[k+1]
        n_H = 0
        for l_a in range(l1_a, l2_a):
            assert i_a[l_a] == k