Sampling

Ensembles can be generated using a Generator. A generator is any object which has a class with a step() method that takes a configuration as a dictionary and returns a new configuration as a dictionary. Some generators can measure observables as the step is taken and return them inside the step, in which case the Configurations need a place to store them, which is created with the inline_observables() method.

class supervillain.generator.Generator[source]
step(configuration)[source]
Parameters

configuration (dictionary) – Contains (at least) the relevant field variables from which to update the Markov chain.

Returns

The next configuration (and any inline observables), also as a dictionary.

Return type

dictionary

inline_observables(steps)[source]
Parameters

steps (int)

Returns

Containing an initialized set of inline observables, each with an outer dimension of size steps. Presumably it makes sense to initialize to 0, but it’s not required.

Return type

dict

For example, a very dumb generator is

class DoNothing(Generator):

    def __init__(self):
        pass

    def step(self, x):
        r'''
        Just copy the incoming configuration and 'measure' the constant 1
        and add it to the configuration.
        '''
        return x.copy() | {'one': 1}

    def inline_observables(self, steps):
        r'''
        Obviously there is no *actual* reason to compute an observable whose value
        is identically 1, but it's useful to show how inline observables work.
        '''
        return {'one': Batch(steps, shape=(), dtype=float)}

Then, if you

L = supervillain.Lattice2D(5)
S = supervillain.Villain(L, kappa=0.1)
e = supervillain.Ensemble(S).generate(17, DoNothing(), start='cold')

the result will be an ensemble e that has 17 identical configurations (because DoNothing makes no updates). Obviously, for good sampling, DoNothing is the worst imaginable algorithm! It has an infinitely bad ergodicity problem!

The Villain Formulation

A less dumb algorithm is a local update, which changes only fields in some small area of the lattice.

As an example, we can formulate an update scheme offering localized changes to the \(\phi\) and \(n\) fields in the Villain formulation.

Picking a site \(x\) at random and proposing a change

\[\begin{split}\begin{aligned} \Delta\phi_x &\sim \text{uniform}(-\texttt{interval\_phi}, +\texttt{interval\_phi}) \\ \Delta n_\ell &\sim [-\texttt{interval\_n}, +\texttt{interval\_n}] \end{aligned}\end{split}\]

for the \(\phi\) on \(x\) and \(n\) on links \(\ell\) which touch \(x\) is ergodic (once swept over the lattice) and satisfies detailed balance so long as we accept the proposal based on the change of action. The NeighborhoodUpdateSlow generator implements this update algorithm but suffers from a variety of defects.

First, its implementation makes a lot of calls. We could ‘unfactor’ the proposal, single-site accept/reject, and sweep of the lattice for some benefits in speed. Moreover, to compute the change in action it evaluates the action for both configurations and takes the difference. But, we know that’s silly, most of the \(\phi\)s and \(n\)s haven’t changed at all and those links will contribute to both actions equally, giving zero difference. We could reduce the arithmetic substantially by computing the difference directly. Finally, for ease of thinking, each proposal reckons locations relative the origin and therefore moves all the fields around in order to update ergodically. All of that data movement adds cost, especially as the lattice gets large.

It was easy to write the NeighborhoodUpdateSlow but we can do better for production.

class supervillain.generator.villain.NeighborhoodUpdate(action, interval_phi=3.141592653589793, interval_n=1)[source]

This performs the same update as NeighborhoodUpdateSlow but is streamlined to eliminate calls, to calculate the change in action directly, and to avoid data movement.

Proposals are drawn according to

\[\begin{split}\begin{aligned} \Delta\phi_x &\sim \text{uniform}(-\texttt{interval\_phi}, +\texttt{interval\_phi}) \\ \Delta n_\ell &\sim W \times [-\texttt{interval\_n}, +\texttt{interval\_n}] \end{aligned}\end{split}\]

