Untitled

def initialize_spinful_random(fc_up, fc_dn, fcdag_up, fcdag_dn, net, n_up, n_down):
    """"
    Initializes a spinful lattice randomly.

    Parameters
    ----------
    fc_up : Annihilation operator for spin up fermions.
    fc_dn : Annihilation operator for spin down fermions.
    fcdag_up : Creation operator for spin up fermions.
    fcdag_dn : Creation operator for spin down fermions.
    net : class Lattice
    n_up, n_down: number of up and down-spin lectrons distributed

    Returns
    -------
    gamma : Peps object
        The post-sampling state tensor network.
    """

    Nx, Ny = net.Nx, net.Ny
    out = {(x, y): 3 for x in range(Nx) for y in range(Ny)}
    available_sites = [(x, y) for x in range(Nx) for y in range(Ny)]
    for _ in range(n_up):
        site = random.choice(available_sites)
        if out[site] == 3:
            out[site] = 0  # Only up-spin
        available_sites.remove(site)
    available_sites = [(x, y) for x in range(Nx) for y in range(Ny)]
    for _ in range(n_down):
        site = random.choice(available_sites)
        if out[site] == 3:
            out[site] = 1  # Only down-spin
        elif out[site] == 0:
            out[site] = 2  # Both up-spin and down-spin
        available_sites.remove(site)
        
    n_up, n_dn, h_up, h_dn = fcdag_up @ fc_up, fcdag_dn @ fc_dn, fc_up @ fcdag_up, fc_dn @ fcdag_dn
    nn_up, nn_dn, nn_do, nn_hole = n_up @ h_dn, h_up @ n_dn, n_up @ n_dn, h_up @ h_dn   # up - 0; down - 1; double occupancy - 2; hole - 3
    nn_up = reduce_operators(nn_up)
    nn_dn = reduce_operators(nn_dn)
    nn_do = reduce_operators(nn_do)
    nn_hole = reduce_operators(nn_hole)

    tt = {0: nn_up, 1: nn_dn, 2: nn_do, 3: nn_hole}
   
    gamma = fpeps.Peps(net.lattice, net.dims, net.boundary)
    for kk in gamma.sites():
        Ga = tt[out[kk]].fuse_legs(axes=[(0, 1)])
        for s in (-1, 1, 1, -1):
            Ga = Ga.add_leg(axis=0, s=s)
        gamma[kk] = Ga.fuse_legs(axes=((0, 1), (2, 3), 4))
        
    return gamma