add invertibility test for local_stout

lib/gpt/qcd/gauge/smear/local_stout.py

@@ -182,6 +182,29 @@ def __call__(self, fields):
         )
         return U_prime
 
+    def inv(self, fields, max_iter=100):
+        C_mu, U, fm = self.get_C(fields)
+        mu = self.params["dimension"]
+        U_prime = g.copy(U)
+        for it in range(max_iter):
+            U_prime_mu_last = g.copy(U_prime[mu])
+            U_prime[mu] @= g(
+                g.matrix.exp(
+                    -g.qcd.gauge.project.traceless_anti_hermitian(C_mu * g.adj(U_prime[mu]))
+                )
+                * U[mu]
+            )
+            eps2 = g.norm2(U_prime_mu_last - U_prime[mu]) / U_prime_mu_last.grid.gsites
+            if eps2 < U_prime_mu_last.grid.precision.eps**2:
+                break
+        if it == max_iter - 1:
+            # indicate failure
+            g.message(
+                f"Warning: local_stout could not be inverted; last eps^2 = {eps2} after {max_iter} iterations"
+            )
+            return None
+        return U_prime
+
     def jacobian(self, fields, fields_prime, src):
         nd = fields[0].grid.nd
         U_prime = fields_prime[0:nd]

tests/qcd/gauge.py

@@ -246,47 +246,18 @@
     assert eps2 < 1e-25
 
 # test reversibility of local_stout
-# ran once, confirmed |rho| < 1/8 bound, disable for now
-# class inv_functional(g.core.group.differentiable_functional):
-#     def __init__(self, U0, fnc):
-#         self.U0 = U0
-#         self.fnc = fnc
-
-#     def __call__(self, U1):
-#         U0test = self.fnc(U1)
-#         vol = U1[0].grid.gsites*3*4
-#         return vol - sum([g.inner_product(self.U0[mu], U0test[mu]).real for mu in range(len(self.U0))])
-
-#     def gradient(self, U1, dU1):
-#         assert U1 == dU1
-
-#         U0test = self.fnc(U1)
-
-#         ret = [g(-1.0 / 2j * g.qcd.gauge.project.traceless_anti_hermitian(self.U0[mu] * g.adj(U0test[mu]))) for mu in range(len(self.U0))]
-#         for r in ret:
-#             r.otype = U1[0].otype.cartesian()
-
-#         dsrc = self.fnc.jacobian(U1, U0test, ret)
-
-#         return dsrc
-
-# for rho in [0.25, 0.05, 0.124]:
-#     lsm = g.qcd.gauge.smear.local_stout(rho=rho, dimension=1, checkerboard=g.even)
-#     for i in range(10):
-#         rng.element(U, scale=100)
-#     for u in U:
-#         g.project(u, "defect")
-#     U0 = g.copy(U)
-#     g.message(rho, g.qcd.gauge.plaquette(U))
-#     U1 = lsm(U)
-#     # now find reverse and check
-
-#     invf = inv_functional(U1, lsm)
-#     opt = g.algorithms.optimize
-#     ls0 = g.algorithms.optimize.line_search_none
-#     gd = g.algorithms.optimize.gradient_descent(maxiter=2000, eps=1e-10, step=0.15, line_search=ls0)
-#     rng.element(U, scale=10)
-#     gd(invf)(U, U)
-
-#     for mu in range(4):
-#         g.message(mu, g.norm2(U[mu] - U0[mu]), g.norm2(U[mu]))
+for rho in [0.05, 0.1, 0.124, 0.25]:
+    for mu in range(4):
+        for cb in [g.even, g.odd]:
+            g.message(f"Testing reversibility for rho={rho}, mu={mu}, cb={cb.__name__}")
+            lsm = g.qcd.gauge.smear.local_stout(rho=rho, dimension=mu, checkerboard=cb)
+            Uprime = lsm(U)
+            U0 = lsm.inv(Uprime)
+            if U0 is None:
+                assert rho > 1 / 8
+            else:
+                eps2 = 0.0
+                for nu in range(4):
+                    eps2 += g.norm2(U[nu] - U0[nu]) / g.norm2(U0[nu])
+                g.message(eps2)
+                assert eps2 < 1e-28