We pick \(\Delta n_\ell\) to be a multiple of the constraint integer \(W\) so that if the adjacent plaquettes satisfy the winding constraint \(dn \equiv 0 \text{ mod }W\) before the update they satisfy it after as well.

Each site update changes \(\phi\) at one site and \(n\) on the \(2D\) links touching that site (one forward and one backward link per direction). A checkerboard sweep is used so that the \(2D\) adjacent links of each same-color site are disjoint, enabling vectorized proposals and a single face_sum() call for \(\Delta S\).

See also

On a small 5×5 example this generator yields about three times as many updates per second than NeighborhoodUpdateSlow on my machine. This ratio should improve for larger lattices because the change in action is computed directly and is of fixed cost, rather than scaling with the volume.

step(cfg)[source]

Make a volume’s worth of random single-site updates.

Uses checkerboarding so that the 2D links adjacent to each same-color site are disjoint, allowing vectorized proposals and a single face_sum() to aggregate \(\Delta S\) per site.

Parameters

cfg (dict) – A dictionary with phi and n field variables.

Returns

Another configuration of fields.

Return type

dict

report()[source]

Returns a string with some summarizing statistics.

One reason the NeighborhoodUpdate suffers is that offering updates to 5 fields at onces (\(\phi\) and the 4 ns that touch it) means that the change in action can be large. We can decouple these proposals.

class supervillain.generator.villain.SiteUpdate(action, interval_phi=3.141592653589793)[source]

This performs the same update to \(\phi\) as NeighborhoodUpdate but leaves \(n\) untouched.

Proposals are drawn according to

\[\Delta \phi_x \sim \text{uniform}(-\texttt{interval\_phi}, +\texttt{interval\_phi})\]
step(cfg)[source]

Make volume’s worth of random single-site updates to \(\phi\).

Parameters

cfg (dict) – A dictionary with phi and n as Forms.

Returns

Updated configuration.

Return type

dict

inline_observables(steps)[source]
Parameters

steps (int)

Returns

Containing an initialized set of inline observables, each with an outer dimension of size steps. Presumably it makes sense to initialize to 0, but it’s not required.

Return type

dict

report()[source]
class supervillain.generator.villain.LinkUpdate(action, interval_n=1)[source]

This performs the same update to \(n\) as NeighborhoodUpdate but leaves \(\phi\) untouched.

Proposals are drawn according to

\[\begin{aligned} \Delta n_\ell &\sim W \times [-\texttt{interval\_n}, +\texttt{interval\_n}] \setminus \{0\} \end{aligned}\]

We pick \(\Delta n_\ell\) to be a multiple of the constraint integer \(W\) so that if the adjacent plaquettes satisfy the winding constraint \(dn \equiv 0 \text{ mod }W\) before the update they satisfy it after as well.

Note

You can run python supervillain/generator/villain/link.py to compare a pure \(W=1\) NeighborhoodUpdate ensemble against an ensemble which also does LinkUpdates. Note that adding the LinkUpdate costs essentially 0 time because all the links are done in parallel and there are no python-level for loops.

step(cfg)[source]

Make volume’s worth of random single-link updates to \(n\).

Parameters

cfg (dict) – A dictionary with phi and n as Forms.

Returns

Updated configuration.

Return type

dict

inline_observables(steps)[source]
Parameters

steps (int)

Returns

Containing an initialized set of inline observables, each with an outer dimension of size steps. Presumably it makes sense to initialize to 0, but it’s not required.

Return type

dict

report()[source]

When \(W=1\) the combination of the SiteUpdate and LinkUpdate are ergodic. But when \(W>1\) the LinkUpdate only offers changes to \(n\) by multiples of \(W\) to preserve the constraint \(dn = 0 \text{ mod }W\). For an ergodic algorithm when \(W>1\) we need to offer ways to change \(n\) by 1 (less than \(W\)) while maintaining the constraint. We need to make closed updates to \(n\), which can be broken up into ExactUpdates (which are automatically closed) and CohomologyUpdates (which change the global winding sector).

