Commit 4d19595d authored by Christoph Conrads's avatar Christoph Conrads

Improve Lüscher trajectory code

* move Python code out of Bash script into its own file
* fit an exponential function to data and compute time to `t` to chaos
  using the fitted function
* use a more stringent condition to check for `chaos`
* plot measured data and fitted function
parent 979fce77
......@@ -6,58 +6,6 @@
# 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/.
read -r -d '' python_plot_code <<EOF
import matplotlib.pyplot as pp
import numpy as np
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)
type = 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(type, w, r, s)
if w < 16 or r <= s or s < 2:
print('Error: w={:d} r={:d} s={:d}'.format(w, r, s), file=sys.stderr)
sys.exit(1)
xs = np.loadtxt(sys.stdin)
if len(xs.shape) != 1:
sys.exit(1)
if min(xs) < 0 or max(xs) >= 1:
print('Error: min={:e} max={:e}'.format(min(xs),max(xs)), file=sys.stderr)
sys.exit(1)
nnz, = np.nonzero(xs > 0.45)
k = -1 if nnz.size == 0 else min(nnz)
k_str = '>{:2d}'.format(xs.size) if k == -1 else '{:3d}'.format(k)
print('{:4s} {:d} {:2d} {:2d} {:s}'.format(type, w, r, s, k_str))
pp.plot(xs)
pp.grid(which="both")
pp.ylim(top=1.0)
pp.ylim(bottom=2.0**-(w+1))
pp.yscale('log', basey=2)
pp.xlabel("t")
pp.ylabel("d(x(t),y(t))")
if type == 'swc':
pp.title("{:s}(2^{:d}, {:d}, {:d})".format('SWB', w, s, r))
else:
pp.title("{:s}(2^{:d}, {:d}, {:d})".format(type.upper(), w, r, s))
pp.savefig(filename)
EOF
set -e
set -u
......@@ -118,8 +66,8 @@ make_luscher_plot()
configure_cxx_code $type $w $r $s >"$cxx_filename"
c++ -Wextra -Wall -std=c++11 -pedantic "$cxx_filename" -o "$binary"
local stats="$("./$binary" | awk '{print $2}' | \
python3 -c "$python_plot_code" $type $w $r $s)"
local stats="$("./$binary" | \
python3 "../plot-luscher-trajectory.py" $type $w $r $s)"
echo -n "$stats" ' '
......
#!/usr/bin/python3
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))
print('{:4s} {:d} {:2d} {:2d} {:3d}'.format(kind, w, r, s, k))
zs = np.exp(c * (xs-1)) / b
pp.plot(xs, ys, 'b.', markersize=10, label='Measurement')
pp.plot(xs, zs, 'g', label='1/b exp({:.2f} (x-1))'.format(c))
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)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment