diff --git a/applications/hmc/cpn.py b/applications/hmc/cpn.py new file mode 100755 index 000000000..59e6599ff --- /dev/null +++ b/applications/hmc/cpn.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python3 +# +# Authors: Christoph Lehner, Mattia Bruno, Gabriele Morandi 2023 +# +# HMC for 2D CP^N-1 theory +# +import gpt as g +import sys, os +import numpy + +g.default.set_verbose("step_size", False) +grid = g.grid([42, 42], g.double) +rng = g.random("hmc-cpn-model") + +# action for conj. momenta: +g.message(f"Lattice = {grid.fdimensions}") +g.message("Actions:") +a0 = g.qcd.scalar.action.mass_term() +g.message(f" - {a0.__name__}") + +# CP^N-1 action +beta = 0.70 +N = 10 +a1 = g.qcd.scalar.action.cpn(N, beta) +g.message(f" - {a1.__name__}") + +# fields +z = g.vcomplex(grid, N) +a1.draw(z, rng) +l = [g.u1(grid) for _ in range(grid.nd)] +rng.element(l) +fields = [z] + l + +# conjugate momenta +mom_z = g.group.cartesian(z) +mom_l = g.group.cartesian(l) +moms = [mom_z] + mom_l + + +def hamiltonian(): + return a0(mom_z) + a0(mom_l) + a1(fields) + + +# molecular dynamics +sympl = g.algorithms.integrator.symplectic + +ip_z = sympl.update_p(mom_z, lambda: a1.gradient(fields, z)) +ip_l = sympl.update_p(mom_l, lambda: a1.gradient(fields, l)) + + +class constrained_iq(sympl.symplectic_base): + def __init__(self, fields, moms): + z = fields[0] + l = fields[1:] + mom_z = moms[0] + mom_l = moms[1:] + + iq_l = sympl.update_q(l, lambda: a0.gradient(mom_l, mom_l)) + + def inner(eps): + a1.constrained_leap_frog(eps, z, mom_z) + iq_l(eps) + + super().__init__(1, [], [inner], None, "constrained iq") + + +iq = constrained_iq(fields, moms) + +# integrator +mdint = sympl.leap_frog(50, [ip_z, ip_l], iq) +g.message(f"Integration scheme:\n{mdint}") + +# metropolis +metro = g.algorithms.markov.metropolis(rng) + +# MD units +tau = 1.0 +g.message(f"tau = {tau} MD units") + + +def hmc(tau, moms): + mom_z = moms[0] + mom_l = moms[1:] + + rng.normal_element(mom_l) + a1.draw(mom_z, rng, z) + + accrej = metro(fields) + h0 = hamiltonian() + mdint(tau) + h1 = hamiltonian() + return [accrej(h1, h0), h1 - h0] + + +# thermalization +for i in range(1, 21): + h = [] + for _ in range(10): + h += [hmc(tau, moms)] + h = numpy.array(h) + g.message(f"{i*5} % of thermalization completed") + g.message( + f"Action = {a1(fields)}, Acceptance = {numpy.mean(h[:,0]):.2f}, |dH| = {numpy.mean(numpy.abs(h[:,1])):.4e}" + ) + +# measure action +def measure(): + return [a1(fields) / (N * beta * grid.fsites * grid.nd)] + + +# production +history = [] +data = [] +for i in range(100): + for k in range(10): + history += [hmc(tau, moms)] + data += [measure()] + g.message(f"Trajectory {i}") + +history = numpy.array(history) +g.message(f"Acceptance rate = {numpy.mean(history[:,0]):.2f}") +g.message(f"<|dH|> = {numpy.mean(numpy.abs(history[:,1])):.4e}") + +data = numpy.array(data) +g.message(f"Energy density = {numpy.mean(data[:,0])}") diff --git a/lib/cgpt/lib/eval/mul_vlat_vlat.h b/lib/cgpt/lib/eval/mul_vlat_vlat.h index 022b88f99..62ff122ae 100644 --- a/lib/cgpt/lib/eval/mul_vlat_vlat.h +++ b/lib/cgpt/lib/eval/mul_vlat_vlat.h @@ -59,9 +59,9 @@ void eval_mul_vlat_vlat(std::vector & dst_vl, // VV -> S/M if (lhs_singlet_rank == 1 && rhs_singlet_rank == 1) { + ASSERT(lhs_singlet_dim == rhs_singlet_dim); if (lhs_unary == 0 && rhs_unary == (BIT_TRANS|BIT_CONJ)) { // outer product -> M - ASSERT(lhs_singlet_dim == rhs_singlet_dim); dst_vl.resize(lhs_singlet_dim*rhs_singlet_dim, 0); for (int i=0;i & dst_vl, } return; } else if (lhs_unary == (BIT_TRANS|BIT_CONJ) && rhs_unary == 0) { - ERR("Not implemented"); + //ERR("Not implemented"); // inner product -> S - /*ASSERT(lhs_singlet_dim == rhs_singlet_dim); dst_vl.resize(1, 0); bool _ac = ac; for (int i=0;i & dst_vl, _ac = true; } } - return;*/ + return; } else { ERR("Invalid combination of two vectors"); } diff --git a/lib/cgpt/lib/expression/unary.h b/lib/cgpt/lib/expression/unary.h index a06d9e606..c80bc7cf5 100644 --- a/lib/cgpt/lib/expression/unary.h +++ b/lib/cgpt/lib/expression/unary.h @@ -78,7 +78,7 @@ cgpt_Lattice_base* lattice_lat(cgpt_Lattice_base* dst, bool ac, const A& lat, Co template cgpt_Lattice_base* lattice_unary_lat(cgpt_Lattice_base* dst, bool ac, const A& la,int unary_expr,ComplexD coef) { if (unary_expr == 0) { - return lattice_lat(dst, ac, la, coef); + return lattice_lat(dst, ac, closure(ToSinglet(la)), coef); } else if (unary_expr == (BIT_SPINTRACE|BIT_COLORTRACE)) { return lattice_expr(dst, ac, coef*ToSinglet(trace(la))); } else if (unary_expr == BIT_SPINTRACE) { diff --git a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc index 5b98f13e5..d139d1e75 100644 --- a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc +++ b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc @@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a typedef vComplexD vtype; _COMPATIBLE_MSR_(iSinglet); _OUTER_PRODUCT_(iVSinglet10); - // _INNER_PRODUCT_(iVSinglet10); + _INNER_PRODUCT_(iVSinglet10); ERR("Not implemented"); } diff --git a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc index 90bc07f7c..492c1d376 100644 --- a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc +++ b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc @@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a typedef vComplexD vtype; _COMPATIBLE_MSR_(iSinglet); _OUTER_PRODUCT_(iVSinglet30); - // _INNER_PRODUCT_(iVSinglet30); + _INNER_PRODUCT_(iVSinglet30); ERR("Not implemented"); } diff --git a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc index 760f10284..7f6325874 100644 --- a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc +++ b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc @@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a typedef vComplexD vtype; _COMPATIBLE_MSR_(iSinglet); _OUTER_PRODUCT_(iVSinglet4); - // _INNER_PRODUCT_(iVSinglet4); + _INNER_PRODUCT_(iVSinglet4); ERR("Not implemented"); } diff --git a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc index ff1b993d1..f2be2da1e 100644 --- a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc +++ b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc @@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a typedef vComplexF vtype; _COMPATIBLE_MSR_(iSinglet); _OUTER_PRODUCT_(iVSinglet10); - // _INNER_PRODUCT_(iVSinglet10); + _INNER_PRODUCT_(iVSinglet10); ERR("Not implemented"); } diff --git a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc index 8dd56db10..9391e0c7d 100644 --- a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc +++ b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc @@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a typedef vComplexF vtype; _COMPATIBLE_MSR_(iSinglet); _OUTER_PRODUCT_(iVSinglet30); - // _INNER_PRODUCT_(iVSinglet30); + _INNER_PRODUCT_(iVSinglet30); ERR("Not implemented"); } diff --git a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc index 7abd5db97..904e7faef 100644 --- a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc +++ b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc @@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a typedef vComplexF vtype; _COMPATIBLE_MSR_(iSinglet); _OUTER_PRODUCT_(iVSinglet4); - // _INNER_PRODUCT_(iVSinglet4); + _INNER_PRODUCT_(iVSinglet4); ERR("Not implemented"); } diff --git a/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template b/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template index 092ee1b81..d7eb6f275 100644 --- a/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template +++ b/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template @@ -40,6 +40,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a typedef {precision_vector} vtype; _COMPATIBLE_MSR_(iSinglet); _OUTER_PRODUCT_(iVSinglet{tensor_arg_1}); - // _INNER_PRODUCT_(iVSinglet{tensor_arg_1}); + _INNER_PRODUCT_(iVSinglet{tensor_arg_1}); ERR("Not implemented"); } diff --git a/lib/gpt/core/domain/sparse.py b/lib/gpt/core/domain/sparse.py index ce330266e..9d5522ca0 100644 --- a/lib/gpt/core/domain/sparse.py +++ b/lib/gpt/core/domain/sparse.py @@ -104,18 +104,30 @@ def slice(self, fields, ortho_dim): length = self.grid.gdimensions[ortho_dim] return gpt.indexed_sum(fields, self.coordinate_lattices()[ortho_dim], length) - def weight(self): - if self.weight_cache is not None: - return self.weight_cache + def global_coordinates(self): + # broadcast below uses grid.processor not gpt.rank() + rank_id = self.grid.processor # this function is mostly used for tests, so it is not performance critical for rank in range(self.grid.Nprocessors): - local_coordinates = np.copy(self.local_coordinates) + nlc = np.array([len(self.local_coordinates)]) + self.grid.broadcast(rank, nlc) + if rank == rank_id: + local_coordinates = np.copy(self.local_coordinates) + else: + local_coordinates = np.zeros((nlc[0], self.grid.nd), dtype=self.local_coordinates.dtype) self.grid.broadcast(rank, local_coordinates) if rank == 0: global_coordinates = local_coordinates else: global_coordinates = np.concatenate((global_coordinates, local_coordinates)) + return global_coordinates + + def weight(self): + if self.weight_cache is not None: + return self.weight_cache + + global_coordinates = self.global_coordinates() unique_coordinates, count = np.unique(global_coordinates, axis=0, return_counts=True) unique_coordinates = unique_coordinates.view(type(global_coordinates)) @@ -208,6 +220,27 @@ def converted(self, precision): mask=mask, ) + # every rank must pass the same list of positions, otherwise code breaks + def restricted(self, positions): + # get emb_coor of wanted positions on each rank to build mask + ec = self.unique_embedded_coordinates(positions) + + # creates full mask + mask = self.lattice(gpt.ot_real_additive_group()) + mask[:] = 0 + mask[ec] = 1 + + cl = self.coordinate_lattices(mark_empty=-1) + local_coordinates = np.hstack(tuple([x[:].real.astype(np.int32) for x in cl])) + mask_local = [m[0]==1 for m in mask[:]] + + return gpt.core.domain.sparse( + self.kernel.grid, + local_coordinates, + dimensions_divisible_by=cl[0].grid.fdimensions, + mask = mask_local + ) + def weight(self): return self.kernel.weight() diff --git a/lib/gpt/qcd/scalar/action/__init__.py b/lib/gpt/qcd/scalar/action/__init__.py index defb512fc..1b3994099 100644 --- a/lib/gpt/qcd/scalar/action/__init__.py +++ b/lib/gpt/qcd/scalar/action/__init__.py @@ -18,6 +18,7 @@ # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # from gpt.qcd.scalar.action.phi4 import phi4 +from gpt.qcd.scalar.action.cpn import cpn from gpt.qcd.scalar.action.mass_term import mass_term, fourier_mass_term, general_mass_term from gpt.qcd.scalar.action.coupling import coupling import gpt.qcd.scalar.action.hermitian_kernel diff --git a/lib/gpt/qcd/scalar/action/cpn.py b/lib/gpt/qcd/scalar/action/cpn.py new file mode 100644 index 000000000..3dd3f01e5 --- /dev/null +++ b/lib/gpt/qcd/scalar/action/cpn.py @@ -0,0 +1,122 @@ +# +# GPT - Grid Python Toolkit +# Copyright (C) 2023 Christoph Lehner (christoph.lehner@ur.de, https://github.com/lehner/gpt) +# 2023 Mattia Bruno, Gabriele Morandi +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +# +import gpt as g +from gpt.core.group import differentiable_functional +import numpy + +# CP^{N-1} model + +# S[z,l] = -N * beta * sum_n,mu (z_{n+mu}^dag * z_n * l_{n,mu} + z_n^dag * z_{n+mu} * l_{n,mu}^dag) +# = -2 * N * beta * sum_n,mu [Re(z_{n+mu}^dag * z_n * l_{n,mu}) - 1] +class cpn(differentiable_functional): + def __init__(self, N, b): + self.N = N + self.beta = b + self.__name__ = f"cpn(N = {self.N}, beta = {self.beta})" + + # z = fields[0], l = fields[1:] + def split(self, fields): + return fields[0], fields[1:] + + def __call__(self, fields): + z, l = self.split(fields) + + J = g.lattice(z) + J[:] = 0.0 + for mu in range(z.grid.nd): + J += g.cshift(z, mu, +1) * g.adj(l[mu]) + + action = -2 * self.N * self.beta * (g.inner_product(J, z).real - z.grid.fsites * z.grid.nd) + return action + + @differentiable_functional.multi_field_gradient + def gradient(self, fields, dfields): + def gradient_l(z, l, mu): + frc = g.lattice(l[0]) + frc @= ( + 2 + * self.beta + * self.N + * g.component.imag(g.trace(z * g.adj(g.cshift(z, mu, +1))) * l[mu]) + ) + frc.otype = l[0].otype.cartesian() + return frc + + def gradient_z(z, l): + J = g.lattice(z) + J[:] = 0.0 + for mu in range(z.grid.nd): + J += g.cshift(z, mu, +1) * g.adj(l[mu]) + J += g.cshift(z * l[mu], mu, -1) + + frc = g.lattice(z) + frc @= -2 * self.beta * self.N * J + + frc -= g.trace(frc * g.adj(z)) * z + return frc + + z, l = self.split(fields) + frc = [] + for df in g.core.util.to_list(dfields): + k = fields.index(df) + if k == 0: + frc.append(gradient_z(z, l)) + else: + frc.append(gradient_l(z, l, k - 1)) + return frc + + # https://arxiv.org/abs/1102.1852 + def constrained_leap_frog(self, eps, z, mom_z): + # TO DO: replace with g.adj(v1) * v2 + def dot(v1, v2): + return g.trace(v2 * g.adj(v1)) + + n = g.real(z.grid) + n @= g.component.sqrt(g.component.real(dot(mom_z, mom_z))) + + # z' = cos(alpha) z + (1/|pi|) sin(alpha) mom_z + # mom_z' = -|pi| sin(alpha) z + cos(alpha) mom_z + # alpha = eps |pi| + _z = g.lattice(z) + _z @= z + + cos = g.real(z.grid) + cos @= g.component.cos(eps * n) + + sin = g.real(z.grid) + sin @= g.component.sin(eps * n) + + z @= cos * _z + g(g.component.inv(n) * sin) * mom_z + mom_z @= -g(n * sin) * _z + cos * mom_z + del _z, cos, sin, n + + # https://arxiv.org/abs/1102.1852 + def draw(self, field, rng, constraint=None): + if constraint is None: + z = field + rng.element(z) + n = g.component.real(g.trace(z * g.adj(z))) + z @= z * g.component.inv(g.component.sqrt(n)) + else: + mom_z = field + z = constraint + rng.normal_element(mom_z) + # TO DO change to z * g(g.adj(z) * mom_z) + mom_z @= mom_z - g(z * g.adj(z)) * mom_z diff --git a/tests/core/sparse_domain.py b/tests/core/sparse_domain.py new file mode 100644 index 000000000..d0100318c --- /dev/null +++ b/tests/core/sparse_domain.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 +# +# Authors: Christoph Lehner 2020 Mattia Bruno 2026 +# +# Desc.: Test small features of sparse domains +# +import gpt as g +import numpy as np +import os + +grid = g.grid([16,16,16,32], g.double) + +position_seed = 'sparse-domain' +rng = g.random(position_seed) + +nsnk = 128 +all_positions = np.array([ + [rng.uniform_int(min=0, max=grid.fdimensions[i] - 1) for i in range(4)] for j in range(nsnk) +], dtype=np.int32) + +nsnk_local = nsnk // grid.Nprocessors + +sdomain = g.domain.sparse(grid, all_positions[g.rank() * nsnk_local: (g.rank()+1) * nsnk_local]) + +assert abs(g.sum(sdomain.weight()) - nsnk) < 1e-16 +assert abs(g.sum(sdomain.kernel.cached_one_mask()) - nsnk) < 1e-16 + +nsrc = 16 +N = nsnk - nsrc +subset = all_positions[nsrc:] + +sdomain2 = sdomain.restricted(subset) + +g.message(f'Restricted sparse domain from {g.sum(sdomain.weight())} points down to {g.sum(sdomain2.weight())}') + +assert abs(g.sum(sdomain2.weight()) - N) < 1e-16 +assert abs(g.sum(sdomain2.kernel.cached_one_mask()) - N) < 1e-16 + +# sparse domain irregularly distributed over nodes + +def gsum(lat): + idx = g.real(lat.grid) + idx[:] = 0 + return g.indexed_sum(lat, idx, 1)[0] + +np.random.seed(46) + +def random_partition_sizes(N, M, spread=2): + base = N // M + sizes = [base] * M + if M==1: + return sizes + + # distribute remainder randomly + for _ in range(N - base * M): + sizes[np.random.randrange(M)] += 1 + + # random perturbations + for _ in range(spread * M): + i, j = np.random.choice(M, size=2, replace=False) + + if sizes[i] > 1: + sizes[i] -= 1 + sizes[j] += 1 + + return sizes + +nsnk_local = random_partition_sizes(nsnk, g.ranks(), spread=int(nsnk*0.1)) +sdomain3 = g.domain.sparse(grid, all_positions[g.rank() * nsnk_local[g.rank()]: (g.rank()+1) * nsnk_local[g.rank()]]) + +assert abs(gsum(sdomain.weight()) - nsnk) < 1e-16 +assert abs(gsum(sdomain.kernel.cached_one_mask()) - nsnk) < 1e-16 \ No newline at end of file diff --git a/tests/core/vcomplex.py b/tests/core/vcomplex.py new file mode 100644 index 000000000..d1b5097f7 --- /dev/null +++ b/tests/core/vcomplex.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python3 +# +# Authors: Christoph Lehner 2020, Mattia Bruno 2023 +# +# Desc.: Illustrate core concepts and features +# +import gpt as g + +L = [8, 4, 4, 4] +grid = g.grid(L, g.double) +rng = g.random("vcomplex") + +def test_local_inner_product(v): + i0 = g(g.adj(v) * v) + i1 = g(g.trace(v * g.adj(v))) + assert g.norm2(i0 - i1) < 1e-15 + +for N in [4, 10, 30]: + v1 = g.vcomplex(grid, N) + test_local_inner_product(v1) + del v1