class supervillain.generator.villain.ExactUpdate(action, interval_z=1)[source]

The LinkUpdate only updates n by multiples of \(W\) in order to preserve the constraint \(dn = 0\; (\text{mod } W\)). Another way to preserve the constraint is to update n around a given site in a coordinated way so that \(dn\) is not changed on any of the neighboring plaquettes.

One way to accomplish this coordinated change is to start with a zero-form \(z\) and update \(n\) by \(dz\). Then, \(dn = d^2z = 0\).

Proposals are drawn according to

\[\begin{aligned} z_x &\sim [-\texttt{interval\_z}, +\texttt{interval\_z}] \setminus \{0\} \end{aligned}\]
step(cfg)[source]

Make a volume’s worth of locally-exact updates to n.

Parameters

cfg (dict) – A dictionary with phi and n as Forms.

Returns

Updated configuration.

Return type

dict

inline_observables(steps)[source]
Parameters

steps (int)

Returns

Containing an initialized set of inline observables, each with an outer dimension of size steps. Presumably it makes sense to initialize to 0, but it’s not required.

Return type

dict

report()[source]
class supervillain.generator.villain.CohomologyUpdate(action, interval_h=1)[source]

Local updates generate only exact changes \(\Delta n = dk\), which have zero cohomology class. No combination of them can change the global winding numbers

\[w_\mu = \sum_{x_\mu} n_\mu(x_\mu, x_\perp)\]

which are \(x_\perp\)-independent whenever \(dn = 0\). These winding numbers label sectors of \(H^1(\mathbb{T}^D, \mathbb{Z}) = \mathbb{Z}^D\).

This update proposes, for each direction \(\mu\), adding a constant \(h_\mu\) to \(n_\mu\) on the single slice \(x_\mu = 0\) (all perpendicular positions). Because $\Delta n_\mu$ is constant on that slice and zero elsewhere, \(d(\Delta n) = 0\) exactly, so the constraint \(dn \equiv 0\; (\text{mod}\; W)\) is preserved automatically for any \(W\). The winding number \(w_\mu\) changes by \(h_\mu\).

Proposals are drawn from

\[h_\mu \sim [-\texttt{interval\_h},\, +\texttt{interval\_h}] \setminus \{0\}\]

One independent scalar \(h_\mu\) is drawn and Metropolized per direction, changing \(N^{D-1}\) links at once. Acceptance is typically low at large \(\kappa\) because the action barrier scales as \(O(\kappa N^{D-1})\) — this is the genuine physical cost of tunneling between winding sectors.

step(configuration)[source]

Offer one winding-sector update per direction.

For each \(\mu\), draws a scalar \(h_\mu\), computes \(\Delta S\) on the \(N^{D-1}\) links of the slice \(x_\mu = 0\), and Metropolizes. The D directions are independent and processed sequentially; accepted changes are applied incrementally to the residual \(r = d\phi - 2\pi n\) so that later directions see the updated residuals.

Parameters

configuration (dict) – A dictionary with \(\phi\) and \(n\) as Forms.

Returns

Updated configuration.

Return type

dict

inline_observables(steps)[source]
Parameters

steps (int)

Returns

Containing an initialized set of inline observables, each with an outer dimension of size steps. Presumably it makes sense to initialize to 0, but it’s not required.

Return type

dict

report()[source]

The combination of the SiteUpdate, LinkUpdate, ExactUpdate, and CohomologyUpdate is ergodic even when \(W>1\). But it can be slow to decorrelate. As mentioned, the CohomologyUpdate often rejects because it touches a macroscopic number of variables (\(N^{D-1}\) links per direction). A major issue is that the route across the torus is very rigid: it’s just a straight shot. Smarter worm algorithms can make high-acceptance updates to many variables across the lattice, which can help overcome critical slowing down.

Remember that in the Villain case we are trying to sample according to

