1"""Small complex cross-ratio lattice utilities."""
 2import numpy as np
 3
 4
 5def complete_d(a, b, c, eps=1e-12):
 6    """Return d such that (a,b,c,d) has cross-ratio -1.
 7
 8    The returned singular mask identifies near-zero denominators; singular
 9    values are represented by complex infinity rather than NaN.
10    """
11    den = (b - c) - (a - b)
12    num = (b - c) * a - (a - b) * c
13    singular = np.abs(den) < eps
14    d = np.empty(np.broadcast(a, b, c).shape, dtype=np.result_type(a, b, c))
15    np.divide(num, den, out=d, where=~singular)
16    d = np.where(singular, np.inf + 0j, d)
17    return d, singular
18
19
20def complete_c(a, b, d, eps=1e-12):
21    """Inverse completion: solve the same equation for c."""
22    # From d=((b-c)a-(a-b)c)/((b-c)-(a-b)), collect c.
23    # c = ((a-b)*d - (b)*a + b*a? use direct Möbius inversion below)
24    # Cross-ratio equation gives (a-b)(c-d)+(b-c)(d-a)=0.
25    # c*((a-b)-(d-a)) = (a-b)d - b*(d-a)
26    den = (a - b) - (d - a)
27    num = (a - b) * d - b * (d - a)
28    singular = np.abs(den) < eps
29    c = np.empty(np.broadcast(a, b, d).shape, dtype=np.result_type(a, b, d))
30    np.divide(num, den, out=c, where=~singular)
31    c = np.where(singular, np.inf + 0j, c)
32    return c, singular
33
34
35def cross_ratio(a, b, c, d):
36    return (a-b)*(c-d)/((b-c)*(d-a))
37
38
39def residual(a, b, c, d, eps=1e-14):
40    den = (b-c)*(d-a)
41    out = np.full(np.broadcast(a,b,c,d).shape, np.inf, dtype=float)
42    good = np.abs(den) > eps
43    out[good] = np.abs(cross_ratio(a,b,c,d)[good] + 1)
44    return out
45
46
47def generate_lattice(n, seed=0, dtype=np.complex128):
48    """Generate a lattice from random first row/column and exact completion."""
49    rng = np.random.default_rng(seed)
50    z = np.zeros((n,n), dtype=dtype)
51    z[0, :] = (rng.uniform(-1,1,n) + 1j*rng.uniform(-1,1,n)).astype(dtype)
52    z[:, 0] = (rng.uniform(-1,1,n) + 1j*rng.uniform(-1,1,n)).astype(dtype)
53    # Avoid an exactly shared random origin while retaining a valid boundary.
54    z[0,0] = 0.2 + 0.15j
55    singular = 0
56    for i in range(n-1):
57        for j in range(n-1):
58            c, bad = complete_c(z[i,j], z[i+1,j], z[i,j+1])
59            if bad: singular += 1
60            z[i+1,j+1] = c
61    return z, singular
62
63
64def plaquette_residuals(z):
65    return residual(z[:-1,:-1], z[1:,:-1], z[1:,1:], z[:-1,1:])
66
67
68def constrained_from_boundary(top, left):
69    n = len(top)
70    z = np.zeros((n,n), dtype=np.result_type(top, left))
71    z[0,:], z[:,0] = top, left
72    for i in range(n-1):
73        for j in range(n-1):
74            z[i+1,j+1], _ = complete_c(z[i,j], z[i+1,j], z[i,j+1])
75    return z
76
77
78def additive_from_boundary(top, left):
79    n = len(top)
80    z = np.zeros((n,n), dtype=np.result_type(top, left))
81    z[0,:], z[:,0] = top, left
82    for i in range(n-1):
83        for j in range(n-1):
84            z[i+1,j+1] = z[i+1,j] + z[i,j+1] - z[i,j]
85    return z