Before changing model to include transitions of SPORE to DEAD

.gitignore

@@ -3,4 +3,6 @@ venv/
 .vscode
 simulation_claudius.py
 OLD*
-.DS_Store
+.DS_Store
+*.backup
+*.dat

extinction_time.py

@@ -5,6 +5,7 @@
 
 
 BASE_PATH = "./data/dormant-life/extinction-time"
+os.makedirs(BASE_PATH, exist_ok=True)
 
 
 def find_extinction_time(ca: CellularAutomaton, t_max: int,

lifetime_distribution.py

@@ -0,0 +1,114 @@
+import sys, os
+import numpy as np
+from collections import Counter
+from gol import CellularAutomaton, DormantLife
+from gol import ALIVE, DEAD, DORM
+from util import random_init_grid
+import matplotlib.pyplot as plt
+
+
+def lifetime_distribution(state: int, dl: DormantLife, t_max: int, t_trans: int,
+                          ignore_transient_dynamics: bool = True) -> Counter:
+    """
+    Generate the lifetime distribution for a given state. The given DormantLife
+    is stepped until t_max (skipping the transient time t_trans) and the
+    distibution how long any cell in the grid stayed in the state state before
+    transitioning is returned as a counter. Optionally, ignore cells that
+    transition to state within the transient time.
+    """
+    lifetime_grid = np.zeros((dl.N, dl.N), dtype=np.intc)
+    dl.step_until(t_trans)
+    # Mask to exclude cells that are state from the get-go
+    if ignore_transient_dynamics:
+        trans_mask = dl.grid != state
+    else:
+        trans_mask = np.ones((dl.N, dl.N), dtype=np.bool_)
+    data = Counter() # initialize distribution
+    while dl.t <= t_max:
+        old_grid = dl.grid.copy()
+        new_grid = dl.step()
+        # Cells that are *still* in state add to lifetime grid
+        lifetime_grid += (trans_mask & 
+                          ((old_grid == state) & (new_grid == state)))
+        # Add lifetimes of cells that left state and reset lifetime grid
+        _mask = (trans_mask & 
+                 ((old_grid == state) & (new_grid != state)))
+        data.update(lifetime_grid[_mask])
+        lifetime_grid[_mask] = 0
+        # Update transient mask: set mask values to true for cells initially in
+        # state that left state.
+        trans_mask |= (np.logical_not(trans_mask) & (new_grid != state)) # FIXME
+    return data
+
+
+def tau_distribution(dl: DormantLife, t_max: int, t_trans: int,
+                     ignore_transient_dynamics: bool = True) -> Counter:
+    tau_grid = np.zeros((dl.N, dl.N), dtype=np.intc)
+    dl.step_until(t_trans)
+    if ignore_transient_dynamics:
+        trans_mask = dl.grid != DORM
+    else:
+        trans_mask = np.ones((dl.N, dl.N), dtype=np.bool_)
+    data = Counter()
+    while dl.t <= t_max:
+        old_grid = dl.grid.copy()
+        new_grid = dl.step()
+        # Count tau
+        tau_grid += (trans_mask & ((old_grid == DORM) & (new_grid == DORM)))
+        # Update distribution
+        data.update(tau_grid[(trans_mask & (new_grid == DORM))])
+        # Reset tau grid
+        # NOTE: Could also add trans_mask & in front, but since tau grid is not updated for cells that are excluded through trans_mask it does not make a difference
+        tau_grid[(old_grid == DORM) & (new_grid != DORM)] = 0
+        # Update trans_mask
+        trans_mask |= (np.logical_not(trans_mask) & (new_grid != DORM))
+    return data
+
+
+if __name__ == "__main__":
+    # # CHANGE N
+    # t_max = 10_000
+    # fig,(ax, axx)=plt.subplots(figsize=(3*7.2, 3*3.2), ncols=2)
+    # alpha = .4
+    # for N in [30, 300]:
+    #     init_grid = random_init_grid(N, seed=100)
+    #     dl = DormantLife(init_grid, alpha, seed=100)
+    #     distr = lifetime_distribution(DORM, dl, t_max, 150, ignore_transient_dynamics=1)
+    #     ax.scatter(distr.keys(), distr.values(), label=N)
+    #     dl = DormantLife(init_grid, alpha, seed=100)
+    #     axx.plot(dl.dorm_count_time_series(t_max))
+    # # ax.set_xscale("log")
+    # ax.set_yscale("log")
+    # ax.legend()
+    # axx.legend()
+
+
+    # # CHANGE ALPHA
+    # t_max = 10_000
+    # fig,ax=plt.subplots(ncols=2)
+    # for alpha in [1, .1]:
+    #     init_grid = random_init_grid(100, seed=100)
+    #     dl = DormantLife(init_grid, alpha, seed=100)
+    #     distr = lifetime_distribution(DORM, dl, t_max, 150, ignore_transient_dynamics=1)
+    #     ax[0].scatter(distr.keys(), distr.values(), label=alpha)
+    #     dl = DormantLife(init_grid, alpha, seed=100)
+    #     ax[1].plot(dl.dorm_count_time_series(t_max))
+    #     print(sum(list(distr.values())))
+    # # ax[0].set_xscale("log")
+    # ax[0].set_yscale("log")
+    # ax[0].legend()
+    # ax[1].legend()
+
+
+    # TAU DISTRIBUTION
+    t_max, t_trans = 10_000, 200
+    fig, ax = plt.subplots()
+    ax.set(yscale="log")
+    for alpha in [1, .8, .5, .2]:
+        init_grid = random_init_grid(100, seed=100)
+        dl = DormantLife(init_grid, alpha, seed=100)
+        distr = tau_distribution(dl, t_max, t_trans)
+        vals, counts = np.fromiter(distr.keys(), dtype=np.intc), np.fromiter(distr.values(), dtype=np.intc)
+        ax.scatter(vals, counts/(t_max-t_trans), label=alpha)
+    ax.legend()
+    plt.show()

paper.mplstyle

@@ -35,7 +35,7 @@ ytick.right : True
 # Set line widths
 axes.linewidth : 0.5
 grid.linewidth : 0.5
-lines.linewidth : 1.
+lines.linewidth : 2.
 
 # Remove legend frame
 legend.frameon : False

simulation.py

@@ -3,16 +3,18 @@
 import numpy as np
 from matplotlib.colors import ListedColormap, BoundaryNorm
 from gol import GameOfLife, DormantLife, DEAD, ALIVE, DORM
+from lifetime_distribution import lifetime_distribution
 
 if __name__ == "__main__":
-    init_grid = np.random.choice([0, 1], p=[0.80, 0.20], size=[30, 30])
-    # init_grid = np.array([
-    #     [0, 0, 0],
-    #     [1, 1, 0],
-    #     [1, 1, 0]
-    # ])
+    # init_grid = np.random.choice([0, 1], p=[0.80, 0.20], size=[30, 30])
+    init_grid = np.array([
+        [0, 0, 1, 0],
+        [0, 1, 0, 0],
+        [1, 0, 0, 0],
+        [1, 0, 0, 0]
+    ])
     gol = GameOfLife(init_grid)
-    dl = DormantLife(init_grid, alpha=1)
+    dl = DormantLife(init_grid, alpha=.3)
     colors = ["white", "tab:blue", "lightblue"]
     cmap = ListedColormap(colors)
     fig, ax = plt.subplots(figsize=(7.2, 3.2), ncols=2)
@@ -21,11 +23,12 @@
     mat_gol = ax[0].matshow(gol.grid, cmap=cmap, vmin=0, vmax=2)
     mat_dl = ax[1].matshow(dl.grid, cmap=cmap, vmin=0, vmax=2)
     def update(frame):
+        print(dl.p_grid)
         mat_gol.set_data(gol.grid)
         mat_dl.set_data(dl.grid)
         ax[0].set(title=r"Game of Life: $N_\text{alive} = %d$"%gol.alive_count)
         ax[1].set(title=r"Dormant Life: $N_\text{alive} = %d$"%dl.alive_count)
         gol.step()
         dl.step()
-    ani = animation.FuncAnimation(fig, update, interval=10, save_count=100)
+    ani = animation.FuncAnimation(fig, update, interval=1000, save_count=100)
     plt.show()

state_transitions.py

@@ -8,6 +8,7 @@
 
 
 BASE_PATH = "data/dormant-life/state-transitions/to-alive-transitions"
+os.makedirs(BASE_PATH, exist_ok=True)
 
 
 def count_transitions(pre_grid: np.ndarray, post_grid: np.ndarray,

tests.py

@@ -1,7 +1,7 @@
 import unittest
 import numpy as np
 from gol import GameOfLife, DormantLife, ALIVE, DORM, DEAD
-
+from lifetime_distribution import lifetime_distribution
 
 class TestDormantLife(unittest.TestCase):
     def test_ALIVE_DORM_conversion(self):
@@ -85,18 +85,42 @@ def test_stochastic_update(self):
         ])
         np.testing.assert_array_equal(grid_step, res)
 
