How to use the netket.operator.GraphOperator function in NetKet

[9, 10],
    [10, 11],
    [11, 12],
    [12, 13],
    [13, 14],
    [14, 15],
    [15, 16],
    [16, 17],
    [17, 18],
    [18, 19],
    [19, 0],
]

g = nk.graph.CustomGraph(edges=edges)
hi = nk.hilbert.CustomHilbert(local_states=[-1, 1], graph=g)
ha = nk.operator.GraphOperator(
    hi, siteops=[sigmax], bondops=[mszsz], bondops_colors=[0]
)
operators["Graph Hamiltonian"] = ha

# Custom Hamiltonian
sx = [[0, 1], [1, 0]]
sy = [[0, 1.0j], [-1.0j, 0]]
sz = [[1, 0], [0, -1]]
g = nk.graph.CustomGraph(edges=[[i, i + 1] for i in range(20)])
hi = nk.hilbert.CustomHilbert(local_states=[1, -1], graph=g)

sx_hat = nk.operator.LocalOperator(hi, [sx] * 3, [[0], [1], [5]])
sy_hat = nk.operator.LocalOperator(hi, [sy] * 4, [[2], [3], [4], [9]])
szsz_hat = nk.operator.LocalOperator(hi, sz, [0]) * nk.operator.LocalOperator(
    hi, sz, [1]
)

else:
        hi = netket.hilbert.Spin(s=0.5,  graph=g)
    sigmaz = [[1, 0], [0, -1]]
    sigmax = [[0, 1], [1, 0]]
    sigmay = [[0, -1j], [1j, 0]]

    interaction = numpy.kron(sigmaz, sigmaz) + numpy.kron(sigmax, sigmax) + numpy.kron(sigmay, sigmay)

    bond_operator = [
        (J[0] * interaction).tolist(),
        (J[1] * interaction).tolist(),
    ]

    bond_color = [1, 2]

    return netket.operator.GraphOperator(hi, bondops=bond_operator, bondops_colors=bond_color)

mszsz = (np.kron(sigmaz, sigmaz)).tolist()

# Notice that the Transverse-Field Ising model as defined here has sign problem
L = 20
site_operator = [sigmax]
bond_operator = [mszsz]

# Hypercube
g = nk.graph.Hypercube(length=L, n_dim=1, pbc=True)

# Custom Hilbert Space
hi = nk.hilbert.Spin(s=0.5, graph=g)

# Graph Operator
op = nk.operator.GraphOperator(hi, siteops=site_operator, bondops=bond_operator)

# Restricted Boltzmann Machine
ma = nk.machine.RbmSpin(hilbert=hi, alpha=1)
ma.init_random_parameters(seed=1234, sigma=0.01)

# Local Metropolis Sampling
sa = nk.sampler.MetropolisLocal(machine=ma)

# Optimizer
opt = nk.optimizer.AdaMax()

# Stochastic reconfiguration
gs = nk.variational.Vmc(
    hamiltonian=op,
    sampler=sa,
    optimizer=opt,

g = netket.graph.CustomGraph(edge_colors)
    hi = netket.hilbert.Spin(s=0.5, graph=g)

    sigmaz = [[1, 0], [0, -1]]
    sigmax = [[0, 1], [1, 0]]
    sigmay = [[0, -1j], [1j, 0]]

    interaction = numpy.kron(sigmaz, sigmaz) + numpy.kron(sigmax, sigmax) + numpy.kron(sigmay, sigmay)

    bond_operator = [
        (interaction).tolist(),
    ]

    bond_color = [1]

    return netket.operator.GraphOperator(hi, bondops=bond_operator, bondops_colors=bond_color)

# Define custom graph
G = nx.Graph()
for i in range(L):
    G.add_edge(i, (i + 1) % L, color=1)
    G.add_edge(i, (i + 2) % L, color=2)

edge_colors = [[u, v, G[u][v]["color"]] for u, v in G.edges]

# Custom Graph
g = nk.graph.CustomGraph(edge_colors)

# Spin based Hilbert Space
hi = nk.hilbert.Spin(s=0.5, total_sz=0.0, graph=g)

# Custom Hamiltonian operator
op = nk.operator.GraphOperator(hi, bondops=bond_operator, bondops_colors=bond_color)

# Restricted Boltzmann Machine
ma = nk.machine.RbmSpin(hi, alpha=1)
ma.init_random_parameters(seed=1234, sigma=0.01)

# Sampler
sa = nk.sampler.MetropolisHamiltonianPt(machine=ma, hamiltonian=op, n_replicas=16)

# Optimizer
opt = nk.optimizer.Sgd(learning_rate=0.01)

# Variational Monte Carlo
gs = nk.variational.Vmc(
    hamiltonian=op,
    sampler=sa,
    optimizer=opt,

Sup = [[0, np.sqrt(2), 0], [0, 0, np.sqrt(2)], [0, 0, 0]]
Sdn = [[0, 0, 0], [np.sqrt(2), 0, 0], [0, np.sqrt(2), 0]]

# Heisenberg term
heisenberg = 0.5 * (np.kron(Sup, Sdn) + np.kron(Sdn, Sup)) + np.kron(Sz, Sz)
# AKLT two-site projector
P2_AKLT = 0.5 * heisenberg + np.dot(heisenberg, heisenberg) / 6.0 + np.identity(9) / 3.0

# 1D Lattice
g = nk.graph.Hypercube(length=10, n_dim=1, pbc=True)

# Hilbert space of spin-1s on the graph
hi = nk.hilbert.Spin(s=1, graph=g)

# AKLT model Hamiltonian as graph
ha = nk.operator.GraphOperator(hilbert=hi, bondops=[P2_AKLT.tolist()])

# Perform Lanczos Exact Diagonalization to get lowest three eigenvalues
res = nk.exact.lanczos_ed(ha, first_n=3, compute_eigenvectors=True)

# Print eigenvalues
print("eigenvalues:", res.eigenvalues)

# Compute energy of ground state
print("ground state energy:", res.mean(ha, 0))

# Compute energy of first excited state
print("first excited energy:", res.mean(ha, 1))