Reference Implementations

Sometimes the obvious implementation can be accelerated dramatically.

It is useful to have a reference implementation where things are done the totally obvious way so that smarter or faster approaches have a fixed target for correctness.

Note

Everything here can be considered an implementation detail and a user should not need the reference implementations at all.

Generators

Villain NeighborhoodUpdate

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

A neighborhood update changes only fields in some small area of the lattice.

In particular, this updating scheme changes the \(\phi\) and \(n\) fields in the Villain formulation.

It works by picking a site \(x\) at random, 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\).

Warning

Because we currently restrict to \(W=1\) for the Villain formulation we do not update \(v\).

Parameters

action (Villain) – The action from which we sample.
interval_phi (float) – A single float used to construct the uniform distribution for \(\phi\).
interval_n (int) – A single integer that gives the biggest allowed changes to \(n\).

proposal(cfg, dx)[source]

Parameters

cfg (dict) – A dictionary with \(\phi\) and \(n\) to update.
dx (Lattice coordinates) – Which site to move to the origin and update.

Returns

A new configuration with updated \(\phi\) and \(n\).

Return type

dict

site(cfg, dx)[source]

Rather than accepting every proposal() we perform importance sampling by doing a Metropolis accept/reject step [5] on every single-site proposal.

Parameters

cfg (dict) – A dictionary with \(\phi\) and \(n\) to update.
dx (Lattice coordinates) – Which site to move to the origin and update.

Returns

dict – A configuration; either the provided one a new one changed by a proposal.
float – The Metropolis-Hastings acceptance probability.
int – 1 if the proposal was accepted, 0 otherwise.

step(cfg)[source]

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

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.

Villain Classic Worm

class supervillain.generator.reference_implementation.villain.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.

Warning

Because the algorithm is about moving a single defect around the lattice, when implemented in pure python the python-level loop can severely impact performance. While this reference implementation was done in pure python, the production-ready generator uses numba for acceleration.

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.

Todo

Generalize to D>2. The current implementation hardcodes D=2: four directions (east/north/west/south) in _neighboring_plaquettes() and _surrounding_links(), 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]

Wordline Classic Worm

class supervillain.generator.reference_implementation.worldline.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.

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]

Observables

Spin Correlations

class supervillain.observable.reference_implementation.spin.Spin_SpinSloppy[source]

Bases: Observable

This performs the same measurement as the non-Sloppy version but does not get all the juice out of every Worldline configuration.

See the Spin_Spin documentation for detailed descriptions.

static Villain(S, phi)[source]: The same as in the Spin_Spin.

static Worldline(S, Links)[source]: See the Spin_Spin documentation.

Note

This taxicab-path measurement is only implemented for \(D=2\) and raises NotImplementedError otherwise.

class supervillain.observable.reference_implementation.spin.Spin_SpinSlow[source]

Bases: Observable

We can deform \(Z_J \rightarrow Z_{J}[x,y]\) to include the creation of a boson at \(y\) and the destruction of a boson at \(x\) in the action. We define the expectation value

\[S_{x,y} = \frac{1}{Z_J} Z_J[x,y]\]

and reduce to a single relative coordinate

\[\texttt{Spin\_Spin}_{\Delta x} = S_{\Delta x} = \frac{1}{\Lambda} \sum_x S_{x,x-\Delta x}\]

Differential Forms

Interlaced Forms

The interlaced picture can be taken literally: a \(p\)-form on an \(N^D\) lattice is stored in the elements of a \((2N)^D\) array with exactly \(p\) odd coordinates. However, this can be very wasteful since only \(\binom{D}{p}/2^D\) of the array elements are nonzero. Nevertheless, it is useful to have an implementation that realizes this picture directly so that the production Form, which is space-efficient, can be tested against it.

The operators of the interlaced format are independent implementations of the exterior calculus with the same sign conventions, giving a fixed target for correctness.

Below, let the unit step \(\varepsilon_k\) advance \(\xi_k\) by one (half a physical lattice step in direction \(k\)).

class supervillain.lattice.interlaced.Lattice(D, N)[source]

A D-dimensional hypercubic lattice of size N, in the interlaced layout.