-    def test_transitions(self):
+    # def test_transitions(self):
+    #     test_grid = np.array([
+    #         [DEAD, ALIVE, DORM],
+    #         [DEAD, ALIVE, DORM],
+    #         [DEAD, DEAD, DEAD]
+    #     ])
+    #     gol = DormantLife(test_grid)
+    #     grid_step = gol.step()
+    #     self.assertEqual(gol.transitions_from(test_grid, ALIVE, DORM), 2)
+    #     self.assertEqual(gol.transitions_from(test_grid, DORM, ALIVE), 2)
+    #     self.assertEqual(gol.transitions_from(test_grid, ALIVE, DEAD), 0)
+    #     self.assertEqual(gol.transitions_from(test_grid, DEAD, ALIVE), 0)
+
+
+class TestLifetimeDistribution(unittest.TestCase):
+    def test_lifetime_measuring(self):
         test_grid = np.array([
-            [DEAD, ALIVE, DORM],
-            [DEAD, ALIVE, DORM],
-            [DEAD, DEAD, DEAD]
+                [DEAD, ALIVE, DEAD],
+                [DEAD, ALIVE, ALIVE],
+                [DEAD, DEAD, DEAD]
         ])
-        gol = DormantLife(test_grid)
-        grid_step = gol.step()
-        self.assertEqual(gol.transitions_from(test_grid, ALIVE, DORM), 2)
-        self.assertEqual(gol.transitions_from(test_grid, DORM, ALIVE), 2)
-        self.assertEqual(gol.transitions_from(test_grid, ALIVE, DEAD), 0)
-        self.assertEqual(gol.transitions_from(test_grid, DEAD, ALIVE), 0)
+        dl = DormantLife(test_grid)
+        self.assertDictEqual(
+            lifetime_distribution(ALIVE, dl, 3, 0, exclude_alive_after_trans=0),
+            {0: 6, 1: 3})
+
+    def test_exclude_alive_after_trans(self):
+        test_grid = np.array([
+                [DEAD, ALIVE, DEAD],
+                [DEAD, ALIVE, ALIVE],
+                [DEAD, DEAD, DEAD]
+        ])
+        dl = DormantLife(test_grid)
+        self.assertDictEqual(
+            lifetime_distribution(ALIVE, dl, 3, 0, exclude_alive_after_trans=1),
+            {0: 6})
 
 
 if __name__ == '__main__':