\[\begin{split}\begin{aligned} Z &= \sum\hspace{-1.33em}\int D\phi\; Dn\; Dv\; e^{-S[\phi, n, v]} \\ S[\phi, n, v] &= \frac{\kappa}{2} \sum_{\ell} (d\phi - 2\pi n)_\ell^2 + 2\pi i \sum_p v_p (dn)_p / W \end{aligned}\end{split}\]

and that we may directly path-integrate out the Lagrange multiplier \(v\) in favor of a constraint

\[Z = \sum\hspace{-1.33em}\int D\phi\; Dn\; e^{-S[\phi, n, v]} \prod_p [dn_p \equiv 0 \text{ mod }{W}]\]

One observable of interest is the Vortex_Vortex correlation function,

\[V_{x,y} = \left\langle e^{2\pi i (v_x - v_y) / W} \right\rangle\]

which poses a tricky problem to evaluate, since if we sample configurations of \(Z\) we integrate \(v\) out first and cannot easily compute the observable. Instead, we can think of adding an insertion of the vortex creation and annihilation operators, and we can absorb them into the action before path-integrating out \(v\) and use the insertions to shift the constraint at the insertions

(1)\[S[\phi, n, v] - 2\pi i (v_x - v_y) / W \rightarrow V_{x,y} = \frac{1}{Z} \sum\hspace{-1.33em}\int D\phi\; Dn\; e^{-S[\phi, n, v]} \prod_p [dn_p \equiv \delta_{px} - \delta_{py} \text{ mod }{W}]\]

where \(x\) and \(y\) label plaquettes. The built-in Vortex_Vortex observable measures this by picking a path between \(x\) and \(y\) on which to change \(n\) to satisfy the new shifted constraint. Depending on the parameters this inevitably hits an overlap problem, where we must change so many links between \(x\) and \(y\) that the change in action is very big and the correlator is naturally very small, except on rare configurations where it is big.

To address overlap problems of this kind Prokof’ev and Svistunov invented worm algorithms [2] which operate by introducing defects where constraints may be broken, propagating those defects, and when they meet and annihilate the configuration again obeys the constraint. Consider the mixed regular+path integral \(G\) with unspecified normalization \(N\) (that will cancel from all interesting quatities that follow) that integrates over all sectors of possible constraints

\[G = \frac{1}{N} \sum\hspace{-1.33em}\int D\phi\; Dn\; dh\; dt\; e^{-S[\phi, n, v]} \prod_p [dn_p \equiv \delta_{ph} - \delta_{pt} \text{ mod }{W}].\]

A configuration of \(G\) consists of two plaquettes \(h\) and \(t\), \(\phi\in\mathbb{R}\) on sites and \(n\in\mathbb{Z}\) on links which satisfy the \((h, t)\) constraint.

Worm algorithms sample from the larger space of configurations of \(G\) and, when the head \(h\) and tail \(t\) coincide, each \(G\) configuration is a valid \(Z\) configuration appearing with the right relative frequency. This is very different from thinking that we need to carefully design updates to always maintain the constraint. No! Violate the constraint as part of the evolution!

If we want to use a worm as an update in our Markov chain, we need to emit configurations that satisfy the original winding constraint. We can imagine constructing a Markov chain to sample configurations of both \(Z\) and \(G\). To go from a \(Z\) configuration insert the head and tail on the same randomly-chosen location (in this case, a plaquette); the change in action is 0 and this is automatically accepted. This configuration now ‘has a worm’ even though the fields have not changed.

Next we evolve the configuration in the larger set of all \(G\) configurations which need not satisfy the original winding constraint but instead violates it in exactly two places (the head and the tail), in opposite ways. Starting from a ‘diagonal’ configuration which has \(h=t\), we can move the head by one plaquette. As the head moves it changes the \(n\) on the link it crosses so that it changes \(dn\) on both adjacent plaquettes. With a clever choice we can ensure that as the head leaves a plaquette it restores the constraint there and breaks it on the destination plaquette. In this way the head moves around the lattice. This change of \(n\) changes the action \(S\), and so the moves need to be Metropolis-tested in some sense.