All forms are (2N,)*D arrays. lat.random(p) returns a p-form with random values at p-form sites and zero everywhere else.

form(p)[source]: Return the p-form filter: 1 at p-form sites, 0 elsewhere.

zeros(p, dtype=<class 'float'>)[source]: Zero p-form array, shape (2N,)*D.

random(p)[source]: Random p-form: uniform [0,1) at p-form sites, 0 elsewhere.

Translation

Forms can be translated around the (periodic) interlaced lattice. Only a completely even shift maps a \(p\)-form to another \(p\)-form.

supervillain.lattice.interlaced.pull(data, shift)[source]

Pull data from shift steps away: result[n] = data[n + shift] (periodic).

Equivalent to push() with the sign of shift reversed.

Parameters

data (np.ndarray) – A \((2N)^D\) interlaced form array.
shift (tuple of int) – One integer per spatial direction.

Returns

Translated array of the same shape.

Return type

np.ndarray

supervillain.lattice.interlaced.push(data, shift)[source]

Translate data by shift steps: result[n + shift] = data[n] (periodic).

Parameters

data (np.ndarray) – A \((2N)^D\) interlaced form array.
shift (tuple of int) – One integer per spatial direction; \(D\) is inferred from data.ndim.

Returns

Translated array of the same shape.

Return type

np.ndarray

Exterior Derivative

At each \((p+1)\)-form output site \(\xi\) with odd directions \(O = (o_0, \ldots, o_p)\):

\[(df)[\xi] = \sum_{j=0}^{p} (-1)^j \bigl(f[\xi + \varepsilon_{o_j}] - f[\xi - \varepsilon_{o_j}]\bigr)\]

Each term is a centered (symmetric) difference in interlaced direction \(o_j\), spanning one physical lattice step.

supervillain.lattice.interlaced.d(data)[source]

Exterior derivative. Infers D from data.ndim.

For each (p+1)-form output site \(x\) with odd directions \(O = (o_0, \ldots, o_p)\), the contribution is

\[\sum_{j} (-1)^j \big( \texttt{data}[x + \varepsilon_{o_j}] - \texttt{data}[x - \varepsilon_{o_j}] \big)\]

where \(\varepsilon_k\) is the unit vector in direction \(k\).

Parameters: data (np.ndarray) – A \((2N)^D\) interlaced \(p\)-form array.
Returns: The \((2N)^D\) interlaced \((p+1)\)-form \(df\).
Return type: np.ndarray

This \(d\) is nilpotent \(d^2 = 0\) holds as checked by test_d_nilpotent in test/test_lattice_interlaced.py.

Codifferential

The codifferential is the formal adjoint of the exterior derivative,

(1)\[\langle da, b \rangle = \langle a, \delta b \rangle\]

tested by test_adjointness in test/test_lattice_interlaced.py.

At each \((p-1)\)-form output site \(\xi\) with odd directions \(M\) and even directions \(E = (e_0, \ldots, e_{D-p})\):

\[(\delta f)_M[\xi] = -\sum_{i=0}^{D-p} (-1)^{e_i - i} \bigl(f_{M \cup \{e_i\}}[\xi + \varepsilon_{e_i}] - f_{M \cup \{e_i\}}[\xi - \varepsilon_{e_i}]\bigr)\]

