Adds scalar field code

Author: gkanwar

Diff for: pybin/scalar_field_aughmc.py

@@ -0,0 +1,76 @@
+"""
+HMC for (complex) scalar fields in arbitrary dimensions.
+"""
+
+import argparse
+
+from scalar_field.scalar_field_hmc import *
+
+
+def main(out_prefix, init_cfg, action, eps, n_leap, iters, skip, therm):
+    ensemble, actions = aug_hmc_ensemble(init_cfg, action, eps, n_leap, iters, skip, therm)
+    print('Ensemble shape = {}.'.format(np.array(ensemble).shape))
+    fname = out_prefix+'.npy'
+    print('Writing ensemble to {}.'.format(fname))
+    np.save(fname, np.array(ensemble))
+    fname = out_prefix+'.S.npy'
+    print('Writing actions to {}.'.format(fname))
+    np.save(fname, np.array(actions))
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description="HMC (with sign flip aug) ensemble for scalar field.")
+    # action
+    parser.add_argument('--m2', type=float, required=True,
+                        help='Lattice M^2 setting inverse spacing')
+    parser.add_argument('--lam', type=float, default=0.0,
+                        help='Coeffient lambda of phi^4 term')
+    parser.add_argument('--complex_type', action='store_true',
+                        help='Use complex phi fields')
+    # logistics
+    parser.add_argument('--out_prefix', type=str, required=True,
+                        help='Output prefix for ensemble')
+    parser.add_argument('--eps', type=float, required=True,
+                        help='HMC traj eps')
+    parser.add_argument('--n_leap', type=int, required=True,
+                        help='HMC traj leaps')
+    parser.add_argument('--iters', type=int, required=True,
+                        help='Number of measurements')
+    parser.add_argument('--therm', type=int, required=True,
+                        help='Number of therm iters')
+    parser.add_argument('--skip', type=int, required=True,
+                        help='Number skipped between meas')
+    parser.add_argument('--seed', type=int)
+    parser.add_argument('--seed_file', type=str)
+    # lattice
+    parser.add_argument('dims', metavar='d', type=int, nargs='+',
+                        help='Dimensions')
+    args = parser.parse_args()
+    print('Running with args = {}'.format(args))
+    # parse args
+    action = Action([ScalarKineticMassTerm(args.m2)])
+    if args.lam != 0.0:
+        action.terms.append(ScalarPhi4Term(args.lam))
+    print('Using action = {}'.format(action))
+    if args.seed is None:
+        args.seed = np.random.randint(np.iinfo('uint32').max)
+        print("Generated random seed = {}".format(args.seed))
+    np.random.seed(args.seed)
+    print("Using seed = {}.".format(args.seed))
+    if args.seed_file is not None:
+        print("Loading seed cfg from {}".format(args.seed_file))
+        dtype = np.complex128 if args.complex_type else np.float64
+        phi = np.fromfile(
+            args.seed_file, dtype=dtype, count=np.prod(args.dims))
+        phi = phi.reshape(args.dims)
+    else:
+        print("Generating hot start cfg.")
+        phi = 0.3*np.random.normal(size=args.dims)
+        if args.complex_type:
+            phi = phi + 1j * 0.3*np.random.normal(size=args.dims)
+    # run
+    main(args.out_prefix, phi, action, args.eps, args.n_leap,
+         args.iters, args.skip, args.therm)
+    # final logs
+    print('Again, args = {}'.format(args))
+

Diff for: pybin/scalar_field_hmc.py