Finally, when the head reaches the tail, the constraint is restored everywhere. We might allow the worm to keep evolving, but since it satisfies the winding constraint it is possible to reach a \(Z\) configuration too. To go from a \(G\) configuration to a \(Z\) configuration we require the head and tail to coincide and satisfy the winding constraint; the change in action is 0 and the acceptance is likewise automatically accepted if this change is proposed. Now we have a \(Z\) configuration and add it to our Markov chain (or update it with different generators first). Inside the worm algorithm we need not provide a direct route from any \(Z\) configuration to any other; the only way between \(Z\) configurations is through \(G\), though nothing prevents us from using other generators to go directly from \(Z\) configuration to \(Z\) configuration.

Notice that

(2)\[V_{x,y} = \frac{\left\langle \delta_{xh} \delta_{yt} \right\rangle_G}{\left\langle \delta_{ht} \right\rangle_G}\]

where the expectation values are with respect to configurations drawn from \(G\) (not \(Z\)!). Let us understand this expression for the vortex correlator. It says that if we draw from the larger space of \(G\) configurations and make a histogram in \(x\) and \(y\), we can normalize that histogram by its value at zero displacement to get the \(V_{x,y}\). We can measure the histogram as we go and report with the \(\phi\) and \(n\). The reason to save the histogram as a sample rather than to just accumulate a histogram for the whole worm’s evolution is that once we reach a \(Z\) configuration we update the other variables; that histogram is conditional on those variables; when they change the histogram will too. So, the worm’s displacement histogram can be saved inline as Vortex_Vortex, as long as we remember to normalize any DerivedQuantity that depends on it by the element of the expectation at the origin.

class supervillain.generator.villain.worm.ClassicWorm(S)[source]

This implements the classic worm of Prokof’ev and Svistunov [2] for the Villain links \(n\in\mathbb{Z}\) which satisfy \(dn \equiv 0 \) (mod W) on every plaquette.

On top of a constraint-satisfying configuration we put down a worm and let the head move, changing the crossed links. We uniformly propose a move in all 4 directions and Metropolize the change.

Additionally, when the head and tail coincide, we allow a fifth possible move, where we remove the worm and emit the updated \(z\) configuration into the Markov chain.

As we evolve the worm we tally the histogram that yields the Vortex_Vortex correlation function.

Warning

This update algorithm is not ergodic on its own. It doesn’t change \(\phi\) at all and even leaves \(dn\) alone (while changing \(n\) itself). It can be used, for example, Sequentially with the SiteUpdate and LinkUpdate for an ergodic method.

Note

Because it doesn’t change \(dn\) at all, this algorithm can play an important role in sampling the \(W=\infty\) sector, where all vortices are completely killed, though updates to \(\phi\) would still be needed.

Note

This class contains kernels accelerated using numba.

Todo

Generalize to D>2. The current implementation is hardcoded for D=2: it uses the _Lattice2D numba struct, four fixed directions (east/north/west/south), and 2D displacement histograms. Raises NotImplementedError for \(D \neq 2\).

inline_observables(steps)[source]

The worm algorithm can measure the Vortex_Vortex correlator. We also store the Worm_Length for each step.

step(configuration)[source]

Given a constraint-satisfying configuration, returns another constraint-satisfying configuration udpated via worm as described above.

report()[source]

The worm is not ergodic on its own—it doesn’t update \(\phi\), for example, and it cannot change a link by ±W. But, in combination of with SiteUpdate and LinkUpdate it is ergodic; the worm can replace the combination of ExactUpdate and CohomologyUpdate. The ExactUpdate can be understood as a very simple worm that takes the tightest nontrivial path, while the CohomologyUpdate can be understood as a worm that goes once around the world in a straight shot. The worm offers dynamically determined constraint-preserving updates and is much more flexible. In can change the holonomy, for example, by finding a route around the torus that isn’t a straight shot but runs through the valley of the action.