The sign \((-1)^{e_i - i}\) equals \((-1)^j\) where \(j = \#\{m \in M : m < e_i\}\) is the insertion position of \(e_i\) into \(M\).

supervillain.lattice.interlaced.delta(data)[source]

Codifferential (formal adjoint of d). Infers D from data.ndim.

For each \((p-1)\)-form output site \(x\) with even directions \(E = (e_0, \ldots, e_{D-p})\), the contribution is

\[-\sum_{i=0}^{D-p} (-1)^{e_i - i} \big( \texttt{data}[x + \varepsilon_{e_i}] - \texttt{data}[x - \varepsilon_{e_i}] \big).\]

Parameters: data (np.ndarray) – A \((2N)^D\) interlaced \(p\)-form array.
Returns: The \((2N)^D\) interlaced \((p-1)\)-form \(\delta f\).
Return type: np.ndarray

The nilpotency \(\delta^2 = 0\) holds as checked by test_delta_nilpotent in test/test_lattice_interlaced.py.

The continuum identity \(\delta = (-1)^{D(k+1)+1}\,\star\,d\,\star\) holds on the interlaced lattice only up to a translation; see the note below.

Hodge Star

For each output \((D-p)\)-form site \(\xi\) with odd directions \(J\), let \(I = \{0,\ldots,D-1\} \setminus J\) be its complement:

\[(\star f)_J[\xi] = \sigma(I \frown J)\; f_I[\xi - \varepsilon_{\mathrm{all}}]\]

where \(\varepsilon_{\mathrm{all}} = (1, \ldots, 1)\). The all-direction shift \(-\varepsilon_{\mathrm{all}}\) flips every coordinate’s parity, mapping the \((D-p)\)-form site \(\xi\) to the \(p\)-form site \(\xi - \varepsilon_{\mathrm{all}}\).

supervillain.lattice.interlaced.star(data)[source]

The Hodge star maps a \(p\)-form to a \((D-p)\)-form. On the interlaced lattice this entails (bijectively) shifting the data by \(\varepsilon_{\mathrm{all}} = (1,\ldots,1)\), flipping every coordinate’s parity.

The sign for each output site \(\xi\) with odd directions \(J\) depends only on \(J\):

\[(\star f)_J[\xi] = \sigma(I \frown J)\; f_I[\xi - \varepsilon_{\mathrm{all}}]\]

where \(I = \{0,\ldots,D-1\} \setminus J\) and \(\sigma(I \frown J)\) is the sign of the permutation sorting \((I \frown J)\).

Parameters: data (np.ndarray) – A \((2N)^D\) interlaced \(p\)-form array.
Returns: The \((2N)^D\) interlaced \((D-p)\)-form \(\star f\).
Return type: np.ndarray

The Hodge inner-product identity \(\langle a, b \rangle = \sum_{\xi} (a \wedge \star b)_{(0,\ldots,D-1)}[\xi]\) holds exactly; it is cross-validated by test_compact_star_matches_interlaced in test/test_lattice.py.

Wedge Product

At each \((n+m)\)-form output site \(\xi\) with odd directions \(O\), summing over all shuffles \(O = A \sqcup B\) (\(A\) the \(n\) a-directions, \(B\) the \(m\) b-directions):

\[(a \wedge b)[\xi] = \sum_{O = A \sqcup B} \sigma(B \frown A)\; a[\xi - \hat{\varepsilon}_A]\; b[\xi + \hat{\varepsilon}_B]\]

where \(\hat{\varepsilon}_A = \sum_{k \in A} \varepsilon_k\) and \(\hat{\varepsilon}_B = \sum_{k \in B} \varepsilon_k\). The sign \(\sigma(B \frown A)\) sorts \((B, A)\) back into \(O\); it equals \((-1)^{nm}\,\sigma(A \frown B)\). In contrast to the compact formula (which evaluates \(b\) at \(x + \hat{e}_A\)), the output site \(\xi\) lies between \(a\) at \(\xi - \hat{\varepsilon}_A\) and \(b\) at \(\xi + \hat{\varepsilon}_B\).

Why σ(B⌢A) and not σ(A⌢B)? That's what appears in the dense format...

The compact wedge evaluates \(a\) at the output site and shifts \(b\) to the far side of the \(a\)-cell:

\[(a \wedge b)_O[x] = \sum_{O = A \sqcup B} \sigma(A \frown B)\; a_A[x]\; b_B\!\left[x + \hat{e}_A\right]\]

The interlaced wedge is instead centered: \(a\) is pulled back to \(\xi - \hat{\varepsilon}_A\) and \(b\) pushed forward to \(\xi + \hat{\varepsilon}_B\). Shifting \(a\)’s evaluation point from \(\xi\) to \(\xi - \hat{\varepsilon}_A\) introduces exactly the factor \((-1)^{nm}\) that converts \(\sigma(A \frown B)\) into \(\sigma(B \frown A)\), so both formulas give the same physical result.

A direct check for \(D=2\), \(n=m=1\), \(O=(0,1)\):

shuffle	compact	interlaced (translated to physical)
\(A=(0),\,B=(1)\)	\(+a_0[x]\,b_1[x+\hat{e}_0]\)	\(-a_1[x]\,b_0[x+\hat{e}_1]\)
\(A=(1),\,B=(0)\)	\(-a_1[x]\,b_0[x+\hat{e}_1]\)	\(+a_0[x]\,b_1[x+\hat{e}_0]\)

The rows swap between the two, but the total \(a_0[x]\,b_1[x+\hat{e}_0] - a_1[x]\,b_0[x+\hat{e}_1]\) is the same.

supervillain.lattice.interlaced.wedge(n, m, a, b)[source]

Wedge product of an \(n\)-form \(a\) and an \(m\)-form \(b\). Infers \(D\) from \texttt{a.ndim}.

At each \((n+m)\)-form output site \(x\) with odd directions \(O\), summing over all shuffles \(O = A \sqcup B\) (\(A\) the \(n\) directions of \(a\), \(B\) the \(m\) directions of \(b\)):

\[(a \wedge b)[x] = \sum_{O = A \sqcup B} \sigma(B \frown A)\; a[x - \hat{\varepsilon}_A]\; b[x + \hat{\varepsilon}_B]\]

where \(\hat{\varepsilon}_A = \sum_{k \in A} \varepsilon_k\), \(\hat{\varepsilon}_B = \sum_{k \in B} \varepsilon_k\), and \(\sigma(B \frown A)\) is the sign of the permutation sorting \((B, A)\) back into \(O\).

Parameters

n (int) – Degree of \(a\).
m (int) – Degree of \(b\).
a (np.ndarray) – A \((2N)^D\) interlaced \(n\)-form array.
b (np.ndarray) – A \((2N)^D\) interlaced \(m\)-form array on the same lattice.

Returns

The \((2N)^D\) interlaced \((n+m)\)-form \(a \wedge b\).

Return type

np.ndarray

Bilinearity, associativity, and the Leibniz rule \(d(a \wedge b) = da \wedge b + (-1)^n\,a \wedge db\) all hold exactly; they are cross-validated by test_compact_wedge_matches_interlaced in test/test_lattice.py.

Unlike the continuum, however, the lattice wedge product is not anti-commutative; see the note below.

Differences from the Continuum

In the continuum, on a Riemannian manifold with Euclidean signature, the Hodge star squares to

\[\star\star a = (-1)^{p(D-p)}\, a \qquad \text{on a } p\text{-form.}\]

Danger

On the lattice the identity picks up a spatial shift. In interlaced coordinates:

\[(\star\star f)[\xi] = (-1)^{p(D-p)}\, f[\xi - 2\varepsilon_{\mathrm{all}}]\]

a backward shift of \(2\varepsilon_{\mathrm{all}} = (2,\ldots,2)\) — one physical lattice step in each direction. The shift drops out of any periodic sum. This is tested by test_star_star in test/test_lattice_interlaced.py.

In the continuum we expect anti/commutativity

\[(a \wedge b) = (-1)^{n m} (b \wedge a)\]

to hold.

Danger

On the lattice this property fails!

In the continuum, on a Riemannian manifold with Euclidean signature, the Hodge star satisfies

\[\delta = (-1)^{D(k+1)+1}\, \star\, d\, \star \qquad \text{on } k\text{-forms}.\]

Danger

On the lattice this identity picks up a spatial shift. In interlaced coordinates the shift is \(2\varepsilon_{\mathrm{all}} = (2, \ldots, 2)\) — one physical step in each direction:

\[\delta = (-1)^{D(k+1)+1}\, T_{2\varepsilon_{\mathrm{all}}}\, \star\, d\, \star\]

where \(T_{\Delta\xi}\) is the translation operator (pull() by \(\Delta\xi\)). The shift drops out of any periodic sum, which is why (1) holds exactly despite it. This is tested by test_star_d_star_equals_shifted_delta in test/test_lattice_interlaced.py.