Latest version

arm61 · Jul 30, 2024 · e690e20 · e690e20
1 parent 9b08308
commit e690e20
Show file tree

Hide file tree

Showing 27 changed files with 798 additions and 1,918 deletions.
diff --git 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
diff --git 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,58 +368,60 @@ 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'
             cache:
                 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"
-
+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
@@ -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
+      - git+https://github.com/bjmorgan/kinisi.git
diff --git 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
@@ -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
@@ -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)