#!/usr/bin/env python
import numpy as np
import supervillain.action
from supervillain.generator import Generator
from supervillain.h5 import ReadWriteable
from supervillain.lattice import delta
import logging
logger = logging.getLogger(__name__)
[docs]class CoexactUpdate(ReadWriteable, Generator):
r'''
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
.. math ::
\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 :class:`~.worldline.VortexUpdate`)
nor can it produce all coclosed changes---only coexact changes.
For a coïnexact coclosed update consider the :class:`~.worldline.WrappingUpdate` or the :class:`worm <supervillain.generator.worldline.worm.Classic>`.
'''
def __init__(self, action, interval_t = 1):
if not isinstance(action, supervillain.action.Worldline):
raise ValueError('Need a Worldline action')
self.Action = action
self.Lattice = action.Lattice
self.kappa = action.kappa
self.interval_t = interval_t
self.ts = tuple(t for t in range(-interval_t, 0)) + tuple(t for t in range(1, interval_t+1))
self.rng = np.random.default_rng()
self.accepted = 0
self.proposed = 0
self.acceptance = 0.
self.sweeps = 0
def __str__(self):
return 'SiteUpdate'
[docs] def step(self, cfg):
r'''
Make a volume's worth of locally-exact updates to m.
Parameters
----------
cfg: dict
A dictionary with m and v field variables.
Returns
-------
dict
Another configuration of fields.
'''
self.sweeps += 1
total_acceptance = 0
total_accepted = 0
v = cfg['v']
delta_v_by_W = delta(v) / self.Action._W
m = cfg['m'].copy()
L = self.Lattice
n_comps = len(L.components[2])
# One independent Metropolis draw per (2-form component, site).
metropolis = self.rng.uniform(0, 1, (n_comps,) + L.dims)
# The idea is to make coordinated changes to m that keep δm=0. We can do that by letting the change in m
# be a coexact form δt with t a two-form so that the change in δm is δ^2t = 0.
# However, rather than a python-level for loop over space, we can accomplish a lot more at the numpy level,
# as in the villain.ExactUpdate.
# Since the coordinated updates of m are derived from t we can think like we think in the ExactUpdate:
# What enters the change in action (per link) is δt, which knows about t on two plaquettes.
#
# That poses a small problem because if we change the action by changing m, we want to be able to track
# that change back to a change in t on ONE particular plaquette, and to accept or reject that change independently
# from other changes in t. Therefore, we use checkerboarding.
#
# In D>2 there are C(D,2) independent 2-form components. Two plaquettes of the same component at
# same-color sites never share boundary links, so we can propose and accept/reject each component
# independently. We process components sequentially to avoid conflicts from shared boundary links
# between different components at the same site.
for color in L.checkerboarding:
for comp_idx in range(n_comps):
# We only offer changes to t[comp_idx] on a single color at once.
t = L.form(2, dtype=int)
t[comp_idx][color] = self.rng.choice(self.ts, len(color[0]))
# To keep δm=0 we let the change in m be given by δt, so that δ(change_m) = δ^2(t) = 0.
change_m = delta(t)
dS_link = 0.5 / self.Action.kappa * change_m * (2*(m - delta_v_by_W) + change_m)
# The change in action originating from each plaquette is the sum of changes on its
# boundary links. coface_sum() accumulates those unsigned boundary contributions,
# giving a 2-form where dS[comp_idx][x] is the ΔS for the (comp_idx) plaquette at x.
dS = dS_link.coface_sum()
# dS is not 0 on off-color plaquettes---those still have links touching the current color.
# Only accept/reject on the current color.
acceptance = np.clip(np.exp(-dS[comp_idx][color]), a_min=0, a_max=1)
accepted = (metropolis[comp_idx][color] < acceptance)
total_accepted += accepted.sum()
total_acceptance += acceptance.sum()
# Update m where the change is accepted.
t[comp_idx][color] *= accepted
m = m + delta(t)
self.proposed += L.cells_of_degree[2]
self.acceptance += total_acceptance / L.cells_of_degree[2]
self.accepted += total_accepted
logger.debug(f'Average proposal acceptance {total_acceptance / L.cells_of_degree[2]:.6f}; Actually accepted {total_accepted} / {L.cells_of_degree[2]} = {total_accepted / L.cells_of_degree[2]}')
return cfg | {'m': m}
[docs] def report(self):
return (
f'There were {self.accepted} coexact proposals accepted of {self.proposed} proposed updates.'
+'\n'+
f' {self.accepted/self.proposed:.6f} acceptance rate'
+'\n'+
f' {self.acceptance / self.sweeps:.6f} average Metropolis acceptance probability.'
)