Finally, we provide a convenience function which provides an ergodic generator.

supervillain.generator.villain.Hammer(S, worms=1)[source]

The Hammer is just syntactic sugar for a Sequentially applied ergodic combination of generators. It may change from version to version as new generators become available or get improved.

Note

When \(W=\infty\) we only include updates that leave \(dn=0\) (NOT \(\text{mod }W\)!).

Note

The ClassicWorm is currently only implemented for \(D=2\). In higher dimensions the Hammer omits it, so the returned combination is not ergodic on its own until a \(D>2\) worm is available.

Parameters
  • S (a Villain action)

  • worms (int) – A positive integer saying how many worms to do per iteration.

Return type

An ergodic generator for updating Villain configurations.

The Worldline Formulation

In the Worldline formulation the constraint \(\delta m = 0\) restricts which kinds of updates we could propose. For example, changing only a single link is guaranteed to break the constraint on both ends. So, we need clever generators to maintain the constraint.

class supervillain.generator.worldline.PlaquetteUpdate(action)[source]

Ref. [3] suggests a simple update scheme where the links surrounding a single plaquette are updated in concert so that the Worldline constraint is maintained.

The links on a plaquette are updated in the same (oriented) way, which guarantees that \(\delta m\) is invariant before and after the update, so that if it started 0 everywhere so it remains. The coordinated change is randomly chosen from ±1.

Warning

HOWEVER this algorithm is not ergodic on its own. The issue is that no proposal can change the worldline TorusWrapping. Instead, if you start cold with \(m=0\), which has global wrapping of (0,0) you stay in the (0,0) sector.

step(cfg)[source]

Performs a sweep over all (site, direction-pair) plaquettes in a randomized order.

For each plaquette at site x with directions (μ, ν), the 4 boundary links of the plaquette are updated in concert (±Δm with orientation) and v[comp, x] is also updated. This preserves δm=0 because the boundary change is exact: δ(∂A)=0. Direction pairs are processed sequentially to avoid shared-link conflicts.

report()[source]

We can decouple the proposals in the PlaquetteUpdate, and update the vortex fields and the worldlines independently.

class supervillain.generator.worldline.VortexUpdate(action, interval_v=1)[source]

This performs the same update to \(v\) as PlaquetteUpdate but leaves \(m\) untouched.

Proposals are drawn according to

\[\begin{split}\begin{aligned} \Delta v_p &\sim [-\texttt{interval\_v}, +\texttt{interval\_v}] \setminus \{0\} &&(W<\infty) \\ \Delta v_p &\sim \text{uniform}(-\texttt{interval\_v}, +\texttt{interval\_v}) &&(W=\infty) \end{aligned}\end{split}\]

on each plaquette \(p\) independently.

Warning

This update is not ergodic on its own, since it does not change \(m\) at all.

step(cfg)[source]

Make a volume’s worth of changes to v.

Parameters

cfg (dict) – A dictionary with m and v field variables.

Returns

Another configuration of fields.

Return type

dict

report()[source]
class supervillain.generator.worldline.CoexactUpdate(action, interval_t=1)[source]

One way to guarantee that \(\delta m = 0\) is to change \(m\) by \(\delta t\) where \(t\) is an integer-valued two-form.

Proposals are drawn according to

\[\begin{aligned} t_p &\sim [-\texttt{interval\_t}, +\texttt{interval\_t}] \setminus \{0\} & \Delta m_\ell &= (\delta t)_\ell \end{aligned}\]

Warning

This algorithm is not ergodic on its own. It does not change \(v\) (see the VortexUpdate) nor can it produce all coclosed changes—only coexact changes. For a coïnexact coclosed update consider the WrappingUpdate or the worm.

step(cfg)[source]

Make a volume’s worth of locally-exact updates to m.

Parameters

cfg (dict) – A dictionary with m and v field variables.

Returns

Another configuration of fields.

Return type

dict

report()[source]

