2
$\begingroup$

I am trying to simulate the spread of mold in a petri dish using a cellular automata based approach. Thanks to the answer in my other question Stochastic cellular automata - algorithm limited by 1 cell per timestep, I am attempting a different approach using a priority queue.

Context

The petri dish can be thought of as a grid of 1mm x 1mm squares, each called a Cell. Each square has a value for the rate of spread (RoS) of mold in that type of food, which is available as empirically measured data. For this question, let's just assume the petri dish is all the same food type with an RoS of 1 mm/hr.

We seed the dish by saying a particular Cell is moldy. Then we kick off an iteration. For every cell, we calculate at what time its 8 neighbors will become moldy. We add those cells to a priority queue, with the priority being the time at which they should become moldy. For example, examining the first step:

enter image description here

At the first step, there is only 1 moldy cell. We can expect that neighbors 2, 4, 5, and 7 will become moldy after 1 hour (it takes 1 hour for mold to spread 1 mm for a RoS of 1mm/hr), and we can expect that neighbors 1, 3, 6, and 8 would become moldy after $\sqrt(2)$ hours, as the distance between diagonal cells is not 1mm, but rather $\sqrt(2)$ mm. We can say that the time at which a neighbor will become moldy can be returned by a function $g(food, distance)$. So, when iterating over this first cell's neighbors, we should expect 2 items in the queue: neighbors 2, 4, 5, and 7 will become moldy after 1 hour, and neighbors 1, 3, 6, and 8 will become moldy after $\sqrt(2)$ hours.

We have thus completed the first iteration over all moldy cells, and we move on to the next item in the queue, which is to iterate over all cells that will become moldy at 1 hour, and calculate their neighbors times to become moldy. A quick pseudocode of the algorithm might look like this

mold_next = PriortyQueue({Initially moldy cells each flagged to t=0})
while mold_next is not empty and t<max_timestep:
  c = pop(mold_next)
  if c is moldy:  # Very important
    continue at next iteration of while loop
  make c moldy
  for all neighbours n of c:
    add n to mold_next at a time drawn from g(food, distance)

The problem

My iteration of this (written in typescript and too complicated to post, or even extract to a minimum reproducible example), has a problem. The shape that emerges is not a circle, as one might expect, but rather, always a square. Let's look at the first few timesteps:

enter image description here

Note the last few steps are not cut off - there are no moldy cells outside the bounds of the image. Take a look at a gif showing more timesteps:

enter image description here

The pattern that emerges is that as the shape grows, it "grows in squares". Meaning from the first cell, moldy cells emerge at the four cadinal direction of the cell. Then, the corners are filled in. Then again, moldy cells emergy from the four cardinal directions at the center of each side, and begin "filling in" towards the edge, until again we have a square. It is not until a square "completes" that the next cell emerges in the cardinal direction.

Searching for clues

I am trying to understand why this is. Looking at the first iteration:

enter image description here

Let's say we are beginning the iteration in this state: The first cell generated its 2 queue items, as described earlier in this question. We are now iterating over the subsequent queue item, which says that cells 2, 4, 5, and 7, will become moldy. We've set them to moldy. Now we iterate over them. The cell in question is the one the arrows are pointing to, say cell $C_2$. The 2 queue items that the first cell generated says that this cell will become moldy at $t_0 + \sqrt(2)$. The cells from which the black arrow originates mold at $t_0 + 1$, and they each add an item to the queue for $C_2$, which would be at $t_0 + 2$, but by $t_0 + 2$, this cell would have already become moldy (at $t_0 + \sqrt(2)$), and as the algorithm says, if c is moldy, continue at next iteration of while loop.

But somehow this pattern has created a squarish emergence, where the corners are always filling in first. Looking a few timesteps down:

enter image description here

Iterating over these cells, I would expect $C_n$ to mold before all $C_1$, $C_2$, and $C_3$, but that is not the case. The order of molding is $C_1$, $C_2$, $C_3$, then $C_n$. I am struggling to understand why this happens.

Fundamentally flawed logic vs implentation bug?

Is there a fundamental problem with my algorithm design? Or do I perhaps have an implementation bug? Should the algorithm I've described not cause a shape that is roughly circular to emerge?

Improve this question
$\endgroup$
1
  • $\begingroup$ Is your priority queue sorting the wrong way? $\endgroup$
    – Richard
    Commented Sep 11, 2021 at 0:58

1 Answer 1

Reset to default
2
$\begingroup$

I think either your priority queue is sorted the wrong direction (largest first) or you're using some kind of non-priority queue.

The Python code below implements the algorithm you describe and exhibits approximately circular growth. I've also included a switch to turn on a more random growth.

#!/usr/bin/env python3

import time
from heapq import heappop, heappush
from typing import Iterator, List, Tuple

import numpy as np
from numpy.random import exponential


def neighbours(x: int, y: int, grid: np.ndarray) -> Iterator[Tuple[int, int, float]]:
  """Generate a list of neighbour coordinates for `x,y`"""
  dxs = [-1, -1,  0,  1, 1, 1, 0, -1]
  dys = [ 0, -1, -1, -1, 0, 1, 1,  1]
  for dx, dy in zip(dxs, dys):
    nx, ny = x + dx, y + dy
    dist = 1 if (dx==0 or dy==0) else np.sqrt(2)
    if nx<0 or ny<0 or ny==grid.shape[0] or nx==grid.shape[1]:
      continue
    yield nx, ny, dist


grid = np.zeros((20, 20))
pq: List[Tuple[float, int, int]] = [] # Time, x, y
heappush(pq, (0.0, 9, 9))
print(pq)

while pq:
  t, x, y = heappop(pq)
  if grid[y,x]==1: # Already moldy
    continue
  grid[y,x] = 1 # Make it moldy
  for nx, ny, dist in neighbours(x, y, grid):
    heappush(pq, (t + dist, nx, ny))  # Circular growth
    # Enable these lines and disable above for random growth
    # Draw from exponential distribution to convert rate to interval
    # rint = exponential(dist)
    # heappush(pq, (t + rint, nx, ny))
  time.sleep(0.05)
  print("\033[2J") # Clear the screen
  print(grid)
```
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you for this answer! - I see that the behavior is as expected. Give me some time to digest what is going on here and I'll get back to this when I can $\endgroup$
    – Seth Lutske
    Commented Sep 13, 2021 at 4:17
  • 1
    $\begingroup$ So, it turns out the issue was a tiny implementation detail. Rather than setting each new $pq_n$'s time ($t_{pq_{n}}$) to the $t_{pq_{n-1}} + dist$, I was setting it to $t_{pq_{n-1}} + k * n$, where $k$ was a constant left over from an old implementation. But your ability to boil the entire algorithm down to this very concise script was extremely helpful is comparing and debugging what was going wrong. Thanks again! $\endgroup$
    – Seth Lutske
    Commented Sep 13, 2021 at 22:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.