#!/usr/bin/env python
# coding: utf-8
#
# compact.py — Differential forms on a D-dimensional hypercubic lattice,
# stored without wasted space.
#
# MOTIVATION
# ----------
# The interlaced representation (see interlaced.py) places p-form data at
# sites with exactly p odd coordinates in a (2N)^D array. That is correct
# but sparse: only C(D,p)/2^D of array elements are used. Here we store a
# p-form compactly as an array
# of shape
#
# (C(D,p), N, N, ..., N) [D spatial axes follow]
#
# The first axis indexes the C(D,p) "components", listed in lexicographic
# order by the sorted tuple of directions that are "form directions".
#
# Example — D=3:
# 0-form components: () → shape (1, N, N, N)
# 1-form components: (0,),(1,),(2,) → shape (3, N, N, N)
# 2-form components: (0,1),(0,2),(1,2) → shape (3, N, N, N)
# 3-form components: (0,1,2) → shape (1, N, N, N)
from functools import cached_property
from itertools import combinations, permutations as _permutations
from math import comb, factorial
import numpy as np
from supervillain.h5 import ReadWriteable
from supervillain.lattice import _dimension
def _hyperoctant_pair_mask(coords, b, D):
"""
Sites in one pair of opposite hyperoctants, indexed by b in 0..2^(D-1)-1.
The representative hyperoctant has coords[0] >= 0; bits of b fix the sign
of coords[1], …, coords[D-1]. The paired hyperoctant flips every sign.
"""
mask_pos = coords[0] >= 0
mask_neg = coords[0] < 0
for k in range(1, D):
bit = (b >> (k - 1)) & 1
if bit == 0:
mask_pos &= coords[k] >= 0
mask_neg &= coords[k] < 0
else:
mask_pos &= coords[k] < 0
mask_neg &= coords[k] >= 0
return mask_pos | mask_neg
# ---------------------------------------------------------------------------
# Lattice
# ---------------------------------------------------------------------------
[docs]class Lattice(ReadWriteable):
"""
A D-dimensional hypercubic lattice with N sites per direction.
The lattice knows how to enumerate the components of any p-form and
can create zero or random p-forms on demand.
"""
def __init__(self, D, N):
self.D = D
self.N = N
# components[p] — ordered list of component tuples for p-forms.
# Each tuple is a sorted sequence of p direction indices (0-based).
# Length = C(D, p). Order is lexicographic.
self.components = {
p: list(combinations(range(D), p))
for p in range(D + 1)
}
# comp_index[p][dirs] — reverse map: tuple of dirs → integer index.
# Use comp_index[p][(i, j)] to find where component (i,j) lives
# along axis 0 of a p-form array.
self.comp_index = {
p: {comp: idx for idx, comp in enumerate(self.components[p])}
for p in range(D + 1)
}
[docs] def to_h5(self, group, _top=True):
"""Write D and N to HDF5; everything else is recomputed on load."""
from supervillain.h5 import Data
Data.write(group, 'D', self.D)
Data.write(group, 'N', self.N)
[docs] @classmethod
def from_h5(cls, group, strict=True, _top=True):
"""Reconstruct from HDF5 by reading D and N reconstructing the rest."""
from supervillain.h5 import Data
D = Data.read(group['D'], strict)
N = Data.read(group['N'], strict)
return cls(D=D, N=N)
@cached_property
def sites(self):
"""Total number of sites $N^D$."""
return self.N ** self.D
@cached_property
def dims(self):
"""Tuple of side lengths (N, N, ..., N), length D."""
return (self.N,) * self.D
@cached_property
def origin(self):
r"""The origin site $(0, \ldots, 0)$, a tuple of length $D$. Use it to index the
zero-displacement element of a scalar field or correlator of shape $(N,)^D$
(e.g. ``correlator[lattice.origin]``) without assuming $D=2$. A
:class:`~supervillain.lattice.Form` has a leading component axis, so index the
component first: ``form[c][lattice.origin]``."""
return (0,) * self.D
@property
def dim(self):
"""Number of dimensions (alias for D)."""
return self.D
@cached_property
def links(self):
"""Total number of links (1-form components) $D N^D$."""
return self.D * self.sites
@cached_property
def cells_of_degree(self):
"""
dict[int, int]: p → number of p-cells = $C(D,p) N^D$."""
return {p: comb(self.D, p) * self.sites for p in range(self.D + 1)}
@cached_property
def cells_of_codegree(self):
"""
dict[int, int]: q → number of (D−q)-cells = $C(D, D−q) N^D$."""
return {q: self.cells_of_degree[self.D - q] for q in range(self.D + 1)}
@cached_property
def coords(self):
"""FFT-convention coordinate for each lattice site, shape (D, N, …, N)."""
return np.stack(
np.meshgrid(*(_dimension(self.N) for _ in range(self.D)), indexing='ij'),
axis=0,
)
@cached_property
def checkerboarding(self):
r"""
Partition lattice sites into colors with no same-color nearest neighbors.
On a lattice of even size both the sites and plaquettes can be bipartitioned so that no
simplex of one color has a neighbor of the same color. On an odd-sized lattice the
periodic boundary conditions make it impossible to accomplish this with only 2 colors, but
a similar construction with more colors is possible.
With even $N$ there are two colors (coordinate-sum parity); with odd $N$ there are $2^{\max(D, 2)}$
to ensure that no nearest-neighbor sites share a color.
The figure below shows both cases for D=2:
.. plot:: example/plot/checkerboarding.py
Returns
-------
tuple of index-array tuples
One ``np.where`` result per color; each element selects that
color's sites from the D spatial axes of a form::
for i, color in enumerate(L.checkerboarding):
form[(slice(None), *color)] = i
.. warning::
No promise is made about the future sizes of the color partitions.
For example, it might be wiser for performance to split the odd-N colors less evenly.
All that is promised is that within each color no site shares a nearest-neighbor edge
with a site of the same color.
"""
D, N = self.D, self.N
coords = self.coords
parity = np.mod(self.coords.sum(axis=0), 2)
if N % 2 == 0:
return tuple(np.where(parity == c) for c in (0, 1))
colors = []
n_pairs = 1 << max(D - 1, 1)
for b in range(n_pairs):
if D == 1:
pair = coords[0] >= 0 if b == 0 else coords[0] < 0
else:
pair = _hyperoctant_pair_mask(list(coords), b, D)
for c in (0, 1):
colors.append(np.where(pair & (parity == c)))
return tuple(colors)
# ------------------------------------------------------------------
# Factory methods
[docs] def zeros(self, p, dtype=float):
"""
Return a zero p-form.
Parameters
----------
p : int
Form degree.
dtype : data-type, optional
np dtype for the underlying array (default ``float``).
Returns
-------
Form
Zero p-form of shape ``(C(D,p), N,...,N)``.
"""
shape = (comb(self.D, p),) + (self.N,) * self.D
return Form(np.zeros(shape, dtype=dtype), degree=p, lattice=self)
[docs] def random(self, p):
"""
Return a p-form with entries drawn uniformly from [0, 1).
Parameters
----------
p : int
Form degree.
Returns
-------
Form
Random p-form of shape ``(C(D,p), N,...,N)``.
"""
shape = (comb(self.D, p),) + (self.N,) * self.D
return Form(np.random.random(shape), degree=p, lattice=self)
@cached_property
def _coord_1d(self):
return _dimension(self.N)
@cached_property
def R_squared(self):
"""Distance-squared from the origin at each site, shape (N,...,N)."""
return np.sum(self.coords**2, axis=0)
[docs] def distance_squared(self, a, b):
r"""Squared lattice distance between coordinate vectors ``a`` and ``b``.
Accounts for periodic boundary conditions: the distance is computed
via ``mod(a - b)``, so it is the shortest-path distance on the torus.
Parameters
----------
a, b : array_like
Coordinate vectors of shape ``(D,)`` for a single pair, or
``(..., D)`` for a batch.
Returns
-------
np.ndarray
Scalar for a single pair, shape ``(...)`` for a batch.
"""
d = self.mod(np.asarray(a) - np.asarray(b))
return np.sum(d**2, axis=-1)
@cached_property
def coordinates(self):
"""Array of shape (sites, D) listing every site's coordinates."""
return np.stack(
[c.flatten() for c in np.meshgrid(*[self._coord_1d] * self.D, indexing='ij')],
axis=1,
)
[docs] def mod(self, x):
"""
Map integer coordinates into FFT-convention lattice values.
Parameters
----------
x : array_like
Integer site coordinates (any shape).
Returns
-------
np.ndarray
Coordinates in the range ``[-(N//2), N//2)`` (FFT convention).
"""
x = np.asarray(x)
modded = np.mod(x, self.N)
return self._coord_1d[modded]
# ------------------------------------------------------------------
# Fourier methods
#
# Convention matches Lattice2D: ortho normalization, spatial axes last.
# For a p-form the spatial axes are the last D axes; the default axes
# argument reflects this so callers seldom need to override it.
def _spatial_axes(self):
"""The last D axes, i.e. the spatial directions of any p-form."""
return tuple(range(-self.D, 0))
[docs] def fft(self, form, axes=None):
r"""
D-dimensional DFT over the spatial axes of ``form``.
.. math::
F_{\boldsymbol{k}} = \frac{1}{N^{D/2}}
\sum_{\boldsymbol{x}} e^{-2\pi i \boldsymbol{k}\cdot\boldsymbol{x}/N} f_{\boldsymbol{x}}
Parameters
----------
form: np.ndarray
The data to transform. Spatial axes are the last D axes.
axes: tuple of int, optional
Override which axes to transform. Defaults to the last D axes.
Returns
-------
np.ndarray
The Fourier-transformed array (complex).
"""
return np.fft.fftn(form, axes=(axes if axes is not None else self._spatial_axes()), norm='ortho')
[docs] def ifft(self, form, axes=None):
r"""
D-dimensional inverse DFT over the spatial axes of ``form``.
Parameters
----------
form : np.ndarray
The data to transform. Spatial axes are the last D axes.
axes : tuple of int, optional
Override which axes to transform. Defaults to the last D axes.
Returns
-------
np.ndarray
The inverse-Fourier-transformed array (complex).
"""
return np.fft.ifftn(form, axes=(axes if axes is not None else self._spatial_axes()), norm='ortho')
[docs] def convolution(self, f, g, axes=None):
r"""
The `convolution <https://en.wikipedia.org/wiki/Convolution>`_ is given by
.. math::
\texttt{convolution}(f,g)(\boldsymbol{r})
= (f * g)(\boldsymbol{r})
= \sum_{\boldsymbol{x}} f(\boldsymbol{x})\,g(\boldsymbol{r}-\boldsymbol{x})
.. collapse:: The convolution is Fourier accelerated.
:class: note
With the ortho-normalized DFT convention
:math:`f(\boldsymbol{x}) = N^{-D/2} \sum_{\boldsymbol{k}} F_{\boldsymbol{k}}\, e^{2\pi i\,\boldsymbol{k}\cdot\boldsymbol{x}/N}`,
.. math::
\begin{aligned}
(f * g)(\boldsymbol{r})
&= \sum_{\boldsymbol{x}}
\Bigl(\tfrac{1}{N^{D/2}}\sum_{\boldsymbol{k}} F_{\boldsymbol{k}}\,
e^{2\pi i\,\boldsymbol{k}\cdot\boldsymbol{x}/N}\Bigr)
\Bigl(\tfrac{1}{N^{D/2}}\sum_{\boldsymbol{k}'} G_{\boldsymbol{k}'}\,
e^{2\pi i\,\boldsymbol{k}'\cdot(\boldsymbol{r}-\boldsymbol{x})/N}\Bigr)
\\
&= \frac{1}{N^D}\sum_{\boldsymbol{k},\boldsymbol{k}'}
F_{\boldsymbol{k}} G_{\boldsymbol{k}'}\,
e^{2\pi i\,\boldsymbol{k}'\cdot\boldsymbol{r}/N}
\underbrace{
\sum_{\boldsymbol{x}} e^{2\pi i\,(\boldsymbol{k}-\boldsymbol{k}')\cdot\boldsymbol{x}/N}
}_{N^D\,\delta_{\boldsymbol{k}\boldsymbol{k}'}}
\\
&= \sum_{\boldsymbol{k}} F_{\boldsymbol{k}} G_{\boldsymbol{k}}\,
e^{2\pi i\,\boldsymbol{k}\cdot\boldsymbol{r}/N}
\\
&= \sqrt{N^D}\times
\underbrace{\tfrac{1}{N^{D/2}}\sum_{\boldsymbol{k}}
\bigl(F_{\boldsymbol{k}} G_{\boldsymbol{k}}\bigr)\,
e^{2\pi i\,\boldsymbol{k}\cdot\boldsymbol{r}/N}
}_{\texttt{ifft}(\hat f \hat g)(\boldsymbol{r})}
\\
\texttt{convolution}(f,g) &= \sqrt{N^D}\;\texttt{ifft}\!\bigl(\texttt{fft}(f)\cdot\texttt{fft}(g)\bigr)
\end{aligned}
Parameters
----------
f, g: np.ndarray
Forms on this lattice; spatial axes are the last D axes.
axes: tuple of int, optional
Returns
-------
np.ndarray
"""
ax = axes if axes is not None else self._spatial_axes()
return np.sqrt(self.sites) * self.ifft(self.fft(f, axes=ax) * self.fft(g, axes=ax), axes=ax)
[docs] def correlation(self, f, g, axes=None):
r"""
The `cross-correlation <https://en.wikipedia.org/wiki/Cross-correlation>`_ is given by
.. math::
\texttt{correlation}(f,g)(\boldsymbol{r})
= (f \star g)(\boldsymbol{r})
= \frac{1}{N^D} \sum_{\boldsymbol{x}} f(\boldsymbol{x})^*\,g(\boldsymbol{x}-\boldsymbol{r})
where :math:`f^*` is the complex conjugate of :math:`f`.
.. collapse:: The cross-correlation is Fourier accelerated.
:class: note
With the ortho-normalized DFT convention
:math:`f(\boldsymbol{x}) = N^{-D/2} \sum_{\boldsymbol{k}} F_{\boldsymbol{k}}\, e^{2\pi i\,\boldsymbol{k}\cdot\boldsymbol{x}/N}`,
.. math::
\begin{aligned}
(f \star g)(\boldsymbol{r})
&= \frac{1}{N^D}\sum_{\boldsymbol{x}}
\Bigl(\tfrac{1}{N^{D/2}}\sum_{\boldsymbol{k}} F_{\boldsymbol{k}}\,
e^{2\pi i\,\boldsymbol{k}\cdot\boldsymbol{x}/N}\Bigr)^*
\Bigl(\tfrac{1}{N^{D/2}}\sum_{\boldsymbol{k}'} G_{\boldsymbol{k}'}\,
e^{2\pi i\,\boldsymbol{k}'\cdot(\boldsymbol{x}-\boldsymbol{r})/N}\Bigr)
\\
&= \frac{1}{N^{2D}}\sum_{\boldsymbol{k},\boldsymbol{k}'}
F_{\boldsymbol{k}}^* G_{\boldsymbol{k}'}\,
e^{-2\pi i\,\boldsymbol{k}'\cdot\boldsymbol{r}/N}
\underbrace{
\sum_{\boldsymbol{x}} e^{2\pi i\,(\boldsymbol{k}'-\boldsymbol{k})\cdot\boldsymbol{x}/N}
}_{N^D\,\delta_{\boldsymbol{k}\boldsymbol{k}'}}
\\
&= \frac{1}{N^D}\sum_{\boldsymbol{k}}
F_{\boldsymbol{k}}^* G_{\boldsymbol{k}}\,
e^{-2\pi i\,\boldsymbol{k}\cdot\boldsymbol{r}/N}
\\
&= \frac{1}{\sqrt{N^D}}\times
\underbrace{\tfrac{1}{N^{D/2}}\sum_{\boldsymbol{k}}
\bigl(F_{\boldsymbol{k}}^* G_{\boldsymbol{k}}\bigr)\,
e^{-2\pi i\,\boldsymbol{k}\cdot\boldsymbol{r}/N}
}_{\texttt{fft}(\hat f^*\!\cdot\hat g)(\boldsymbol{r})}
\\
\texttt{correlation}(f,g) &= \texttt{fft}\!\bigl(\texttt{fft}(f)^*\cdot\texttt{fft}(g)\bigr) / \sqrt{N^D}
\end{aligned}
.. warning::
We have :math:`g(\boldsymbol{x}-\boldsymbol{r})` whereas
`Wikipedia <https://en.wikipedia.org/wiki/Cross-correlation>`_ has
:math:`g(\boldsymbol{x}+\boldsymbol{r})`.
The difference is just the sign of the relative coordinate.
.. warning::
We normalize by the spacetime volume :math:`N^D`;
`Wikipedia <https://en.wikipedia.org/wiki/Cross-correlation>`_ does not.
Parameters
----------
f, g: np.ndarray
Forms on this lattice; spatial axes are the last D axes.
axes: tuple of int, optional
Returns
-------
np.ndarray
"""
ax = axes if axes is not None else self._spatial_axes()
return self.fft(self.fft(f, axes=ax).conj() * self.fft(g, axes=ax), axes=ax) / np.sqrt(self.sites)
[docs] def linearize(self, v, dims=(-1,)):
r"""Flatten the D spatial dims of ``v`` into a single axis of size ``sites``.
Parameters
----------
v : np.ndarray
Array whose D adjacent spatial dims will be collapsed into one.
dims : tuple of int
Axes of the *result* that come from flattening. Pass the same value
to :meth:`coordinatize` for a round-trip.
Returns
-------
np.ndarray
"""
shape = v.shape
v_dims = len(shape)
dm = set(dims)
future_dims = v_dims - (self.D - 1) * len(dm)
dm = set(d % future_dims for d in dm)
new_shape = []
idx = 0
for i in range(future_dims):
if i not in dm:
new_shape.append(shape[idx])
idx += 1
else:
new_shape.append(self.sites)
idx += self.D
return v.reshape(new_shape)
[docs] def coordinatize(self, v, dims=(-1,), center_origin=False):
r"""Unflatten a linear site-index axis back into D spatial axes.
Parameters
----------
v : np.ndarray
Array with a ``sites``-sized axis to expand.
dims : tuple of int
Axes of ``v`` to unflatten. Each must have size ``sites``.
center_origin : bool
If True, roll each spatial block so that the origin sits in the
centre of the array. Useful for visualisation; not invertible.
Returns
-------
np.ndarray
"""
v_dims = len(v.shape)
to_reshape = np.sort(np.remainder(np.array(dims), v_dims))
new_shape = ()
for i, s in enumerate(v.shape):
new_shape += ((s,) if i not in to_reshape else self.dims)
reshaped = v.reshape(new_shape)
if not center_origin:
return reshaped
# Each original axis at position a expands to D axes; subsequent blocks
# shift right by (D-1) per earlier expansion.
axes = to_reshape + np.arange(len(to_reshape)) * (self.D - 1)
for a in axes:
for d in range(self.D):
reshaped = np.roll(reshaped, self.N // 2, axis=int(a) + d)
return reshaped
@cached_property
def _hyperoctahedral_permutations(self):
"""Site-index permutations for each of the D!·2^D hyperoctahedral group elements."""
coords = self.coordinates # (sites, D)
coord_to_idx = {tuple(c): k for k, c in enumerate(coords)}
result = []
# The hyperoctahedral group B_D = (Z/2Z)^D ⋊ S_D combines the D! permutations
# of coordinate axes with 2^D independent per-axis sign flips.
for perm in _permutations(range(self.D)):
# signs is NOT the sign (parity) of the permutation — it is D independent
# bits, one per axis, encoding which coordinates get negated.
for signs in np.ndindex(*([2] * self.D)):
# Map the 0/1 bits to +1/−1: a signed permutation x → sign_vec * x[perm].
sign_vec = np.array([1 - 2 * s for s in signs])
idx_perm = np.array([
coord_to_idx[tuple(self.mod(sign_vec * coords[i][list(perm)]))]
for i in range(self.sites)
])
result.append(idx_perm)
return result
[docs] def symmetrize(self, correlator, dims=(-1,)):
r"""Average ``correlator`` over the hyperoctahedral group (D!·2^D signed permutations).
Projects onto the totally-symmetric (A₁/trivial) irrep of the lattice point
group.
.. plot:: example/plot/symmetrize.py
Parameters
----------
correlator : np.ndarray
Spatial shape ``(N,...,N)`` (or any array whose ``dims`` axes span the
sites of this lattice after linearization).
dims : tuple of int
Axes to symmetrize (same convention as :meth:`linearize`).
Returns
-------
np.ndarray
Same shape as ``correlator``.
"""
C = self.linearize(correlator, dims=dims)
v_dims = len(C.shape)
sites_axis = list(dims)[0] % v_dims
perms = self._hyperoctahedral_permutations
result = np.sum([np.take(C, p, axis=sites_axis) for p in perms], axis=0)
return self.coordinatize(result / len(perms), dims=dims)
def __repr__(self):
return f"Lattice(D={self.D}, N={self.N})"
# ---------------------------------------------------------------------------
# Form (numpy.ndarray subclass)
# ---------------------------------------------------------------------------
# ---------------------------------------------------------------------------
# Translation operators
# ---------------------------------------------------------------------------
[docs]def push(form, shift):
r"""Translate the form forward by $\Delta x$.
.. math ::
\texttt{push}(f, \Delta x)[\ldots, x] = f[\ldots, x - \Delta x] \quad \text{(periodic)}
Parameters
----------
form : np.ndarray
Array whose last ``len(shift)`` axes are the spatial directions.
shift : sequence of int
One integer per spatial direction.
Returns
-------
np.ndarray
"""
result = form
for i, s in enumerate(shift):
if s:
result = np.roll(result, s, axis=i - len(shift))
return result
[docs]def pull(form, shift):
r"""Translation operator $T_{\Delta x}$: pull content from position $x + \Delta x$ to $x$.
.. math ::
T_{\Delta x} f[\ldots, x] = \texttt{pull}(f, \Delta x)[\ldots, x] = f[\ldots, x + \Delta x] \quad \text{(periodic)}
Parameters
----------
form : np.ndarray
Array whose last ``len(shift)`` axes are the spatial directions.
shift : sequence of int
One integer per spatial direction.
Returns
-------
np.ndarray
"""
return push(form, tuple(-s for s in shift))
# ---------------------------------------------------------------------------
# Exterior derivative d : Ω^p → Ω^{p+1}
# ---------------------------------------------------------------------------
[docs]def d(f):
r"""
The exterior derivative of a p-form, a (p+1)-form.
For each output component $O = (o_0, \ldots, o_p)$:
.. math::
(df)_O[x] = \sum_{j=0}^{p} (-1)^j \, \Delta_{o_j} f_{O \setminus \{o_j\}}[x]
where $\Delta_k A[x] = A[x + \hat{e}_k] - A[x]$ is the forward finite difference.
The sign $(-1)^j$ is the signature of the permutation sorting $o_j$
into the remaining directions (see :ref:`sign-conventions`).
Parameters
----------
f : Form
A p-form on a :class:`Lattice`.
Returns
-------
Form
The (p+1)-form $df$, or the scalar ``0`` if $f$ is a D-form.
"""
lat = f.lattice
p = f.degree
if p == lat.D:
return 0
# d is an exact ±1 integer combination of finite differences, so it preserves the
# input dtype (an integer form stays integer).
result = lat.zeros(p + 1, dtype=f.dtype)
for out_comp in lat.components[p + 1]: # each (p+1)-form component
out_idx = lat.comp_index[p + 1][out_comp]
for j, k_j in enumerate(out_comp):
# Remove direction k_j from out_comp to get the source p-form component.
in_comp = tuple(k for k in out_comp if k != k_j)
in_idx = lat.comp_index[p][in_comp]
sign = (-1) ** j
# Forward finite difference of f[in_idx] in direction k_j.
# np.roll(A, -1, axis=k) shifts data so result[n] = A[n+1 mod N].
spatial = f[in_idx] # shape (N,...,N)
fwd_diff = np.roll(spatial, -1, axis=k_j) - spatial
result[out_idx] += sign * fwd_diff
return result
# ---------------------------------------------------------------------------
# Codifferential δ : Ω^p → Ω^{p-1}
# ---------------------------------------------------------------------------
[docs]def delta(f):
r"""
Codifferential (formal adjoint of d) of a p-form, returning a (p-1)-form.
For each output component $M = (m_0, \ldots, m_{p-1})$ and each direction
$e \notin M$, let $j = \#\{m \in M : m < e\}$ (the position where $e$
would be inserted to keep $M \cup \{e\}$ sorted):
.. math::
(\delta f)_M[x] = - \sum_{e \notin M} (-1)^j \, \nabla^*_e f_{M \cup \{e\}}[x]
where $\nabla^*_e A[x] = A[x] - A[x - \hat{e}_e]$ is the backward finite difference.
With the overall minus, $\delta$ is the formal adjoint of :func:`d` under the componentwise inner product (:eq:`d-delta-formal-adjoint`).
Parameters
----------
f : Form
A p-form on a :class:`Lattice`.
Returns
-------
Form
The (p-1)-form $\delta f$, or the scalar ``0`` if $f$ is a 0-form.
"""
lat = f.lattice
p = f.degree
if p == 0:
return 0
# δ is an exact ±1 integer combination of finite differences, so it preserves the
# input dtype (an integer form stays integer).
result = lat.zeros(p - 1, dtype=f.dtype)
all_dirs = set(range(lat.D))
for out_comp in lat.components[p - 1]: # each (p-1)-form component
out_idx = lat.comp_index[p - 1][out_comp]
M_set = set(out_comp)
for e in sorted(all_dirs - M_set): # directions not in M
# Position where e is inserted into sorted(M ∪ {e}).
j = sum(1 for m in out_comp if m < e)
sign = (-1) ** j
in_comp = tuple(sorted(M_set | {e}))
in_idx = lat.comp_index[p][in_comp]
# Backward finite difference of F[in_idx] in direction e.
# np.roll(A, +1, axis=e) shifts so result[n] = A[n-1 mod N].
spatial = f[in_idx]
bwd_diff = spatial - np.roll(spatial, +1, axis=e)
result[out_idx] -= sign * bwd_diff
return result
δ = delta
# ---------------------------------------------------------------------------
# Hodge–de Rham Laplacian Δ : Ω^p → Ω^p
# ---------------------------------------------------------------------------
[docs]def laplacian(f):
r"""
Hodge–de Rham Laplacian of a p-form, a p-form of the same degree.
The Laplacian (or Laplace–de Rham operator) is
.. math::
\Delta = d\delta + \delta d,
the symmetric combination of the :func:`exterior derivative <d>` and its
formal adjoint the :func:`codifferential <delta>`. Because :func:`delta`
is the adjoint of :func:`d`, the Laplacian is self-adjoint and positive
semidefinite under the componentwise inner product,
.. math::
\langle \Delta f, f \rangle
= \langle d f, d f \rangle + \langle \delta f, \delta f \rangle \geq 0.
Parameters
----------
f : Form
A p-form on a :class:`Lattice`.
Returns
-------
Form
The p-form $\Delta f = (d\delta + \delta d) f$.
"""
lat = f.lattice
D = lat.D
# Rather than composing d and δ, evaluate Δ directly: on a flat periodic
# lattice d and δ are constant-coefficient combinations of the commuting
# shift operators T_k A[x] = A[x + ê_k]. With d = Σ_k (ê_k∧) ∂_k for the
# forward difference ∂_k = T_k − 1 and δ = −Σ_k ι_k ∂*_k for the backward
# difference ∂*_k = 1 − T_k⁻¹, the cross terms cancel through the
# anticommutator {ê_k∧, ι_l} = δ_kl, leaving
# dδ + δd = −Σ_k ∂_k ∂*_k = −Σ_k (T_k − 2 + T_k⁻¹).
# So Δ acts diagonally on every component I as the negative of the ordinary
# nearest-neighbor scalar Laplacian — no mixing between the C(D,p) components:
# (Δf)_I[x] = Σ_k (2 f_I[x] − f_I[x + ê_k] − f_I[x − ê_k]).
# This is 2D array shifts, independent of p, and agrees with dδ + δd
# (test_laplacian_matches_d_delta).
#
# Δ is an exact integer combination of shifts, so it preserves the input dtype.
# Spatial axes are the last D axes, so direction k is axis k - D.
result = (2 * D) * f
for k in range(D):
axis = k - D
result = result - np.roll(f, -1, axis=axis) - np.roll(f, +1, axis=axis)
return result
# ---------------------------------------------------------------------------
# Hodge star ★ : Ω^p → Ω^{D-p}
# ---------------------------------------------------------------------------
def _perm_sign(seq):
"""Sign of the permutation that sorts a sequence of distinct integers."""
return (-1) ** sum(
1 for i in range(len(seq)) for j in range(i + 1, len(seq))
if seq[i] > seq[j]
)
[docs]def star(f):
r"""
Hodge star of a p-form, a (D-p)-form.
For each output component $J$ (a sorted $(D-p)$-tuple of directions),
let $I$ be the complement of $J$ in $\{0, \ldots, D-1\}$:
.. math::
(\star f)_J[x] = \sigma(I \frown J) \; f_I[x - \hat{e}_I]
where $\sigma(I \frown J) = (-1)^{\#\{(i, j) \in I \times J \,:\, i > j\}}$ is
the sign of the permutation sorting the concatenation $(I \frown J)$ and
$\hat{e}_I = \sum_{k \in I} \hat{e}_k$.
In the discrete :ref:`interlaced <interlaced>` geometry, a p-form and its
Hodge dual are centered at different lattice positions. The shift aligns
them so that the inner-product identity holds after summing over the lattice:
.. math::
\sum_{x, I} a_I[x] \, b_I[x] = \sum_x (a \wedge \star b)_{(0, \ldots, D-1)}[x]
For $p = 0$ and $p = D$ the shift is trivial ($\hat{e}_I = 0$ or the
shifts cancel in the wedge).
Parameters
----------
f : Form
A p-form on a :class:`Lattice`.
Returns
-------
Form
The (D-p)-form $\star f$.
"""
lat = f.lattice
p = f.degree
D = lat.D
# ★ applies only ±1 signs, so it preserves the input dtype.
result = lat.zeros(D - p, dtype=f.dtype)
for J_comp in lat.components[D - p]:
J_set = set(J_comp)
I_comp = tuple(k for k in range(D) if k not in J_set)
sign = _perm_sign(I_comp + J_comp)
# Roll f[I_comp] by +1 in each direction k ∈ I_comp to
# evaluate f_I at site n − ê_I.
spatial = f[lat.comp_index[p][I_comp]]
for k in I_comp:
spatial = np.roll(spatial, +1, axis=k)
result[lat.comp_index[D - p][J_comp]] = sign * spatial
return result
# ---------------------------------------------------------------------------
# Wedge product ∧ : Ω^n × Ω^m → Ω^{n+m}
# ---------------------------------------------------------------------------
[docs]def wedge(a, b):
r"""
Wedge product of an n-form a and an m-form b, an (n+m)-form.
For each output component $O = (o_0, \ldots, o_{n+m-1})$, the sum over all
shuffles $O = A \sqcup B$ ($A$ the $n$ $a$-directions, $B$ the $m$
$b$-directions):
.. math::
(a \wedge b)_O[x] = \sum_{O = A \sqcup B} \sigma(A \frown B) \; a_A[x] \; b_B[x + \hat{e}_A]
where $\hat{e}_A = \sum_{k \in A} \hat{e}_k$ and
$\sigma(A \frown B) = (-1)^{\#\{(k, j) \in A \times B \,:\, j < k\}}$ is the
sign of the permutation sorting $(A \frown B)$ back into $O$.
Parameters
----------
a : Form
An n-form on a :class:`Lattice`.
b : Form
An m-form on the same :class:`Lattice`.
Returns
-------
Form
The (n+m)-form $a \wedge b$.
Raises
------
ValueError
If $n + m > D$.
"""
lat = a.lattice
n, m = a.degree, b.degree
if n + m > lat.D:
raise ValueError(
f"Cannot wedge a {n}-form and a {m}-form in D={lat.D}: {n}+{m} > {lat.D}"
)
# The wedge product is bilinear, so the result dtype is the promotion of the inputs'.
result = lat.zeros(n + m, dtype=np.result_type(a.dtype, b.dtype))
for out_comp in lat.components[n + m]:
out_idx = lat.comp_index[n + m][out_comp]
# Enumerate all n-element (a-direction) subsets of out_comp as A_dirs;
# the complementary m-element (b-direction) subset becomes B_dirs.
for A_dirs in combinations(out_comp, n):
B_dirs = tuple(k for k in out_comp if k not in A_dirs)
# Sign σ(A⌢B): inversions in (A, B) relative to sorted O,
# i.e. pairs (a, b) with a ∈ A, b ∈ B, b < a.
inversions = sum(
1
for k in A_dirs
for j in B_dirs
if j < k
)
sign = (-1) ** inversions
a_idx = lat.comp_index[n][A_dirs]
a_spatial = a[a_idx] # shape (N,...,N)
# b is evaluated at x + hat_e_A: roll by -1 in each direction k ∈ A_dirs.
b_idx = lat.comp_index[m][B_dirs]
b_spatial = b[b_idx] # shape (N,...,N)
for k in A_dirs:
b_spatial = np.roll(b_spatial, -1, axis=k)
result[out_idx] += sign * a_spatial * b_spatial
return result
# ---------------------------------------------------------------------------
# Tests
# ---------------------------------------------------------------------------
if __name__ == '__main__':
import argparse
from itertools import product as iproduct
parser = argparse.ArgumentParser(description="Test compact differential forms.")
parser.add_argument('--D', type=int, default=4, help='Spacetime dimension')
parser.add_argument('--N', type=int, default=6, help='Lattice size per direction')
args = parser.parse_args()
D = args.D
N = args.N
lat = Lattice(D, N)
print(f"Testing on {lat}\n")
# --- Sparsity report ---------------------------------------------------
# How many numbers do we actually store per lattice site?
print("STORAGE PER SITE")
for p in range(D + 1):
n_comps = comb(D, p)
old_fraction = n_comps / 2**D
print(f" {p}-form: {n_comps} component(s) per site "
f"(was {old_fraction:.4f} × 2^D = {n_comps} nonzero in old layout)")
# --- Nilpotency of d ---------------------------------------------------
print("\nNILPOTENCY OF d")
for p in range(D + 1):
a = lat.random(p)
da = d(a)
if p == D:
assert da == 0, f"d({p}-form) should be scalar 0"
print(f' ✅ d({p}-form) = 0')
else:
assert not np.isclose(da, 0).all(), f"d({p}-form) is unexpectedly zero"
dda = d(da)
# np.asarray handles both the scalar-0 case (p == D-1) and the
# Form case (p < D-1) uniformly.
assert np.isclose(np.asarray(dda), 0).all(), f"d² ≠ 0 on {p}-form"
print(f' ✅ d² = 0 on {p}-forms')
# --- Nilpotency of δ ---------------------------------------------------
print("\nNILPOTENCY OF δ")
for p in range(D + 1):
a = lat.random(p)
da = delta(a)
if p == 0:
assert da == 0, f"δ({p}-form) should be scalar 0"
print(f' ✅ δ(0-form) = 0')
else:
assert not np.isclose(da, 0).all(), f"δ({p}-form) is unexpectedly zero"
dda = delta(da)
assert np.isclose(np.asarray(dda), 0).all(), f"δ² ≠ 0 on {p}-form"
print(f' ✅ δ² = 0 on {p}-forms')
# --- Adjointness Σ (da)·b = -Σ a·(δb) --------------------------------
# Inner product of two p-forms: sum over all components and all sites.
# In compact layout every array element is a genuine DOF, so np.sum()
# is the correct lattice L² inner product.
print("\nADJOINTNESS ⟨da, b⟩ = ⟨a, δb⟩")
for p in range(D):
a = lat.random(p)
b = lat.random(p + 1)
LHS = +( d(a) * b ).sum()
RHS = +( a * delta(b) ).sum()
assert np.isclose(LHS, RHS), f"Adjointness failed for p={p}: {LHS} vs {RHS}"
print(f" ✅ ⟨da, b⟩ = ⟨a, δb⟩ for a ∈ Ω^{p}, b ∈ Ω^{p+1}")
# --- Wedge: bilinearity ------------------------------------------------
print("\nBILINEARITY")
for n, m in iproduct(range(D + 1), repeat=2):
if n + m > D:
continue
a = lat.random(n); b = lat.random(n); c = lat.random(m); e = lat.random(m)
assert np.isclose(wedge(a + b, c), wedge(a, c) + wedge(b, c)).all()
assert np.isclose(wedge(a, c + e), wedge(a, c) + wedge(a, e)).all()
print(f" ✅ ∧ bilinear for {n} ∧ {m}")
# --- Wedge: Leibniz rule d(a∧b) = da∧b + (-1)^n a∧db ----------------
print("\nLEIBNIZ RULE")
for n, m in iproduct(range(D + 1), repeat=2):
if n + m + 1 > D:
continue
a = lat.random(n); b = lat.random(m)
LHS = d(wedge(a, b))
RHS = wedge(d(a), b) + (-1)**n * wedge(a, d(b))
assert not np.isclose(LHS, 0).all()
assert np.isclose(LHS, RHS).all(), f"Leibniz failed for {n}∧{m}"
print(f" ✅ d(a∧b) = da∧b + (-1)^{n} a∧db for a ∈ Ω^{n}, b ∈ Ω^{m}")
# --- Wedge: associativity (a∧b)∧c = a∧(b∧c) -------------------------
print("\nASSOCIATIVITY")
for n, m, p in iproduct(range(D), repeat=3):
if n + m + p > D:
continue
a = lat.random(n); b = lat.random(m); c = lat.random(p)
LHS = wedge(wedge(a, b), c)
RHS = wedge(a, wedge(b, c))
assert not np.isclose(LHS, 0).all()
assert np.isclose(LHS, RHS).all(), f"Associativity failed for {n},{m},{p}"
print(f" ✅ (a∧b)∧c = a∧(b∧c) for degrees {n}, {m}, {p}")
# --- Laplacian: direct stencil == dδ + δd ---------------------------
# The performant component-diagonal Laplacian must reproduce the explicit
# composition dδ + δd, including the degenerate ends of the complex where
# δf (p=0) or df (p=D) is the scalar 0.
print("\nLAPLACIAN Δ = dδ + δd")
for p in range(D + 1):
f = lat.random(p)
direct = np.asarray(laplacian(f))
ref = np.zeros_like(direct)
if p > 0:
ref = ref + np.asarray(d(delta(f)))
if p < D:
ref = ref + np.asarray(delta(d(f)))
assert np.allclose(direct, ref), f"Δ ≠ dδ+δd on {p}-forms"
# Self-adjoint, and the Weitzenböck identity ⟨Δf,f⟩ = ⟨df,df⟩+⟨δf,δf⟩ ≥ 0.
b = lat.random(p)
assert np.isclose((laplacian(f) * b).sum(), (f * laplacian(b)).sum()), \
f"Δ not self-adjoint on {p}-forms"
df2 = (np.asarray(d(f)) ** 2).sum()
δf2 = (np.asarray(delta(f)) ** 2).sum()
assert np.isclose((laplacian(f) * f).sum(), df2 + δf2), \
f"⟨Δf,f⟩ ≠ ⟨df,df⟩+⟨δf,δf⟩ on {p}-forms"
assert df2 + δf2 >= -1e-9, f"Δ not positive on {p}-forms"
print(f" ✅ Δ = dδ+δd, self-adjoint, ⟨Δf,f⟩ = ⟨df,df⟩+⟨δf,δf⟩ ≥ 0 on {p}-forms")
# --- Hodge star: inner product identity --------------------------------
# Σ_{n,I} a_I[n] b_I[n] = Σ_n (a ∧ ★b)_{0,…,D-1}[n]
print("\nHODGE STAR Σ a·b = Σ (a∧★b)")
for p in range(D + 1):
a = lat.random(p)
b = lat.random(p)
LHS = (a * b).sum()
RHS = wedge(a, star(b)).sum()
assert np.isclose(LHS, RHS), \
f"Hodge identity failed for p={p}: {LHS} vs {RHS}"
print(f" ✅ Σ a·b = Σ (a∧★b) for p={p}")
# --- Discrete δ = (-1)^{n(k+1)+1} pull(★d★, ê_all) ---------------------
# The continuum identity δ = (-1)^{n(k+1)+1} ★d★ acquires a spatial shift
# on the compact lattice because the staggered Hodge star accumulates one
# step back per application. The exact discrete relation is:
# ★d★ f = (-1)^{n(k+1)+1} · push(δf, (1,...,1))
# equivalently:
# δf = (-1)^{n(k+1)+1} · pull(★d★ f, (1,...,1))
print("\nCODIFFERENTIAL ★d★ = (-1)^{n(k+1)+1} push(δ, ê_all)")
for k in range(1, D + 1):
f = lat.random(k)
sign = (-1) ** (D * (k + 1) + 1)
lhs = np.asarray(star(d(star(f))))
rhs = sign * np.asarray(push(delta(f), (1,) * D))
assert np.isclose(lhs, rhs).all(), \
f"★d★ ≠ (-1)^{{n(k+1)+1}} push(δ, ê_all) on {k}-forms (sign={sign:+d})"
print(f" ✅ ★d★ = (-1)^{{n(k+1)+1}} push(δ, ê_all) on {k}-forms")
# --- to_interlaced round-trip ------------------------------------------
# Verify that embedding a compact form into the (2N)^D array places
# nonzero values only at sites with exactly p odd coordinates.
print("\nTO_INTERLACED")
TXYZ = np.mgrid[(slice(0, 2*N),) * D]
odds = np.mod(TXYZ, 2).sum(axis=0) # number of odd coords per site
for p in range(D + 1):
a = lat.random(p)
big = a.to_interlaced() # shape (2N,...,2N)
# Nonzero only where exactly p coordinates are odd.
wrong_sites = big[odds != p]
assert np.all(wrong_sites == 0), f"to_interlaced put values at non-{p}-form sites"
# All p-form sites are filled (random data is generically nonzero).
right_sites = big[odds == p]
assert not np.all(right_sites == 0), f"to_interlaced left {p}-form sites empty"
print(f" ✅ {p}-form embeds correctly into (2N)^D array")
# --- Round-trip tests for from_interlaced / to_interlaced --------------
print("\nROUND-TRIP: compact → interlaced → compact")
for p in range(D + 1):
f = lat.random(p)
g = Form.from_interlaced(p, f.to_interlaced())
assert (f == g).all(), f"compact→interlaced→compact failed for p={p}"
print(f" ✅ from_interlaced(to_interlaced(f)) == f for p={p}")
print("\nROUND-TRIP: interlaced → compact → interlaced")
from supervillain.lattice import interlaced as il
il_lat = il.Lattice(D, N)
for p in range(D + 1):
data = il_lat.random(p)
assert (Form.from_interlaced(p, data).to_interlaced() == data).all(), \
f"interlaced→compact→interlaced failed for p={p}"
print(f" ✅ to_interlaced(from_interlaced(data)) == data for p={p}")
# --- Cross-validation against interlaced.py ----------------------------
print("\nCROSS-VALIDATION: compact d vs interlaced d")
for p in range(D):
a = lat.random(p)
path_A = d(a).to_interlaced()
path_B = il.d(a.to_interlaced())
assert (path_A == path_B).all(), f"compact d ≠ interlaced d on {p}-forms"
print(f" ✅ compact d = interlaced d on {p}-forms")
print("\nCROSS-VALIDATION: compact δ vs interlaced δ")
for p in range(1, D + 1):
a = lat.random(p)
path_A = delta(a).to_interlaced()
path_B = il.delta(a.to_interlaced())
assert (path_A == path_B).all(), f"compact δ ≠ interlaced δ on {p}-forms"
print(f" ✅ compact δ = interlaced δ on {p}-forms")
print("\nCROSS-VALIDATION: compact wedge vs interlaced wedge")
for p, q in iproduct(range(D + 1), repeat=2):
if p + q > D:
continue
a = lat.random(p)
b = lat.random(q)
path_A = wedge(a, b).to_interlaced()
path_B = il.wedge(p, q, a.to_interlaced(), b.to_interlaced())
# The two paths can take different floating-point arithmetic orderings,
# so compare with np.isclose. Other cross-validations should be bit-exact matches.
assert np.isclose(path_A, path_B).all(), f"compact ∧ ≠ interlaced ∧ for {p}∧{q}"
print(f" ✅ compact wedge = interlaced wedge for {p}∧{q}")
print("\nCROSS-VALIDATION: compact ★ vs interlaced ★")
for p in range(D + 1):
a = lat.random(p)
path_A = star(a).to_interlaced()
path_B = il.star(a.to_interlaced())
assert (path_A == path_B).all(), f"compact ★ ≠ interlaced ★ on {p}-forms"
print(f" ✅ compact ★ = interlaced ★ on {p}-forms")
print("\nAll tests passed.")