To have a fully ergodic algorithm we will also need to update the TorusWrapping of the worldlines.

class supervillain.generator.worldline.WrappingUpdate(action, interval_w=1)[source]

Because CoexactUpdate fails to change the wrapping, we should separately offer wrapping-changing proposals.

We propose coordinated changes on all the x-direction links on a single timeslice and coordinated changes on all the t-direction links on a single spatial slice.

That is, on all the \(m\) on a cycle around the torus we propose to change \(m\) according to

\[\Delta m \sim [-\texttt{interval\_w}, +\texttt{interval\_w}] \setminus \{0\}\]

Warning

This algorithm is not ergodic on its own. The issue is that no proposal can generate coexact changes.

step(cfg)[source]

Propose independent updates of \(m\) on every timeslice and on every spatial slice.

In principle all the proposals may be made in parallel but we just do them sequentially.

report()[source]

The combination of the VortexUpdate, CoexactUpdate, and WrappingUpdate are ergodic. However, the may be suboptimal. In particular the WrappingUpdate touches a large number of links which often leads to very large changes in action and rejection. Just as in the Villain case we can make smarter updates to a dynamically-determined set of variables with high acceptance by using a worm algorithm.

Unlike the Villain formulation, the Worldline formulation has a constraint even when \(W=1\), \(\delta m = 0\) everywhere, from path-integrating \(\phi\)

\[\begin{split}\begin{aligned} Z &= (2\pi\kappa)^{-|\ell|/2}\sum\hspace{-1.33em}\int D\phi\; Dm\; Dv\; e^{-S[\phi, m, v]} \\ S[\phi, m, v] &= \frac{1}{2\kappa} \sum_\ell \left(m - \frac{\delta v}{W} \right)_\ell^2 - i \sum_x \left(\delta m\right)_x \phi_x \end{aligned}\end{split}\]

Thinking about the Spin_Spin correlation function, we want to insert \(e^{i(\phi_x-\phi_y)}\), which shifts the constraint to

\[(\delta m)_s = \delta_{sx} - \delta_{sy}\]

and we can perform updates in the mixed regular+path integral \(G\) that integrates over all sectors of possible constraints

\[G = \frac{1}{N} (2\pi\kappa)^{-|\ell|/2} \sum\hspace{-1.33em}\int Dm\; Dv\; dh\; dt\; e^{-S[m, v]} \prod_s [\delta m = \delta_{sh} - \delta_{st}]\]

As in the Villain case, when the head \(h\) and tail \(t\) coincide a configuration of \(G\) is a valid \(Z\) configuration and each appears with the correct relative frequency. Again, the philosophy is to evolve in a larger space with the constraint lifted and celebrate when we receive a constraint-satisfying configuration.

Just like the Villain case we can measure a two-point correlator inline, but in this case the constraint is \(\delta m = 0\) everywhere and constraint-violating insertions are of \(e^{\pm i \phi}\) (as in Spin_Spin),

\[S_{x,y} = \frac{\left\langle \delta_{xh} \delta_{yt} \right\rangle}{\left\langle \delta_{ht} \right\rangle}\]

which amounts to constructing a normalized histogram after ensemble averaging.

class supervillain.generator.worldline.worm.ClassicWorm(S)[source]

This implements the classic worm of Prokof’ev and Svistunov [2] for the worldline links \(m\in\mathbb{Z}\) which satisfy \(\delta m = 0\) on every site.

On top of a constraint-satisfying configuration we put down a worm and let the head move, changing the crossed links. We uniformly propose a move in all \(2D\) directions and Metropolize the change.

Additionally, when the head and tail coincide, we allow a \((2D+1)\)-th possible move, where we remove the worm and emit the updated configuration into the Markov chain.

As we evolve the worm we tally the histogram that yields the Spin_Spin correlation function.

Warning

When \(W>1\) this update algorithm is not ergodic on its own. It doesn’t change \(v\) at all. However, when \(W=1\) we can always pick \(v=0\) (any other choice may be absorbed into \(m\)), and this generator can stand alone.

