plot-luscher-trajectory.py 2.11 KB
Newer Older
1 2
#!/usr/bin/python3

3 4 5 6 7 8
# Copyright 2019 Christoph Conrads
#
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
import matplotlib.pyplot as pp
import numpy as np
import scipy.optimize
import sys

if len(sys.argv) < 5:
	msg = '5 command line arguments needed, {:d} given'.format(len(sys.argv))
	print(msg, file=sys.stderr)
	sys.exit(1)

kind = sys.argv[1]
w = int(sys.argv[2])
r = int(sys.argv[3])
s = int(sys.argv[4])
filename = "luscher-trajectory-{:s}-{:d}-{:d}-{:d}.png".format(kind, w, r, s)

if w < 8 or r <= s or s < 2:
	print('Error: w={:d} r={:d} s={:d}'.format(w, r, s), file=sys.stderr)
	sys.exit(1)

data = np.loadtxt(sys.stdin)
xs = data[:,0]
ys = data[:,1]

if len(xs.shape) != 1:
	sys.exit(1)

if min(ys) < 0 or max(ys) >= 1:
	print('Error: min={:e} max={:e}'.format(min(ys),max(ys)), file=sys.stderr)
	sys.exit(1)


# fit function
nnz, = np.nonzero(ys > 2**-3)
k0 = None if nnz.size == 0 else min(nnz)

b = 2**w
t = (xs > 0) & (xs > 0 if k0 is None else xs <= k0)

assert np.any(t)

result = scipy.optimize.lsq_linear(np.expand_dims(xs[t]-1,1), np.log(1.0*b*ys[t]))
c = np.asscalar(result.x)

# The constant below is the expected distance of uniformly distributed points in
# a hypercube with word size b, cf. Lüscher (1993), Section 4.1.
# This expression is hidden behind the magic value  12/25 on p. 9.
y_min = (b/2) / (b+1)
# Compute the smallest x rounded up such that 1/b * e^(c(x-1)) >= y_min
k = int(np.ceil(np.log(y_min*b) / c + 1))


61
print('{:4s} {:2d} {:2d} {:2d}  {:3d}'.format(kind, w, r, s, k))
62 63 64 65

zs = np.exp(c * (xs-1)) / b

pp.plot(xs, ys, 'b.', markersize=10, label='Measurement')
66
pp.plot(xs, zs, 'g', label='1/b exp({:.2f} (t-1))'.format(c))
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
pp.grid(which="both")

pp.ylim(top=1.0)
pp.ylim(bottom=2.0**-(w+1))
pp.yscale('log', basey=2)

pp.legend(loc='lower right')
pp.xlabel("t")
pp.ylabel("d(x(t),y(t))")
if kind == 'swc':
	pp.title("{:s}(2^{:d}, {:d}, {:d})".format('SWB', w, s, r))
else:
	pp.title("{:s}(2^{:d}, {:d}, {:d})".format(kind.upper(), w, r, s))

pp.savefig(filename)