From 98e6f3cb7a62c491bfabce3a81f2f2f272932b64 Mon Sep 17 00:00:00 2001 From: Christoph Lehner Date: Mon, 8 Jul 2024 20:13:34 +0200 Subject: [PATCH] switch to automatic differentiation framework for dQ and dE --- lib/gpt/ad/reverse/node.py | 22 +- lib/gpt/core/foundation/tensor.py | 5 + lib/gpt/qcd/gauge/__init__.py | 9 +- lib/gpt/qcd/gauge/loops.py | 30 ++ lib/gpt/qcd/gauge/stencil/__init__.py | 2 +- lib/gpt/qcd/gauge/stencil/topology.py | 560 -------------------------- tests/ad/ad.py | 2 +- tests/qcd/gauge.py | 18 +- 8 files changed, 67 insertions(+), 581 deletions(-) diff --git a/lib/gpt/ad/reverse/node.py b/lib/gpt/ad/reverse/node.py index 6371b0f3..a0b0ac24 100644 --- a/lib/gpt/ad/reverse/node.py +++ b/lib/gpt/ad/reverse/node.py @@ -69,27 +69,31 @@ def gradient(self, fields, dfields): for i in indices: self.arguments[i].gradient = None self.arguments[i].with_gradient = True + for i in range(len(fields)): + assert fields[i] not in [a.value for a in self.arguments] + self.arguments[i].value @= fields[i] self.node() return [self.arguments[i].gradient for i in indices] def str_traverse(node, indent=0): if not callable(node._forward): - return (" "*indent) + "leaf(" + str(node._container) + ")" + return (" " * indent) + "leaf(" + str(node._container) + ")" else: - pre = " "*indent + pre = " " * indent if node._tag is not None: tag = node._tag else: tag = str(node._forward) ret = pre + "(" + tag + "):" for x in node._children: - ret = ret + "\n" + str_traverse(x, indent+1) + ret = ret + "\n" + str_traverse(x, indent + 1) return ret - + # gctr = 0 + class node_base(base): foundation = foundation @@ -102,7 +106,7 @@ def __init__( with_gradient=True, infinitesimal_to_cartesian=True, _container=None, - _tag=None + _tag=None, ): # global gctr # gctr+=1 @@ -189,8 +193,8 @@ def setter(y, z): return x.project(getter, setter) def project(x, getter, setter): - assert False # for future use - + assert False # for future use + def _forward(): return getter(x.value) @@ -282,9 +286,9 @@ def backward(self, nodes, first_gradient, initial_gradient): raise Exception( "Expression evaluates to a field. Gradient calculation is not unique." ) - #if isinstance(self._container[0], complex) and abs(self.value.imag) > 1e-12 * abs( + # if isinstance(self._container[0], complex) and abs(self.value.imag) > 1e-12 * abs( # self.value.real - #): + # ): # raise Exception( # f"Expression does not evaluate to a real number ({self.value}). Gradient calculation is not unique." # ) diff --git a/lib/gpt/core/foundation/tensor.py b/lib/gpt/core/foundation/tensor.py index 6a395568..42bd4519 100644 --- a/lib/gpt/core/foundation/tensor.py +++ b/lib/gpt/core/foundation/tensor.py @@ -62,3 +62,8 @@ def component_multiply(a, b): res = a.new() res.array = numpy.multiply(a.array, b.array) return res + + +def copy(a, b): + for i in range(len(a)): + a[i].array[:] = b[i].array[:] diff --git a/lib/gpt/qcd/gauge/__init__.py b/lib/gpt/qcd/gauge/__init__.py index f39a9191..12e52c65 100644 --- a/lib/gpt/qcd/gauge/__init__.py +++ b/lib/gpt/qcd/gauge/__init__.py @@ -17,7 +17,12 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # from gpt.qcd.gauge.create import random, unit -from gpt.qcd.gauge.loops import rectangle, field_strength +from gpt.qcd.gauge.loops import ( + rectangle, + field_strength, + differentiable_topology, + differentiable_energy_density, +) from gpt.qcd.gauge.topology import topological_charge_5LI from gpt.qcd.gauge.staples import staple from gpt.qcd.gauge.transformation import transformed @@ -26,8 +31,6 @@ staple_sum, energy_density, topological_charge, - differentiable_topology, - projected_topology_gradient, algebra_laplace, ) import gpt.qcd.gauge.project diff --git a/lib/gpt/qcd/gauge/loops.py b/lib/gpt/qcd/gauge/loops.py index 4da3520b..50f0f1d4 100644 --- a/lib/gpt/qcd/gauge/loops.py +++ b/lib/gpt/qcd/gauge/loops.py @@ -170,3 +170,33 @@ def field_strength(U, mu, nu): F = g.eval(U[mu] * v + g.cshift(v * U[mu], mu, -1)) F = 0.125 * (F - g.adj(F)) return F + + +def differentiable_topology(aU): + Bx = field_strength(aU, 1, 2) + By = field_strength(aU, 2, 0) + Bz = field_strength(aU, 0, 1) + + Ex = field_strength(aU, 3, 0) + Ey = field_strength(aU, 3, 1) + Ez = field_strength(aU, 3, 2) + + coeff = 8.0 / (32.0 * np.pi**2) + + Q = g.sum(g.trace(Bx * Ex) + g.trace(By * Ey) + g.trace(Bz * Ez)) * coeff + + return Q + + +def differentiable_energy_density(aU): + Nd = len(aU) + grid = aU[0].grid + res = None + for mu in range(Nd): + for nu in range(mu): + Fmunu = field_strength(aU, mu, nu) + if res is None: + res = Fmunu * Fmunu + else: + res += Fmunu * Fmunu + return (-1.0 / grid.gsites) * g.sum(g.trace(res)) diff --git a/lib/gpt/qcd/gauge/stencil/__init__.py b/lib/gpt/qcd/gauge/stencil/__init__.py index 96cb8327..ebabbd20 100644 --- a/lib/gpt/qcd/gauge/stencil/__init__.py +++ b/lib/gpt/qcd/gauge/stencil/__init__.py @@ -19,5 +19,5 @@ from gpt.qcd.gauge.stencil.plaquette import plaquette from gpt.qcd.gauge.stencil.staple import staple_sum from gpt.qcd.gauge.stencil.energy_density import energy_density -from gpt.qcd.gauge.stencil.topology import topological_charge, differentiable_topology, projected_topology_gradient +from gpt.qcd.gauge.stencil.topology import topological_charge from gpt.qcd.gauge.stencil.algebra_laplace import algebra_laplace diff --git a/lib/gpt/qcd/gauge/stencil/topology.py b/lib/gpt/qcd/gauge/stencil/topology.py index 35e76e2b..29435155 100644 --- a/lib/gpt/qcd/gauge/stencil/topology.py +++ b/lib/gpt/qcd/gauge/stencil/topology.py @@ -18,7 +18,6 @@ # import gpt as g import numpy as np -from gpt.core.group import differentiable_functional default_topological_charge_cache = {} @@ -87,562 +86,3 @@ def topological_charge( T = g.sum(T).real / T.grid.gsites return T - - -def v_projected_gradient(U, mu, nu, right_left): - Nd = len(U) - assert Nd == 4 - - grad = [g.group.cartesian(U[0]) for i in range(Nd)] - for gr in grad: - gr[:] = 0 - - grad[mu] += g.cshift(g.adj(U[nu]) * right_left * g.cshift(U[nu], mu, 1), nu, -1) * g.adj( - U[mu] - ) - grad[mu] -= U[nu] * g.cshift(right_left, nu, 1) * g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]) - - grad[nu] -= U[nu] * g.cshift( - g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]) * right_left, mu, -1 - ) - grad[nu] -= g.cshift(g.adj(U[mu]) * U[nu] * g.cshift(right_left, nu, 1), mu, -1) * g.adj( - U[nu] - ) - - grad[nu] += right_left * g.cshift(U[nu], mu, 1) * g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]) - grad[nu] += U[nu] * g.cshift(right_left, nu, 1) * g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]) - - grad[mu] @= g.qcd.gauge.project.traceless_hermitian(grad[mu] / 2j) - grad[nu] @= g.qcd.gauge.project.traceless_hermitian(grad[nu] / 2j) - - return grad - - -def v(U, mu, nu): - return g.eval( - g.cshift(U[nu], mu, 1) * g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]) - - g.cshift(g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]) * U[nu], nu, -1) - ) - - -def F_projected_gradient(U, mu, nu, right_left): - Nd = len(U) - assert Nd == 4 - - right_left_s = g.cshift(right_left, mu, 1) - v_val = v(U, mu, nu) - grad = v_projected_gradient(U, mu, nu, g(right_left * U[mu])) - grad2 = v_projected_gradient(U, mu, nu, g(U[mu] * right_left_s)) - for rho in range(Nd): - grad[rho] += grad2[rho] - - grad[mu] -= g.qcd.gauge.project.traceless_hermitian(U[mu] * v_val * right_left / 2j) - grad[mu] -= g.qcd.gauge.project.traceless_hermitian(U[mu] * right_left_s * v_val / 2j) - - # left * U[mu] * v * right + left * g.cshift(v * U[mu], mu, -1) * right - return grad - - -def field_strength_projected_gradient(U, mu, nu, right_left): - Nd = len(U) - assert Nd == 4 - - fg1 = F_projected_gradient(U, mu, nu, g(right_left)) - fg2 = F_projected_gradient(U, mu, nu, g(g.adj(right_left))) - grad = [] - for mu in range(Nd): - grad.append(g(0.125 * (fg1[mu] - g.adj(fg2[mu])))) - return grad - - -def topological_charge_gradient(U, mask=None): - field_strength = g.qcd.gauge.field_strength - - Bx = field_strength(U, 1, 2) - By = field_strength(U, 2, 0) - Bz = field_strength(U, 0, 1) - - Ex = field_strength(U, 3, 0) - Ey = field_strength(U, 3, 1) - Ez = field_strength(U, 3, 2) - - if mask is None: - mask = g.complex(U[0].grid) - mask[:] = 1.0 - - delta = g( - g.expr(field_strength_projected_gradient(U, 1, 2, g(mask * Ex))) - + g.expr(field_strength_projected_gradient(U, 3, 0, g(mask * Bx))) - + g.expr(field_strength_projected_gradient(U, 2, 0, g(mask * Ey))) - + g.expr(field_strength_projected_gradient(U, 3, 1, g(mask * By))) - + g.expr(field_strength_projected_gradient(U, 0, 1, g(mask * Ez))) - + g.expr(field_strength_projected_gradient(U, 3, 2, g(mask * Bz))) - ) - - coeff = 8.0 / (32.0 * np.pi**2) - - ret = [] - for d in delta: - dQ = g(coeff * d) - dQ.otype = U[0].otype.cartesian() - ret.append(dQ) - - return ret - - -class differentiable_topology(differentiable_functional): - def __init__(self, mask=None): - self.mask = mask - - def __call__(self, U): - return topological_charge(U, mask=self.mask) - - def gradient(self, U, dU): - gradient = topological_charge_gradient(U, mask=self.mask) - ret = [] - for u in dU: - ret.append(gradient[U.index(u)]) - return ret - - -def projected_gradient_F_projected_gradient(U, mu, nu, right_left, outer_right_left): - Nd = len(U) - assert Nd == 4 - - right_left_s = g.cshift(right_left, mu, 1) - v_val = v(U, mu, nu) - - one = g.identity(g.group.cartesian(U[0])) - grad = projected_gradient_v_projected_gradient(U, mu, nu, right_left, one, outer_right_left) - grad2 = projected_gradient_v_projected_gradient(U, mu, nu, one, right_left_s, outer_right_left) - for rho in range(Nd): - grad[rho] += grad2[rho] - - one = g.identity(g.lattice(grad[0])) - otype = grad[0].otype - - #grad[mu] -= g.qcd.gauge.project.traceless_hermitian(U[mu] * v_val * right_left / 2j) - # v_val_1 - grad[mu] -= projected_gradient_traceless_hermitian( - U[mu] * g.cshift(U[nu], mu, 1) * g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]) * right_left / 2j, - one, - outer_right_left[mu], - otype - ) - - grad[mu] -= projected_gradient_traceless_hermitian( - - U[mu] * g.adj(g.cshift(U[mu] * g.cshift(U[nu], mu, 1), nu, -1)) / 2j, - g.cshift(g.adj(right_left) * U[nu], nu, -1), - g.cshift(outer_right_left[mu], nu, -1), - otype - ) - - grad[nu] -= projected_gradient_traceless_hermitian( - U[nu] * g.cshift(g.cshift(g.adj(U[mu]), nu, 1) * g.adj(U[nu]) * right_left, mu, -1) / 2j, - g.cshift(U[mu], mu, -1), - g.cshift(outer_right_left[mu], mu, -1), - otype - ) - - grad[nu] -= projected_gradient_traceless_hermitian( - - U[nu] * g.cshift(U[mu], nu, 1) * g.cshift(g.adj(U[nu]), mu, 1) * g.adj(U[mu]) / 2j, - g.adj(right_left), - outer_right_left[mu], - otype - ) - - # v_val_2 - grad[mu] += projected_gradient_traceless_hermitian( - U[mu] * g.cshift(g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]) * U[nu], nu, -1) * right_left / 2j, - one, - outer_right_left[mu], - otype - ) - - grad[mu] += projected_gradient_traceless_hermitian( - - U[mu] * g.adj(g.cshift(U[mu], nu, 1) * g.adj(g.cshift(U[nu], mu, 1)))/ 2j, - g.adj(U[nu] * g.cshift(right_left, nu, 1)), - g.cshift(outer_right_left[mu], nu, 1), - otype - ) - - grad[nu] += projected_gradient_traceless_hermitian( - - U[nu] * g.cshift(g.cshift(g.adj(U[mu]), nu, 1), mu, -1) / 2j, - g.adj(g.cshift(g.adj(U[mu]) * U[nu] * g.cshift(right_left, nu, 1), mu, -1)), - g.cshift(g.cshift(outer_right_left[mu], nu, 1), mu, -1), - otype - ) - - grad[nu] += projected_gradient_traceless_hermitian( - U[nu] * g.cshift(right_left, nu, 1) / 2j, - g.cshift(U[mu], nu, 1) * g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]), - g.cshift(outer_right_left[mu], nu, 1), - otype - ) - - #grad[mu] -= g.qcd.gauge.project.traceless_hermitian(g(U[mu] * right_left_s * v_val) / 2j) - grad[mu] -= projected_gradient_traceless_hermitian( - U[mu] * right_left_s * g.cshift(U[nu], mu, 1) * g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]) / 2j, - one, - outer_right_left[mu], - otype - ) - - grad[mu] -= projected_gradient_traceless_hermitian( - - U[mu] * g.adj(g.cshift(U[mu] * right_left_s * g.cshift(U[nu], mu, 1), nu, -1)) / 2j, - g.adj(g.cshift(g.adj(U[nu]), nu, -1)), - g.cshift(outer_right_left[mu], nu, -1), - otype - ) - - grad[nu] -= projected_gradient_traceless_hermitian( - U[nu] * g.cshift(g.cshift(g.adj(U[mu]), nu, 1) * g.adj(U[nu]), mu, -1) / 2j, - g.cshift(U[mu] * right_left_s, mu, -1), - g.cshift(outer_right_left[mu], mu, -1), - otype - ) - - grad[nu] -= projected_gradient_traceless_hermitian( - - U[nu] * g.cshift(U[mu], nu, 1) * g.cshift(g.adj(U[nu]), mu, 1) * g.adj(U[mu] * right_left_s) / 2j, - one, - outer_right_left[mu], - otype - ) - - #v_val_2 - grad[mu] += projected_gradient_traceless_hermitian( - U[mu] * right_left_s * g.cshift(g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]) * U[nu], nu, -1) / 2j, - one, - outer_right_left[mu], - otype - ) - - grad[mu] += projected_gradient_traceless_hermitian( - - U[mu] * g.adj(g.cshift(U[mu] * right_left_s, nu, 1) * g.adj(g.cshift(U[nu], mu, 1)))/ 2j, - g.adj(U[nu]), - g.cshift(outer_right_left[mu], nu, 1), - otype - ) - - grad[nu] += projected_gradient_traceless_hermitian( - - U[nu] * g.cshift(g.cshift(g.adj(U[mu] * right_left_s), nu, 1), mu, -1) / 2j, - g.adj(g.cshift(g.adj(U[mu]) * U[nu], mu, -1)), - g.cshift(g.cshift(outer_right_left[mu], nu, 1), mu, -1), - otype - ) - - grad[nu] += projected_gradient_traceless_hermitian( - U[nu] / 2j, - g.cshift(U[mu] * right_left_s, nu, 1) * g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]), - g.cshift(outer_right_left[mu], nu, 1), - otype - ) - - return grad - - -def projected_gradient_field_strength_projected_gradient(U, mu, nu, right_left, outer_right_left): - Nd = len(U) - assert Nd == 4 - - fg1 = projected_gradient_F_projected_gradient(U, mu, nu, g(right_left), outer_right_left) - fg2 = projected_gradient_F_projected_gradient(U, mu, nu, g(g.adj(right_left)), outer_right_left) - grad = [] - for mu in range(Nd): - grad.append(g(0.125 * (fg1[mu] - g.adj(fg2[mu])))) - - return grad - - -def topological_charge_gradient2(U): - field_strength = g.qcd.gauge.field_strength - - Bx = field_strength(U, 1, 2) - - delta = g( - + g.expr(field_strength_projected_gradient(U, 3, 0, Bx)) - ) - - coeff = 8.0 / (32.0 * np.pi**2) - - ret = [] - for d in delta: - dQ = g(coeff * d) - dQ.otype = U[0].otype.cartesian() - ret.append(dQ) - - return ret - - -def projected_gradient_topological_charge_gradient(U, outer): - field_strength = g.qcd.gauge.field_strength - - Bx = field_strength(U, 1, 2) - By = field_strength(U, 2, 0) - Bz = field_strength(U, 0, 1) - - Ex = field_strength(U, 3, 0) - Ey = field_strength(U, 3, 1) - Ez = field_strength(U, 3, 2) - - delta = g( - #g.expr(field_strength_projected_gradient(U, 1, 2, g(mask * Ex))) - + g.expr(projected_gradient_field_strength_projected_gradient(U, 3, 0, Bx, outer)) - #+ g.expr(field_strength_projected_gradient(U, 3, 0, field_strength_projected_gradient_dinput(U, 3, 0, outer))) - - #+ g.expr(field_strength_projected_gradient(U, 2, 0, g(mask * Ey))) - #+ g.expr(field_strength_projected_gradient(U, 3, 1, g(mask * By))) - #+ g.expr(field_strength_projected_gradient(U, 0, 1, g(mask * Ez))) - #+ g.expr(field_strength_projected_gradient(U, 3, 2, g(mask * Bz))) - ) - - coeff = 8.0 / (32.0 * np.pi**2) - - ret = [] - for d in delta: - dQ = g(coeff * d) - dQ.otype = U[0].otype.cartesian() - ret.append(dQ) - assert False - return ret - - -def projected_gradient_traceless_hermitian(U_L, U_R, R_L, otype): - U_L = g(U_L) - U_R = g(U_R) - U = g(U_L * U_R) - R_L = g(R_L) - Ndim = otype.shape[0] - A = g(g.qcd.gauge.project.traceless_anti_hermitian(g(U_L * R_L * U_R * 0.5 - 0.5 * g.adj(U_R) * R_L * g.adj(U_L))) * 0.5j) - A -= 1.0 / Ndim / 2.0 * g.qcd.gauge.project.traceless_anti_hermitian(U - g.adj(U)) * g.trace(R_L + g.adj(R_L)) * 0.25j - A.otype = otype - return A - - -def projected_gradient_v_projected_gradient(U, mu, nu, right, left, outer_right_left): - Nd = len(U) - assert Nd == 4 - - right_left = g(right * U[mu] * left) - - grad = [g.group.cartesian(U[0]) for i in range(Nd)] - for gr in grad: - gr[:] = 0 - - otype = grad[0].otype - - one = g.identity(g.lattice(grad[0])) - - #grad[mu] += g.cshift(g.adj(U[nu]) * right * left * g.cshift(U[nu], mu, 1), nu, -1) * g.adj( - # U[mu] - #) - - grad[mu] += projected_gradient_traceless_hermitian( - - U[mu] * g.cshift(g.cshift(g.adj(U[nu]), mu, 1) * g.adj(right_left) * U[nu], nu, -1) / 2j, - one, - outer_right_left[mu], - otype - ) - - grad[mu] += projected_gradient_traceless_hermitian( - U[mu] * left * g.cshift(U[nu], mu, 1) * g.cshift(g.adj(U[mu]), nu, 1) / 2j, - g.adj(U[nu]) * right, - g.cshift(outer_right_left[mu], nu, 1), - otype - ) - - grad[nu] += projected_gradient_traceless_hermitian( - - U[nu] / 2j, - g.cshift(U[mu], nu, 1) * g.cshift(g.adj(U[nu]), mu, 1) * g.adj(right_left), - g.cshift(outer_right_left[mu], nu, 1), - otype - ) - - grad[nu] += projected_gradient_traceless_hermitian( - U[nu] * g.cshift(g.cshift(g.adj(U[mu]), nu, 1), mu, -1) / 2j, - g.cshift(g.adj(U[nu]) * right_left, mu, -1), - g.cshift(g.cshift(outer_right_left[mu], nu, 1), mu, -1), - otype - ) - - # grad[mu] -= U[nu] * g.cshift(right * left, nu, 1) * g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]) - - grad[mu] -= projected_gradient_traceless_hermitian( - - U[mu] * g.cshift(U[nu], mu, 1) * g.cshift(g.adj(right_left), nu, 1) * g.adj(U[nu]) / 2j, - one, - outer_right_left[mu], - otype - ) - - grad[mu] -= projected_gradient_traceless_hermitian( - U[mu] * left * g.cshift(g.cshift(g.adj(U[nu]), mu, 1) * g.adj(U[mu]), nu, -1) / 2j, - g.cshift(U[nu], nu, -1) * right, - g.cshift(outer_right_left[mu], nu, -1), - otype - ) - - grad[nu] -= projected_gradient_traceless_hermitian( - U[nu] * g.cshift(right_left, nu, 1) * g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]) / 2j, - one, - outer_right_left[mu], - otype - ) - - grad[nu] -= projected_gradient_traceless_hermitian( - - U[nu] * g.cshift(g.cshift(g.adj(right_left), nu, 1) * g.adj(U[nu]), mu, -1) / 2j, - g.cshift(U[mu], mu, -1), - g.cshift(outer_right_left[mu], mu, -1), - otype - ) - - #grad[nu] -= U[nu] * g.cshift(g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]) * right * left, mu, -1) - grad[mu] -= projected_gradient_traceless_hermitian( - - U[mu] * g.cshift(g.cshift(g.adj(U[nu]), mu, 1), nu, -1) / 2j, - g.cshift(g.adj(right_left) * U[nu], nu, -1), - g.cshift(g.cshift(outer_right_left[nu], mu, 1), nu, -1), - otype - ) - - grad[mu] -= projected_gradient_traceless_hermitian( - U[mu] * left / 2j, - g.cshift(U[nu], mu, 1) * g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]) * right, - g.cshift(outer_right_left[nu], mu, 1), - otype - ) - - grad[nu] -= projected_gradient_traceless_hermitian( - U[nu] * g.cshift(g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]) * right_left, mu, -1) / 2j, - one, - outer_right_left[nu], - otype - ) - - grad[nu] -= projected_gradient_traceless_hermitian( - -U[nu] * g.cshift(U[mu], nu, 1) * g.cshift(g.adj(U[nu]), mu, 1) / 2j, - g.adj(right_left), - g.cshift(outer_right_left[nu], mu, 1), - otype - ) - - #grad[nu] -= g.cshift(g.adj(U[mu]) * U[nu] * g.cshift(right_left, nu, 1), mu, -1) * g.adj(U[nu]) - grad[mu] -= projected_gradient_traceless_hermitian( - - U[mu] / 2j, - g.cshift(U[nu], mu, 1) * g.cshift(g.adj(right_left), nu, 1) * g.adj(U[nu]), - g.cshift(outer_right_left[nu], mu, 1), - otype - ) - - grad[mu] -= projected_gradient_traceless_hermitian( - U[mu] * left * g.cshift(g.cshift(g.adj(U[nu]), mu, 1), nu, -1) / 2j, - g.cshift(g.adj(U[mu]) * U[nu], nu, -1) * right, - g.cshift(g.cshift(outer_right_left[nu], mu, 1), nu, -1), - otype - ) - - grad[nu] -= projected_gradient_traceless_hermitian( - U[nu] * g.cshift(right_left, nu, 1) * g.cshift(g.adj(U[nu]), mu, 1) / 2j, - g.adj(U[mu]), - g.cshift(outer_right_left[nu], mu, 1), - otype - ) - - grad[nu] -= projected_gradient_traceless_hermitian( - -U[nu] * g.cshift(g.cshift(g.adj(right_left), nu, 1) * g.adj(U[nu]) * U[mu], mu, -1) / 2j, - one, - outer_right_left[nu], - otype - ) - - # grad[nu] += right_left * g.cshift(U[nu], mu, 1) * g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]) - grad[mu] += projected_gradient_traceless_hermitian( - - U[mu] * g.cshift(g.cshift(g.adj(U[nu]), mu, 1) * g.adj(right_left), nu, -1) / 2j, - g.cshift(U[nu], nu, -1), - g.cshift(outer_right_left[nu], nu, -1), - otype - ) - - grad[mu] += projected_gradient_traceless_hermitian( - U[mu] * left * g.cshift(U[nu], mu, 1) * g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]) / 2j, - right, - outer_right_left[nu], - otype - ) - - grad[nu] += projected_gradient_traceless_hermitian( - U[nu] * g.cshift(g.adj(g.cshift(U[mu], nu, 1)) * g.adj(U[nu]), mu, -1) / 2j, - g.cshift(right_left, mu, -1), - g.cshift(outer_right_left[nu], mu, -1), - otype - ) - - grad[nu] += projected_gradient_traceless_hermitian( - - U[nu] * g.cshift(U[mu], nu, 1) * g.cshift(g.adj(U[nu]), mu, 1) * g.adj(right_left) / 2j, - one, - outer_right_left[nu], - otype - ) - - # grad[nu] += U[nu] * g.cshift(right_left, nu, 1) * g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]) - grad[mu] += projected_gradient_traceless_hermitian( - - U[mu] * g.cshift(U[nu], mu, 1) * g.cshift(g.adj(right_left), nu, 1) * g.adj(U[nu]) / 2j, - one, - outer_right_left[nu], - otype - ) - - grad[mu] += projected_gradient_traceless_hermitian( - U[mu] * left * g.cshift(g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]), nu, -1) / 2j, - g.cshift(U[nu], nu, -1) * right, - g.cshift(outer_right_left[nu], nu, -1), - otype - ) - - grad[nu] += projected_gradient_traceless_hermitian( - U[nu] * g.cshift(right_left, nu, 1) * g.adj(g.cshift(U[nu], mu, 1)) * g.adj(U[mu]) / 2j, - one, - outer_right_left[nu], - otype - ) - - grad[nu] += projected_gradient_traceless_hermitian( - - U[nu] * g.adj(g.cshift(U[nu] * g.cshift(right_left, nu, 1), mu, -1)) / 2j, - g.cshift(U[mu], mu, -1), - g.cshift(outer_right_left[nu], mu, -1), - otype - ) - - return grad - - -# tr(U) -> tr(Ta U) Ta = 0.5 * (U + U^dag - 1/3 tr(U+U^dag)) -# now: 0.5 * tr(L (U + U^dag - 1/3 tr(U + Udag)) R) -# -> project( U R L - R L U^dag ) - 1/3 Ta tr(Ta (U + Udag)) tr(R L) -# = project( U R L + R L U^dag ) - 1/3 project(U + Udag) tr(R L) - -class projected_topology_gradient(differentiable_functional): - def __init__(self, vector): - self.vector = vector - self.left = g.mcolor(vector[0].grid) ##g.group.cartesian(g.mcolor(vector[0].grid)) - self.right = g.copy(self.left) - g.random("test").cnormal([self.left, self.right]) - - def __call__(self, U): - # return g.inner_product(self.left, g.qcd.gauge.project.traceless_hermitian(U[0])).real - gradient = v_projected_gradient(U, 0, 1, g(self.right * U[0] * self.left)) - #gradient = F_projected_gradient(U, 0, 1, g(self.right * self.left)) - #gradient = field_strength_projected_gradient(U, 0, 1, g(self.right * self.left)) - #gradient = topological_charge_gradient2(U) - return sum([g.inner_product(self.vector[mu], gradient[mu]).real for mu in range(len(gradient))]) - - def gradient(self, U, dU): - #gradient = projected_gradient_traceless_hermitian(U[0], g.adj(self.left)) - #ret = [gradient] - - gradient = projected_gradient_v_projected_gradient(U, 0, 1, self.right, self.left, g(g.adj(self.vector))) - #gradient = projected_gradient_F_projected_gradient(U, 0, 1, g(self.right * self.left), g(g.adj(self.vector))) - #gradient = projected_gradient_field_strength_projected_gradient(U, 0, 1, g(self.right * self.left), g(g.adj(self.vector))) - #gradient = projected_gradient_topological_charge_gradient(U, g(g.adj(self.vector))) - ret = [] - for u in dU: - ret.append(gradient[U.index(u)]) - return ret diff --git a/tests/ad/ad.py b/tests/ad/ad.py index fbe60820..9ea38311 100755 --- a/tests/ad/ad.py +++ b/tests/ad/ad.py @@ -122,7 +122,7 @@ def real(x): # create something to minimize f = real(c).functional(*args) - ff = [vv.value for vv in args] + ff = [g.copy(vv.value) for vv in args] v0 = f(ff) opt = g.algorithms.optimize.adam(maxiter=40, eps=1e-7, alpha=learn_rate) opt(f)(ff, ff) diff --git a/tests/qcd/gauge.py b/tests/qcd/gauge.py index 8072d73e..2bf80a46 100755 --- a/tests/qcd/gauge.py +++ b/tests/qcd/gauge.py @@ -134,21 +134,25 @@ assert eps < 1e-13 g.message("Test diff top") -diff_Q = g.qcd.gauge.differentiable_topology() +adU = [g.ad.reverse.node(g.copy(u)) for u in U] +dQ = g.qcd.gauge.differentiable_topology(adU) +diff_Q = dQ.functional(*adU) diff_Q.assert_gradient_error(rng, U, U, 1e-3, 1e-8) assert abs(Q - diff_Q(U)) < 1e-13 -g.message("Test masked diff top") -msk = g.complex(grid) -rng.normal(msk) -diff_Q = g.qcd.gauge.differentiable_topology(mask=msk) -diff_Q.assert_gradient_error(rng, U, U, 1e-3, 1e-8) - Q = g.qcd.gauge.topological_charge_5LI(U, cache={}) eps = abs(Q - 0.32270083147744544) g.message(f"Test 5LI Q definition: {eps}") assert eps < 1e-13 +# Test differentiable energy_density +g.message("Test diff E") +E = g.qcd.gauge.energy_density(U) +dE = g.qcd.gauge.differentiable_energy_density(adU) +diff_E = dE.functional(*adU) +diff_E.assert_gradient_error(rng, U, U, 1e-3, 1e-8) +assert abs(E - diff_E(U)) < 1e-13 + # Test gauge actions for action in [g.qcd.gauge.action.wilson(5.43), g.qcd.gauge.action.iwasaki(5.41)]: # test action double precision versus quadruple precision