Note

This class contains kernels accelerated using numba.

inline_observables(steps)[source]

The worm algorithm can measure the Spin_Spin correlator. We also store the Worm_Length for each step.

step(configuration)[source]

Given a constraint-satisfying configuration, returns another constraint-satisfying configuration updated via worm as described above.

report()[source]

Finally, we provide a convenience function which provides an ergodic generator.

supervillain.generator.worldline.Hammer(S, worms=1)[source]

The Hammer is just syntactic sugar for a Sequentially applied ergodic combination of generators. It may change from version to version as new generators become available or get improved.

Parameters
  • S (a Worldline action)

  • worms (int) – A positive integer saying how many worms to do per iteration.

Return type

An ergodic generator for updating Villain configurations.

Combining Generators

You can combine generators. One simple combination is just the sequential application.

class supervillain.generator.combining.Sequentially(generators)[source]

Sequentially applies the list of generators as a single step.

For example we can get an ergodic Worldline update scheme by combining the PlaquetteUpdate and WrappingUpdate.

>>> p = supervillain.generator.worldline.PlaquetteUpdate(S)
>>> h = supervillain.generator.worldline.WrappingUpdate(S)
>>> g = supervillain.generator.combining.Sequentially((p, h))
Parameters

generators (iterable of generators) – The ordered iterable of generators to apply one after the next.

step(cfg)[source]

Apply each generator’s step one after the next, passing each generator’s full output as the input to the next, and return the final configuration.

inline_observables(steps)[source]
Parameters

steps (int)

Returns

Containing an initialized set of inline observables, each with an outer dimension of size steps. Presumably it makes sense to initialize to 0, but it’s not required.

Return type

dict

report()[source]

Returns a string with some summarizing statistics.

Another (trivial) combination is a decorrelating application.

class supervillain.generator.combining.KeepEvery(n, generator, blocked_inline=True)[source]

To decorrelate and get honest error estimates it can be helpful to do MCMC but then only analyze evenly-spaced configurations. Rather than keep every single generated configuration and then throw a bunch away, we can keep only those we might analyze.

Note

The number of updates per second will decrease by a factor of n, but the autocorrelation time should be n less. Generating a fixed number of configurations will take n times longer.

>>> p = supervillain.generator.worldline.PlaquetteUpdate(S)
>>> g = supervillain.generator.combining.KeepEvery(10, p)
Parameters
  • n (int) – How many generator updates to make before emitting a configuration.

  • generator – The generator to use for updates.

  • blocked_inline (bool) – Rather than just keeping the inline_observables() from the last update, all of the inline measurements are averaged across all n updates. This helps capture rare-but-important measurements that would otherwise be missed.

step(cfg)[source]

Applies the generator n times and returns the resulting configuration.

inline_observables(steps)[source]
Parameters

steps (int)

Returns

Containing an initialized set of inline observables, each with an outer dimension of size steps. Presumably it makes sense to initialize to 0, but it’s not required.

Return type

dict

report()[source]

Returns a string with some summarizing statistics.

Monitors

We can make pass-through generators to perform specific tasks, like logging.

class supervillain.generator.monitor.Logger(generator, channel)[source]

Mostly good for debugging.

import logging
logger = logging.getLogger(__name__)

...

L = supervillain.lattice.Lattice2D(N)
S = supervillain.action.Worldline(L, kappa, W)
G = supervillain.generator.worldline.Hammer(S, L.sites)
M = supervillain.generator.monitor.Logger(G, logger.info)
E = supervillain.Ensemble(S).generate(100, M)
Parameters
  • generator (any generator) – The generator to wrap.

  • channel (some function) – Will be applied to the resulting configuration.

step(cfg)[source]

Apply the wrapped generator and channel the result.

inline_observables(steps)[source]

Forwards to the wrapped generator.

report()[source]

Forwards to the wrapped generator.

One could construct monitors that, for example, measured observables inline.