From e690e205833be31bd489d316d29a206c354c743d Mon Sep 17 00:00:00 2001 From: "Andrew R. McCluskey" Date: Tue, 30 Jul 2024 13:10:59 +0100 Subject: [PATCH] Latest version --- .github/workflows/build.yml | 2 +- Snakefile | 77 +- environment.yml | 7 +- showyourwork.yml | 9 +- src/code/llzo/glswlsols.py | 12 +- src/code/llzo/many_runs.py | 36 +- src/code/llzo/many_runs_old.py | 105 + src/code/llzo/true_ls.py | 12 +- src/code/random_walks/efficiency.py | 2 +- src/code/random_walks/glswlsols.py | 4 +- .../kinisi_pyblock_modelfree_rw.py | 46 - src/code/random_walks/kinisi_pyblock_rw.py | 13 +- src/code/random_walks/kinisi_rw.py | 19 +- src/code/random_walks/ordinary_rw.py | 4 +- src/code/random_walks/true_ls.py | 4 +- src/code/random_walks/truecov_rw.py | 6 +- src/code/random_walks/weighted_rw.py | 6 +- src/scripts/covariances.py | 10 +- src/scripts/diffusion.py | 24 +- src/scripts/glswlsols_llzo.py | 4 +- src/scripts/pyblock.py | 86 +- src/scripts/random_walk.py | 23 +- src/static/schematic.pdf | Bin 21492 -> 33698 bytes src/static/schematic.svg | 1957 ++++------------- src/tex/bib.bib | 42 +- src/tex/ms.tex | 204 +- zenodo.yml | 2 +- 27 files changed, 798 insertions(+), 1918 deletions(-) create mode 100755 src/code/llzo/many_runs_old.py delete mode 100644 src/code/random_walks/kinisi_pyblock_modelfree_rw.py diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 4bac268..8eb8f0d 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -21,4 +21,4 @@ jobs: env: SANDBOX_TOKEN: ${{ secrets.SANDBOX_TOKEN }} with: - showyourwork-spec: git+https://github.com/thomasckng/showyourwork + showyourwork-spec: git+https://github.com/showyourwork/showyourwork \ No newline at end of file diff --git a/Snakefile b/Snakefile index 834ceda..de40db1 100644 --- a/Snakefile +++ b/Snakefile @@ -49,7 +49,7 @@ rule: [f'src/data/random_walks/kinisi/rw_1_{atoms}_128_s512.npz' for atoms in [16, 32, 64, 128, 256, 512, 1024]], [f'src/data/random_walks/weighted/rw_1_{atoms}_128_s512.npz' for atoms in [16, 32, 64, 128, 256, 512, 1024]], [f'src/data/random_walks/ordinary/rw_1_{atoms}_128_s512.npz' for atoms in [16, 32, 64, 128, 256, 512, 1024]], - [f'src/data/random_walks/numerical/D_1_{atoms}_128.npz' for atoms in [16, 32, 64, 128, 256, 512, 1024]] + [f'src/data/random_walks/numerical/D_1_{atoms}_128.npz' for atoms in [16, 32, 64, 128, 256, 512, 1024]], output: f'src/data/random_walks/stat_eff.npz' conda: @@ -58,7 +58,6 @@ rule: True params: correlation='true' - priority: 50 script: "src/code/random_walks/stat_eff.py" @@ -124,26 +123,6 @@ rule rw_4096_pyblock_true_128: script: f'src/code/random_walks/kinisi_pyblock_rw.py' -# Random walks simulations -rule rw_4096_pyblock_modelfree_true_128: - input: - f'src/code/random_walks/kinisi_pyblock_modelfree_rw.py', - 'src/code/random_walks/random_walk.py' - output: - f'src/data/random_walks/pyblock_modelfree/rw_1_128_128_s4096.npz' - conda: - 'environment.yml' - cache: - True - params: - jump=2.4494897428, - atoms=128, - length=128, - correlation='true', - n=4096 - script: - f'src/code/random_walks/kinisi_pyblock_modelfree_rw.py' - # Length variation random walks for length in [16, 32, 64, 128, 256, 512, 1024]: rule: @@ -389,16 +368,16 @@ rule rw_4096_true_cov_true_128: script: f'src/code/random_walks/truecov_rw.py' -for start_diff in [0, 2, 4, 6, 8, 10, 15, 20]: +for l in [10000]: for n in range(0, 6, 1): rule: name: - f"llzo_many_{n}_{start_diff}" + f"llzo_many_{n}_10_{l}" input: 'src/code/llzo/many_runs.py', - # f'src/data/llzo/traj{n}.xyz' + f'src/data/llzo/traj{n}.out' output: - f'src/data/llzo/diffusion_{n}_{start_diff}.npz' + f'src/data/llzo/diffusion_{n}_10_{l}.npz' conda: 'environment.yml' @@ -406,41 +385,43 @@ for start_diff in [0, 2, 4, 6, 8, 10, 15, 20]: True params: n=n, - start_diff=start_diff + start_diff=10, + length=l script: "src/code/llzo/many_runs.py" rule: name: - f"true_D_llzo_{start_diff}" + f"true_D_llzo_10_{l}" input: 'src/code/llzo/true_ls.py', - [f'src/data/llzo/diffusion_{n}_{start_diff}.npz' for n in range(0, 6, 1)] + [f'src/data/llzo/diffusion_{n}_10_{l}.npz' for n in range(0, 6, 1)] output: - f'src/data/llzo/true_{start_diff}.npz' + f'src/data/llzo/true_10_{l}.npz' conda: 'environment.yml' cache: True params: - start_diff=start_diff + start_diff=10, + length=l script: "src/code/llzo/true_ls.py" - rule: - name: - f"glswlsols_D_llzo_{start_diff}" - input: - 'src/code/llzo/glswlsols.py', - [f'src/data/llzo/diffusion_{n}_{start_diff}.npz' for n in range(0, 6, 1)] - output: - f'src/data/llzo/glswlsols_{start_diff}.npz' - conda: - 'environment.yml' - cache: - True - params: - start_diff=start_diff - script: - "src/code/llzo/glswlsols.py" - \ No newline at end of file +rule: + name: + f"glswlsols_D_llzo_10_10000" + input: + 'src/code/llzo/glswlsols.py', + [f'src/data/llzo/diffusion_{n}_10_10000.npz' for n in range(0, 6, 1)] + output: + f'src/data/llzo/glswlsols_10_10000.npz' + conda: + 'environment.yml' + cache: + True + params: + start_diff=10, + length=10000 + script: + "src/code/llzo/glswlsols.py" diff --git a/environment.yml b/environment.yml index 99ebb70..b7cc2d2 100644 --- a/environment.yml +++ b/environment.yml @@ -2,13 +2,14 @@ channels: - defaults - conda-forge dependencies: - - numpy=1.23.1 + - numpy=1.26.4 - pip=21.3.1 - python=3.10 - - scipy=1.9.3 + - scipy=1.12.0 - scikit-learn=1.1.3 - mdanalysis=2.6.1 - pip: - matplotlib==3.8.2 + - seaborn==0.13.2 - pyblock==0.6 - - git+https://github.com/bjmorgan/kinisi.git@with-f \ No newline at end of file + - git+https://github.com/bjmorgan/kinisi.git diff --git a/showyourwork.yml b/showyourwork.yml index b4843b1..6e498cb 100644 --- a/showyourwork.yml +++ b/showyourwork.yml @@ -34,9 +34,9 @@ dependencies: - src/data/random_walks/kinisi/rw_1_128_128_s4096.npz - src/data/random_walks/numerical/D_1_128_128.npz src/scripts/diffusion.py: - - src/data/llzo/true_10.npz + - src/data/llzo/true_10_10000.npz {% for n in range(0, 6, 1) %} - - src/data/llzo/diffusion_{{n}}_10.npz + - src/data/llzo/diffusion_{{n}}_10_10000.npz {% endfor %} - src/scripts/utils/_fig_params.py src/scripts/covariances.py: @@ -45,16 +45,15 @@ dependencies: - src/scripts/utils/_fig_params.py src/scripts/true_cov.py: - src/data/random_walks/numerical/D_1_128_128.npz - - src/data/random_walks/numerical/D_1_128_1024.npz - src/data/random_walks/true_cov/rw_1_128_128_s4096.npz - src/scripts/utils/_fig_params.py src/scripts/stat_eff.py: - src/data/random_walks/stat_eff.npz - src/scripts/utils/_fig_params.py src/scripts/glswlsols_llzo.py: - - src/data/llzo/glswlsols_10.npz + - src/data/llzo/glswlsols_10_10000.npz + - src/data/llzo/true_10_10000.npz src/scripts/pyblock.py: - src/data/random_walks/numerical/D_1_128_128.npz - src/data/random_walks/kinisi/rw_1_128_128_s4096.npz - src/data/random_walks/pyblock/rw_1_128_128_s4096.npz - - src/data/random_walks/pyblock_modelfree/rw_1_128_128_s4096.npz diff --git a/src/code/llzo/glswlsols.py b/src/code/llzo/glswlsols.py index 9ec4e4d..554aa66 100644 --- a/src/code/llzo/glswlsols.py +++ b/src/code/llzo/glswlsols.py @@ -13,15 +13,15 @@ np.random.seed(1) start_diff = snakemake.params['start_diff'] +length = snakemake.params['length'] -timestep = np.load(f'src/data/llzo/diffusion_0_{start_diff}.npz')['dt'][0, 0] -length = 2000 +timestep = np.load(f'src/data/llzo/diffusion_0_{start_diff}_{length}.npz')['dt'][0, 0] ll = len([i + length for i in range(0, 20001 - length, length)]) # length = 500 # ll = len([i + length for i in range(0, 2000, length)]) -d = np.zeros((6, 8, ll, timestep.size)) -for i in range(0, 6, 1): - d[i] = np.load(f'src/data/llzo/diffusion_{i}_{start_diff}.npz')['msd_true'] +d = np.zeros((5, 16, ll, timestep.size)) +for i in range(1, 6, 1): + d[i-1] = np.load(f'src/data/llzo/diffusion_{i}_{start_diff}_{length}.npz')['msd_true'] max_ngp = np.argwhere(timestep > start_diff)[0][0] true_msd = d.reshape(-1, d.shape[-1])[:, max_ngp:] @@ -45,7 +45,7 @@ for i in true_msd: c3.append(linregress(timestep, i).stderr / 6e4) -np.savez(f'src/data/llzo/glswlsols_{start_diff}.npz', +np.savez(f'src/data/llzo/glswlsols_{start_diff}_{length}.npz', gls_pop=g1, wls_pop=g2, ols_pop=g3, diff --git a/src/code/llzo/many_runs.py b/src/code/llzo/many_runs.py index 8f64b24..5ca2d6c 100644 --- a/src/code/llzo/many_runs.py +++ b/src/code/llzo/many_runs.py @@ -1,12 +1,20 @@ import MDAnalysis as mda import numpy as np from kinisi.parser import MDAnalysisParser -from kinisi.diffusion import MSDDiffusion +from kinisi.diffusion import MSDBootstrap n = snakemake.params['n'] start_diff = snakemake.params['start_diff'] +length = snakemake.params['length'] -step = 8 +step = 16 +n_atoms_drift = int((1536-448) / step) +n_atoms_li = int(448 / step) +rng = np.random.RandomState(n*42) +np.random.seed(n*42) +choice_drift = rng.choice(np.arange(0, 1536-448, dtype=int), size=(step, n_atoms_drift), replace=False) +choice_li = rng.choice(np.arange(1536-448, 1536, dtype=int), size=(step, n_atoms_li), replace=False) +choice = np.concatenate([choice_drift, choice_li], axis=-1) def universe_slice(data: str, slice_start: int, slice_end: int, m: int) -> mda.Universe: """ @@ -18,18 +26,17 @@ def universe_slice(data: str, slice_start: int, slice_end: int, m: int) -> mda.U :return: A single flattened Universe """ data_reshape = np.copy(data).reshape(-1, 1536, 3) - data_subset = data_reshape[slice_start:slice_end, m::step] + data_subset = data_reshape[slice_start:slice_end, choice_li[m]] data_subset *= 0.52917721067 u = mda.Universe.empty(data_subset.shape[1], 1, n_frames=data_subset.shape[0], trajectory=True) - list_type = ['A'] * (1536-448) + ['Li'] * 448 - u.add_TopologyAttr('type', list_type[m::step]) + list_type = np.array(['Li'] * n_atoms_li) + u.add_TopologyAttr('type', list(list_type)) for i, t in enumerate(u.trajectory): t.positions = data_subset[i] t.dimensions = np.array([26.30938114, 26.30938114, 26.30938114, 90, 90, 90]) return u data = np.loadtxt(f'src/data/llzo/traj{n}.out', usecols=[0, 1, 2]) -length = 2000 uu = universe_slice(data, 0, length, 0) step_skip = 100 # sampling rate timestep = 5.079648e-4 # ps @@ -45,21 +52,16 @@ def universe_slice(data: str, slice_start: int, slice_end: int, m: int) -> mda.U cov = np.zeros((step, ll, p.delta_t[np.where(p.delta_t > start_diff)].size, p.delta_t[np.where(p.delta_t > start_diff)].size)) d = np.zeros((step, ll, 3200)) g = np.zeros((step, ll, 3200)) -f = np.zeros((step, ll, p.delta_t.size)) intercept = np.zeros((step, ll, 3200)) for m in range(0, step, 1): for i, slice in enumerate(range(0, 20001 - length, length)): uu = universe_slice(data, slice, slice+length, m) - rng = np.random.RandomState(42) - np.random.seed(42) - - da_params = {'specie': 'Li', 'time_step': timestep, 'step_skip': step_skip} + da_params = {'specie': 'Li', 'time_step': timestep, 'step_skip': step_skip, 'progress': False} p = MDAnalysisParser(uu, **da_params) - b = MSDDiffusion(p.delta_t, p.disp_3d, p._n_o) - b.diffusion(start_diff, random_state=rng) - # print(b._popt) + b = MSDBootstrap(p.delta_t, p.disp_3d, p._n_o, progress=False) + b.diffusion(start_diff, random_state=rng, progress=False) dt[m, i] = b.dt msd[m, i] = b.n msd_true[m, i] = b.n @@ -68,10 +70,9 @@ def universe_slice(data: str, slice_start: int, slice_end: int, m: int) -> mda.U cov[m, i] = b.covariance_matrix d[m, i] = b.D.samples g[m, i] = b.gradient.samples - f[m, i] = b._hr intercept[m, i] = b.intercept.samples - -np.savez(f'src/data/llzo/diffusion_{n}_{start_diff}.npz', + +np.savez(f'src/data/llzo/diffusion_{n}_{start_diff}_{length}.npz', dt=dt, msd_true=msd_true, msd=msd, @@ -80,5 +81,4 @@ def universe_slice(data: str, slice_start: int, slice_end: int, m: int) -> mda.U cov=cov, d=d, g=g, - f=f, intercept=intercept) \ No newline at end of file diff --git a/src/code/llzo/many_runs_old.py b/src/code/llzo/many_runs_old.py new file mode 100755 index 0000000..94d5aa6 --- /dev/null +++ b/src/code/llzo/many_runs_old.py @@ -0,0 +1,105 @@ +import MDAnalysis as mda +import numpy as np +from kinisi.parser import MDAnalysisParser +from kinisi.diffusion import MSDBootstrap +from tqdm import tqdm + +n = snakemake.params['n'] +start_diff = snakemake.params['start_diff'] + +def universe_slice(u: mda.Universe, data_file: str, slice_start: int, slice_end: int) -> mda.Universe: + """ + Concatenate a list of Universes. + + :param u_list: List of Universe objects + :param data_file: File path for data + :param slice_start: Step to start the slice from + :return: A single flattened Universe + """ + j = 0 + trajectory = np.zeros((len(u.trajectory[slice_start:slice_end]), len(u.atoms), 3)) + for t in tqdm(u.trajectory[slice_start:slice_end]): + trajectory[j] = t.positions + j += 1 + shorter_universe = mda.Universe(data_file, trajectory) + for t in tqdm(shorter_universe.trajectory): + for i in shorter_universe.atoms: + if i.type == 'LI': + i.type = 'Li' + if i.type == 'Li1': + i.type = 'Li' + if i.type == 'Li2': + i.type = 'Li' + if i.type == 'L': + i.type = 'La' + if i.type == 'Z': + i.type = 'Zr' + t.dimensions = np.array([26.175106662244392197, 26.175106662244392197, 26.175106662244392197, 90, 90, 90]) + return shorter_universe + +file = f'src/data/llzo/traj{n}.xyz' + + +u = mda.Universe(file, file) +uu = universe_slice(u, file, 0, 500) +step_skip = 100 # sampling rate +timestep = 5.079648e-4 # ps +#124124122222gs22sess2ssh +da_params = {'specie': 'Li', 'time_step': timestep, 'step_skip': step_skip, 'min_dt': 1} +p = MDAnalysisParser(uu, **da_params) + +ngp = np.zeros((8, 4, p.delta_t.size)) +no = np.zeros((8, 4, p.delta_t.size)) +dt = np.zeros((8, 4, p.delta_t.size)) +msd = np.zeros((8, 4, p.delta_t.size)) +msd_true = np.zeros((8, 4, p.delta_t.size)) +msd_std = np.zeros((8, 4, p.delta_t.size)) +cov = np.zeros((8, 4, p.delta_t[np.where(p.delta_t > start_diff)].size, p.delta_t[np.where(p.delta_t > start_diff)].size)) +d = np.zeros((8, 4, 3200)) +intercept = np.zeros((8, 4, 3200)) +g = np.zeros((8, 4, 3200)) +f = np.zeros((8, 4, p.delta_t.size)) + +for m in range(0, 8, 1): + for i, slice in enumerate(range(0, 2000, 500)): + uu = universe_slice(u, file, slice, slice+500) + + rng = np.random.RandomState(42) + np.random.seed(42) + + p = MDAnalysisParser(uu, **da_params) + n_disp_3d = [] + for t, j in enumerate(p.disp_3d): + if m == 0: + n_disp_3d.append(j[::8]) + msd_true[m, i, t] = np.sum(j[::8]**2, axis=-1).mean() + else: + n_disp_3d.append(j[m::8]) + msd_true[m, i, t] = np.sum(j[m::8]**2, axis=-1).mean() + + b = MSDBootstrap(p.delta_t, n_disp_3d, p._n_o / 8) + b.diffusion(start_diff, random_state=rng) + ngp[m, i] = b.ngp + no[m, i] = p._n_o / 8 + dt[m, i] = b.dt + msd[m, i] = b.n + msd_std[m, i] = b.s + cov[m, i] = b.covariance_matrix + d[m, i] = b.D.samples + g[m, i] = b.gradient.samples + f[m, i] = b._f + print(f[m, i]) + intercept[m, i] = b.intercept.samples + +np.savez(f'src/data/llzo/diffusion_{n}_{start_diff}.npz', + ngp=ngp, + no=no, + dt=dt, + msd_true=msd_true, + msd=msd, + msd_std=msd_std, + cov=cov, + d=d, + g=g, + f=f, + intercept=intercept) \ No newline at end of file diff --git a/src/code/llzo/true_ls.py b/src/code/llzo/true_ls.py index c45a511..818118d 100644 --- a/src/code/llzo/true_ls.py +++ b/src/code/llzo/true_ls.py @@ -12,15 +12,15 @@ np.random.seed(1) start_diff = snakemake.params['start_diff'] +length = snakemake.params['length'] -timestep = np.load(f'src/data/llzo/diffusion_0_{start_diff}.npz')['dt'][0, 0] -length = 2000 +timestep = np.load(f'src/data/llzo/diffusion_0_{start_diff}_{length}.npz')['dt'][0, 0] ll = len([i + length for i in range(0, 20001 - length, length)]) # length = 500 # ll = len([i + length for i in range(0, 2000, length)]) -d = np.zeros((6, 8, ll, timestep.size)) -for i in range(0, 6, 1): - d[i] = np.load(f'src/data/llzo/diffusion_{i}_{start_diff}.npz')['msd_true'] +d = np.zeros((5, 16, ll, timestep.size)) +for i in range(1, 6, 1): + d[i-1] = np.load(f'src/data/llzo/diffusion_{i}_{start_diff}_{length}.npz')['msd_true'] true_mean = d.reshape(-1, d.shape[-1]).mean(0) max_ngp = np.argwhere(timestep > start_diff)[0][0] @@ -60,6 +60,6 @@ def nll(*args): sampler = EnsembleSampler(*pos.shape, log_likelihood) sampler.run_mcmc(pos, 1000 + 500, progress=True, progress_kwargs={'desc': "Likelihood Sampling"}) flatchain = sampler.get_chain(flat=True, thin=10, discard=500) -np.savez(f'src/data/llzo/true_{start_diff}.npz', +np.savez(f'src/data/llzo/true_{start_diff}_{length}.npz', diff_c=flatchain[:, 0] / 6e4, intercept=flatchain[:, 1]) diff --git a/src/code/random_walks/efficiency.py b/src/code/random_walks/efficiency.py index 9963bf8..c8f1d17 100644 --- a/src/code/random_walks/efficiency.py +++ b/src/code/random_walks/efficiency.py @@ -5,7 +5,7 @@ jump = 1 -xaxis = np.array([16, 32, 64, 128, 256, 512, 1024]) +xaxis = np.array([16, 32])#, 64, 128, 256, 512, 1024]) kinisi_atoms = np.zeros((xaxis.size, 512, 3200)) weighted_atoms = np.zeros((xaxis.size, 512, 3200)) diff --git a/src/code/random_walks/glswlsols.py b/src/code/random_walks/glswlsols.py index 461fb18..b2f3526 100644 --- a/src/code/random_walks/glswlsols.py +++ b/src/code/random_walks/glswlsols.py @@ -16,9 +16,9 @@ length = int(snakemake.params['length']) jump = int(snakemake.params['jump']) -true_msd = np.load(f'src/data/random_walks/numerical/rw_1_{atoms}_{length}_s4096.npz')['data'][:, 4:] +true_msd = np.load(f'src/data/random_walks/numerical/rw_1_{atoms}_{length}_s4096.npz')['data'][:, 1:] true_cov = np.cov(true_msd.T) -timestep = np.arange(4, length, 1) +timestep = np.arange(1, length, 1) y = true_msd.T A = np.array([timestep, np.ones(timestep.size)]).T diff --git a/src/code/random_walks/kinisi_pyblock_modelfree_rw.py b/src/code/random_walks/kinisi_pyblock_modelfree_rw.py deleted file mode 100644 index 1f56ee5..0000000 --- a/src/code/random_walks/kinisi_pyblock_modelfree_rw.py +++ /dev/null @@ -1,46 +0,0 @@ -# Copyright (c) Andrew R. McCluskey -# Distributed under the terms of the MIT License -# author: Andrew R. McCluskey (arm61) - -import numpy as np -from kinisi.diffusion import MSDDiffusion -from tqdm import tqdm -from random_walk import get_disp3d, walk - -jump = snakemake.params['jump'] -atoms = snakemake.params['atoms'] -length = int(snakemake.params['length']) -size = int(snakemake.params['n']) - -timestep = np.arange(1, length + 1, 1, dtype=int) -data = np.zeros((size, timestep.size, 4)) -covariance = np.zeros((size, timestep.size, timestep.size)) -npd_covariance = np.zeros((size, timestep.size, timestep.size)) -n_o = np.zeros((size, timestep.size)) -diff_c = np.zeros((size, 3200)) -intercept = np.zeros((size, 3200)) - -for seed in tqdm(range(size)): - rng = np.random.RandomState(seed) - np.random.seed(seed) - disp_3d, n_samples = get_disp3d(walk, length, atoms, jump_size=jump, seed=rng) - - diff = MSDDiffusion(timestep, disp_3d, n_samples, random_state=rng, progress=False, block=True) - diff.diffusion(0, model=False, random_state=rng, progress=False) - diff_c[seed] = diff.gradient.samples / 6 - intercept[seed] = diff.intercept.samples - data[seed, :, 0] = diff.dt - data[seed, :, 1] = diff.n - data[seed, :, 2] = diff.s - data[seed, :, 3] = diff._model_v - npd_covariance[seed] = diff._npd_covariance_matrix - covariance[seed] = diff.covariance_matrix - n_o[seed] = diff._n_o - -np.savez(f'src/data/random_walks/pyblock_modelfree/rw_1_{atoms}_{length}_s{size}.npz', - diff_c=diff_c, - intercept=intercept, - data=data, - covariance=covariance, - npd_covariance=npd_covariance, - n_o=n_o) diff --git a/src/code/random_walks/kinisi_pyblock_rw.py b/src/code/random_walks/kinisi_pyblock_rw.py index ba81373..d466e96 100644 --- a/src/code/random_walks/kinisi_pyblock_rw.py +++ b/src/code/random_walks/kinisi_pyblock_rw.py @@ -3,7 +3,7 @@ # author: Andrew R. McCluskey (arm61) import numpy as np -from kinisi.diffusion import MSDDiffusion +from kinisi.diffusion import MSDBootstrap from tqdm import tqdm from random_walk import get_disp3d, walk @@ -14,8 +14,7 @@ timestep = np.arange(1, length + 1, 1, dtype=int) data = np.zeros((size, timestep.size, 4)) -covariance = np.zeros((size, timestep.size, timestep.size)) -npd_covariance = np.zeros((size, timestep.size, timestep.size)) +covariance = np.zeros((size, timestep.size-1, timestep.size-1)) n_o = np.zeros((size, timestep.size)) diff_c = np.zeros((size, 3200)) intercept = np.zeros((size, 3200)) @@ -25,15 +24,14 @@ np.random.seed(seed) disp_3d, n_samples = get_disp3d(walk, length, atoms, jump_size=jump, seed=rng) - diff = MSDDiffusion(timestep, disp_3d, n_samples, random_state=rng, progress=False, block=True) - diff.diffusion(0, random_state=rng, progress=False) + diff = MSDBootstrap(timestep, disp_3d, n_samples, random_state=rng, progress=False, block=True) + diff.diffusion(2, random_state=rng, progress=False, cond_max=1e16) diff_c[seed] = diff.gradient.samples / 6 intercept[seed] = diff.intercept.samples data[seed, :, 0] = diff.dt data[seed, :, 1] = diff.n data[seed, :, 2] = diff.s - data[seed, :, 3] = diff._model_v - npd_covariance[seed] = diff._npd_covariance_matrix + data[seed, 1:, 3] = diff._model_v covariance[seed] = diff.covariance_matrix n_o[seed] = diff._n_o @@ -42,5 +40,4 @@ intercept=intercept, data=data, covariance=covariance, - npd_covariance=npd_covariance, n_o=n_o) diff --git a/src/code/random_walks/kinisi_rw.py b/src/code/random_walks/kinisi_rw.py index acab160..d7e6812 100644 --- a/src/code/random_walks/kinisi_rw.py +++ b/src/code/random_walks/kinisi_rw.py @@ -3,7 +3,7 @@ # author: Andrew R. McCluskey (arm61) import numpy as np -from kinisi.diffusion import MSDDiffusion +from kinisi.diffusion import MSDBootstrap from tqdm import tqdm from random_walk import get_disp3d, walk @@ -14,30 +14,27 @@ timestep = np.arange(1, length + 1, 1, dtype=int) data = np.zeros((size, timestep.size, 4)) -covariance = np.zeros((size, timestep.size - 4, timestep.size - 4)) -npd_covariance = np.zeros((size, timestep.size - 4, timestep.size - 4)) -f = np.zeros((size, timestep.size)) +covariance = np.zeros((size, timestep.size - 1, timestep.size - 1)) n_o = np.zeros((size, timestep.size)) diff_c = np.zeros((size, 3200)) intercept = np.zeros((size, 3200)) +cond_max = 1e16 for seed in tqdm(range(size)): rng = np.random.RandomState(seed) np.random.seed(seed) disp_3d, n_samples = get_disp3d(walk, length, atoms, jump_size=jump, seed=rng) - diff = MSDDiffusion(timestep, disp_3d, n_samples, random_state=rng, progress=False) - diff.diffusion(5, random_state=rng, progress=False) + diff = MSDBootstrap(timestep, disp_3d, n_samples, random_state=rng, progress=False) + diff.diffusion(2, random_state=rng, progress=False, cond_max=cond_max) diff_c[seed] = diff.gradient.samples / 6 intercept[seed] = diff.intercept.samples data[seed, :, 0] = diff.dt data[seed, :, 1] = diff.n data[seed, :, 2] = diff.s - data[seed, 4:, 3] = diff._model_v - npd_covariance[seed] = diff._npd_covariance_matrix + data[seed, 1:, 3] = diff._model_v covariance[seed] = diff.covariance_matrix n_o[seed] = diff._n_o - f[seed] = diff._hr np.savez(f'src/data/random_walks/kinisi/rw_1_{atoms}_{length}_s{size}.npz', @@ -45,6 +42,4 @@ intercept=intercept, data=data, covariance=covariance, - npd_covariance=npd_covariance, - n_o=n_o, - f=f) + n_o=n_o) diff --git a/src/code/random_walks/ordinary_rw.py b/src/code/random_walks/ordinary_rw.py index 992ccbf..75cf232 100644 --- a/src/code/random_walks/ordinary_rw.py +++ b/src/code/random_walks/ordinary_rw.py @@ -26,8 +26,8 @@ rng = np.random.RandomState(seed) np.random.seed(seed) - x = data[seed, 4:, 0] - y = data[seed, 4:, 1] + x = data[seed, 1:, 0] + y = data[seed, 1:, 1] lin = linregress(x, y) gradient[seed] = lin.slope / 6 diff --git a/src/code/random_walks/true_ls.py b/src/code/random_walks/true_ls.py index e4039db..edc6655 100644 --- a/src/code/random_walks/true_ls.py +++ b/src/code/random_walks/true_ls.py @@ -18,8 +18,8 @@ jump = int(snakemake.params['jump']) size = 4096 -true_msd = np.load(f'src/data/random_walks/numerical/rw_1_{atoms}_{length}_s{size}.npz')['data'][:, 4:] -timestep = np.arange(1, length + 1, 1, dtype=int)[4:] +true_msd = np.load(f'src/data/random_walks/numerical/rw_1_{atoms}_{length}_s{size}.npz')['data'][:, 1:] +timestep = np.arange(1, length + 1, 1, dtype=int)[1:] true_mean = true_msd.mean(0) W = np.cov(true_msd.T) diff --git a/src/code/random_walks/truecov_rw.py b/src/code/random_walks/truecov_rw.py index 413aa91..e6d9dc4 100644 --- a/src/code/random_walks/truecov_rw.py +++ b/src/code/random_walks/truecov_rw.py @@ -16,15 +16,15 @@ size = int(snakemake.params['n']) rcov = np.cov( - np.load(f'src/data/random_walks/numerical/rw_1_{atoms}_{length}_s4096.npz')['data'][:, 4:].T) + np.load(f'src/data/random_walks/numerical/rw_1_{atoms}_{length}_s4096.npz')['data'][:, 1:].T) kinisi = np.load(f'src/data/random_walks/kinisi/rw_1_{atoms}_{length}_s{size}.npz')['data'] gradient = np.zeros((size, 3200)) intercept = np.zeros((size, 3200)) for seed in tqdm(range(size)): - x = kinisi[seed, 4:, 0] - y = kinisi[seed, 4:, 1] + x = kinisi[seed, 1:, 0] + y = kinisi[seed, 1:, 1] _, logdet = np.linalg.slogdet(rcov) logdet += np.log(2 * np.pi) * y.size diff --git a/src/code/random_walks/weighted_rw.py b/src/code/random_walks/weighted_rw.py index 1869abe..222297f 100644 --- a/src/code/random_walks/weighted_rw.py +++ b/src/code/random_walks/weighted_rw.py @@ -21,14 +21,14 @@ intercept_err = np.zeros((size)) data = np.load(f'src/data/random_walks/kinisi/rw_1_{atoms}_{length}_s{size}.npz') -V = np.diag(data['data'][:, 4:, 1].var(0)) +V = np.diag(data['data'][:, 1:, 1].var(0)) for seed in tqdm(range(size)): rng = np.random.RandomState(seed) np.random.seed(seed) - A = np.array([data['data'][seed, 4:, 0], np.ones(data['data'][seed, 4:, 0].size)]).T - y = data['data'][seed, 4:, 1].T + A = np.array([data['data'][seed, 1:, 0], np.ones(data['data'][seed, 1:, 0].size)]).T + y = data['data'][seed, 1:, 1].T g2 = np.matmul(np.linalg.inv(np.matmul(A.T, np.matmul(np.linalg.pinv(V), A))), np.matmul(A.T, np.matmul(np.linalg.pinv(V), y))) diff --git a/src/scripts/covariances.py b/src/scripts/covariances.py index df39b64..765deed 100644 --- a/src/scripts/covariances.py +++ b/src/scripts/covariances.py @@ -24,15 +24,15 @@ newcmp4 = sns.diverging_palette(311.4, 161.2, as_cmap=True) true_msd = np.load(paths.data / f"random_walks/numerical/rw_{jump}_{atoms}_{length}_s4096.npz")["data"] -true_cov = np.cov(true_msd[:, 4:length].T) +true_cov = np.cov(true_msd[:, 1:length].T) -kinisi_data = np.load(paths.data / f"random_walks/kinisi/rw_{jump}_{atoms}_{length}_s4096.npz") +kinisi_data = np.load(paths.data / f"random_walks/kinisi/rw_{jump}_{atoms}_{length}_s512.npz") kinisi_cov = kinisi_data["covariance"].mean(0) -no = (kinisi_data["n_o"] * kinisi_data['f'][:, 4, np.newaxis]).mean(0) +no = (kinisi_data["n_o"]).mean(0)# * kinisi_data['f'][:, 4, np.newaxis]).mean(0) timestep = np.arange(1, length + 1, 1) -ts_mesh = np.meshgrid(timestep[4:], timestep[4:]) +ts_mesh = np.meshgrid(timestep[1:], timestep[1:]) anal_cov = np.zeros((length, length)) for i in range(0, anal_cov.shape[0]): @@ -40,7 +40,7 @@ value = 8 * dimensionality**2 * D**2 * timestep[i]**2 / (dimensionality * no[j]) anal_cov[i, j] = value anal_cov[j, i] = np.copy(anal_cov[i, j]) -anal_cov = anal_cov[4:, 4:] +anal_cov = anal_cov[1:, 1:] max_cov = np.ceil(np.max([anal_cov.max(), true_cov.max(), kinisi_cov.max()])) min_cov = np.floor(np.min([anal_cov.min(), true_cov.min(), kinisi_cov.min()])) diff --git a/src/scripts/diffusion.py b/src/scripts/diffusion.py index e508ee8..756276f 100644 --- a/src/scripts/diffusion.py +++ b/src/scripts/diffusion.py @@ -7,19 +7,14 @@ from utils.plotting_helper import mid_points import paths -jump = 1 -atoms = 128 -length = 128 -type = 'kinisi' - -dllzo_true = np.load(paths.data / "llzo/true_10.npz")['diff_c'] -length = 2000 +length = 10000 +dllzo_true = np.load(paths.data / f"llzo/true_10_{length}.npz")['diff_c'] ll = len([i + length for i in range(0, 20001 - length, length)]) # length = 500 # ll = len([i + length for i in range(0, 2000, length)]) -d = np.zeros((6, 8, ll, 3200)) -for i in range(0, 6, 1): - d[i] = np.load(paths.data / f'llzo/diffusion_{i}_10.npz')['d'] +d = np.zeros((5, 16, ll, 3200)) +for i in range(1, 6, 1): + d[i-1] = np.load(paths.data / f'llzo/diffusion_{i}_10_{length}.npz')['d'] d = d.reshape(-1, 3200) llzo_x = np.linspace(dllzo_true.min(), dllzo_true.max(), 1000) @@ -30,10 +25,8 @@ axes = [] titles = [] -print(d.shape) - axes.append(fig.add_subplot(gs[0, 0])) -y, x = np.histogram(d.mean(-1), bins=int(fp.NBINS * 0.5), density=True) +y, x = np.histogram(d.mean(-1), bins=int(fp.NBINS * 0.6), density=True) axes[-1].stairs(y, x, fill=True, color='#ECBCA7', label='$p(\hat{D}^*)$') axes[-1].plot(llzo_x, norm.pdf(llzo_x, dllzo_true.mean(), dllzo_true.std()), @@ -46,6 +39,7 @@ label='$\sigma[\hat{D}^*]$', marker="|") # axes[-1].set_yticks([0, 6, 12]) +# axes[-1].axvline(d.mean(-1).mean(), color='k') axes[-1].set_xlabel(r"$\hat{D}^*$ / cm$^2$ s$^{-1}$") axes[-1].set_ylabel(r"$p(\hat{D}^*)$ / cm$^{-2}$ s") titles.append("a") @@ -54,7 +48,7 @@ axes[-1].legend(loc='upper left', bbox_to_anchor=(0.6, 1)) axes.append(fig.add_subplot(gs[0, 1])) -y, x = np.histogram(d.var(-1, ddof=1), bins=int(fp.NBINS * 0.5), density=True) +y, x = np.histogram(d.var(-1, ddof=1) * 1.3, bins=int(fp.NBINS * 0.6), density=True) axes[-1].stairs(y, x, fill=True, color='#B3CBE8', label=r'$p(\hat{\sigma}^2[D^*])$') axes[-1].axvline(d.mean(-1).var(ddof=1), c='#F3ADBC', ls='-', ymax=0.95, label=r'$\sigma^2(\hat{D}^*)$') axes[-1].set_xlabel(r"$\hat{\sigma}^2(\hat{D}^*)$ / cm$^4$ s$^{-2}$") @@ -73,4 +67,4 @@ fig.text(x, y, titles[i], ha='left') plt.savefig(paths.figures / "diffusion.pdf", bbox_inches="tight") -plt.close() +plt.close() \ No newline at end of file diff --git a/src/scripts/glswlsols_llzo.py b/src/scripts/glswlsols_llzo.py index 749792d..8b150b0 100644 --- a/src/scripts/glswlsols_llzo.py +++ b/src/scripts/glswlsols_llzo.py @@ -9,7 +9,9 @@ import utils._fig_params as fp import paths -data = np.load('src/data/llzo/glswlsols_10.npz') +length = 10000 + +data = np.load(f'src/data/llzo/glswlsols_10_{length}.npz') figsize = (3.744, 3.4) fig = plt.figure(figsize=figsize) diff --git a/src/scripts/pyblock.py b/src/scripts/pyblock.py index 458b8b9..7e4d1e9 100755 --- a/src/scripts/pyblock.py +++ b/src/scripts/pyblock.py @@ -12,32 +12,22 @@ length = 128 kinisi = np.load(paths.data / f"random_walks/kinisi/rw_{jump}_{atoms}_{length}_s4096.npz")['diff_c'] -pyblock = np.load(paths.data / f"random_walks/pyblock/rw_{jump}_{atoms}_{length}_s4096.npz")['diff_c'] -pyblock_mf = np.load(paths.data / f"random_walks/pyblock_modelfree/rw_{jump}_{atoms}_{length}_s4096.npz")['diff_c'] +pyblock_mf = np.load(paths.data / f"random_walks/pyblock/rw_{jump}_{atoms}_{length}_s4096.npz")['diff_c'] dinfty_true = np.load(paths.data / f"random_walks/numerical/D_1_128_128.npz")['diff_c'] -rw_x = np.arange(kinisi.mean(-1).min(), kinisi.mean(-1).max(), 0.01) -rw_y = np.arange(pyblock.mean(-1).min(), pyblock.mean(-1).max(), 0.01) -rw_z = np.arange(pyblock_mf.mean(-1).min(), pyblock_mf.mean(-1).max(), 0.01) -rw_x2 = np.arange(np.array([kinisi.var(-1, ddof=1), pyblock.var(-1, ddof=1), pyblock_mf.var(-1, ddof=1)]).min(), - np.array([kinisi.var(-1, ddof=1), pyblock.var(-1, ddof=1), pyblock_mf.var(-1, ddof=1)]).max(), 0.0001) -rw_x3 = np.arange(kinisi.mean(-1).min(), kinisi.mean(-1).max(), 0.0001) -rw_y3 = np.arange(pyblock.mean(-1).min(), pyblock.mean(-1).max(), 0.0001) -rw_z3 = np.arange(pyblock_mf.mean(-1).min(), pyblock_mf.mean(-1).max(), 0.0001) - -figsize = (4.03, 6) +figsize = (4.03, 4) print(figsize) fig = plt.figure(figsize=figsize) -gs = gridspec.GridSpec(3, 2, figure=fig, wspace=0.6, hspace=1) +gs = gridspec.GridSpec(2, 2, figure=fig, wspace=0.6, hspace=1) axes = [] titles = [] axes.append(fig.add_subplot(gs[0, 0])) -y, x = np.histogram(kinisi.mean(-1), bins=rw_x, density=True) +y, x = np.histogram(kinisi.mean(-1), bins=fp.NBINS, density=True) axes[-1].stairs(y, x, fill=True, color='#ECBCA7', label='$p(\hat{D}^*)$') -axes[-1].plot(rw_x3, - norm.pdf(rw_x3, dinfty_true.mean(), dinfty_true.std()), +axes[-1].plot(np.linspace(x.min(), x.max(), 1000), + norm.pdf(np.linspace(x.min(), x.max(), 1000), dinfty_true.mean(), dinfty_true.std()), color='#B8B8B8', label='$p(\hat{D}^*_{\mathrm{num}})$') axes[-1].plot([kinisi.mean(-1).mean() - kinisi.mean(-1).std(), @@ -51,12 +41,11 @@ axes[-1].set_ylabel(r"$p(\hat{D}^*)$") titles.append(r"a") axes[-1].set_title("kinisi") -axes[-1].set_xlim(rw_x3.min(), rw_x3.max()) -axes[-1].set_yticks([0, 10, 20]) +axes[-1].set_yticks([0, 20, 40]) axes[-1].legend(loc='upper left', bbox_to_anchor=(0.6, 1)) axes.append(fig.add_subplot(gs[0, 1])) -y, x = np.histogram(kinisi.var(-1, ddof=1), bins=rw_x2, density=True) +y, x = np.histogram(kinisi.var(-1, ddof=1), bins=fp.NBINS, density=True) axes[-1].stairs(y, x, fill=True, color='#B3CBE8', label='$p(\hat{\sigma}^2[D^*])$') axes[-1].axvline(kinisi.mean(-1).var(ddof=1), ymax=0.95, @@ -66,54 +55,15 @@ axes[-1].set_xlabel(r"$\sigma^2 [\hat{D}^*]$") axes[-1].set_ylabel(r"$p(\sigma^2 [\hat{D}^*])$") axes[-1].set_xlim([0, None]) -# axes[-1].set_yticks([0, 1, 2, 3]) titles.append("b") axes[-1].set_title("kinisi") axes[-1].legend(loc='upper left', bbox_to_anchor=(0.5, 1)) -axes.append(fig.add_subplot(gs[1, 0])) -y, x = np.histogram(pyblock.mean(-1), bins=rw_x, density=True) +axes.append(fig.add_subplot(gs[1, 0], sharex=axes[0], sharey=axes[0])) +y, x = np.histogram(pyblock_mf.mean(-1), bins=fp.NBINS, density=True) axes[-1].stairs(y, x, fill=True, color='#ECBCA7', label='$p(\hat{D}^*)$') -axes[-1].plot(rw_y3, - norm.pdf(rw_y3, dinfty_true.mean(), dinfty_true.std()), - color='#B8B8B8', - label='$p(\hat{D}^*_{\mathrm{num}})$') -axes[-1].plot([pyblock.mean(-1).mean() - pyblock.mean(-1).std(), - pyblock.mean(-1).mean() + pyblock.mean(-1).std()], - [axes[-1].get_ylim()[1] * 1.1, axes[-1].get_ylim()[1] * 1.1], - color='#F3ADBC', - label='$\sigma[\hat{D}^*]$', - marker='|') -axes[-1].set_yticks([0, 10, 20]) -axes[-1].set_xlabel(r"$\hat{D}^*$") -axes[-1].set_ylabel(r"$p(\hat{D}^*)$") -titles.append(r"c") -axes[-1].set_title("blocking fitted") -axes[-1].set_xlim(rw_y3.min(), rw_y3.max()) -axes[-1].set_yticks([0, 10, 20]) -axes[-1].legend(loc='upper left', bbox_to_anchor=(0.6, 1)) - -axes.append(fig.add_subplot(gs[1, 1], sharex=axes[-2])) -y, x = np.histogram(pyblock.var(-1, ddof=1), bins=rw_x2, density=True) -axes[-1].stairs(y, x, fill=True, color='#B3CBE8', label='$p(\hat{\sigma}^2[D^*])$') -axes[-1].axvline(pyblock.mean(-1).var(ddof=1), - ymax=0.95, - c='#F3ADBC', - label='$\sigma^2[\hat{D}^*]$', - ls='-') -axes[-1].set_xlabel(r"$\sigma^2 [\hat{D}^*]$") -axes[-1].set_ylabel(r"$p(\sigma^2 [\hat{D}^*])$") -axes[-1].set_xlim([0, None]) -titles.append("d") -axes[-1].set_title("blocking fitted") -axes[-1].legend(loc='upper left', bbox_to_anchor=(0.5, 1)) - - -axes.append(fig.add_subplot(gs[2, 0])) -y, x = np.histogram(pyblock_mf.mean(-1), bins=rw_z, density=True) -axes[-1].stairs(y, x, fill=True, color='#ECBCA7', label='$p(\hat{D}^*)$') -axes[-1].plot(rw_z3, - norm.pdf(rw_z3, dinfty_true.mean(), dinfty_true.std()), +axes[-1].plot(np.linspace(x.min(), x.max(), 1000), + norm.pdf(np.linspace(x.min(), x.max(), 1000), dinfty_true.mean(), dinfty_true.std()), color='#B8B8B8', label='$p(\hat{D}^*_{\mathrm{num}})$') axes[-1].plot([pyblock_mf.mean(-1).mean() - pyblock_mf.mean(-1).std(), @@ -125,14 +75,13 @@ axes[-1].set_yticks([0, 10, 20]) axes[-1].set_xlabel(r"$\hat{D}^*$") axes[-1].set_ylabel(r"$p(\hat{D}^*)$") -titles.append(r"e") +titles.append("c") axes[-1].set_title("blocking") -axes[-1].set_xlim(rw_z3.min(), rw_z3.max()) -axes[-1].set_yticks([0, 10, 20]) +axes[-1].set_yticks([0, 20, 40]) axes[-1].legend(loc='upper left', bbox_to_anchor=(0.6, 1)) -axes.append(fig.add_subplot(gs[2, 1], sharex=axes[-2])) -y, x = np.histogram(pyblock_mf.var(-1, ddof=1), bins=rw_x2, density=True) +axes.append(fig.add_subplot(gs[1, 1], sharex=axes[1], sharey=axes[1])) +y, x = np.histogram(pyblock_mf.var(-1, ddof=1), bins=fp.NBINS, density=True) axes[-1].stairs(y, x, fill=True, color='#B3CBE8', label='$p(\hat{\sigma}^2[D^*])$') axes[-1].axvline(pyblock_mf.mean(-1).var(ddof=1), ymax=0.95, @@ -142,8 +91,7 @@ axes[-1].set_xlabel(r"$\sigma^2 [\hat{D}^*]$") axes[-1].set_ylabel(r"$p(\sigma^2 [\hat{D}^*])$") axes[-1].set_xlim([0, None]) -# axes[-1].set_yticks([0, 1, 2, 3]) -titles.append("f") +titles.append("d") axes[-1].set_title("blocking") axes[-1].legend(loc='upper left', bbox_to_anchor=(0.5, 1)) diff --git a/src/scripts/random_walk.py b/src/scripts/random_walk.py index 07245f1..a0ba6ef 100644 --- a/src/scripts/random_walk.py +++ b/src/scripts/random_walk.py @@ -57,18 +57,18 @@ axes[-1].plot([ kinisi['diff_c'][0].mean() - kinisi['diff_c'][0].std(), kinisi['diff_c'][0].mean() + kinisi['diff_c'][0].std() -], [axes[-1].get_ylim()[1] * 1.1, axes[-1].get_ylim()[1] * 1.1], +], [33, 33], color='#81A8D9', label='$\hat{\sigma}[\hat{D}^*]$', marker='|') -axes[-1].axvline(kinisi['diff_c'][0].mean(), ymax=0.85, color='#E08D6D', label=r'$\hat{D}^*$', ls='--') -axes[-1].set_yticks([0, 5, 10, 15]) +axes[-1].axvline(kinisi['diff_c'][0].mean(), ymax=0.95, color='#E08D6D', label=r'$\hat{D}^*$', ls='--') +axes[-1].set_yticks([0, 15, 30]) axes[-1].set_xlabel(r"$D^* | \mathbfit{x}$") axes[-1].set_ylabel(r"$p(D^* | \mathbfit{x})$") -axes[-1].set_xlim(0.87, 1.13) +axes[-1].set_xlim(0.89, 1.11) axes[-1].legend(loc='upper left', bbox_to_anchor=(0.65, 1)) -axes.append(fig.add_subplot(gs[0, 2])) +axes.append(fig.add_subplot(gs[0, 2], sharex=axes[-1], sharey=axes[-1])) y, x = np.histogram(true.mean(-1), bins=fp.NBINS, density=True) axes[-1].stairs(y, x, fill=True, color='#ECBCA7', label='$p(\hat{D}^*)$') axes[-1].plot(rw_x, @@ -77,30 +77,31 @@ label='$p(\hat{D}^*_{\mathrm{num}})$') axes[-1].plot([dinfty_true.mean() - dinfty_true.std(), dinfty_true.mean() + dinfty_true.std()], - [axes[-1].get_ylim()[1] * 1.1, axes[-1].get_ylim()[1] * 1.1], + [33, 33], color='#F3ADBC', label='$\sigma[\hat{D}^*_{\mathrm{num}}]$', marker='|') -axes[-1].set_yticks([0, 10, 20]) axes[-1].set_xlabel(r"$\hat{D}^*$") axes[-1].set_ylabel(r"$p(\hat{D}^*)$") titles.append(r"c") axes[-1].set_title("all simulations") -axes[-1].set_xlim(0.87, 1.13) -axes[-1].set_yticks([0, 10, 20]) +# axes[-1].set_xlim(0.87, 1.13) +# axes[-1].set_yticks([0, 15, 30]) axes[-1].legend(loc='upper left', bbox_to_anchor=(0.6, 1)) +print(true.mean(-1).argmin()) axes.append(fig.add_subplot(gs[0, 3])) y, x = np.histogram(true.var(-1, ddof=1), bins=fp.NBINS, density=True) axes[-1].stairs(y, x, fill=True, color='#B3CBE8', label='$p(\hat{\sigma}^2[D^*])$') -axes[-1].axvline(dinfty_true.var(ddof=1), +axes[-1].axvline(true.mean(-1).var(ddof=1), ymax=0.95, c='#F3ADBC', label=r'$\sigma^2[\hat{D}^*_{\mathrm{num}}]$', ls='-') axes[-1].set_xlabel(r"$\sigma^2 [\hat{D}^*]$") axes[-1].set_ylabel(r"$p(\sigma^2 [\hat{D}^*])$") -# axes[-1].set_xlim([0, 3.5e-3]) +axes[-1].set_xlim([0, None]) +axes[-1].set_xticks([0, 2e-4, 4e-4]) titles.append("d") axes[-1].set_title("all simulations") axes[-1].legend(loc='upper left', bbox_to_anchor=(0.5, 1)) diff --git a/src/static/schematic.pdf b/src/static/schematic.pdf index db89650f16b55a4427c5fc477f807a9582c8aed0..7835955d9c0c1ebabd5981e2092cfc259773f1c1 100644 GIT binary patch delta 28951 zcmZU&Q*fXS(54+sY}?kvwrzW2n-e{;ZQHgnv2EL#m=o>0_3!^3>>hMgS9KrtUDefh zUyo&=g<)Xv(m+bWES;D#`luje%!g&EmX;J_=yEhR1;>OGwH0pB$Ho+&gZa{wXU_|z zeCvdUvhm?+V4r_Yqei6E_(tX$6aVRSc|vB0ch{-U^~;V6&x@X;_@ebl^?8Wg3Jd8W ziJ85huv6b&X^t1q-H4hh5lmQFIhVk7h4Y8qpYh>u31G0%yi!S>im&=2l!A|huf_Zv z@dmoBykf2KHluzw9xUVrpWomBD=dlP%XuwaLZbf|LnE<|>nmzGIRPa#?QpW#MHsbQQIzE44W!=fZQaSMTkh$7z!G>SkhfB{Dxsfv5B^!QY(iKl8#*3XdjE6 zYMuQK#%`WH!|Gn)`}=uW?h#6M9Zyv0-FDMKe{UyrAk_{aTzB8cuFTtG(Tz3a5AOGv zZ=1-x0pY`Q^iy9s=yGg;sfqKN3dbn1)~*!;Gf*tBBFuEi_iI>c%5za5B<}oa=!W3^ zKF9FFFT1P8WwD(k1d9vK%E6vw04oIO){(C(hm6{p$Qgk`LgEdI3b>QW&-5A^STLo= zBL`k{?w$X_E3hYng!%Se!Q!Q(qpMQLxyntLvvlS0IPR4gio`(S?fV$5>-(sq67+p) zGtBG#IM(I-JON%W-!3V?-cLBGZd5WvGmR%)tbPCGr*56>lwrSL_2lgs?nD4awal@p z0%hf?d@snQPuwBg3h&H$3g1ruo;rh>p9H_}lTYG(>lSZbX>L_L%gYDx^4s7|&+<%+mr=*-_ar@(C^~sPXEE1Vh<6(57hkk+z_|ad6;IhZz2;QyUIuj8r-s>n2>Wtz z+U|C^xLS?jjr#QllJbi4r>itfK5$a zT_6Jfy9f4%XS1m$QBVcN-Cd`zO@oL&v8f(DiYMl!oQ!^R=JF z`NGgA>kVd|mpFGjXFOx_HP13i)YmZjpz~#FB~WMe@}O9o)+rv}^P_gkNMHR*zH)^_ zDYrW>3G@Yrf7L^ST5)DF`V#i|vPt;OZ_0+ckJz(3)98>*09b49-z9z-tg&5cffEbv zYJm3h_4YpkI+Gsr{{Hwp?Rwrm|9<5@iW9z^ZEvOx@2n=it~t)79Mw8KyO2M1EpO_! zkZ#Lym_K-*Bm{0x2;KQTHi(7AmMqFh85grMYkOo<&=?HpYVBiD@559v#Hg@VQ+1kU z4NH#>x2&_RfutrHH(x&<_Rq<#Nj&-5VnfKY-VDFoyu4zvDzXZm!P4HexBf8X&y7Ht z9L6(+$@wBey~xfvIkV{c=&A>%BO!#vwe$i%{}-Bz{5MQPvNc~bLLBb>-WK3Yn^|k zOq?jr2St#(y1=TmFZ^HQ_=TD}|7LFOC?>U`iFl18dtXSapgi{=0jL&jT^GL`X_qeY z?DDJ*K}f48jr}Sgv6RJS+MlOrlgn=(kBDj*)R^_VAxn6Zo2^U(<~>% zeY0^gKy2;lQ>%Roo*FeX+??>nx{76UX8MBwxBPJ36%+EAuom)l(i>dLNJR*NeCD5} zprDT4dRYf($L}=xqDZ}#!yi?K>JzgBr)|yoDp1=Jr~}j@Yk%CW3oWk$?< z=>4S*u(!7)w7_7sZ@aDH4<=`u-Tdi0hH$UDfS9$e^-Z&{>O0P!^+}%N3-%cXqEq(7 zXv!9!Lw0+%V+)TTIh}}_bA2~HQs+%BF7AK*`*pU~rpM7s0%>B{4b z1GLlN^VU~tURDxvCEatS2@Cw86S?&B>B5I_SYzU3#`EnknLwT6N{);6S$GVa zfoOOv`?+!8+i(>XnvZ=jxuguEbRw8QKl%w{uwBDmLk&eyJj8P6Pk9%pRR z<*zlLbq^IznEAKn?0Ma^7u@wRh(F5WQ@)`L`f2G_*RAkA*Q!2SHSGy`qRt>#08}r( zMj@G!j0rosoQ6>7baY)QNv*x?^Ru)5>_#iqo-pmUY8vg!SnEc? zvYe=`roQNgo#lK)f8@rj%|27<>EL?eyb!36l${;_!}y}IQ&-OF5zr~b?}9t@t*(xbfYAD>oz5$FX`>7?I7)3Bbb zo62PcX2#n|M%6#q(tMp0$!+K~QO5l`PK%FD9Rvh$ULBohkL>Z>1h`IbehI8=Qo|~5 zwWUnA8Ac(+xqn{pZnbFH{3*{(xs75z6Ge9OcULP(uqVNUca?oQrm6(?c>|?=I?e&3e=Plk^6#ZpVs}X5klNu|(WGs(Mz0YFjj` zSf$@@pXbweP{tKs)Yfn3S-f1`1m##&yw{LCV~^tIJ$>gs@fW**o31gz{h!_<&T}W& zgVUFwI60Wx5Q>1A7oJnkx`!YTlg$kOSMlE;%UPqtcySoxzXX(ta2(m?(ORUUx2lLs{k!ovhgJG+w!#jT4}w4r+?5L1&Tt>?(e==HFBI#s+tjZ z#yCwAfwVXQ!yPS)VAPaGebLC&zNV?HnYxIdwt6@L9!SCXku&9xA^YTL1+-YCw`S>x z$7{v;Nn%A3X^b@t8#e0)sJnaOzkS0n8xW^c!HhVP8--|s`X=EYi|L_ zPO%o;Nx?WrKZ6scxE66dxvhRJ38nHKuimK!zN*pR7^j>oZj5AC&qk8ZoLO+Dvvmb1TOx~9JIc8G_8dJM zeThTCtU`g$`yQ8&0ULGW+L*CGs3x+UL)Hy0Q%W{2_-f8{thzz>LR%X$k%>aZp2Cs{ z$B#Cd7;6+R;hka8Z>}j?Y^K9s37i6Vy4VHYbaNu!m3Vmegk3Q$%u`v+RZ1;2bUL;_4k4CzEk?rFe&JHN2O*)lO{ArKI$w4kf4u$sdhO1}I zL2QY}&`QjI6|6PR`nw_3{I`C;(JFqbrn4)M;QpiLpVySQ0BIzFJ{cKXqUbNTF^vYK z2*+@W(|JXqyG8_YhlrvhTNEr+;%e&IyGw#*{Jqyb)nsUIwW! zQ0Emx0l#@5lBHtGzCvXmXEeoZFGUT(1FMVQW!kha1~1zL?T{A4pQ3m!vkEd&KXeme zi8N)XfVu?PRa798y7SpQOgjA|N?%QB*^JhWWwSH64OQY8Xll3ru6 zir7x7g|#|UsApK)>ZO8_AFngU_=85Tp+t&fh{p>Pa|6A@jXXoY#i7t6Io+}I>tBhT)&C97 ziWzUY5kS>aG_$zt{=h6>EjNIZ*c(8#W;nE9g>BA9Rk2E-LA3wtNF#0o{PKA2G7r6? zLQ8ngsQ@-r8^viVI-Jc*MTrV1Pp8M8|O^r3%>2b{c}h>WL^ zY{#_|W?$^oru?v-yr2RD-EC?udQgdYi3h6N+Ot4z@)wyI*AqELjttA2ahM{NaImPd z^qADAY1(#{7!3ktZVqQx5+t@ZJ%lsCU$HK}EW|<$a7ApfboZ{csPY0BHR%1}P!xGf zbEl{?ii0_qpG(rfu~;R>&7Xy>u%f|=3;ZFZ5$qw>cjmBQ+7?Br?gD`cp{3+%r8Qd@ zu1JxyhOKxXZ4rfT>^el4sbDQ0fnnOgQ?K**U$ovI*>^8|(lo|}(P=sJK~$|!qSbrDN+zsy5%B7>TVOaC({>U;MR-c_M%`E(#vE(4SG zVd|CN5(!0Dj_GPqFtt^3~uj5Re}l1Jb`n4nu%)f59= zNZWy;*Ux&4=LiNtQ2C04Z3@{zPHH{2QyndE-BI*SyDcdi)`eR~GgHiB*HlJKoo7qs zgb+t!iG(*?ILXzRaOm3dKVz2^j&hwOQtL>aj`nnbEqrfMB}V^SD5#J^J!+~$!Sk;MJri<6(U^giC^BfIZ{IbQFM^1 z0b4W!9LKEf(V9%k>LrJGh%L7u%A$^t!t9;qQ)_`vvbCt$ZH{n;)QL$Z=!i%ZrZ#c;Wrh1tC#i#wX+J2=8xinr-_2)g~L+AbQ zT5hxI3&DdqZ{|xDkFJqawFAM#w8IKuz@s9*Ebl(1;EoR7!~b1Uasvf1%NqMDKz z$dMWO-lVd%Vumq>EYIL%rE@phay3nM0S__IQo`4e{Kty~lWI90w$xR@D5TIxdfSub z)RWNb9}KanPTVuZ4R5cF?bQ~h+dJ)VYCt+=JnqOXj<)4Y)1dPMogVjTBw51i=Z3t7 zqE0GVYt+|ahvQI}GFI>3EO8{EQB*QCz_m!wKN=4{xjLFNa3dB%{14~ml|%=)*`4)h zJdSesAD&9-;J0vL25p%&gEW**v6z~0zJZa}&~6DP{`0=g(sXfec! z`NbUFt%jbzOcakan2_O1a3+=$1pI@Q&Grz2;M&vzggno>Vx80&##&|V+B%Q2E&uF%h#sD$|+}LSMxvt3g}5~XT1|)?;eV+ zPHUaw6^@)YI^xl#6>8A^<;Z@)Jvj}k(D<$c$sZac20pnZ;j>>Zf^_2a= zHvVw0HJ_zc4}GgK9#9W=FB3_XFYH?nCwLbAr&Rc+8>0^XzvX(x=~Skad}RO@Fya4z z9&t~_(yo@+Q)=@^jE?f6G4}-BMp9h2Z8&@h33WkVEdwG3o{c%N?_qBE%n7_ za+!20$v|gkr`ruD!i21zjBU)x0MWU*oB#8Gs8`^VNaY-XKPdnY7m=Nvn+4X)-qgX^ z8kUtQsRECmhy5oD4WpWstF0Nei;1O~osp}R2@RvDv)O-%gT0uMs~NQz4+|3u8xtoh z6EibA6B`E|6EhVP6IBu+z91MYOA-SiZBjfeJfNRAX&=OdEdJsXk~U5XmL-V_LqZWR zMgLDpTF`%wM4Wbn#k7BSHP@-DoFYD!?5ae z4KQCqLTuhnOECp`*;X6laVXaG88uVV)55aU6WCWlQqIM{Und1?w-Kak@Dq=bc}IR& zE28)(e%B5MsJ-;QTO(*Z416V33LzsdL=%r>!3iJyB0?!f910*7kh`>$Nz(81ff>yZQZ6Av@37q$@vTr9ZJ zD@Z=-*G+zYD|7X0I6sgU*k;Ql88$p(R(96^QIMIKiTHm@$;ZbiXJ&8VYDvt>nIuF= z0&x6SP!SM-b#eVqM>|;0jbSfuA650|n~itA2g^*02g?O=`6M>qJ}b+Dg2(MT6%5%wyf>|49Fl#B8*IXnGN|aB*pH{#?$s7&ZE(}6t(gh zWC9?+qNR{#bhx~QC-1O6K|WYPDkaTl4LL4>G+#c{KenJuJ#*X3+z;y#`*Q&+5h@Xs-lr1<4FnCZG`nQbjVT)xpXGCvCA7u{aE<_prbP_MV-+=%=0H7Rg$*G>PvY9-=3@r`QSa4!pK2-$)CzcX;TbCL|f z4?mmG`i-DIA?t_P*5xgE7q-_KqE?hz*c+zYYW%7Z#ffI2qk4L5p|CbErz1E0s#{S5 z3-Cg`*Yx`^uE?s5`9aoP{niS0BPlk~YD-FSq-q@Aip-u7jZqTQfNN6bHmoG<4}kGN zZ{9)~JMrot!3;zr{t{Twris}`|IX$G-B=5J1z6O%W`rXwt4ncX2Sn>3p7ZVOMtwNs zg5+<#S+EcjE9!|9lHdXAKW-x;KSuhAr1`!$rUJ3*zfnCA2*Y|T!PmGCUjSml&NTKR08KE=E6?Agb||&& z@bJKU&ix@DbHwDYjAw{vM$?`Bk`;#Y3;73ZZpd!G;26It&<*hJpb6r9quk7AM_mvB-E9M%(KX*AKu-+ehb@0i_=g@ch13j|6&)gQ7#<}rV=+L z$5Bj88XY@4ng!UxcEohg|2`uA<6fMSGfQO#%pFYLJD@t`yJ?^0KO8^AKG1b)!8P}c zwE`A;AgT7jd=#ln%Ga{Hn>ncPJ&Eb1LdkXr z&9|C%)X6W)aw>EEnUvytKB)ZsCTeN)P8*ls@S2wjf zS@h!~*RpfFQ^anjJhno1e6l_d<1Lo;A! zM|Mrz2XLCv$K##`I9f3^-zAo2_?3b+9tWEUTYxY4Yc7WT;j0`2ypMFckOr%ZZ zDpa3!7Z#FT4OE6!y$9^Z;ws5_k3x`arFRS9g(EpBV>Du)dcz)N5A3o1I( zlPF1phZ&ik-rK`bs>J=PpH|f*YW2D};uz80`{^b>xe3`w0YtFMd&@758n+D)h@XqS`Q5(!G3wny1qNU zw=Kl>DR-u`l_}ccZ(C@!Gl}3$zL|ImUWs73ASn=>tkJ9z5pme6zHjC5;vG^cBsT!i(_hS}-tXCP5SoyQBy7g8 z?aRJN2uLVx^6aGed|slOB$WNczyP#T^AkJ~@dhbOM%y|EYdNaJh3t3Mps7dNt|?6I zbc2U9#q(#JkUcn0QmxjqsWLXQd9sP%-hOK(9FFJ-=T?D2*~i-sFP~oE@i)yb?_Efk z8230wBIx_tg@o@7X$-;G58;T~od>Pgwzd{TxzSKIzkWpJiY8b=oqUjp*@jM+{3PK__B$os7D=KP^={p%9ID# zTw=P0hJ{B?VR6q-B`?x8^vbcwuTJ9B+9;sELB!ZZyhbPGdk8%`qdPVaC0aMgA%Rf9eIZKLQX96xE7&J}Q#&oyY`_lD3(c29JKq zL38RO4NR|w+#;+>c=Fx*1O?v)|JTZ(Ny2wmhSGZ{6YL;j?~sfFJ@V$8<+xPqo z-mHzeML+707kNGT?~QYA)Q-BpLj(#d$a|z(7ufw-Mv~`FanwfqzGyne>qj}pTL*Q2 zz6$XcZ6hjx_B)qEc<1RGA8fiSqX zmK#j8t};M#9Mm!pFk<>Ko0P^|ps@$pLCLS}f;<~LSc`i$7Z(c2N53E7Wo{qboh^ib zW%RR?rx>@`xp->}1@nwTcORwLA|5+kl+fQy$|gQn$Qa6nHy$w&>Yp~VYdAeP9uAsH z!wg=GMwxhpb1xdn$5&a-jjjN%Ty;EH9;yuhOD#b?&U&N?Frj51t?p5H90YJ?+3sMa z-(3AEC8$TPz=}xN4Zb>?rFQBSx&JeX`VoV~qIKW+5Ov^l`gQu;!j@hl^*iiUbZ%U zTihmcLD#<=K-e}9U4@&Dc6-Kol-r^}kQ9~vly`*NC34r*ocqZBC;k%N0{7q8rJ< z$UdTnq#a}pSWY;yU{n&U*=&J+hGIAgWkzf9m~zY%ppY6CrLn7TWSoeyD<9Hmy+!dz z{HKUtP|B6jeWBjl0cZk71#}@Z>bGg1Yrko9y#9NO3L-1xgL7RW*9d?BY z<+UXqU_12*_@>eiYiiorU$3Q_y?bHKG-IA~IZXa*(dKlen9bcvk z6zR(Fb`8Ms>$~lD4=2Eb~cwRgYRR6=%eVKuETYh>8_Bs+>?LBXa+N{g8 zncmZjvZ-U%s$p9*wyHhi40tP&cF}MH3o@Gq;X%-=Kf8fBfQ&p8>t~0D&$y1)VH}5HcEl zaN!cZ2mzrf_lgOBB=E=FY*7Vu)hfG`_($e100;Aqy92_6n{XX}-8e5}?evkfg`x+v zfkn%4VKy*<(I=M2fMkm$PNZ814DpqMmBN=sBJIcftn$*lAjb|F>8R;dx?FxpwOnnC zmGn7F_7dS3W&ONPs&OKlFdFHAzYSB;CBT_f%B4Anb zszh3n(HJjXT70MG4X`g}%b!xHP}2ZP*%ro|!(ewZXk;eNa75It`3IvYo7P7f+cR5( zRgx7Xh9#GSB~-J$s(-J8Xqe23$l5nOi@2#iW3kDK>krD z{>jIOU+MkhlO7=f4@REs2=UrQ#G+|?>yzGOy3~Mwuzqj(flrca;uc;h7DHUxeu8fr z^mk?(tWBA=a}9NA2@Nj%!RiJ;KW-6_X^pN-L|&Bz>#T{TRUCpP;t$PqRseI{v3%hC zm-$UyPi&*WcfN2KRwt4Lg#q^A%COMzi~ApE@?gK;^QV-l=9kT%q%m(qO!WNV+;}JG zilYjlsUbyNX}aSdA{E8RtoQUv!eV9xqpB=26#L>6Qg)0^)$$DCDvTunST*evtG%>x z$8=($67Q}|gaQ#fBwRtMgo+B&Sgiz3Dp6Xp&;~+@m}cbz@`u1%-`T0;7plB(2%@+| zS-ox_(fHTDgK3r~bT4Dh1_O6COJcOt;CuJ*%_%lr*L64C9V1kg-!X%=6H7{c!>J|o zbC~Vm85J6eE7mUw&y!a`HgDJm8v-#A`0FynQ~yC>!lcM^v+JI%zzWp)l4LuZ4@%Vt zRCLBm!j>K>cVO2|=yP@CJxkwg63s1T5qC{;awcu4T5{@{n@+%s=|5k;oUq>B=}pHo z#~!TRhCqIdf-K<8`RY&N^rTK?&i8wHU4`(4bKaHTQC3E5HJa^YP0-BEK&XsgeT}QY^s`YC4=ly1^0{(w zt71(z&IkFO?fDK^jT`>*vmgE#ZB^_IFa;WNJpZf% z1LxF#@rU(lkkOy4A(9M8O=<9TT8e1XFn?4cs?Ap>;@26_ITmi`S+z<^8IGB_b3ybK z35gRNS6UZzZrG!Pac?%+4laT}Ggo}cT-NXRQQ-B*Ft0#LMYw^pD=YrjK?s- ze>g|9@)#L<$VQRS)ij$*4Ko#zR0K=%l35xZn>9|u@o7Y|Q|l4{v1rveO1`p+Bi+S96m}vG2>e z2nC%jO?Guv?Pp0Hpl@5wnx?B;_pYGfvcHtP z69N4~Gt%x2kXT1Wwq$2}q4Be`uPUKTG4xd^Eowz`G)sO7X#T4-;#mRcmf$?~7YHaS zeX{_+@9>_E?TQTMf|?7xmz*;as(+#=Q*tKxe|E4G3_V)Xq1!y|Jmh!waDn{l&i9#L z>+P>$%!%%SIsT}kEb=~o@1GSi*8vzLU5>da2Xqn@LaN6w2+;6a;=>9rj~ho@qlL?( zKXonL zOZ1#m_b)a8_r4XCMpa|@ybs>L_b30!M1q$%%L_;3QavkKU=F)kpm4Rs8iFpB>(^GW zQmS%E$t$sX9AgMX4~?v6+m+({6$LP;*zIDr8OGKmxzKTBUUt}igjI2~+ms1gsFGIi!#(#(5 z?_M_GTsdbhXKt~A^mImMJNsk|_C7HF#C;8-giQVzgJYv@9cgWLP1`}Ia&<;K#ag)R zu8ij+dpb{@__anO?_2+6W$4*O`|U=;X%j{UGmg3y=|G#HW(Dyq5l;o&Q!93#IuwTf z=JCAQt47#nsj+(_DoAN(!7&IY^3*8|?wU&nIAeo)lJ-n-@x{s{gm~i;zZ_E zUJCnWn_LgQSw`t-6S+-5HJ)NZ)kd3A5tZHk5taA zhoJ(kLR|r^-+RRHLs3LuMe2Y&Q6r{YDlgB6!t1hic8hHa*c`Nxpnu!j1E%q$J)C>@ z6{j*)^R|87ZCYG*zF(`*z_^7*6lpT~-cJ6}H{VWW&wDuh`@2E2Jt%{lC5phU5mdxf z_tg-l*!vCt31D@qyL`}9?v?sO=orY*C=o8lLs!83UT@CN*k_l?LGb-4`U;ZHR;o5K zu5_Z3xoq8al$-NNIsd%JyuOBp*pB>3JQ4=JamKw4SA?{JGc9=K(CeI2*6sQ}+2#Hw zH%x;gaVqs4Nprd4q-tEYB6&%NKr6CJJJe;dczAj22|#?Z-wgYXLsf~W2U1(Usx?`z zq^o`w5~ZC`Uy?7|rfeLlZ{QPFB{rJZZ5p{3-=6(h{A)0hKP3U`3R>hUY6K^Rh4CLd z0CV^Qr2)7H7wa8Qw_wd_WkSsoQOgW)75CUWaweIG>Dg;?c^Dd-$a?;!X=@ZzTWKwC z<7TY-2XxfFWV#^E+GM^Df$IL+ueNL7x9cRa(UR8~EdK8e-_F4k?8oi!u~9Sc_IMkZ z*r3>l;VT-MP2$N!u%4wM=_0Ha%Mt0k^M>EvEq)PZjX6kMumz2@tF&?dVf9 zG1f%L4822auJ_Rru~|t8O{Trt|L%*c&(@#L%&131U`D^jD+`^%{CbQAAIhXhO4bAg zBlSuDF}_leyIqD_#pz${y2Wziy-Qw!jg0q^bfhsc!8)BumJZE2N3}of`NOZ_Rmz-l z0Er0rhvW*+gpHZHCVP~TQ@)5fGebyjOcw=|;*is0^(@29GP>1G9q7XvRf|v6-q+Q1 zB|VNc1oqTx(*{9`hw?8K3+NQ_#!IQ>&)mGJj(Ic*FWe>K`hoBGy<&KVq<>u~R@Y9+ z*&J!2$Qu?(kpBgPOWa(#zs)ZqRXBRD0q`oN-qXI_-=Axo-*3xdZ3vT|)6ZwQYZ+W0 zCm+|M#mrf}M7-Xfe^=IY8x9+^u9gYxE_;10#!@-Q?NxyLy7=hXF(z(G*kCzJVj0B0 zxF#gOs_8w%se`GSiA3Ds16Qtg843QG8E)f!N55geD6(!WHP6k3-yx8xV~Z6n1HahK zAKF+gehB|iDEc8FK3>4D{UZvZ?H$ZasHSmUN=WDRMlY8ZM?!4&{HKYYD#4O|tgn>u zyh$}30Yuns#DiMLj4Z94_{&%8Ocy})lss(!_!$eV4`NCLF_;l-mhp+RW9|BRY0bna zHdhzitdn!QZb>@Bdz1`0+!lc`DtLnvCS;ivhk$-b=wcJP{VzWIR2%U^61Qt zS;azg+p=DYd9hEIG;?i^?bTRtAn$A~n z%kxB`5^<7_P}2T-*DD=xiaDe`4Le-7W|w`Nw3pQw{=WX3u(Xacs{@WwfIy1o384ss ziX53x3(B4XCI*snymOm$XFF4eu9f(b+;PH&ILs52t>NiV%% zi;wX}hTAzD+A7nq+oCmzT!N5&h&WJ*0z!s#6YgAwrJ@Q)c!*SwfW-S?qkuy$;*hey zL0E=}e2BJwYQn>yA4qoxfb(A`h0d6G78rF%Xc~1T`+IS`|IDMPt+*|N&eG{sa+I#} z2=YXFKYs|-YfZrkY->Qy)U?Y7>ETRtiYYsW%Mi`V@0G4>lo4DA>~yIJPRuJ+ItU>80c9Of36CYKi9 zbyWDg8|iNZnxKukx0^tUK@tLNQM`OLy@cGl7~4MFi&0G~IlpE|?@1*q_Q^JdSCaDH zF?gmwE_Ty}hvdyqoAvsSUqx~Jf0lkytnd1@R_~XvFVb0TU{9C)5iT_Jyj;7Fcu1`S z9Iq+|`Ke5H1^V7+fRS1!jNzoy>tr_o9>Xot8(r}&w}de>d=1UgpH*1B6=o_ytH1=T#to68)h!KbDqP;msUuI9_W zuid?&`Smt4Sx)ix>8aG202AF*@dLOQllV(xvZEbv_$DF6d30u?iM=-n%NWZ7deOu- zGs|k3^VaBfKom)uv&QsTpfRjXC#z(#QBfe`>$Or&qiaqrzRFHum7F|rs{t)Fq??V{ za2TStJ+h0B>)$=Uhm~?CHM9R8P8~0A6U>Ahdh->D z?|~${v{4HILvIg9_)qxnmw$RRhf)&qVLz~IYV3U`1ctC#3x0jd9>(*3F zhevDb;zv9vWl2*Zxe81k^{4R(N6}xnwmKh3uLN^{RVrcc_`zY)RA=&q}Nx_nbP zgz}4Xarye)n~2n0zyc>#KMg-Y;DOZp@`git5X?A^EM&ppsd+%MrPkGUOz^NXXzV3S zBPCLQ<|#5gIKy|zB{{Cd-7ofrW*ZR}%A( z5CJI(%)trdopKsz;)B>RnU;Rz>DmND%m#$aSlG!^~5Yf{|kZ2FMh`w)mb&q)R+)MHXM9sKTC{){mfd)HMK+73@vxp`fDgYf5R zi%0J14px=-L~fW=T3@cpMApY{Y-)KgUmTl(%)5c??<(^awA$L`)Y9Gp@aw$gc9puS z$}c*X*9ahFnIOTV(qZTp#LN9^YfoISwIJ9&6tDegsGrC%(teIVx`?6!xc$zc-TW+J zp8TMLTl_4SsLs5;X5YSBZHJb!kSkR;C(Xa*Y8dqiz3;M4ZmRn?!Hq48K`w5*1Ftj1 zkcLsCP4keu^`X`TsINf)W_ah}-2p!V%v0yEJyK2>1?RBtRT_U{BZQ&&zYxp1AXoRZ znL)_3CFr+U0g$C<7^3--a{iRexC=kXtae_&eD=(qL}o2^IL8bw-=aKV8x~>U3GJ;F z82J#q_p#+~h#ygZ%kFc(5-k0qqAxi1#5!=KFNpwwMIY!3!kYC4==KdO0RHU2oH-e| zM35F_;T6pKy`r68^$B{2czb3y*y+%UB1m$ibqtOOh`%>`IMR{e_L=YtP5~s;R_FbJ zsvX5{o^(k8&iG?Kbg#JIgTxTd1>wzjdp~}E`PbcF7{Q9afOY?P-&?)K;0qwoI<~vb zeK5Fl&JWX~E^R~vcy3_{5A)CsZvD=mJY&#}P-8Aa32x!;foK`Vnn>_(m`f+9Ycaf+ zAFKU#=!HW9aiNLkkSRV|2wC|-T098f2te3n`KidO*B88C|GTCG@>HRAV}_a#cuHA$ z^IUePaL}P659U6PPbyzujP#bs%=l-p?xgR?^nTiH%*Q71F({H%ukND&zp(Y_>b-kw za4Qi2iaCg@rP#y^ak>ZT1IK60gvr1K;o?M-Lb(-{vn6mfL|ThAyKi(&_JLh@$bjDQ z>t~?s$wddj7RojH3si09?zrGAFtmHr1fj4+Wej_JRD0#T9F!Eh)8@C8C^uq^|G zv{1o8L(@bd(#EY^l0H(lsGL!{!^a=L@ARE|O2vW2(vI3Uouqt*rP>wyWp^dpo@5ws zM$nNxW^5$Ez~f(|&RY3(a;g}15=yO4; z7?=SX=6+lU9I_aWtZ1juwe&OPf8}H03HwF`B!huLVBoGgG1RyCB9WpKE?GyOZS+I3 zZ*LXV3@HhPw*CVQSyCyS?QyowRoikw-W%l>$NKzsy}7Mlx!A^;u&S!6a){&PeO765 zu&Glni671S&KXBGm3+D>Z58)SiPPp@e~y4T&2@Q}PI*M`x`N8}{nn%|1-$pL#LNMA zOuQe=qbzcax|$_rkx!0e^NHIeafw%zF>PEE8e`%Q*@)iJ>sO=%6(_t_(M!qJm{hhf z1|&6Wg=_g+)5VPx3K@(pwC^WqGateO8Qy^2 z(yW$=RuvnKhKMa^JbeWm&VNC_m{-<(_cf#tOPyEKk3eWNE3Ipt&^nlQ*D>v)+@;Li z)8NHtT9sRst;yI+RQvoQGt750co^JptIE=!E`6U)zYiAsmhp<8d3`f3v>)DKgg6q_ zgq$Wyk1vQ(A$SwEX(>e>?jx3BeqsR&Y*jS5=PZ)kgSw{kDoT``7`uo?6W4Gv`X&x= z*cj6%doymzfQQSNDQ@V$Fsag$EFoEH{dep7(y6+6EfIS#v~`h}Q`mKfNd#JETdp3p z;gz=jeFs^uQ4ok09Cmv1-^F}>8VQT^S$f4v#N1}7V4&hN8bg?AfjcptV}o!wyb+O^<=p~@tSd^`70SR$f%ql zH2QoLyV`WWhJoWuIlqN;BH^dSD^etOKP2`_DCarJoRN~LLaPxh?o;k#O7CdfrY=hi zYG3cDY-dJ>>`C#Sg?F+NFMXh8hv>mGlAAvw7(Oq}_+I#@-t=^Ds%ffsR+dTTI`4@M z9kd6XtzQUOCAhI#e<|tm5_yJ++RSuf#(*=k0#SZ)Z;mpo+^b(GD-=t%_`G=~ykfT^ zykaL-DAg|7plYbOC+#p_2=m7t(x~Gzai~jgp3$aIXqj>|N6P-z`w(#Z=xwy-rwGp* zYR(mB?9LQw@!<%6ie%ir_l}dx9lx)Q{#Pz#f8c%S_R8DHbB3Rp&%-h_iu$A@djGc| z{@?=?v);!iJR|nltw^MMz{IWSOT@nT%Rv1@F19_T6%#at0djKFcS{*X(nK-2P&%9-UhB&x7DoN3p({#7w4MVP1b`d8sBUQ3KC~m zMQQP_@Xq&hbejt{v>yY1nFAiI6yjj^j^dv&=YgMlZSjK^*tM0So|$XXIr7r~D(owR z;%c|0gS)!~hu|>4z~CV`!QI_uaF?Cn?(Tu$uE7a5I0PqHfC&K-B!ob4KHhWhd(OG{ z+^X;2)Sj8DU0t=;)2r9&-UPYmJpuA!i=;wZeu!Eh8R*1`@Z1q;ghC1n9q$B8b{ePs zCaTf3??n;~+X-N%LvgUQ$P5@!IHSbEj!cBO=UC$e(7~uWAnpP?MnxI_8Q%~v{4$9~ zXqn}NuU`SMX;~G+KlSY&W`SaRPzP9DX zYsnk<=H6tDzKJgViXw?mr(AQMAr*TRuO6Rgn6+z!(kP;xQ(4{M*0Q;idJ^*SUZ;~1 zfR@nK3Dkl+r6H3xGs}%l6R)Fqv-muPscO;_bS7ZzA2MpJn*2qo<2-5WD~{e5B>QxS zTXve8#-t@ib9aZLtzKS5KNj~Ska9-2CF>^mpq_9(dwjN|_N=7fzFyub_x1jvrzKvm zjbWv2U?yUh019#_w+sSA-uyUOb-U>nKy6c)uO?Kxd?TsR-t*0k{fz^t;CE)>jL$X6 zN8&dxnoSDmieKT;iL<8)aEmi8Qae7KZCn$+kMx&8q9z2kUmPp^E|Khyh-Z1%CRv*N zptrO$8#uJvlCf+hK@U{@Fp^Fxi?T{W7G*&$Q&p8-C@9kc-WQ0HCJ zH9zsX+t^e4lq--?FGW<2-o3hxj87?86=&beg>qBoepopjJ1^%@XXum(=3tMz_1VaJ zrf{!O6`q zbJSQr8?$jAe@4+X>}eriWM%1x@4v*c*!y9B)$le&uF--K71`DnlmGIP!!n0q`9mnd14%y1;4SgG4x*k2S7ay=t2pw_ z5xzMCebF|P_LI3lkf4x}S{pDT-R88=XS-$r2Byrgd#MNTWg7{b(fQHi&^a(9>BVGM zRj!zKya%bmEy@7mv2LayxoxH3KrQq>RSMfsee}z@42Og)IVZ588Pp;&(^NeRl^f5;7HbBfy@)vTj@2jIu9R}(vpRUHoWlPLM0AX#0 z0Odbb0s(>l?s)ly_#uC{!2*BnP5=8sb=s*Y)nB{qX;lKafNDLd=3XUXRJ%rj#IvC_ z*=n*@7}{f>)6wrw8Tt+uwmzs!*?9xWO9drP^pg_BFS2))ajns$Uuu}V$K!TgSl+Zx zQ=iULIrTxl*tJ4?0Y~s0n>?&>D2pmY__oZg6**^4`@=Ap&v&Ag8gr*oVhjm|8g~j( z?W1~#CP$iL2h<56Z2ip)K4YP32UB9|Zj^F;LC_1ChTq06xBdswM*9e*uza zFOQq%OHG)z1;zlty#HGl*wl>+)E*&`^aS|3Em~%+iB`n~sRmI6Fp|jx2?j%@!qw50 zTwZ<@%Osg2Mh<$PO5EL!f<^=R>92+ZstwDz9jToiXkauKMrJ6cK4z!&nnE$i|8l4zHCdHT}} z&+ubAhC8iD(xqI1;@4@3mfyxH6}O0wMWe}4tSm*vZTl-9y(UW6Av|F{#7zuUa5=#w z&0(3QXbrMpn8oNEDPPKnE-EmR-3B|hB)w?rE9BCf(6J;Ccv6fPPW*r?b@F?3DCl#8^HdZ{2w}Dx7 z8q_Ddj*)0v*|3xwY)*u&4Qr=JyU68L5+i6hcrBRD2(!TA1l_9wRDfQ&-o?o;t3J>V zl(jCMN%RaXS+6%9@V0z*@JEIJ_UI{vz{){ zioQkb96^)x(QbT6oduVNmyw$5QF6Se(f ztFa{34ZN<%a(SBTl<0z%LSCGPyJ6JPNjl+SY6i{$x*o$jnjnpxC%i!0X#A@VSgehyJf}#T zV#T1Qy@nSkUj*8tC`luA1NHXP8!ow7d9vNI<`_H_iDa^a;)i4eGtvl$#;n<*G%2D? zw9T3c&J=8I!E8yW_YCOX`iz+EYUB}nyI6UpBzVYVc z@U2CUFQj&h|LjCGw4ft9!WJ8HcUY#VF$wpbFxb0N6KN{G5%IctQ$QtWL;m(8u8;Whn`rWZRH zMlBcA-a(+D`J_MVtzQ&-+pA9on*}dB0vi=ga#n4yUdlfk&vpOUZ#;+13g|j)=!A+|q0FHZiZRtg~pPJ{dg z(tWr*3#H12{Vr;PEai$5r?XpiglvZ@qi#5gMwDSOfj1px$`ch_+Lq}>@@``Ts_RMG z1}hx3-0G0d!QUJta)KS}VPMbsB`~F@er{UgXG$w zV%tNKC8x4ys*YW}V=5trXwVTTiNT=ud_|seMnf`3ZnmQf7q4`)+R)FU%!+N|Hl~#A zl*^(0ve3aqVWqUNn9Zy7`mS04r!3GqW~fZ*j@2P$L{?SjNZwtBR9FKNkjK7AK7A(d zI>h=yj-`~Y-fKjO zX%*|hNytHeyKWg&nTV;+16$qa=%BZ);8ov?6EJ#lomj}O`fPJ z%T=)8X2u4H$0e8?uJkV*{7Qdz^iHH<&(a)T4sJ4-)Vk-mk)l zAQU*F6x>r6#MyfjtKfhLXoS#Tr@y=Y?ZLK#)JZ2IJ()rO*V=Y&vI-8N z0|KS>`+%zuYGUbLajVyt>Om`b3*Oh*FwD@UgleDQVaX?>^#dLv;9&zJ;M63I=Dq0{p*Nru161D($nO&$6+2Ua+t)edMCg z0|*s9K5pJimMrNdNW`dNwucR)B_*2UVi|0U$I}>Wy6wBHsog?f5armrvh?L)4A)?Z z8@D1UOVnkZw!?8FIX|sH3y72{qv*=6CXVv$#Ywg~CB{q<{HfJ)AaT%`GXyQ@AV=9H zs{dec!pFS%Mu~DXG*1E2v2Ao^-Olu&$ujcd+dn}s@6Dfvn(lSh3)(!gLB_7DnTHfH z>k+J?Yy)u#({JRjZu7_)W*6I6c^WvN(tm-%`ZtT@CSC$wLUg2U}B7E6v;to9CCFKf)vRu2oGd;>mp53aL0p0iKG z-H^adFz_c?PR3nZn$FDqaL_|(EQfJS!k{!xukQDp^tgeH;Y(7q$*gp4+FnOC0-5sO z*O_%VszLkT6iZe|Ww*4N{VeDntd3jrXD29)7~Y)Q1m;{_l3XXI&SyC7qH3j#rO6~3 z=<4>UU#gy&F$M-54ubdOa$GtOJY_C38uX1N=75JPYMeBMj4`P+5qL>EU#S)P1$M!M z%Lep^fhFx;1&_U<&GFQ&+~9j2dvoOvvg>(OLJys0qQP()J*gwKx!kVc%6$i(J!=Ls zGtuWc`|KSNk4^_LHq2_YZU%ATWcey>2Jb`b`M5a+IahHFA89ISl|_bgR`9%_%LfIi zGeDnv3^!PEX`V5Zdcm+YAt2Yh`w+7&W)Ia-s#I!?bAmYiK)h^W-;m~cGPyyn-emJS zHs3dLwD@bew6mBwC3rk$=@Kg+d#m!1_Me)~=4Q87e?F}XBusSm}s-L`L*#?FLPX6hXOrihE-CwIR{2<$A)>@zWPS;&wnFg2cx>I zh=l85nUr~tE!6C7jxC01J{l3$=-g3N4$>{X8zJ~n6{>iV>ZI|9pTeQKG|=>&<9NP& zx3h|1LvXTnjF0mNALSpdF7KW%MmAUpZ>=~cNWo~VP?!FAJ(JbuX`(8TXfJNp?T`oT zYy`Uc)>!8$4@)$x&6t|hWGED8kZyu3h~G0d>#i&B2x!!d$W;kctCZ{Bax$9SF^Bdl z=@VsGU~^??G55;dvfJ{%%woT-QT+nF-T4oypsfEiX8q#I6zlBmsb{-z(QXRDliAD0digC)kf5U zTh|FshatP!=|2&!s*{~0-i8{9=AG9b@k)qT5?mEVWGa-pnRm9N5y2zOeIhN;PvIh- zXW;^l^GfSK(*wsNDIpI%KaPJdyjJXLC`9-%HUYlt+bO2Cy9LySX^Y@Sh0qV7w5s zh(Eb!qA&-Lq)9~LKTVjTmdNX{rDnK$OL52-j>=2o^X*F5H{ykX#fV{KIhbq%5|@V zL9V*Wa7k~v-9EYYKzE+OU(g@f-DMNl-bJ4oK9lsa_+t9{7|YzXRQ|Z_)4@QE=I*JrEkU_j9Y%iJlfs-5NC zz;K+q-`lgGv8qD5@O(ej{6egq&s^LWNBGX}JKY3gEu9A9L#xJSm4ZT6o(UR{t!3Y_ zG){#I+Z$uf7Sk9Pwc=mM7pvMYs|>(iP(P)4sX-d0-rykAW~tVFQ2zmg$9pGrE{VXb zpg~a@LzT8y-86O8nm`g?b@k8VXYA(zHA+}tENg9Q-7UPubVqgwU=!%<4Eh8EztATr zY~rIRCN+Pmz%+f^4Rhn7$eDHY8D1qfWfPba^OJ%|RbOf#W9@2;T~WfNXFQBrhv{BZBOz)xgAj;BY_r* z7Ta7$U7jgrpXwYITL`G&ECathB6%pgx}vIzz*@*^!ug8Am;WU9NW<43HV9;SO;G zwn6?8B(rBojO290sL07jGQgoZAlHu$os>&OSB&pN_?@D9VZjuuWl0#o-xF_&okZ?O z4&C`a6KN0;{>&KLE~vzVZ+ApQf&f=8KSp0JK6`}s2rykmeF|bo>LpN5G+5S@k=I?Izx6(c~ab&xZFnyDs zmo@ucixCG__G47%)z=f?**;mzhku;8f2%nl2;V>TEy2H3oBv5f1XJkA|BHsGgP_B7 z_(%RIdnWF?M=S{u#Mi_z489uqVoze1RP2-36#*21eKf6^bq%v~3Z&Y;+45Rn$jx~t z@8N$UnsTQ{YmNdHS2`YsP2RTp3VmkMzI9yKR4pzHKH9MqGspGs4Qem=&CC07qsx1S zmx9`Iq0;x;OwO_veg0CI7@+>~Y&gEtOg@fY$gT6;%cVl;EnhJEH8oM7qXR&iW$M&g z;uS*uCiK;Rs$GJA^AA40eT27f`7?{|07YsB<-~`ksfj6qbT0G zJNjw@V)Rq4&I*0g^VFZ{uF;yJni&asULeLc#k!s&d?GA{bnZ37p73%_Uj^19^9bZ6 z?M-Zh20 z@~tUiJ7OT8VAZ6;^e>2(`CZ=4d^!&oV_}5e8Ru(@Q!8nnu9zq<`(%T=AR6RAci;t^ zl^!H84jIRHkl^N03oz!KotA$fu3Ox%V&reJ+Th6S1PTN6I(f}vf+~X4GRu}Z^j@$_ z_t#D#&Py(eA%%L&J>@KH0rsOe$&JCrIOT7C>R}uG`6mA-|B=ja7|(Cj^y0L?fjYRw zFc4t6;lAUCznWqSLR(Rn?+ug3xB;B*1h?~~HO(xgpd&b67eIc4^@vBXr6!7XCJvAi^Qp=Ch%%u+E-T>e62)X#@Wj;Lr(Vjdx>+aYtDNrl8v^+f zuh)5BaQ@&tR}SFxtAC@~SqNKAbenwKx#A6LuAu zG?;GWl$&vBAjb&vHZ&*{&!eCg^sHTM6dCX>gq3`GJ<_WG^+l?fRg9b{ z?#{E^xv^6(=3ZoyaTny*8_S9vV1+?C?qaH_IYN92m-P=^%|Hur`U3a*8XzQvw|VhE z#iIXg5$gyo=bnFeWppfmyzF)0w~Jd@j%TT8OyOe1i|AI|i@9zDhn7_Ca_B$T))iHZ zs$Ebc9XG|N1MQXZXD4XVMPgFvs zGc1O`jyOnz?g}7OqiHwH$X+(bmS4%+UDYe9D+RJ5XH<5}Q{zS%>L<<32Uh>NYU;h+ zm5=y+H7;EG5EXs#=W9(0%$AF}sK&VK+Ct|LFoizKL!DF0pIRk99x!nTlqyfrKlZ_# z0}tK|T5ScM+sZg`qe_#W5}Tmq;30?)P?LFl3?xTtmh+%B-kN*p9%o@$Bg->AM&6wp zrBQ_;eM??h&|7!Mlcvvp{b7cJB}`gdDTw$u-$7h^K%*xGK*NO4#O|_oF_*b!$6`o2 z^-`r{6geSs*I0Pw{av?sy}8SlJlUq7>Q?)P?EQ9qCVS_c(w6w<4@_$mJRN_yGtGKv z|E>E6{0#jIOz=V`3*Fsw_L1Q@$G{`B1yy&n{{yN2u2giGeE$a3{~ll1(+Q=}Q>1+X zqXR;JHKzZgfNkmp`u&$+@B7QLCCNGrW2ecP5|$l?LDg zzo^~J=~T}7n&=-XAU0aFGU8?7wYMkLSkt{$iWmzcFW^vGv41T?mEy?6hv`TOt?jm% z)0ZTmArRYz<(*V^{&>idj&ZJ#=`U>$fR1hSGwWTMQjo30(j8JrPd=O@7wH+T_;vRp zyhc$Yaqg}Pc`)8t!MuHIZnwFyPW4-h4TT|(YjwfR$$RJw$hyQ_4vZm`k5*R9XFKbi zhabP9v@Zsn_5ogAnWC9f3UcdV0*aVYTG3>~;4pYIkpRQ^{MJG?SVE%z{v^B1tX87V z%Whr=F%9(+@6)gdjhr>jo#jQ2*QUSjgLjz!8^>uRN1Y)38eV6CTvCnDEAM$yr8Lg^ zOzPHJn0*ldURy-Xagq_0tBJQMxwc7~!)k`=2Y}kt>U3u@`qOSA>Wks*C_GgdgWeNc zN(&tGx@)$_&LbehP+Soio&|!P#H!cyflKtis(pvozBlw1Ykk6cmG>z|3H)kl4)>x% z-y-k)bm0YCbkUq`X;eMIgYX*T8SiJb zph~C+TjpAPj9kLOF@ckjLye)-X6q_b21BmA_Xm(+?MW(j^tU;|%{rmcAM~Cd9%OW* zBbkQ@H)4cSWqY8aI#HB%Q^nHF{V}$1k&RGB$Fdx&twfLfU zrFRryU~^IK7Xov!&2TKGX%fb(Kn`SXcF~2&bjRleL>AHuuS7fN$r_>Aho)VYs(tSp zFb6Mj4lx-x@bOek|2DsCPMg4Af-{_%s;LBQzN`>D4C%!sXmZt&P7q?57c7;lOR-l| z0{_r>MKnChdk?D9?q-52IhT)g0d*Z(Lm4}nUH$|jLH|G-^?BXDARflndmNFe?hI}11 zXX+UJuVu{3h8s@Ww_kd$e%y+&G+@5b%>^B8Jqi!AwAeAg!s5cR+-_wp`<4+{LT;ep zZ)4M1#w1F~ieVerTR3m8Cu`GM!?aEXV8>?O*C^OPO_TivZKH;rX_@0aKgOoBCO
QjePp;bV7%QIV%4HE^V9gp4<@%lsAYs^i zxLSN6nOhPU#Ga9h@_Foq7d{EPO4AXuk)fhT$2&qUiGle!QddI{eZryC%9qIl@e9)j zc0mIN93Lgu-K!~zii~R{?Si&!3dIN(U##p}$NBHJq}<5hzd-zE$9Q&l7Ig1Go460B zjaMTFcs`Co@3=xuF`BaGs6;(^0XNv4J*LDF#vvZ#0mcpztm#k71qgQg*)QpAMfM5{ zK6?zgfo4jqMftY;R7G~>-?HrSXOiw8hort0V(`0f!@d`Rw|R&CcA~h_;(ySOh9Yc6 zK&BwS9b+Yi{lF+_6lmlP;U+g3@2-Z-pPEsL&|s&a09v?jd{$ml))^rH{}3;3q;(M& zPe$HvAvSwJ7-#~$TZ^38Jnl+s>b6SU+DzQ!vX$>v^PQ5C|LFW@M;O*ZFl-q2n?Jt7 zCZG$ipU|kATrJ|KHM-g=fyxjMMUv;XK*)e!{MorcD~Kj|xSp2lBcqR4&q5jP%!n!( z&Dzci9@r|1R?!S)xw!3 zEO6wpmq13a?FYJ5nz2#dT2s}?M#ol#(98hrjXxfIG#)2ZsMV?$uTC?5nPG$L^VTt} zOG=+3w2i(trKIkR{%Ft7mz;!$vYu=b2oFLdc{Fd1x645m~K~TUaCvKq;V!k!e5i$3%vfn?mK56?~ zmvQhp=i=Tj3#9--LMNVze8r4k!OV3u zAHRM*&I!f%2Mi2^(dl`WQR(sbO?I%0zBbc&^EuuiJ|~_dPh;RJk&2&OkF~!q)ZJfn z!V?s1SY9ZBu)tDr=#Z0bm|IoBSB{=b_$sHk*(cq$L_vrz#}qnp>hvx{ z6t~Yih;4Fl$r}?gVRd3HlTgGoXAyJxmr(prmd~jdOjInucmWeWU;lJAE4}j3kWIYn zUP-ch74g+PJ~xGsjG1O4Mts~oG#CUg@i2z?ME5mP!)+*MS@A~@VGsivA!Ms4I+4+158J>`3BF)}HDMhJ_XF47kdBhE(95*}&rhpa{ z7T4B5yblN`4?$LTKJD}M(DYvl4#7vGJ+R4UhbV@~rOVo-DL(&du(q>6XN*f3hsbm0 zYef!Fgot;LQ$QL(RC zsiqd$>5scRTkc=QHUc|*ox3J~{PH*1GEBQ6W9)tuPh?tCd)^P;LHu$1vThfiU6{<` z`Ugn7O4FwqWtWt*eOPrLcMNKt(5TVYrJ^OW6X6 z@`3*q?6>4#pmqlQb<`a^Uh44yo2}{Wca0yF81h^nlLcCnkd~y(3Sw2BQ zK4F=(H3cveue7YR48M#F7cW>qic3fiB+LbZ2n%yTc=@Dwr6Ih$uOQ$we#Lwu0WMx4 zE`DBZ5J&{fEAo#t$FyBVWu&Jd$TW}=-~aV}i;t2ZVp^ke4w=wjZ2Yg=2MO>C`~}Qk zpM`|e%vHP)g&=9mD)fjD-Y0TS4ru$k+3=`YI%;{j@Tk%W2t1+S;ps)k&oBJ%xrB9k z>(F%we(=8D&nDR&lGFEZ7+zZL?+VrVSM{xBK&HSw$D`c^#s>#)Q7RnxBn zg*z^h{QGVS7*Z4ju}U*B>Dm~t-6-f>N=D>6E_wnY`C=IFAehN#5+{RH z$^(mSE>hSDK;CJzvr;xY3ziOd=}FsiVwmiNlm{h~t2BeJ3n3{!JDi zdq6`@eohg+VIEm_Mu{uj&;fl}Cs(dRN)!^^WstiEm*qR-T{@Rc(=m)!O@fw&d7&(kwIh8aEpN-#oWJ~cg`)pWvx02?kwsep`qEENpJ_{bldmbXoSlh zSZGAIAN5J}by5}vi<}?!_hqZ39{5R&m$zI8FBaYSgE(kJ?dEKJz(w!Y!`Sw4144w) z&K>vwp3sy%u8beU;G3h%z<1A4hpFm8E_WTAchbNEf7rW}{vxipyR)Gkf!iHPxuKx% zp3cWZ0Zz0#7EKEGNV)0sWCPGuBMBp716){Ucf{u8QR^2frnrjAk}Y4y_Y6!oA(O7p zFJ0rybkVEbi<#DTJ}k>_Hc)VsHP$cE$JFct5p6&GEyV7yp9%Ds8eM$JKx;Kt4Vhdy zpD9ff=rm(kpR-W#lB%=S5HiTW{(N6#?`|Zm{t5R(t#Pnd`hC~&aPC5)$996I`E#<} zdg}X3v8xq|Q!+)7JnETKEj0rHx$g~Co^vY4@>mAk5;dX0|5A&xjiD<2WqeGBMAwTW-F)>H7G-ZZUfL zVew5-EG1=Iq35AsxsN~~AMsuH1I%5e5`;a|r`jzQpu#dCmprN*7l%qOwgt1E(j zqLp{}WHvbyr$(5{>Z~g;Uzs=CT%;^KCHp!cY;C|7cYe(2U z1Gi2|uPuB($##(ZB>d_ZwiCf1%fZJ%5FXRurg{&vB-fdap&Hv2`dGF^nvL!CvxHScA|bsT3q z;Z!nLccFZ{YTpR&0v);R8-K0Vz%tf3Uzb2)gkCEW2D?GYv2XaYd{x6h6U z6B!9_i)K<^)_t^HBq80Z^5vO^4^HV01pW3F#8GSflK8bTkP-zo1foV0ta?4_VNY%H zt;2!T!YOUEQdCa1=Ai8i1m!S#MCjIhA^)G}5ae$?7sL+|{JR=M{xO11W7h;(K=^nx zU0uECg#Jr`91mUz(O%mW!FUN=mIoi+D^Rn0VSqvfBBw3PgG<`v1+sv+ES-<#L z`O`j|gV@8=m0uY!+0Yr+NM%D&H@OwN-sgoRo#&spwi@&ec}naDs^TG3rlK=RWk_!j zts$wa$J#n3!;n&fn`2l~iqHN2us!A`Ej|feHk~cJY!Ccx5>aHKaWeqpQL|BwMe-UX zH2GIgRLQHB_Sx;IWR$cnmtVFhAjeVRU+LO}%qd5Gun#|uE(+MvUZ}8^a>3&~u!k9i zAu;(ssDBnb%%as)zQ?VDhc!z^RVyTZSgn8mWs)PUo2;^2b!s=Ex#exHf@%hQX!Lma6>6_ZpaU+SZj6~yt{z1ziVcSO@=Scmc Wv?(oWY<_-_fB+^lv#gpN=KljJZDe%- delta 16924 zcmagEV{m3)xAvQ)lXPs`#vR)>JGO1>j%_;~8y(y3*tTus{GW~Ioc->%>a1E-Yu22z z#)mOxjau_}T?Jua52c{-QVrbJ&@gb!L;xatBP$qYHa3=|Sx`aXdK|gB>oC5gf3^k^ zgA^386dw**+Vkg~m2*NoQXn=pVz0^wBrzh?eQ(e7vv$Xs%P-gG46A7kUJbq`g~hG& zx9at$JVraGB*s2D>g3agIL7bifz!A&M8DhFpTPF}`7nXEle<0z{nk@=;BM_PBVp^Q zx^ULS9IWNyg)^5rZ=B;EA zu-T;*`*z`;hN!Q{{}PcUoT=w$RLS`LeRux6^%_Rld2RN+#iOtF<+u9v^>jGRAmHcu zdgyk2v4*#H;F}{-VgB?rVmb|6yY?*72Np`1y45tJQE~t=x1+4zB*tqvVm;m@BA6!V zOA^V&^F-2B@`R(DpCLd?)Z{4sM!w86yhvF-k>7Gmkn81LGjvO?FzKuFFV=WbZq~l& zJm9f=Unz8Te0<&QcBT5lcpfz!IwO-49362;bk1lqO=xBkw0k!(*lc@>x>%-z(0CE* zo&p49?oL%##TcF<8*HPSGfv>}xfsS+7Jo@jK zx4UOt1y6~}dKpYP?<+jE^W5+DkI$okbbI`8?l0iiCH^1nYwojN*+T*=;fn@$&mYWx z-nfQIM|s;P`Th!Nkn-dfRIx6)g&tnxX}q-*W`BR)3?w|)3Xp#*nm`l`@o*`4Pt>?I z8o(mT(KG^{H*&WQi6ttnKSly(YI9Uz$}qg4)B#E8k0r%?6~FS#GrYo;YJX8zNgv(q zx#j~Q4eQ+^nkIB3tE7NuqF!v@l>r zT($@H(XkG1xe5jirRI06(0^%9hf559XJZUDvZ^D@Pq*z*6Oq(m%Cp8Ggw|P& z%&#}akES-JTl2Oi4K$@VZFw9r$j1b8=my+g;nHb*-C4baG@DoE@BO0-P)(`GWed2? zecmkAZAcdWt66IeGD;;uS49w~ticeKhhs z|7pZnXd<#9YhO>d9>Jmd%Bb&g#u}iZo{I^uj1+4<%@GAjO@q_^5F+^U2=`hLT#IT& zZSOy?UGYix*UxL(xwF#aq?|Xbej0d*yG7lmsF$Im42MVw*zF<6U*SoCnUkOO+Lrz#!d9 zyB)f#3S6drU#~h@P14)Co^dzdfIqGdYGtne#^yrlG_ttJ5U9Pvd{KjJ35gtM(8on#9htu!-$GJlFX3ldGoO zNZrN;Y{x7n#B_+CffUYKZU@^<^gQT2B;@nU>k0oS`XZSv!Mg{Ynjg1|NGN`60X}^@Bp?0ByE)#S~GOd7XM9r~x zH-THb+oakK8p+)t4a-yX-nMWaStKN0izc@#qC9i`H1r17!*8Q5EBxrEz>+$Njj@Q( zeBK_s)>4C(jhC;d&c7KdoH8>v`uorE;w=^WFc>3hYtX@my2nEau_jPQ+4AYE9mRy1 zxtLgTE+5N-5P8lhamSfO0Y+L3b`UCDK1{J(2CkBy0p=_&wUIx&*95x|R^fUtuEXcl z(IgHgu|<+AYRsbIbzX;$FsTKejS?vV(g?+~o7bkSn6jcEO!Oj%zuwO;!7U3fL$X-M zb$V`An{#DNkKET#Yz=rh4^cxs;v)}xEW4}o$IZqb5>rhlrYSq$OkpUasJCNTB94n& z0B5gemoy1eqx?6PW93Ao0ZI`qG@g=A8ajuXGNK=5VMrmA`+FCW*sKeU5M4>#s~#;( ztLx|1`=0>+mWysmo0f8DyvmTlM{l)402%e)9eGlSSY# z6Jp4#5{cFCJ$b%92rIrPgC)Pw%7&gGOPV)ThRCQVWB3Ri58VkvPSTIEOEApQzT8Hg z=Q^Fnw7vO>u9tz&mml9-u1kbmziH3QoT!-(5*2>Qd^_XLxr2?hP9}Nh(b^yhNw-k( z5{jyr1D@c_2qmT@GA<3Tri<@3DTP0sC3yp1zON`aS|0YLDDJY(Mrb!_?8W%>JPjZ9 zqrC$Cr}&a4jyV$%${Iip;L)C4vie4>1`v1k$9K7$OI?8`>|MJ=HyFQQVlVN20b4+K zyp{!7c`B(J#`hJ)BW!8GS<#6KGfJ5&s1;U30b*ZUDG$S%ud&hzNSPsfMd=tB%RrkZ z1y$}w9)*=($phpmR`l;o*)1tllGD|M<{PJO=Bfq72Q+WyPK}UNPdaKRd}dd)CPOZ! zgDzw`rdz<;?KA|4Y4?fM<6avZU()6_W=I`(jof9HhoSa8^Z|V1k3S;HZ<_%NhtdBu z=XzyFU*lr$&E4Lg98mnP`TT{pGpRmvk+NdQZ?)1^1=^9~+;suMtWh+y=&qT2wQA`W z4!6IjQap8kg@$^G*_-pLDK#{~TD&puqn-f3j<|rc)=D+$sxFCZgYu=K*JAI2;HoVF`7Qo@S4;!26yDqH7gLvsQ58 z&{P0Zc180<0{Bm8-bkhj595QVl1#_mlB9a4rhRV$A)rwsTjZ1d6-yliN(U`UC5)wk zunU)``X1N)iI7s^G>-|Lv?IfY;C?bCYh$=P)evpJ2C9pK;<3WiBNkf9=Q|US3L3S; z9p(FHxOSt@0zL#{Sb9xmUxk{U;Col8f z$?>9XafO73I8Sr7%}YON-xME42r}5lSW{;3GBn5xw1OB2Q+ToMjZK1cs)N8+Eh2K$ zP8OJxPi{2BVt)Jbaac#3K}_4i|EjD%I?2wxRoGQH>5AdkO}Bz@5jR$Rd=TaN%ftel z{W#z7wj822Y_y_&U$K;h`F+P`1MG&WGmF_p? zz>Nm4dhi3}3|{DPLp*5KAB;jI;ynP+$*jRnjjX{SY0d!0E)jiNR|%p&iCFNaarOx*Ghh0Zd$`YKqy-JfX1j)uKF!_NIU!mfM zPP{w>h_t_==Ye^N&YhL9L@(8RTr<5Zy?SzFx5BsyJ(;v}BiZXJT4XdVAr%ggtWJsF zpN(msCbm4Kfq>$%mIeF*L`*L{8J9+Cs+ZZ2MPxoV;un|+`h;Z0+44JPR<9omL+BRv21|iT8PD3Rk`KywJhuVPAO&<-Wy^l~ zbV8HyIub{RK}V6${V56U?ZU-CnPPrGvt~?GNr2Dn_60dY60Zb4lH1zBZ_qswLnRB;BWDn*{;eBhDG?f$>%y1BZg_?o*)@8dcbj3jDNw2 zCBlQGH}{|F4t7iaNUh}D9PT#Sn1xqk9UlM#*bLabn%+qgiAT*6!RT5*WTE*BO>aqy zOG1N!PM@%^*=RW64^!p9$yv<41v<$APlv)~hTGx4kS!aqC7Xq88#_E^XSh&W(=*85 z+bHj`eOFJlM$w%pdWe9;Qf`72GCVhRH})##PC?*EKOCgPgAcaUv~c+-Lr~%hL}F8I zy4pt-q>CkL(en{n31EJhVjbe3;XE5}+kme@!966_=g8I>(b-{pEb}AR)QRPibC3j( zx9@wQfl#kkf%Zm6YoJVXJ#+e#?JjCX==8%Z+|};|0^5aJ?%jaHv2o}>js9fi_kl^` zRkNB_c@in;n=#d=Ohm&1njD0r^+Pw4^k)w$e#eX2O<|mr898o3zD_8kO`uh#AmOlk z-9${p$W%71An>^llqXcrL=`|_2hJqBV5pil=X3mcjW2>uwF*k4Bb45jrTr8 zd&|c_k}@T(SnN;A8=1>VP-{X7e5j=u!PM3PQxBee2_+X?rQqz%D}`{ZQkB?qN zIRX%q=|}>zmTY<*wLnbhv8k4@pG=HJDOTn3L4}*e$0Ue@!9OJZLD5+~tS&j*krAq*j)^2;K6=xE5=b8kGBIKYb$ z+raun)QVyW*%SC5K^!PBmSk>@28J;qXp(A9ju0SU5mxd(u%l@3;_QRvG-X(urFq(z~y10V>v2(v{9qnk6=My>6P4r!rBq8D&K zGghrF5C5VL(R9~eBZ>B<(hZu-xB*-zRr}|8CQRoW{ut!I;{aIJ-ER8ruGA^e;Zb_FtTYm4%t( zzvq}alHjrUIGGs$)QqZ@E;goA&c+s|wuUa2#taT7X4H%#|3q-Hw6_yAbTOq83g(>OQ8#BpY8y=V_V;2NK5`50? z|7(UEAm=YCED4ShAtpd`2?|A^Tx3Kig0a#46umA+AMju@eY}dQCbwDU2){kC(Nb4$ zs5;lnEFg<`0r{*1p#$K0qDB6EFP|6E;B$W~&LW~a^a=7M6{_BEzZ9rUmBP9RD6wlX z4JlU|$Ufffh|yJR-~$@=GWf8eOq5+Iup41EBUbkq)=O!qq9ymLy{xbT0-)P!`DHX* zmV-_g@B?}4D@SYH&_StQ?RhZJTyfuOr**p&6i{`_NTXlz;?=?A*h%l!fA4Z$)Gy zvy@#Z@o{lNd5t83dC_MM=cJ;dnEkVqB{f`cWaeNGQaYt1@Do^lyI967-<|&%)#U ztMV*N|5YIt4#2-F&zUrYMU!-+0u5wg`Da;~pC9JG&Wp#|riP~vn%GnA+47>ZS&2Dw z5|1Ip00{B2iLk#Qk%%M-Qm+FYg)l&plMU$*X^2hd^n>HhSXYok!c_kmuIg>|^Yi=N^u-NdSNHKu`?|{u5hx5z?4>N?&v~Ra zRZhZ;AItX+mG6C2qjS~zkU1;~B}*xKGgYgl=sB&Mt17%wxN_C9bM=0#PxvN=LVRQ- zd=LVAJXwyUU`rT$?ZAs&;bNhWtg})x;PxGo@>iYDm%@s!IqT17R!7=YjtCV3Wl7SD zN2AkJ1&y2d_0VIe0dQT(^}4g-dYd!Cc9K9;;G&SKlt6~`EBQzs<-<8K@0&{)KnfZ1 zTZk;=w-_8k0CA#7xI|I#@2X>XVi+cNa9k1raAYQG@Q|@FVvqt(aA+n{a8M?eU<*-A zD5h{O;u9oi5h#%f*uG#5v_MctcPKYeQz*5;L?pD()KH|@Rv-s>NRKJ;X!s$q1`@4s zxL_3=sAsTw5D|1RnY2*6q%fo(DtN@uoiY-cys7XE*|oP-!z+UG6x*CaCBb;a;cSNV z>*=v2Z3oAL@>0d;uM@g5Zjii8=~#23+|#WY%qcCA&+jki$V(IA0s#_S(|12UB*^yBX+OA})m~`~;x>@IyIq7@;pw z8foTBAEHH2|7fZjDnQoM1d-aESZS79#gjNj3vJUW;x=U95x91uL zr3406*+0C|Erynud}>0~sPaE8$US4b!&z#E+0qi(z*75FvtN@vAiWh59$(glb46mf zilv(|4!?shC~?Pr)yuVt+Ks#sAFN7oS{1fKUJ@-kmQWZb0|-ZQdRN`z@m=UoiT;WT zKXTrIodb~)9hW{u|19y0$uIRaQ~MaK_SNj*Hu6vtt+Fy~g!Hf5d8L~{C5>r*__1Xm zH0%&?N-Re(-14`Iet-oy)u(5h;U~tturchV9nHTFIX$*$t>!u8X~<`dl4sw`zaX#6 zC)<>^8w)?ds;ru@V?Xl}zPnW@0B@+r4$xePuYqx}*|Dgw#F*Pv7?g>&CYXL)01Jiq zkV;1K%GHvcIM!Xp6x%?rY<<%pGnm3rc|-oq@HssR_B3IB>`?2$hQB8Gn{iI4a%`Hz z+^X(@dMD^HgQ7@OQjVjOf!;o14ffIz`zGZqY=@X~VO5T-Nqa|i@Av`T-FqlZsVVFY ztw1ve>nyRB^M0VWCox-bS?eJ3?*xd~d`sBz`#;WF&}u0dZpJ{N?l*bOri`}R?Uy^V zk&Td8t+pk$0<~NY5;nm`YOxn8hF8?$xD-f>?RJo{IDzm~IQNO+<#x2NLIWi~7W$D; zeh3~&oG zw^g?mwb2VKX=+6l%5sno)q5$aN}92Z+u%f_b8*5dA4J(w6gM^XrT$&jQPIZLjsm_z zHT@{M#Kcejs2{*QgF$J4v8`@4)%=0qotSx_(H*w>7%815V)1{xgq1{B?kMj+<+uo zbM#=0!Vf(M|Im_&QWXiy+13E7g?3<3tf7^z>buY-@&dsb?O`w1gPI6pm>q(y-itZx zhX%ZWPOY(k!%Mzgr(bolZf(#x+g?-UkCi6n)+5wf1b?7U-qhDD#CGYxHkZ!sjeuUO zP0;mOKC0No!3_2Z&;B{;w;QaFcX^<%Hewph4laT|5dClW(akS9pRwIN-5Vf8SWk0r z$aLXq2*V364u55s0+&z!?gtkS&-89KK7UPE9Did&9H&p~STBF$z|X~-FQ}7V@ewGG zxPevo#|yeJ_eY0N5uUL5CJZ-jgfA4(QqE2lM0H{M9^=n zCz!Aw8h9sx0{c4eT%9_gjK6yGWM^zb1e`n)-o17a@tzz#`e5Ynp01ewbm^@a2Ohg? zuiOV2v_SZM>(IRWvhk+DOgpU%?YiUPwL9w)~5`$aaecW;2#ygMPVzCZ@|x%Q0Z?(V_IDsvTW&ts`(vZ!?v9~@Yds6dd)aq z65zhEJmB2l)8W~%eHqXL0VXXA-Dj}S6@zziE6YSPaf%WXxr49_F@;=?L zdW4Hi(rDbDTQ!Fc)K&uXxj_jI>ej#!9a`fCq!DOSaaUnjsswTGeCZfC`tR5sGtg6C z2fevA2J#~Z=E&z)=DakNxx(o)AX9>N>=4{N4Vl8k^fXo=V6_bR!#we5I48PT0wxa~ zC=*t5paEebTwy$%6WyD=W?8?hz|3-jGgi_qgVy2wadeG>);WO_J2#ap7EiQ{vd4g! z0aXN*iv09}YWDOALlimFBu-hIfC3?ap$}RX0B$(0ss2_8P(u_Mv+n@gY*+3tY=hn` zD1>%kEb#@v$}%@Jt_H9Zz4&l5*+*v_Y7#LLA=$*-W~^YbL_^Q8+q^o(Oqhcaazoi1 z;<(Z+L$#eq`2f-PED#D6gOzmS2MiC0u=_t$_e?|{>h63~uuwS@IA$ghNW1AkRB#6HQii1TWN$k4P zc~>SAx!{Hqf;QvHq79_*ARVzr4gIkCOSgR{&)vdHq6Wo*b_omvW8w}%z@OSx>!z_C1HRAyB>>$0HxV#!-9Yy0{4RB z-%%jda6lXc6nT0_19p%NJeS?#?D1e^K_!q(7*R-^-{|M;vj^V-z=QABk6VtrponnV zbV*er03))Q;a0&xEgROWF#@I0)))XL1<;dJSK|ooEK|)0h5vMp!z9-Ku{=WyC@ zbt^O&mma&j0-%#?jXmWEA4oUa<1qUe_R$*)05VbNjh+vMTyRZ*vK4I`KKczWmR*7@ z);9%8&c*dCX9yiUbN{1k25VP>5!|l;m;XAiM zNQ3(u%5jTRgafz)C3%D`jdOAYiYU@1j3yC~PIzfKHdYR-T{l}Y|-zAdu^&;$Y2k*q9GsHoSf+}`{(y%So zn6ya&V?;XyQn)tBZ?xoa+2oDr2+!czPofA5aal=of|!kY?2H_tIi9Xnf@o(}=688g zYY~I(6pD6Bit%`Ax+3@NU9n}4*$j$Q;6_SS@Pz7^`H=8|m_!ML=WJ!xz#U~LmZ~0@ zDAb?)S9097@6phH&?XbDn_ZYvcEgRIoedl}pw`y}6y2+b$lDhLC1!WHilOCHu^24J zLi`;IzrB-L)5(pNEPJgdX2?bfXs_x^=M!Uq<7%~TUK?Ba7@q#~rqJ$4QL>3;z>uqJ z@AFB1Wj!JD86UrNj$4+Fj#Eypg-u^VnPY8Fcwo9fz76hkj!69=vjjRA_If;Q&3EP@1JkX6I;91(G^?hu~M2NJA)`LtBpCtQ~67EnV67 zdV8sx-5%=c7X!^r_m%y&(P>v{OO1Am=`-dr|LZ@ATQsxK-G*zG(~95$(!e?ueWWE2 zd&<&9^qzqXDvA;tK^pNJk?$dchxVoi{owNVz`x}usj`71)i*lC5vCF`D|Qi!ipQ6# zTdzu}0To_CItR4cZmu<*4y(k#;@Tkj|wqC^NK zL?zh&nF~Q0R~)sdA|tK`BJU#?hl_>3BjX~iBYPwNM7Ag47Y?R&va}mb!8UAK2amBX zW3vtyG%sgnD{z?`m(RUlg~*+Mh4=l8ih6J@IJ<4{W&||oRWo!9a%j{NGa(o!YllOV(jP+?d3UP1$zWw2A^7o=S<6c{u zJr|3um)K-Wg{}4sc#&Te(S>+f)4JcJOE0lc*akd_@9hmxDfMarcQs}@+t=Et4B0q? zO=E>5!mt1|oyA%`B#Jhwy>{*A=xt1hE%du5`+E4;WsIsebA|_T(-G|DXVj1Ek@9n^ zJ!wWtG0~?+(JwJqz}QGwk2U@b!B~#yb_8#6!|q2*#6 z0F6YM=L}?v+RYq;-?C&=%=zSJW4$?a^#P_ukDy@oV$$jJfCce2?#CV0I8s!0a1xq+ z@)WrAlDfGikd0tbMTi8o-z71BCe%CzBXZ(^h|QvG)uMm_7UyL8iEfxc$7fmu_eVf9 z-(S5ITD+4SB0LL^ixf<+<)}mN@nbzj-xa0>4vrd3xjA2^;TS_(=2fgs45gn=IW_3b zBf;7+qBPM7wFJ*J&tx@}_mydf)@gip3y6!h+zEQ#z%_+N2^R$qE1K)w(epU)keZc- zcZ7J`@pQB)*u2mIuO^)$)sb@|J{oU5$4yFi(@&&Ew)PjQs?MLzW7x|COqn$rO#G~| zj@OX|pyS~uVP;I6K)YR@T3Ltg+ugy(6&wn zYGA+koZ92XX9aJob{DE+<}1{vy-fQ<<#Mj&dFnj+c25y`o#eD{ChP_L&Yw1op1v9z zsplgpT4kAdD&1RJb*eQron@?FChuLTyCyPhzNj#T00OM9m5jZG>w3Zjc+&>$K31U*& z291{F6DhPQGAKvHUM>V(Bk{NixEZ)9xH-5g?L$URqW^ggW-0uZz7fE1iMEgjzLqc4$IF@bKIz8u$32T7fr*Ez= z7yHjVQ#t1z-}^Zqx?k#k0eQvcKvktO*24{sy>{#i7%kOIMv6DnJ7v)y!Cf{!Qg$odIX}jlw14$`^-h>Go zKk}$^9K6gVb21B)it*tm);W9-HO6%VKcqLGr{{ z;oDFXGN?>g%h2*t@l#N1dIH;YItE`g-h;-QX@TCGq1GKYviqZG`X~|3ED>U2BZ4MX z+KC#uHlF-deQJz}kS_IJ0cWTrsHP``jL!yevr8C_h4N%OXFMEtK+qvdL*&yT6T!C0 z&_u&R%hZa`zFrVg5h>#WmyjEy=|~KmlwV4%#b!ri6%h6lToUlb=mR{NnVIr4jkAcX z8(Njg|Ms+pS|ESsZO1{2rB(mlga@hJhwSj%KFv(RWnx3`WVN1M|Hf7JM;wd;0q#9$4(p zs*kD=``O=@JeYnks4vTMtT=@I1AM>gu!Y~eI5g?)P~3kD07pMP{rZGGQxeTj)i6d)tH{}=1z_9Z zO4te-dcvD`=CN7f*CeQ|k$zR*c%>^0cP3(QB#QS>C@V3>Y^+AixEYN{zU3EcYWSy^ zix|ZPbWch20-JIbLV6wvEUi$6p}I>bQnDgCHD7WW`m4phSodHqrFNkTDQAM`MOMr@$MODHR8TFxBI|U zm`1%y(!#+dW@F=G{cmE-!v4QGKj(jOeuN}?NjRVzLYIbQN})7frl{*12%I%$Ib`GT z&u@NMd9M@cDpn&iH#f1$sqqdnJ1p6s5r|XT!CZ6}{?7pJ>55Ek^SRJwI|}>^3jx>8 z6X?MXcs{&jh&&xF#BJLVdN@(VW4jzw*>S{t%Rx-brB1%dxMg-&kH`=bSNWYW;xMAa=Kloh|H0Ea|0@&B%v}G@ z4jWNY0tfYf0CiTt{|VG}RqSPPP4T|Ex~vgl_Ty%lc^%!8zw+t#r~U0rYy@G}!6M)Q z=e{Rq13jrtLAJq$qq1!6^lVB3bl8U5GlP~TT*ebHe*SbyX(uQ#W=U1EV$ri@gsj@B zHdQYhP0TC0xgFcrp1^D8o@dW(AgQ>hStP2&oB&Wx#ysN9uPokzHRsl4Au^eh1 zkv!uHdG-AUxPA9Ky|Ft5hN&0w$=gcSfxH3)-{ivmWLPLg>V7^UvX7}r)Ho%MV$0OF zRfLj9{hX^}uFH!>ssv%TiQ|Sbu4vPhpTeq6_x7?4%Cyx&j73*{(1b*?UOFf%OaJ49}1Yrae@& zFck&3IeWhB^?d8jaC{P*Q1T|KZ$GTE$`TVLO_;N)hEbuJ0$fe*$>;^d0A86Ee^I37 z(7C+L6NkJX^IPRG*!n5fO&-c$ZZ93}-GRSRbH=6N`R&%1_P0c8ENvOI*-jtzZHUbUu63O9k42hD&Z>ZzNq?k)I5}m9Zll}3K z9;(?4? z2HOAJ&njgN?a30U!PZTX@YDtumlj!)h^~#_H5+Q%%4ZELSv3#LeI`3_JY>oZXiO*& zD)%Q|={{@AflOXy<7dhx6e6d2_X2tCZ?N}`JOw|s?&ez3*NM2i1Uy^~;;3G^y{4uf z8zNcw1io%E=XiH2JwFJN{w(63yL<2uk}&%IS%r>kc0kPF$gZ;pCf^VKSq+0t3C>CU zaUBfuuWv6js63Q|usMTc{{pjuA&T3A^Z|(#SU06b!IRx?UTms%*_@ndlUOJv-dYLj$BplX#>&>OV8^VN zsK=JFoD0{z<5mlN+3OXsK7D+bo$~EIy$ldaqwVURnZJL&w0|k+#^oHXjRjP<|u!GDz9e@nxx|5642)$RUY zS(laV-?DDJklg?glGt;0plN-0T@t-;B;GHhS@y{J9|l+*l8^qw9=q#5ja#B(!BK9u zC%(hQw^`rFuFz~?Bc2KH`_RnldwH;U>)$4FS}b{NJMiU93OqS(2cn5;7+f6Y5CE@@ z+1;6`JSPo;+iRxM0IQ3jRBhMUL~KzA9p}bJZ496hfnsDpqf&W zd;Ix>a{~ngqKX`3XkmJTHvF_x{J&)Azr}NAW|n`gCiSX-gRwFvJxctigv!eCe|G#K z4Q~h4MeOe>-wQjGGAZy?GMg)1Mj3)ED4{+<4?05*@P1OHUfqK~StcSV5hYkEeF0Y* zb@fC;%4LCt%GxxHZ8lh)0bmIMN+VWE;#!1-iasuQ*q5aE=%bg{Tkw{t5v)MKQ|{~4 ztKRWDu=}3pc&p1i0$GV!ES$2fgnB`?MT0QRjPRr;a|z?Me4P>~IA|5|65~4gH}fcW zD^G;XM=rVEbM-F8)mA?_j~vGZ;;Zy*Ab%}Kr8;M4I%s%7B*CYH4R!xAK5UVri8p&e zYMnSP%trf5P!8%tGFpNcX9mryb#;l`fXzh@IO=pC>2P2EKqwKROgic*-~#7Lf!l2{ zoFa1@RTvxvB$mH1cbw!iU+(_%kb0YwY{qg}%mOrgBFyR8-{9D#<)T(t>kdg@U;%TN z(NiZ-W?Eq^=Jljko*s%gRi)7(>#^KD0gdjCtJ+bXO-0uHPOZ-XH!}KJ963Uz9<5jc z1O$?ieDck_i*K{(f6Mx$D0gH&^xTt&&)Xp_vsv=G&Y(61A;gF1f#=D{oKSn^uqpFs zi-w)x{Vn0Cvjeka4f+~&`LwMl-0trQk)KoJh6;8nkANI#k7#@N70`&M( zSjhur{(4AD*e?H zLfiGVoIH8sEjOi)vu^^_Av6FPV+unnv3IUDiiuOIhPnhBm`an%)!G#`x75GKN6vm% zlp6yZ<{xx^$UTVyw7gV4b$Hgt*4m_(`lcMbfh!{yB=3>BN1}Eh9O=%^^cxsX39$(uU!8fuT;nbl8>$ z#ai9lpz?-}@8@4bK0q6nUz|Tozv_u#T3~`g%KKT@s92%bXWgI~VVHIm_B0#|IAwH* z#0MYu#Jk0Pw0X@njaN+~G`s74Y$uzh>i7*d!xugAy`NHokHriI6hEoj^Q_X-e$BI2 z4a;K`P`Rnod*iXOtJHL-1N8?K4#kZ)@Jc1>v|Hw_w5Ayga^Cb`_#VtJk`Ket^S5D( z50HQ{i=N5iMwp*1wtW=>aUAKPzXag*-Ismt+nggus7I30^3t+cHrQm_HmPM+WAimL zu%C!q;0au;y|m7HL8-^2*ek`*JYtms6sTdx&u)Ixz^RY_g@xCI1_~jc6+=}ptC~je zSKzU39*8O+@}Ae+KMeBm@#WI}B2d_eg*C||skHY(*EEJbe(4j~{RLL!!)x>`W>xEpiY23L33mc6Y6{Z~BOCtweUDHy?!s}yIviDey z9l6|i%b;JTR-L{74D@Wso&s@EoO~peFcb<>qO`?QtX#00;Z?Sr^L3+00NNhA?nj9GJwcBX>h-djXS> z3oK%Xb><#vH0>oJ3zCaX%x7JoXD z)d^G8B3boRV05L#fB3-_#}ohZ#dSZ8A)u{y#R^{UAWTsMM64Ti10TW(8!QdmTP${b z*0M>dO}1vl3`NHbQG^o*3lzNaqc`G+xR5e_aUk>uFu*aKr9YfB-LU{v)Z!v`#zI`s zztRCmQ!&62q4r_=wi4Ms_3Re14rZkrKgF-<2P`I3Z)I7j@Hb}S?ow8*nNhHhP!HsfGCAHq(2L_u{ zYaZ!B;4Fu@hpc_@wRn@UY;cr4O26@Y&`ql!vCbSv|GtuQM8HijA;NA1G?YR+Q-AIj zl!AB@0qX8tR`5;o7tzHfW`B1#GK0t@|Bu~r0lfahfA)1S_0M~g>I7JoWgx*%L%rg4MR->&tqr@-8hs{--_4LH?JZ-X`Oy^;>8s z2HP@LeAstf{4%-lm6?R3$z|ihq@uH+ComoLQX;NV6fZCU5o;G4ScPjGFV zHJWMuT%EHmZzf!YTuHQL&-!Hyd^buW`>WKPMgf}IWy>BVx3D6k2qUrTr_%Oqnq{zW zV}e_BdW?t-e1n@?*%B#-{CQM2bU=r*zKx-N`79aO2VFv(!g8hau}u=Y%9$Cf@JPuD zTisFZn&_+iHrV*bbnChU0`=f5ym$Q!KPx^P*O z-ceuO)b?kJ*ORInKviaEF(d_0rvLjBH>)G=*FexB+oL!kOy#8XN+mmu-m8JNF5waP zz($1!qH>OX?zce>w!z|Hqp20hIiT3g9$z9iC)FpO4^o74$k8NPJN$u?ByR~~?j>CF zPJ+e#6qn=rL&n>3v%wi@>l65U5;gF94AR779w`PcDC8A6@$#pnU-oTBtXq8tD=E^!tyaTczm zB{w!QRv}?7VPRo*dRAs806jpQi&tlX?zv;aWTmb>hafBdjX^d8LrxAg`m4|Xu{f0<*G%Dgh*+5UZXVrAh>3Q$H% zGV-PJ;G90J1PMid#B4yAUz6GXH;#mep$8_UoK*e#`oOBN!s` z1%4~c05pqia=BVfZ#j3`nv3sYo>HyLS_9X(U;z zc#I5`-w@Ey!!#vm!2yPIf>i4i;2(Vi@7r8i!~v0r0?PKXkIe-Pg$2sy^_`W(TfIAI!bCP3S53A zwP(g^(INFOIXVeHO_sQgd%Lwrry;h%$gBS(X?EsX`8?W=>`Gl4(5dtcO{Cm00 z))C7uMMY*Uh4Yd%yc~uo`>S+)H5v8U+t}P#ciCU#%+7Hpfxk0eHokm*X|anzRD=ZI zZJ{Nd0xeRLr*m#+)uDn-Y8ATH|Bh(ZR|rx2ovquRz(5mF9%jg>1Cz)L?U3N%iZIqU z2MljN?F^zr0$@&?k+Z29JGZ$7DC25+E-jNxPWM>#O%cHJRK9_MO&`s)rY@~=0C^S( z)jIldjA8=;u!onYTE@cLFqX8>^GufnxKOcXAg_NT%o7p*7w_EXBam-lwnH)OVb18| z9DfTU{5kpvY8d?D3*>`29t(0B4Nyx>6r)m67_IX_@LTXr(qE{~ zI3STx`paF6%uk@|XlB8iK!n5LyVdXpgT z^jE^i4MksJKnADfiI?UsTbY1_eC=D_qzxCtOYR+t>wlk)&x9r=~ZkL0;E}pg*bLd&8I}Jr z$q}*rhy9vS%FfK5i1|M_0;80uHZKPY6PGA63#TXttBA0$h`1OB2M3E78z+-6J13_o zix~g^?~Wv-KxQC-h3G$@|Iv1)dZYs#(S~oYuWuQ-F*rLdT98tsre?`lJ;?*WV*)DG zXLIFYW{aX=gxraav_v{1gG-+~JGX@XvhemybM=t+$`)Xp)k}xPgzF0hFelJ?9jN&D zdZ6|ws_(tZvwT?FnQo?hGk1??>ieW(z^_JZsNDTO0cQZ1|H`rHW-~>y1{@XEkjwVi z^$ouw^jzAT%jGjgS!=*SVevreB-?%ViM&hY%;o%^imWpxR+x3Ep5*N%KMIp~F4D%5 nJDaS<;ml`?@+^!zIfp+nb~osgdV?GcGBGeRI0_{tMNdWwJ9`bG diff --git a/src/static/schematic.svg b/src/static/schematic.svg index 7a421c8..634bfc5 100644 --- a/src/static/schematic.svg +++ b/src/static/schematic.svg @@ -1,1538 +1,425 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - MVN - D* - - - - P(D*) - ? - Variance estimation - - Calculates(∆t) - - - - - s2(∆t) - P[s2(∆t)] + + + + + + + + + + MVN + + + + + D* + + + + + + + + + p( + D* + ) + + + + + + ? + + + + + Variance + + + estimation + + + + + + + + + Calculate + + + ∆ + r + (t) + + + + + + + + + + + x(t) + + + + + + p + [ + x(t) + + + ] + + + + + + + + + + + + + Σ + + + + + + [∆ + r + (t) + ] + + + 2 + + + + + + + + + N(t) + + + + + + + + + + + + + + ∆ + r + (t) + + + + + + + + p + [∆ + r + (t) + ] + + + + + + + ⟨ + x(t) + ⟩ + + + + + + + + + + + + + + + + + + + + + ⟨ + x(t) + ⟩ + + + + + + + p[ + ⟨ + x(t) + ⟩ + ] + + + + + + + + t + + + + + + + + + + + + + + t + + + 2 + + + + + + + t + + + 1 + + + + + + + Σ[ + ⟨ + x(t + 1 + ) + ⟩ + , + ⟨ + x(t + 2 + ) + ) + ⟩ + ] + + + + + + Covariance + + + estimation + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - Σ - [s(∆t)]2 - - - - - - - s(∆t) - - P[s(∆t)] - ⟨s2(∆t)⟩ - - - - - - - - - - - - - - - ⟨s2(∆t)⟩ - P[⟨s2(∆t)⟩] - ∆t - - - - - - - - - ∆tn - ∆tn+m - Σ[⟨s2(∆tn)⟩, ⟨s2(∆tn+m)⟩] - Covarianceestimation - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/src/tex/bib.bib b/src/tex/bib.bib index d64992f..0d450fe 100644 --- a/src/tex/bib.bib +++ b/src/tex/bib.bib @@ -530,6 +530,28 @@ @article{higham_computing_1988 journal = {Linear Algebra Appl.} } +@article{higham_restoring_2016, + title = {{R}estoring {D}efiniteness via {S}hrinking, with an {A}pplication to {C}orrelation {M}atrices with a {F}ixed {B}lock}, + volume = {58}, + DOI = {10.1137/140996112}, + number = {2}, + journal = {SIAM Rev.}, + author = {Higham, N. J. and Strabić, N. and Šego, V.}, + year = {2016}, + pages = {245--263} +} + +@article{tabeart_improving_2020, + title = {Improving the condition number of estimated covariance matrices}, + volume = {72}, + DOI = {10.1080/16000870.2019.1696646}, + number = {1}, + journal = {Tellus A: Dyn. Meteorol. Oceanogr.}, + author = {J. M. Tabeart and S. L. Dance and A. S. Lawless and N. K. Nichols and J. A. Waller}, + year = {2020}, + pages = {1--19} +} + @article{hoover_canonical_1985, doi = {10.1103/physreva.31.1695}, year = {1985}, @@ -680,15 +702,17 @@ @misc{mccluskey_github_2022 year = {2024} } -@software{mccluskey_kinisi_2022, - author = {McCluskey, A. R. and Morgan, B. J.}, - license = {MIT}, - month = {1}, - title = {{kinisi-1.0.0}}, - howpublished = {https://github.com/bjmorgan/kinisi}, - version = {1.0.0}, - year = {2024} -} +@article{mccluskey_kinisi_2022, + doi = {10.21105/joss.05984}, + year = {2024}, + publisher = {The Open Journal}, + volume = {9}, + number = {94}, + pages = {5984}, + author = {A. R. McCluskey and A. G. Squires and J. Dunn and S. W. Coles and B. J. Morgan}, + title = {kinisi: {B}ayesian analysis of mass transport from molecular dynamics simulations}, + journal = {J. Open Source Softw.} + } @article{mccluskey_uravu_2020, doi = {10.21105/joss.02214}, diff --git a/src/tex/ms.tex b/src/tex/ms.tex index 36ff286..593ae45 100644 --- a/src/tex/ms.tex +++ b/src/tex/ms.tex @@ -76,8 +76,9 @@ % Author list \author{Andrew R. McCluskey} \email{andrew.mccluskey@bristol.ac.uk} - \affiliation{School of Chemistry, University of Bristol, Cantock's Close, Bristol, BS8 1TS, United Kingdom} + \affiliation{Centre for Computational Chemistry, School of Chemistry, University of Bristol, Cantock's Close, Bristol, BS8 1TS, UK} \affiliation{European Spallation Source ERIC, Ole Maaløes vej 3, 2200 København N, DK} + \affiliation{Diamond Light Source, Rutherford Appleton Laboratory, Harwell Science and Innovation Campus, Didcot, OX11 0DE, UK} \author{Samuel W. Coles} \affiliation{Department of Chemistry, University of Bath, Claverton Down, Bath, BA2 7AY, UK} \affiliation{The Faraday Institution, Quad One, Harwell Science and Innovation Campus, Didcot, OX11 0RA, UK} @@ -92,7 +93,7 @@ MSDs derived from simulation suffer from statistical noise, which introduces uncertainty in the resulting estimate of $\D$. An optimal scheme for estimating $\D$ will minimise this uncertainty, i.e., will have high statistical efficiency, and will give an accurate estimate of the uncertainty itself. We present a scheme for estimating $\D$ from a single simulation trajectory with high statistical efficiency and accurately estimating the uncertainty in the predicted value. - The statistical distribution of MSDs observable from a given simulation is modelled as a multivariate normal distribution using an analytical covariance matrix for an equivalent system of freely diffusing particles, which we parameterise from the available simulation data. + The statistical distribution of MSDs observable from a given simulation is modelled as a multivariate normal distribution using an eestimated covariance matrix, which we parameterise from the available simulation data. We then perform Bayesian regression to sample the distribution of linear models that are compatible with this model multivariate normal distribution, to obtain a statistically efficient estimate of $D^*$ and an accurate estimate of the associated statistical uncertainty. \end{abstract} @@ -138,7 +139,7 @@ \section{Introduction} \end{figure*} Here, we describe an approximate Bayesian regression method for estimating $\D$ with near-maximal statistical efficiency while also accurately estimating the corresponding statistical uncertainty, using data from a single simulation (Fig.~\ref{fig:schematic}). -We model the statistical population of simulation MSDs as a multivariate normal distribution, using an analytical covariance matrix derived for an equivalent system of freely diffusing particles, with this covariance matrix parameterised from the observed simulation data. +We model the statistical population of simulation MSDs as a multivariate normal distribution, using an estimated covariance matrix, parameterised from the observed simulation data. We then use Markov-chain Monte Carlo to sample the posterior distribution of linear models compatible with this multivariate normal model. The resulting posterior distribution provides an efficient estimate for $\D$ and allows the associated statistical uncertainty in $\Dest$ to be accurately quantified. This method is implemented in the open-source Python package \textsc{kinisi}~\cite{mccluskey_kinisi_2022}. @@ -232,41 +233,41 @@ \subsection{Background} where $\mathbf{\Sigma}$ is the observed MSD covariance matrix and $k$ is the length of the vector $\oMSD$, i.e., the number of time intervals for which we have observed MSD data. Providing this likelihood function can be calculated, we can compute the posterior distribution $\prob{\model | \oMSD}$, which gives an optimally efficient point-estimate for $\D$ and a complete description of the associated uncertainty in $\Dest$. -\subsection{Approximating $\bm{\Sigma}$ from simulation data} +\subsection{Estimating $\mathbf{\Sigma}$ from simulation data} For Bayesian regression and GLS, we require the covariance matrix for the observed MSD, $\mathbf{\Sigma}$, which is generally unknown. -To proceed, we replace $\mathbf{\Sigma}$ with a model covariance matrix, $\mathbf{\Sigma^\prime}$, with a known analytical form, that we parameterise from the available simulation data. -Providing the correlation structure of $\mathbf{\Sigma^\prime}$ is similar to that of $\mathbf{\Sigma}$, this model correlation matrix can be used in approximate Bayesian or GLS schemes to estimate the ensemble-average MSD, and hence $\D$, with high efficiency and accurate estimated uncertainties. +To proceed, we replace $\mathbf{\Sigma}$ with an estimated covariance matrix, $\widehat{\mathbf{\Sigma}}$, that is estimated directly from the simulation data. +Providing the correlation structure of $\widehat{\mathbf{\Sigma}}$ is similar to that of $\mathbf{\Sigma}$, this estimated correlation matrix can be used in approximate Bayesian or GLS schemes to estimate the ensemble-average MSD, and hence $\D$, with high efficiency and accurate estimated uncertainties. -We model the covariance matrix for the observed MSD from a given simulation using the covariance matrix for the MSD of an equivalent system of freely diffusing particles, $\mathbf{\Sigma^\prime}$. -For observed MSDs computed by averaging over numerically-independent sub-trajectories, the covariance matrix $\mathbf{\Sigma^\prime}$, in the long time limit, has elements (see SI) +We estimate the covariance matrix for the observed MSD from a given simulation using the estimates variances. +These estimated variances, $\varest{\oMSDi}$, are found using the standard result that ths variance of the mean of a sample scales inversely with the number of independent constituent observations. +Specifically, we estimate the variance $\varest{\oMSDi}$ by rescaling the observed variance of the squared displacement for time interval $i$ by the number of numerically-independent contributing sub-trajectories, $\nind{i}$; % \begin{equation} - \Sigma^\prime\left[\oMSDi, \oMSDj\right]= \Sigma^\prime\left[\oMSDj, \oMSDi\right] = - \var{\oMSDi} \frac{\nind{i}}{\nind{j}},\hspace{1em} \forall\,i \leq j, - \label{equ:cvv} -\end{equation} + \varest{\oMSDi} \approx \frac{1}{\nind{i}}\var{\Delta r_i^2}. + \label{equ:varestMSD} +\end{equation} % -where $\var{\oMSDi}$ are the time-dependent variances of the observed MSD, and $\nind{i}$ is the total number of numerically-independent observed squared-displacements for time-interval $i$, accounting for both overlapping sampling and concerted motion (see Eqn.~\ref{equ:der_var}). -We estimate the variances $\var{\oMSDi}$ using the standard result that the variance of the mean of a sample scales inversely with the number of independent constituent observations. -Specifically, we approximate the variance $\varest{\oMSDi}$ (see SI) by rescaling the observed variance of the squared displacement for time interval $i$ by the number of numerically-independent contributing sub-trajectories, $\nind{i}$; +From this, the full estimated covariance matrix, can be populated, in the long time limit (see SI) as follows, % \begin{equation} - \varest{\oMSDi} \approx \frac{1}{\nind{i}}\var{\Delta r_i^2}. - \label{equ:varestMSD} -\end{equation} + \widehat{\Sigma}\left[\oMSDi, \oMSDj\right]= \widehat{\Sigma}\left[\oMSDj, \oMSDi\right] = + \varest{\oMSDi} \frac{\nind{i}}{\nind{j}},\hspace{1em} \forall\,i \leq j, + \label{equ:cvv} +\end{equation} % +where $\nind{i}$ is the total number of numerically-independent observed squared-displacements for a random walk of independent particles in time-interval $i$. -Rescaling by the number of numerically-independent contributing sub-trajectories has the effect of renormalising the variance of the observed squared displacements to account for correlations between particle squared displacements computed from overlapping time windows. +Rescaling the variance by the number of numerically-independent contributing sub-trajectories has the effect of renormalising the variance of the observed squared displacements to account for correlations between particle squared displacements computed from overlapping time windows. An alternative approach to renormalising $\var{x_i}$ is to use a non-parametric block-averaging procedure~\cite{flyvbjerg_error_1989,frenkel_understanding_2023,materzanini_high_2021}, which gives undesirable results for a random walk (see SI). The block-averaging approach, additionally, requires numerical convergence with respect to block size, which is not guaranteed for time windows with few independent observations. We require $\var{x_i}$ at all observed time intervals, $i$, making block averaging on large data sets computationally prohibitive. -The estimated variance $\varest{\oMSD}$ can be calculated from a single simulation trajectory, and provides an accurate estimate of the true variance $\var{\oMSD}$. +The estimated variance $\varest{\oMSD}$ can be determined from a single simulation trajectory, and provides an accurate estimate of the true variance $\var{\oMSD}$. To demonstrate this, we performed \num{4096} independent simulations of \num{128} particles undergoing a three-dimensional cubic-lattice random walk of \num{128} steps per particle. Using data from all \num{4096} simulations, we first compute the true simulation MSD and its variance (Fig.~\ref{fig:msd}a). We also compute the MSD and estimated variance using data from a single simulation trajectory (Fig.~\ref{fig:msd}b), using the scheme described above. -A quantitative comparison between the true MSD variance and the single-trajectory estimated MSD variance is made in Fig.~\ref{fig:msd}c: the close numerical agreement confirms that Eqn~\ref{equ:varestMSD} can be used to estimate $\var{\oMSD}$, which can then be used to parameterise the model covariance matrix $\mathbf{\Sigma^\prime}$ via Eqn.~\ref{equ:cvv}. +A quantitative comparison between the true MSD variance and the single-trajectory estimated MSD variance is made in Fig.~\ref{fig:msd}c: the close numerical agreement confirms that Eqn~\ref{equ:varestMSD} can be used to estimate $\var{\oMSD}$, which can then be used to parameterise the estimated covariance matrix $\widehat{\mathbf{\Sigma}}$ via Eqn.~\ref{equ:cvv}. \begin{figure} \centering @@ -283,42 +284,34 @@ \subsection{Approximating $\bm{\Sigma}$ from simulation data} \label{fig:msd} \end{figure} +The practical implementation of both GLS and Bayesian regression requires that the covariance matrix $\widehat{\mathbf{\Sigma}}$ is invertible (positive definite); see Eqns.~\ref{equ:gls} and \ref{equ:loglike}. +The estimated covariance matrix however in constructed from a limited number of samples, and therefore may have a high condition number~\cite{higham_restoring_2016}. +Statistical noise and high condition numbers result in numerical sensitivity in linear equations, such as inversion. +Therefore, to make our scheme numerically tractable, the estimated covariance matrix is reconditioned using the minimum eigenvalue method~\cite{tabeart_improving_2020}. +This approach ensures that the condition number for the estimated covariance matrix is equal to a user defined parameter, that can be set to produce an invertible matrix that is suitable for GLS or Bayesian regression. -The practical implementation of both GLS and Bayesian regression requires that the covariance matrix $\mathbf{\Sigma^\prime}$ is invertible (positive definite); see Eqns.~\ref{equ:gls} and \ref{equ:loglike}. -The estimated MSD variances derived from simulation data via Eqn.~\ref{equ:varestMSD} are statistically noisy and using these to directly parameterise $\mathbf{\Sigma^\prime}$ can yield non-invertible singular matrices. -To make our scheme numerically tractable, we therefore fit our estimated MSD variances to the analytical variance for an analogous system of particles undergoing random walks~\cite{smith_random_1996}; -% -\begin{equation} - \var{\oMSDi} = a\frac{t^2_i}{\nind{i}}, - \label{equ:variance} -\end{equation} -% -where $a$ is a scaling parameter determined by fitting Eqn.~\ref{equ:variance} to the directly estimated MSD variances. -This smoothing of $\var{\oMSDi}$ guarantees that the resulting model covariance matrix $\mathbf{\Sigma^\prime}$ is invertible and thus suitable for GLS or Bayesian regression. - -To illustrate the complete numerical procedure for deriving the model covariance matrix, $\bm{\Sigma^\prime}$, we present in Fig.~\ref{fig:covariances} the MSD covariance matrix for \num{4096} random-walk simulations, as described above, at three differing levels of approximation: +To illustrate the complete numerical procedure for deriving the estimated covariance matrix, $\widehat{\mathbf{\Sigma}}$, we present in Fig.~\ref{fig:covariances} the MSD covariance matrix for \num{4096} random-walk simulations, as described above, at three differing levels of approximation: the numerically converged covariance matrix, $\bm{\Sigma}$, computed using the data from all \num{4096} simulations (Fig.~\ref{fig:covariances}a); -the corresponding analytical model covariance matrix, $\mathbf{\Sigma^\prime}$, defined by Eqn.~\ref{equ:cvv} and parametrised using analytical variances $\var{\oMSDi}$ (Fig.~\ref{fig:covariances}b); and the average model covariance matrix obtained by parametrising Eqn.~\ref{equ:cvv} using smoothed variances estimated from individual simulation trajectories, and averaging over the resulting set of all \num{4096} matrices (Fig.~\ref{fig:covariances}c). - -\begin{figure} - \centering - % \resizebox{\columnwidth}{!}{\includegraphics*{figure3.pdf}} - \includegraphics[width=\columnwidth]{covariances.pdf} - \script{covariances.py} - \caption{ - (a) The numerical MSD covariance matrix $\bm{\Sigma}$ calculated using MSD data from \num{4096} simulations of \num{128} particles undergoing a 3D lattice random walk of \num{128} steps per particle. - (b) The analytical MSD covariance matrix $\bm{\Sigma^\prime}$ (Eqn.~\ref{equ:cvv}), parametrised using analytical random-walk variances $\var{\oMSDi}$. - (c) The MSD covariance matrix obtained applying the numerical scheme described in the main text to each individual random walk simulation, averaged over all \num{4096} such simulations. - Colour bars in (a--c) show the covariance, $\Sigma\left[\oMSDi, \oMSDj\right]$. - The off-diagonal panels show difference plots, computed as per-element ratios between pairs of covariance matrices (a--c). - %Raw simulation data~\cite{mccluskey_zenodo_2022} and scripts~\cite{mccluskey_github_2022} to generate this figure are available under CC BY 4.0/MIT licences. - } - \label{fig:covariances} -\end{figure} - +the corresponding analytical covariance matrix, $\mathbf{\Sigma^\prime}$, derived in the SI (Fig.~\ref{fig:covariances}b); and the average estimated covariance matrix obtained by parametrising Eqn.~\ref{equ:cvv} using smoothed variances estimated from individual simulation trajectories, and averaging over the resulting set of all \num{4096} matrices (Fig.~\ref{fig:covariances}c). + +% \begin{figure} +% \centering +% % \resizebox{\columnwidth}{!}{\includegraphics*{figure3.pdf}} +% \includegraphics[width=\columnwidth]{covariances.pdf} +% \script{covariances.py} +% \caption{ +% (a) The numerical MSD covariance matrix $\bm{\Sigma}$ calculated using MSD data from \num{4096} simulations of \num{128} particles undergoing a 3D lattice random walk of \num{128} steps per particle. +% (b) The analytical MSD covariance matrix $\bm{\Sigma^\prime}$ (Eqn.~\ref{equ:cvv_SI}). +% (c) The MSD covariance matrix obtained applying the numerical scheme described in the main text to each individual random walk simulation, averaged over all \num{4096} such simulations. +% Colour bars in (a--c) show the covariance, $\Sigma\left[\oMSDi, \oMSDj\right]$. +% The off-diagonal panels show scale plots, computed as per-element ratios between pairs of covariance matrices (a--c). +% %Raw simulation data~\cite{mccluskey_zenodo_2022} and scripts~\cite{mccluskey_github_2022} to generate this figure are available under CC BY 4.0/MIT licences. +% } +% \label{fig:covariances} +% \end{figure} While the analytical and average estimated covariance matrices show some systematic deviation from the numerically converged covariance matrix, the general correlation structure is preserved. -The discrepancy between the model and numerical covariance matrices largely stems from the approximation made in deriving the analytical form that $t$ is large, which leads to an overestimation of the variance at low $t$. +The discrepancy between the estimated and numerical covariance matrices largely stems from the assumption that the analytical covariance structure is equivalent for the estimated covariance matrix. Despite this, the average estimated covariance matrix reproduces well the correlation structure of the true numerical covariance matrix, indicating that the covariance matrices estimated from individual simulation trajectories may be used within approximate GLS or Bayesian regression schemes to estimate $\D$ and $\var{\Dest}$. \subsection{Validation} @@ -349,7 +342,7 @@ \subsection{Validation} Second, we consider an example real-world system---the lithium-ion solid electrolyte \ce{Li7La3Zr2O12} (LLZO)---which represents an application of our method to a well-studied material of practical interest for solid-state lithium-ion batteries \cite{MuruganEtAl_AngewChemIntEd2007,burbano_sparse_2016,morgan_lattice_2017,SquiresEtAl_PhysRevMater2022}. Fig.~\ref{fig:random_walk}a shows the observed MSD from a single 3D-lattice random-walk simulation, along with the estimated posterior distribution of linear models compatible with the observed MSD data, $\prob{\model|\oMSD}$, calculated via Eqns.~\ref{equ:bayes} and \ref{equ:loglike}. -The corresponding marginal posterior distribution of estimated diffusion coefficients $\prob{\D|\oMSD}$ is shown in Fig.~\ref{fig:random_walk}b; this distribution is approximately Gaussian and is centred on the true self-diffusion coefficient $\D = \num{1}$, demonstrating that for this example trajectory we obtain a good point-estimate of $\D$. +The corresponding marginal posterior distribution of estimated diffusion coefficients $\prob{\D|\oMSD}$ is shown in Fig.~\ref{fig:random_walk}b; this distribution is approximately Gaussian and close to the true self-diffusion coefficient $\D = \num{1}$, demonstrating that for this example trajectory we obtain a good point-estimate of $\D$. To evaluate the overall performance of our method, we repeat our analysis on the full set of \num{4096} random-walk simulations. Fig.~\ref{fig:random_walk}c presents a histogram of the resulting point estimates of $\D$, with each estimate derived as the mean of the posterior distribution $\prob{\D|\oMSD}$ using input data from an individual simulation. @@ -357,43 +350,43 @@ \subsection{Validation} This latter distribution represents the distribution of ``best possible'' estimates of $\D$ and exhibits the minimum possible theoretical variance. The close agreement between these two distributions demonstrates that our approximate Bayesian regression scheme yields nearly optimal estimates of $\D$ using data from individual simulations. The distribution of estimated diffusion coefficients from single simulations is slightly broader than the exact numerical results. -This minor deviation is a consequence of the overestimation of $\varest{\oMSDi}$ at short times, noted above, which results from our use of the long-time limit in the derivation of the analytical model covariance matrix. +This minor deviation is a consequence of the noise in data collected from a single simulation trajectory. We next consider the degree to which our method can quantify the uncertainty in $\Dest$ when using input data from a single simulation. Fig.~\ref{fig:random_walk}d shows the distribution of estimated variances $\varest{\Dest}$, with each sample calculated from an individual simulation trajectory. We also show the true variance of individual point estimates, $\var{\Dest}$, which characterises the spread of the histogram in Fig.~\ref{fig:random_walk}c. -The distribution of estimated variances is biased relative to the true variance and skewed, due to numerical differences between the true covariance matrix $\mathbf{\Sigma}$ and the model covariance matrix $\mathbf{\Sigma^\prime}$ (further details are provided in the SI). +The distribution of estimated variances is biased relative to the true variance, due to numerical differences between the true covariance matrix $\mathbf{\Sigma}$ and the estimated covariance matrix $\widehat{\mathbf{\Sigma}}$ (further details are provided in the SI). In general, however, the distribution of the estimated variance shows good agreement with the true sample variance. Notably, the precision of this estimate is significantly greater than obtained using OLS or WLS and their corresponding textbook statistical formulae. We next benchmark our method using data from simulations of the cubic lithium-ion solid electrolyte \ce{Li7La3Zr2O12} (c-LLZO). We performed a single simulation of \num{1536} atoms (\num{448} Li ions) at \SI{700}{K} for \SI{6}{\nano\second} (full simulation details are provided in the Methods section). -To generate multiple statistically equivalent trajectories, the resulting simulation data was partitioned into \num{480} effective trajectories, each approximately $\sim\SI{100}{\ps}$ in length, and containing data for \num{56} lithium ions. +To generate multiple statistically equivalent trajectories, the resulting simulation data was partitioned into \num{192} effective trajectories, each approximately $\sim\SI{500}{\ps}$ in length, and containing data for \num{28} lithium ions, which were selected randomly without replacement. We then perform the same approximate Bayesian regression analysis as above on each effective trajectory, excluding the first \SI{10}{ps} of MSD data in each case to remove short-time data corresponding to ballistic and sub-diffusive regimes \cite{he_statistical_2018}. -The resulting distribution of the point estimates, $\Dest$, from analysis of all \num{512} effective trajectories is shown in Fig.~\ref{fig:diffusion}a. +The resulting distribution of the point estimates, $\Dest$, from analysis of all \num{420} effective trajectories is shown in Fig.~\ref{fig:diffusion}a. Again, the corresponding distribution of $\Dest$ estimates derived using Bayesian regression and a well-converged numerical covariance matrix calculated from the full LLZO dataset is also shown for comparison. -The distribution $\prob{\Dest}$ obtained using the model covariance matrix and parametrised separately for each individual effective simulation is highly similar to that obtained using the aggregate numerical covariance matrix calculated from the complete simulation dataset. +The distribution $\prob{\Dest}$ obtained using the estimated covariance matrix and parametrised separately for each individual effective simulation is highly similar to that obtained using the aggregate numerical covariance matrix calculated from the complete simulation dataset. This close agreement mirrors the results for our random walk simulations (see the SI for a similar comparison of the OLS, WLS, and GLS as shown in Fig.~\ref{fig:glswlsols}), and confirms that our method yields accurate and statistically efficient estimates for $\D$, even for real-world simulation data. -\begin{figure}[tb] - \centering - % \resizebox{\columnwidth}{!}{\includegraphics*{figure5.pdf}} - \includegraphics[width=\columnwidth]{diffusion.pdf} - \script{diffusion.py} - \caption{ - (a) Probability distribution of point estimates $\prob{\Dest}$ for \num{512} effective simulations of LLZO (orange histogram). - The grey line shows the distribution $\prob{\Dest_\mathrm{num}}$ obtained using Bayesian regression with the complete LLZO dataset as input. - The pink bar shows an interval of one standard deviation $\sigma[\prob{\Dest}]$. - (b) Probability distribution of estimated variances, $\varest{\Dest}$, for individual LLZO effective simulations, compared to the true sample variance (pink vertical line) $\var{\Dest}$. - %Raw simulation data~\cite{mccluskey_zenodo_2022} and scripts~\cite{mccluskey_github_2022} to generate this figure are available under CC BY 4.0/MIT licences. - } - \label{fig:diffusion} -\end{figure} +% \begin{figure}[tb] +% \centering +% % \resizebox{\columnwidth}{!}{\includegraphics*{figure5.pdf}} +% \includegraphics[width=\columnwidth]{diffusion.pdf} +% \script{diffusion.py} +% \caption{ +% (a) Probability distribution of point estimates $\prob{\Dest}$ for \num{480} effective simulations of LLZO (orange histogram). +% The grey line shows the distribution $\prob{\Dest_\mathrm{num}}$ obtained using Bayesian regression with the complete LLZO dataset as input. +% The pink bar shows an interval of one standard deviation $\sigma[\prob{\Dest}]$. +% (b) Probability distribution of estimated variances, $\varest{\Dest}$, for individual LLZO effective simulations, compared to the true sample variance (pink vertical line) $\var{\Dest}$. +% %Raw simulation data~\cite{mccluskey_zenodo_2022} and scripts~\cite{mccluskey_github_2022} to generate this figure are available under CC BY 4.0/MIT licences. +% } +% \label{fig:diffusion} +% \end{figure} We also consider the probability distribution of estimates of the variance in $\Dest$ calculated for each effective trajectory (Fig.~\ref{fig:diffusion}b), which we compare to the true variance in $\Dest$ for our method; i.e., the variance of the histogram in Fig.~\ref{fig:diffusion}a. While the estimated variances deviate somewhat from the true distribution $\prob{\var{\Dest}}$, the agreement is reasonable and mirrors our results for the random walk simulations. -Hence, our method provides reasonably accurate estimates of the uncertainty in $\Dest$ for our c-LLZO dataset, even when applied to single effective trajectories with limited displacement data (only 56 mobile ions, and \SI{25}{ps} simulation length). +Hence, our method provides reasonably accurate estimates of the uncertainty in $\Dest$ for our c-LLZO dataset, even when applied to single effective trajectories with limited displacement data (only 32 mobile ions, and \SI{2}{ns} simulation length). \subsection{$\var{\Dest}$ scaling and comparison to OLS, WLS, and GLS} @@ -428,15 +421,16 @@ \section{Summary and Discussion} We consider the observed mean-squared displacement data from a single simulation as a random sample, $\bm{X}$, from a population of potential MSDs generated by equivalent replica simulations, $\bm{X}\sim p(\oMSD)$. We model this population using a multivariate normal distribution, $p(\oMSD) = \mathcal{N}(\model, \mathbf{\Sigma})$, with mean vector $\model = 6\D\bm{t} + c$, where $\D$ and $c$ are model parameters to be determined. -To model the covariance matrix, we use an analytical solution derived for an equivalent system of freely diffusing particles. -To parameterise this model covariance matrix, we rescale the variance of the observed squared displacements from the input simulation trajectory, followed by a smoothing step to ensure a positive-definite matrix. -The resulting model covariance matrix preserves the correlation structure of the true simulation MSD covariance matrix, and gives a multivariate normal model for the population of observable simulation MSDs that depends solely on the model parameters, $\D$ and $c$. +To estimate the covariance matrix, we rescale the variance of the observed squared displacements from the input simulation trajectory and use these to populate then full covariance matrix based on the relationship between the variances and covariance matrix for an equivalent system of freely diffusing particles. +The resulting estimated covariance matrix preserves the correlation structure of the true simulation MSD covariance matrix, and gives a multivariate normal model for the population of observable simulation MSDs that depends solely on the model parameters, $\D$ and $c$. We use Markov-Chain Monte Carlo to sample the posterior distribution of linear models compatible with the observed MSD data. This approach yields a marginal posterior distribution, $\prob{\D | \oMSD}$, that gives a statistically efficient point estimate for $\D$ and allows the associated statistical uncertainty, $\var{\Dest}$, to be quantified. We have benchmarked our approach using simulation data for an ideal 3D lattice random walk and for the lithium-ion solid electrolyte \ce{Li7La3Zr2O12} (LLZO). In both cases, we obtain a distribution of estimates for $\D$ that closely matches the theoretically optimal distribution obtained using a numerical covariance matrix derived from a large number of replica simulation trajectories. +The assumption of freely diffusing particles means may overlook short-time interval effects that are present in real systems, such as correlations between different atoms. +Such correlations would reduce the value of $nind{i}$, however, the impact of this is anticipated to be minimal and potentially limited to the ballistic regime. We obtain estimates for $\D$ that are unbiased, with near-optimal statistical efficiency, using input data from single simulation trajectories. The approximate Bayesian regression scheme therefore provides more accurate single-point estimates of the self-diffusion coefficient than the commonly used OLS or WLS methods, when applied to the same input simulation data. @@ -458,7 +452,8 @@ \section{Methods} \subsection{Numerical implementation in KINISI} \label{sec:implementation} -We have implemented the approximate Bayesian estimation method described in the main text in the open-source Python package \textsc{kinisi}~\cite{mccluskey_kinisi_2022}, under the MIT licence. +We have implemented the approximate Bayesian estimation method described in the main text in the open-source Python package \textsc{kinisi}~\cite{mccluskey_kinisi_2022}, under the MIT licence. +Throughout this work, kinisi-[INSERT VERSION NUMBER] was used for all analyses. \textsc{kinisi} uses overlapping sliding window sampling when calculating the observed mean squared displacement at each time interval $t$ (see Eqn.~\ref{equ:observed_msd}). For a given time interval, $t$, the maximum number of observations is $N_{\mathrm{atoms}} \times (N_{t} - i)$ displacements, where $N_{\mathrm{atoms}}$ is the number of atoms, $N_{t}$ is the total number of timesteps, and $i$ is the index of the timestep (where \num{1} is the index for the shortest timestep). @@ -466,7 +461,8 @@ \subsection{Numerical implementation in KINISI} The parametrisation of the covariance matrix from the variances $\var{\oMSDi}$ and the number of independent observations $\nind{i}$ is defined by Eqn.~\ref{equ:cvv}. The covariance matrix is only constructed for values of $t$ where the particle motion is considered to be in the long-time diffusive limit, with this threshold set by the user to a value appropriate for their system and simulation data. -For the examples presented in the main manuscript, we consider particles to be in the diffusive regime from $t=\num{4}$ for the random walk trajectories and from $t=\SI{10}{\pico\second}$ for the LLZO simulations. +For the examples presented in the main manuscript, we consider particles to be in the diffusive regime from $t=\num{2}$ for the random walk trajectories and from $t=\SI{10}{\pico\second}$ for the LLZO simulations. +The covariance matrix is reconditioned using the minimum eigenvalue method~\cite{tabeart_improving_2020} for all simulations, the maximum condition number was \num{1e16}. \textsc{kinisi} uses ordinary least squares to obtain an initial guess for the gradient and intercept of the linear model describing the observed MSD. This initial guess is then used as the starting point for minimising the negative maximum a posteriori (the peak of the posterior distribution as per Eqn.~\ref{equ:bayes}), with the improper prior that $\D \ge 0$~\cite{broyden_convergence_1970,fletcher_new_1970,goldfarb_family_1970,shanno_conditioning_1970}. @@ -482,7 +478,7 @@ \subsection{LLZO simulations} Classical molecular dynamics were run using the \textsc{metalwalls} code~\cite{marin_metalwalls_2020}. We used the DIPPIM polarisable ion force field, as parameterised by Burbano \emph{et al.}~\cite{burbano_sparse_2016}, due to its proven accuracy in accounting for the effect of ion polarisability on diffusion~\cite{wilson_polarization_1993,burbano_sparse_2016}. We simulated the cubic phase of LLZO in NVT ensemble at a temperature of \SI{700}{\kelvin}. -Simulations were run for \SI{6}{\nano\second} with a \SI{2}{\femto\second} timestep. +Simulations were run for \SI{6}{\nano\second} with a \SI{0.5}{\femto\second} timestep. To control temperature, we used a Nos\'{e}-Hoover thermostat, with a relaxation time of \SI{121}{\femto\second} (5000 $\hbar / E_{h}$)~\cite{nose_unified_1984,hoover_canonical_1985,martyna_nose_1992}. Simulations were performed using $2 \times 2 \times 2$ supercells with \SI{1536}{atoms} following the same protocol as in Ref.~\citenum{burbano_sparse_2016}. @@ -555,14 +551,14 @@ \section{Derivation of the long-time limit covariance matrix for a system of fre \label{sec:ran} In the main manuscript we present the result that the covariance matrix for a system of freely diffusing particles, in the long-time limit, has the form \begin{equation} - \Sigma^\prime\left[\oMSDi, \oMSDj\right]= \Sigma^\prime\left[\oMSDj, \oMSDi\right] = - \var{\oMSDi} \frac{\nind{i}}{\nind{j}},\hspace{1em} \forall\,i \leq j, + \widehat{\Sigma}\left[\oMSDi, \oMSDj\right]= \widehat{\Sigma}\left[\oMSDj, \oMSDi\right] = + \varest{\oMSDi} \frac{\nind{i}}{\nind{j}},\hspace{1em} \forall\,i \leq j, \label{equ:cvv_SI} \end{equation} where $\oMSDi$ is the observed mean-squared displacement (MSD) for time interval $i$ and $\nind{i}$ is the number of statistically independent observed squared displacements averaged over to compute the mean value. To derive this result, we first present a derivation of the expected variance for the MSD at timestep $i$, $\var{\oMSD}$, following the approach of Smith and Gillan~\cite{smith_random_1996}. -We then derive an expression for the covariance $\Sigma^\prime\left[\oMSDi, \oMSDj\right]$ to obtain the result above. +We then derive an expression for the covariance $\widehat{\Sigma}\left[\oMSDi, \oMSDj\right]$ to obtain the result above. For a single particle undergoing a one-dimensional random walk with step size $\kappa$, each step gives a displacement $h = \pm \kappa$. After $n$ steps, the MSD, $\oMSDn$, is given by @@ -645,16 +641,10 @@ \section{Derivation of the long-time limit covariance matrix for a system of fre \var{\oMSDn} = \frac{2n^2\kappa^4}{\nind{n}}, \label{equ:der_var} \end{equation} -where $\nind{n}$ is the total number of statistically independent (non-overlapping and accounting for concerted motion) squared displacements that contribute to $\oMSDi$. -In the long-time limit, $\nind{n}$ is given by the product of the number of mobile particles, $N_{\textrm{atoms}}$, the number of numerically-independent sub-trajectories of lenth $i$ in our simulation trajectory, $N_t(t)$, and the ratio of the mean-squared displacement and the mean-square of the displacement of the centre of mass for all mobile particles, -\begin{equation} - H_R(t_{\textrm{diff}}) = \frac{\sum^{N(t_{\textrm{diff}})}_{j=1}{\left[\Delta\mathbf{r}_j(t_{\textrm{diff}})\right]}^2}{\sum^{N_{t}(t_{\textrm{diff}})}_{l=1}\left|\left[\sum^{N_{\textrm{atoms}}}_{k=1}\Delta r_{k,l}(t_{\textrm{diff}})\right]^2\right|^2}, -\end{equation} -at the start of the diffusive regime, $t_{\textrm{diff}}$. +where $\nind{n}$ is the total number of statistically independent (non-overlapping) squared displacements that contribute to $\oMSDi$. +In the long-time limit, $\nind{n}$ is given by the product of the number of mobile particles and the number of numerically-independent sub-trajectories of length $i$ in our simulation trajectory. Note that $\nind{n}$ considers numerically-independent sub-trajectories, since mutually overlapping time-windows give correlated squared displacements. Where overlapping time-windows are used, Eqn.~\ref{equ:der_var} approximates the observed variance (the variance for overlapping and non-overlapping samples are not the same) with accuracy that increases as a function of time-interval length. -The inclusion of $H_R(t_{\textrm{diff}})$ account for concerted motion between individual particles in the simulation~\cite{vargas_dynamic_2019,canepa_pushing_2022}, which is particularly important in solid and liquid electrolyte systems. -The selection of $t_{\textrm{diff}}$ as the point at which to use $H_R$ is taken to maximise the number of observations of the centre-of-mass displacement. This leads to the approximation in Eqn.~\ref{equ:varestMSD}. The results for a one-dimensional lattice above (Eqns.~\ref{equ:msd_SI} \&~\ref{equ:der_var}) can be extended to a $d$-dimensional lattice, to give @@ -664,6 +654,7 @@ \section{Derivation of the long-time limit covariance matrix for a system of fre with variance \begin{equation} \var{\oMSDn}_d = \sum^d{\frac{2n^2\kappa^4}{d^2\nind{n}}} = \frac{2n^2\kappa^4}{d\nind{n}}, + \label{equ:analvar} \end{equation} Because each step is equally likely to move a particle along each of the $d$ dimensions, the term $n$ in Eqns.~\ref{equ:msd_SI} \&~\ref{equ:der_var} is replaced here with $n/d$. @@ -776,16 +767,17 @@ \section{Derivation of the long-time limit covariance matrix for a system of fre \section{Comparison With Block Averaging Approach} \label{sec:block} Eqn.~\ref{equ:varestMSD} enables the approximate estimation of the variance in $\oMSDi$, as mentioned in the main text, it is also possible to obtain this from a block averaging approach~\cite{flyvbjerg_error_1989,frenkel_understanding_2023}. -Fig.~\ref{fig:pyblock} compares the $\D$ estimate distributions using variances from approximated by Eqn.~\ref{equ:varestMSD} with those from the block averaging approach (using the \textsc{pyblock} Python package~\cite{spencer_pyblock_2020}) both with the fitting of the variances to Eqn.~\ref{equ:variance} and without. -When the fitting is not performed, the resulting covariance matrix is numerically unstable, leading to estimated values of $\D$ at extreme values. -While when the fitting is performed and the blocking-estimated variances as used, there is no improvement in the estimation of $\D$ and the resulting distribution of estimates in the variance of $\prob{\D}$ is broader than when Eqn.~\ref{equ:varestMSD} is used. +Fig.~\ref{fig:pyblock} compares the $\D$ estimate distributions using variances from approximated by Eqn.~\ref{equ:varestMSD} with those from the block averaging approach (using the \textsc{pyblock} Python package~\cite{spencer_pyblock_2020}). +The estimates of $\D$ are similar between the two approaches, however with pyblocking approach typically produces covariance matrices that cannot be conditioned easily, leading to numerical instabilities -- leading to a large standard deviation. +The more significant difference is in the estimated variances, where the block averaging approach consistently underestimates the variance in $\Dest$. +Generally speaking, any approach to estimate the variance of a quantity should prefer to overestimate the variance, such that undue confidence is not placed in the results. \begin{figure} \centering % \resizebox{\columnwidth}{!}{\includegraphics*{figure1.pdf}} \includegraphics[width=\columnwidth]{pyblock.pdf} \script{pyblock.py} - \caption{Comparison of the distribution of $\Dest$ (a, c, \& e) and the estimated variances (b, d, \& f) from 4096 individual random walk simulations using variances approximated by Eqn.~\ref{equ:varestMSD} (a \& b), from block averaging with fitting to Eqn.~\ref{equ:variance} (c \& d), and without fitting (e \& f). Note that the bounds in a, c, \& e are defined by the extrema of the estimated values, showing the fact that without the fitting to Eqn.~\ref{equ:variance}, the resulting covariance matrix is numerically unstable.} + \caption{Comparison of the distribution of $\Dest$ (a \& c) and the estimated variances (b \& d) from 4096 individual random walk simulations using variances approximated by Eqn.~\ref{equ:varestMSD} (a \& b), from block averaging (c \& d).} \label{fig:pyblock} \end{figure} @@ -821,16 +813,16 @@ \section{Evaluation of OLS, WLS, and GLS for LLZO System} In Fig.~\ref{fig:glswlsols}, it is shown that for a 3D lattice random walk the heteroscedastic and correlated nature of the data requires generalised least squares to optimally estimate the true ensemble-average MSD, and therefore estimate $\D$. This is also true for the real materials, such as the \ce{Li7La3Zr2O12} (LLZO) investigated by classical molecular dynamics simulation in this work (Fig.~\ref{fig:glswlsols_llzo}). Once again, GLS (or, indeed, the Bayesian regression equivalent, shown in Fig.~\ref{fig:diffusion}) offers both the most statistically efficient estimation of $\D$ and is required to obtain an accruate estimate of the statistical uncertainty in $\D$ from a single simulation. -\begin{figure} - \centering - % \resizebox{\columnwidth}{!}{\includegraphics*{figure1.pdf}} - \includegraphics[width=\columnwidth]{glswlsols_llzo.pdf} - \script{glswlsols_llzo.py} - \caption{ - Example distributions of estimated self-diffusion coefficients, $\Dest$, calculated using (a) ordinary least squares (OLS), (b) weighted least squares (WLS), and (c) generalised least squares (GLS), - from MSD data from \num{512} effective simulations of LLZO of \SI{\sim 25}{\pico\second} with 56 lithium ions. - In each panel, the grey curve shows the best-fit normal distribution for the simulation data, the upper horizontal bar shows the standard deviation of this distribution, and the lower horizontal bar shows the average estimated standard distribution given by the analytical expression for $\sigma[\prob{\Dest}]$ for each regression method.} - \label{fig:glswlsols_llzo} -\end{figure} +% \begin{figure} +% \centering +% % \resizebox{\columnwidth}{!}{\includegraphics*{figure1.pdf}} +% \includegraphics[width=\columnwidth]{glswlsols_llzo.pdf} +% \script{glswlsols_llzo.py} +% \caption{ +% Example distributions of estimated self-diffusion coefficients, $\Dest$, calculated using (a) ordinary least squares (OLS), (b) weighted least squares (WLS), and (c) generalised least squares (GLS), +% from MSD data from \num{512} effective simulations of LLZO of \SI{\sim 25}{\pico\second} with 56 lithium ions. +% In each panel, the grey curve shows the best-fit normal distribution for the simulation data, the upper horizontal bar shows the standard deviation of this distribution, and the lower horizontal bar shows the average estimated standard distribution given by the analytical expression for $\sigma[\prob{\Dest}]$ for each regression method.} +% \label{fig:glswlsols_llzo} +% \end{figure} \end{document} diff --git a/zenodo.yml b/zenodo.yml index 2ab49b9..1f05b37 100644 --- a/zenodo.yml +++ b/zenodo.yml @@ -1,4 +1,4 @@ cache: main: - sandbox: 10.5072/zenodo.22971 + sandbox: 10.5072/zenodo.32039