improve

lib/gpt/qcd/gauge/stencil/algebra_laplace.py

@@ -20,15 +20,6 @@
 import numpy as np
 
 
-class fixed_gauge_algebra_laplace:
-    def __init__(self, parent, U):
-        self.parent = parent
-        self.U = U
-
-    def __call__(self, dst, src):
-        self.parent(self.U + dst, self.U + src)
-
-
 class algebra_laplace:
     def __init__(self, U):
         self.U = U
@@ -167,5 +158,22 @@ def projected_gradient(self, left, U, right):
             ret[mu] @= g.qcd.gauge.project.traceless_hermitian(ret[mu])
         return ret
 
-    def fixed_gauge(self, U):
-        return fixed_gauge_algebra_laplace(self, U)
+    def inverse(self, inverter):
+        def _mat(dst, src):
+            U = src[0 : self.nd]
+
+            def _inv_mat(_dst, _src):
+                self(U + _dst, U + _src)
+
+            inv_mat = g.matrix_operator(mat=_inv_mat, accept_list=True, accept_guess=(True, False))
+
+            im = inverter(inv_mat)
+
+            g.eval(dst[self.nd :], im * src[self.nd :])
+
+            for d in dst[self.nd :]:
+                d.otype = src[self.nd].otype
+
+            g.copy(dst[0 : self.nd], U)
+
+        return g.matrix_operator(mat=_mat, accept_list=True)

lib/gpt/qcd/scalar/action/mass_term.py

@@ -162,7 +162,9 @@ def draw(self, fields, rng):
 
         value = g.group.inner_product(pi, pi) * 0.5
 
-        pi_prime = self.inv_sqrt_M(pi)
+        fields_prime = self.inv_sqrt_M(fields)
+        pi_prime = fields_prime[n:]
+
         for mu in range(n):
             pi_prime[mu].otype = pi[mu].otype
             pi[mu] @= g.project(pi_prime[mu], "defect")
@@ -180,13 +182,18 @@ def gradient(self, fields, dfields):
 
         dS = []
 
-        pi_prime = self.M(fields)[n:]
         dU = [f for f in dfields if f in U]
         if len(dU) > 0:
             grad_U = self.M_projected_gradient(U, pi)
         else:
             grad_U = None
 
+        dpi = [f for f in dfields if f in pi]
+        if len(dpi) > 0:
+            pi_prime = self.M(fields)[n:]
+        else:
+            pi_prime = None
+
         for _field in dfields:
             mu = fields.index(_field)
 

tests/qcd/gauge.py

@@ -244,3 +244,49 @@
     eps2 = g.norm2(Up[3][0, 0, 0, t] - ref)
     g.message(f"Test temporal gauge at t={t}: {eps2}")
     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]))

tests/qcd/scalar.py

@@ -185,7 +185,7 @@ def __slap(dst, src):
 
 cg = g.algorithms.inverter.block_cg({"eps": 1e-12, "maxiter": 100})
 slap = g.matrix_operator(
-    mat=lap, inv_mat=cg(lap.fixed_gauge(U)), accept_list=True, accept_guess=(False, True)
+    mat=lap, inv_mat=lap.inverse(cg), accept_list=True, accept_guess=(False, True)
 )
 slap2 = slap * slap