time_series.py

@@ -1,11 +1,12 @@
 import sys, os
 import numpy as np
 import multiprocessing
-from gol import CellularAutomaton, DormantLife
+from gol import DormantLife
 
 
-GRID_SIZE = 100
+GRID_SIZE = 20
 BASE_PATH = f"data/dormant-life/time-series/grid-size-{GRID_SIZE}"
+os.makedirs(BASE_PATH, exist_ok=True)
 
 
 def alive_dorm_time_series(dl: DormantLife,
@@ -49,8 +50,10 @@ def dormant_life_time_series(grid_size: int,
 
 def save_data(alive_data, dorm_data, alpha, header: str):
     fname = f"alpha-{str(alpha)[0]}p{str(alpha)[2:]}.dat"
+    os.makedirs(os.path.join(BASE_PATH, "alive"), exist_ok=True)
     np.savetxt(os.path.join(BASE_PATH, "alive", fname),
                (alive_data), header=header)
+    os.makedirs(os.path.join(BASE_PATH, "dorm"), exist_ok=True)
     np.savetxt(os.path.join(BASE_PATH, "dorm", fname),
                (dorm_data), header=header)
 
@@ -71,7 +74,7 @@ def compute(alpha):
     np.savetxt(os.path.join(BASE_PATH, "alpha-range.dat"),
                (alphas), header="Alpha values for which data is stored.")
 
-    with multiprocessing.Pool(processes=3) as pool:
+    with multiprocessing.Pool(processes=4) as pool:
         # Use imap_unordered to apply the function to each value of alpha in 
         # parallel and yield the results as they are ready, regardless of the
         # order