@@ -0,0 +1,75 @@
+"""
+HMC for (complex) scalar fields in arbitrary dimensions.
+"""
+
+import argparse
+
+from scalar_field.scalar_field_hmc import *
+
+
+def main(out_prefix, init_cfg, action, eps, n_leap, iters, skip, therm):
+    ensemble, actions = hmc_ensemble(init_cfg, action, eps, n_leap, iters, skip, therm)
+    print('Ensemble shape = {}.'.format(np.array(ensemble).shape))
+    fname = out_prefix+'.npy'
+    print('Writing ensemble to {}.'.format(fname))
+    np.save(fname, np.array(ensemble))
+    fname = out_prefix+'.S.npy'
+    print('Writing actions to {}.'.format(fname))
+    np.save(fname, np.array(actions))
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="HMC ensemble for scalar field.")
+    # action
+    parser.add_argument('--m2', type=float, required=True,
+                        help='Lattice M^2 setting inverse spacing')
+    parser.add_argument('--lam', type=float, default=0.0,
+                        help='Coeffient lambda of phi^4 term')
+    parser.add_argument('--complex_type', action='store_true',
+                        help='Use complex phi fields')
+    # logistics
+    parser.add_argument('--out_prefix', type=str, required=True,
+                        help='Output prefix for ensemble')
+    parser.add_argument('--eps', type=float, required=True,
+                        help='HMC traj eps')
+    parser.add_argument('--n_leap', type=int, required=True,
+                        help='HMC traj leaps')
+    parser.add_argument('--iters', type=int, required=True,
+                        help='Number of measurements')
+    parser.add_argument('--therm', type=int, required=True,
+                        help='Number of therm iters')
+    parser.add_argument('--skip', type=int, required=True,
+                        help='Number skipped between meas')
+    parser.add_argument('--seed', type=int)
+    parser.add_argument('--seed_file', type=str)
+    # lattice
+    parser.add_argument('dims', metavar='d', type=int, nargs='+',
+                        help='Dimensions')
+    args = parser.parse_args()
+    print('Running with args = {}'.format(args))
+    # parse args
+    action = Action([ScalarKineticMassTerm(args.m2)])
+    if args.lam != 0.0:
+        action.terms.append(ScalarPhi4Term(args.lam))
+    print('Using action = {}'.format(action))
+    if args.seed is None:
+        args.seed = np.random.randint(np.iinfo('uint32').max)
+        print("Generated random seed = {}".format(args.seed))
+    np.random.seed(args.seed)
+    print("Using seed = {}.".format(args.seed))
+    if args.seed_file is not None:
+        print("Loading seed cfg from {}".format(args.seed_file))
+        dtype = np.complex128 if args.complex_type else np.float64
+        phi = np.fromfile(
+            args.seed_file, dtype=dtype, count=np.prod(args.dims))
+        phi = phi.reshape(args.dims)
+    else:
+        print("Generating hot start cfg.")
+        phi = 0.3*np.random.normal(size=args.dims)
+        if args.complex_type:
+            phi = phi + 1j * 0.3*np.random.normal(size=args.dims)
+    # run
+    main(args.out_prefix, phi, action, args.eps, args.n_leap,
+         args.iters, args.skip, args.therm)
+    # final logs
+    print('Again, args = {}'.format(args))
+

Diff for: pybin/scalar_field_twopt.py

