  An essential aspect of topological phases of matter is the existence of
Wilson loop operators that keep the ground state subspace invariant. Here we
present and implement an \it unbiased \rm numerical optimization scheme to
systematically find the Wilson loop operators given a single ground state wave
function of a gapped Hamiltonian on a disk. We then show how these Wilson loop
operators can be cut and glued through further optimization to give operators
that can create, move, and annihilate anyon excitations. We subsequently use
these operators to determine the braiding statistics and topological twists of
the anyons, yielding a way to fully extract topological order from a single
wave function. We apply our method to the ground state of the perturbed toric
code and doubled semion models with a magnetic field that is up to a half of
the critical value. From a contemporary perspective, this can be thought of as
a machine learning approach to discover emergent 1-form symmetries of a ground
state wave function. From an application perspective, our approach can be
relevant to find Wilson loop operators in current quantum simulators.