@@ -0,0 +1,69 @@
+"""
+Measure two-point function given scalar field ensemble.
+"""
+
+from scalar_field.scalar_field import *
+
+import argparse
+import numpy as np
+import os
+import sys
+import time
+
+def compute_twopts(cfgs, verbose):
+    start = time.time()
+    corrs = []
+    Nd = len(cfgs.shape[1:])
+    spatial_axes = tuple(range(0,Nd-1))
+    for i,c in enumerate(cfgs):
+        if verbose and i % 10 == 0:
+            print('Cfg {} / {} ({:.2f}s)'.format(
+                i+1, cfgs.shape[0], time.time()-start))
+        corr = np.array(all_corrs(c, xspace, tspace))
+        # corr = np.sum(np.mean(corr, axis=0), axis=spatial_axes) # mom proj
+        corr = np.mean(corr, axis=0)
+        corrs.append(corr)
+    print('Two-points done. '
+          'Total time = {:.2f}s.'.format(time.time()-start))
+    return corrs
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description='Measure two-pt function for scalar field theory')
+    # ensemble path
+    parser.add_argument('--path', type=str, help='Path to .npy format ensemble')
+    parser.add_argument('--tag', type=str, help='Prefix for .dat format ensemble')
+    parser.add_argument('--Ncfg', type=int)
+    parser.add_argument('--use_real64', action='store_true')
+    parser.add_argument('--use_real32', action='store_true')
+    # corr params
+    parser.add_argument('--quiet', action='store_true')
+    parser.add_argument('dims', metavar='d', type=int, nargs='+')
+    args = parser.parse_args()
+    print('Running with args = {}'.format(args))
+
+    if args.path is not None: # for .npy ensembles
+        print('Reading ensemble from {}.'.format(args.path))
+        cfgs = np.load(args.path)
+        dtype = cfgs.dtype
+        args.tag = os.path.splitext(args.path)[0]
+    else: # for (older format) .dat ensembles
+        assert args.tag is not None
+        assert args.Ncfg is not None
+        assert args.dims is not None
+        if args.use_real32:
+            fname = args.tag + '.cfgs.dat'
+            dtype = np.float32
+        else:
+            fname = args.tag + '.dat'
+            dtype = np.float64 if args.use_real64 else np.complex128
+        print('Reading ensemble from {}.'.format(fname))
+        cfgs = np.fromfile(fname, dtype=dtype)
+        cfgs = cfgs.reshape(args.Ncfg, *args.dims)
+    corrs = compute_twopts(cfgs, verbose=not args.quiet)
+    twopt = np.array(corrs)
+    print('twopt shape {}'.format(twopt.shape))
+    fname = args.tag + '.twopt.dat'
+    print('Writing two-pts to {}.'.format(fname))
+    twopt.tofile(fname)
+            

Diff for: scalar_field/__init__.py

@@ -0,0 +1,161 @@
+"""
+Library for lattice scalar field functionality.
+"""
+
+import itertools
+import math
+import numpy as np
+
+### UTILITIES
+def tup_add(t1, t2):
+    assert len(t1) == len(t2)
+    return tuple([t1[i] + t2[i] for i in range(len(t1))])
+def tup_sub(t1, t2):
+    assert len(t1) == len(t2)
+    return tuple([t1[i] - t2[i] for i in range(len(t1))])
+
+class ActionTerm(object):
+    def local_action(self, x, phi):
+        raise NotImplementedError
+    def action(self, phi):
+        raise NotImplementedError
+    def force(self, phi):
+        raise NotImplementedError
+    def __str__(self):
+        raise NotImplementedError
+    def __repr__(self):
+        raise NotImplementedError
+
+class Action(object):
+    def __init__(self, terms):
+        self.terms = terms
+    def local_action(self, x, phi):
+        out = 0.0
+        for term in self.terms:
+            out += term.local_action(x, phi)
+        return out
+    def action(self, phi, verbose=False):
+        out = np.zeros(phi.shape, dtype=phi.dtype)
+        for term in self.terms:
+            act = term.action(phi)
+            if verbose: print('Action from {} = {}'.format(term, np.sum(act)))
+            out += act
+        return out
+    def force(self, phi, verbose=False):
+        out = np.zeros(phi.shape, dtype=phi.dtype)
+        for term in self.terms:
+            F = term.force(phi)
+            if verbose: print('Force from {} = {:.8g}'.format(term, np.linalg.norm(F)))
+            out += F
+        return out
+    def __str__(self):
+        return str(self.terms)
+    __repr__ = __str__
+
+continuum_param = True
+
+class ScalarKineticMassTerm(ActionTerm):
+    def __init__(self, m2):
+        self.m2 = m2
+    def local_action(self, x, phi):
+        dims = phi.shape
+        Nd = len(dims)
+        coords = [x]
+        for mu in range(Nd):
+            muhat = [0]*Nd
+            muhat[mu] = 1
+            coords.append(tup_add(x, muhat))
+            coords.append(tup_sub(x, muhat))
+        coords_arrs = tuple([[] for mu in range(Nd)])
+        for tup in coords:
+            for mu in range(Nd):
+                coords_arrs[mu].append(tup[mu])
+        lin_inds = np.ravel_multi_index(coords_arrs, dims, mode='wrap')
+        SK = 0.0
+        diag = phi[x]
+        for i in range(Nd):
+            fwd_coord = list(x)
+            fwd_coord[i] = (fwd_coord[i] + 1) % dims[i]
+            bwd_coord = list(x)
+            bwd_coord[i] = (bwd_coord[i] - 1) % dims[i]
+            fwd = phi[tuple(fwd_coord)]
+            bwd = phi[tuple(bwd_coord)]
+            SK += 2*np.conj(diag) * (diag - fwd - bwd) # factor of 2?
+        SK = np.real(SK)
+        SM = self.m2 * np.abs(diag)**2
+        pre = 0.5 if continuum_param else 1.0
+        return pre * (SK + SM)
+    def action(self, phi):
+        SK = np.zeros(phi.shape, dtype=np.float64)
+        for mu in range(len(phi.shape)):
+            fwd = np.roll(phi, -1, axis=mu)
+            bwd = np.roll(phi, 1, axis=mu)
+            SK += np.real(np.conj(phi) * (2*phi - fwd - bwd))
+        SM = self.m2 * np.abs(phi)**2
+        pre = 0.5 if continuum_param else 1.0
+        return pre * (SK + SM)
+    def force(self, phi):
+        Nd = len(phi.shape)
+        out = 2 * self.m2 * np.conj(phi)
+        for mu in range(Nd):
+            out += 2 * (
+                2*np.conj(phi)
+                - np.conj(np.roll(phi, 1, axis=mu))
+                - np.conj(np.roll(phi, -1, axis=mu)))
+        # F = -dS/dphi            
+        pre = 0.5 if continuum_param else 1.0
+        return -np.conj(out) * pre
+    def __str__(self):
+        return "KineticMassTerm(M2={:.4f})".format(self.m2)
+    __repr__ = __str__
+    
+class ScalarPhi4Term(ActionTerm):
+    def __init__(self, lam):
+        self.lam = lam
+    def local_action(self, x, phi):
+        return self.lam * np.abs(phi[x])**4
+    def action(self, phi):
+        return self.lam * np.abs(phi)**4
+    def force(self, phi):
+        return -np.conj(4 * self.lam * np.conj(phi) * np.abs(phi)**2)
+    def __str__(self):
+        return "Phi4Term(lambda={:.8f})".format(self.lam)
+    __repr__ = __str__
+
+def wrap(p):
+    return (p + math.pi) % (2*math.pi) - math.pi
+
+def make_smear_corr(smear_x, smear_t):
+    def smear_corr(corr):
+        L = corr.shape
+        Nd = len(L)
+        corr_fft = np.fft.fftn(corr)
+        ks = [ np.linspace(0, 2*math.pi, num=L_mu, endpoint=False)
+               for L_mu in L ]
+        weights = []
+        for mu in range(Nd):
+            assert(L[mu] % 2 == 0)
+            ks[mu] = (ks[mu] + math.pi) % (2*math.pi) - math.pi
+            radius = smear_t if mu == Nd-1 else smear_x
+            weights.append(np.exp(-radius**2 * np.square(ks[mu]) / Nd))
+        ein_inds = ",".join(map(chr, range(ord('a'), ord('a') + Nd)))
+        weights = np.einsum(ein_inds, *weights) # big outer product
+        assert(corr_fft.shape == weights.shape)
+        corr_fft *= weights
+        return np.fft.ifftn(corr_fft)
+    return smear_corr
+
+def all_corrs(phi, xspace, tspace):
+    corrs = []
+    coord_ranges = []
+    for Lx in phi.shape[:-1]:
+        assert(Lx % xspace == 0)
+        coord_ranges.append(range(0, Lx, xspace))
+    assert(phi.shape[-1] % tspace == 0)
+    coord_ranges.append(range(0, phi.shape[-1], tspace))
+    all_axes = tuple(range(len(phi.shape)))
+    for src in itertools.product(*coord_ranges):
+        shift = tuple(-np.array(src))
+        corr = np.conj(phi[src]) * np.roll(phi, shift, axis=all_axes)
+        corrs.append(corr)
+    return np.array(corrs)