Julia / Cellular Automata: efficient way to get neighborhood

Question

I'd like to implement a cellular automaton (CA) in Julia. Dimensions should be wrapped, this means: the left neighbor of the leftmost cell is the rightmost cell etc.

One crucial question is: how to get the neighbors of one cell to compute it's state in the next generation? As dimensions should be wrapped and Julia does not allow negative indices (as in Python) i had this idea:

Considered a 1D CA, one generation is a one-dimensional array:

0 0 1 0 0

What if we create a two dimensional Array, where the first row is shifted right and the third is shifted left, like this:

0 0 0 1 0
0 0 1 0 0
0 1 0 0 0

Now, the first column contain the states of the first cell and it's neighbors etc.

i think this can easily be generalized for two and more dimensions.

First question: do you think this is a good idea, or is this a wrong track?

EDIT: Answer to first question was no, second Question and code example discarded.

Second question: If the approach is basically ok, please have a look at the following sketch:

EDIT: Other approach, here is a stripped down version of a 1D CA, using mod1() for getting neighborhood-indices, as Bogumił Kamiński suggested.

for any cell: - A array of all indices - B array of all neighborhood states - C states converted to one integer - D lookup next state

function digits2int(digits, base=10)
   int = 0
   for digit in digits
      int = int * base + digit
   end
   return int
end

gen = [0,0,0,0,0,1,0,0,0,0,0]
rule = [0,1,1,1,1,0,0,0]

function nextgen(gen, rule)
  values = [mod1.(x .+ [-1,0,1], size(gen)) for x in 1:length(gen)] # A
  values = [gen[value] for value in values]                         # B
  values = [digits2int(value, 2) for value in values]               # C
  values = [rule[value+1] for value in values]                      # D
  return values
end

for _ in 1:100
  global gen
  println(gen)
  gen = nextgen(gen, rule)
end

Next step should be to extend it to two dimensions, will try it now...

score 4 · Accepted Answer · answered Sep 26 '19 at 11:41

4

The way I typically do it is to use mod1 function for wrapped indexing.

In this approach, no matter what dimensionality of your array a is then when you want to move from position x by delta dx it is enough to write mod1(x+dx, size(a, 1)) if x is the first dimension of an array.

Here is a simple example of a random walk on a 2D torus counting the number of times a given cell was visited (here I additionally use broadcasting to handle all dimensions in one expression):

function randomwalk()
    a = zeros(Int, 8, 8)
    pos = (1,1)
    for _ in 1:10^6
        # Von Neumann neighborhood
        dpos = rand(((1,0), (-1,0), (0,1), (0,-1)))
        pos = mod1.(pos .+ dpos, size(a))
        a[pos...] += 1
    end
    a
end

answered Sep 26 '19 at 11:41

Bogumił Kamiński

66,844
3
80
107

Thank you! I just implemented a one dimensional CA using mod1 for getting the indices, and it is roundabout three times faster than my approach... Still, i have some difficulties understanding how to get from one to two dimensions. Will clarify my question and check back. – universal_amateur Sep 26 '19 at 14:28
In my example I have a two dimensional matrix. – Bogumił Kamiński Sep 26 '19 at 14:31
yeah i see, but i tried to build an example in one dimension where i get *all* the neighbor-indices from *all* cells at once like this: [mod1.(x .+ [-1,0,1], size(gen)) for x in 1:length(gen)] ... this works well, but cannot transfer it to two dimensions... Going to prepare the question. – universal_amateur Sep 26 '19 at 14:42
It would be easier if you asked about more atomic issues in separate questions, as it is hard to answer the question if it is updated several times. – Bogumił Kamiński Sep 26 '19 at 20:39
No problem, I am just hinting what is simpler to answer :). – Bogumił Kamiński Sep 27 '19 at 05:45
Basically, you answered my question, mod1() works perfect for getting neighborhood indices in any "wrapped" discret space. Will give a more detailed example sooner or later :) – universal_amateur Sep 28 '19 at 20:11

score 2 · Answer 2 · answered Sep 26 '19 at 11:31

Usually, if the CA has cells that are only dependent on the cells next to them, it's simpler just to "wrap" the vector by adding the last element to the front and the first element to the back, doing the simulation, and then "unwrap" by taking the first and last elements away again to get the result length the same as the starting array length. For the 1-D case:

const lines = 10
const start = ".........#........."
const rules = [90, 30, 14]

rule2poss(rule) = [rule & (1 << (i - 1)) != 0 for i in 1:8]

cells2bools(cells) = [cells[i] == '#' for i in 1:length(cells)]

bools2cells(bset) = prod([bset[i] ? "#" : "." for i in 1:length(bset)])

function transform(bset, ruleposs)
    newbset = map(x->ruleposs[x],
        [bset[i + 1] * 4 + bset[i] * 2 + bset[i - 1] + 1
        for i in 2:length(bset)-1])
    vcat(newbset[end], newbset, newbset[1])
end

const startset = cells2bools(start)

for rul in rules
    println("\nUsing Rule $rul:")
    bset = vcat(startset[end], startset, startset[1]) # wrap ends
    rp = rule2poss(rul)
    for _ in 1:lines
        println(bools2cells(bset[2:end-1]))  # unwrap ends
        bset = transform(bset, rp)
    end
end

As long as only the adjacent cells are used in the simulation for any given cell, this is correct.

If you extend this to a 2D matrix, you would also "wrap" the first and last rows as well as the first and last columns, and so forth.

From my experience this approach is typically used e.g. when you distribute the computations out-of-core (so you have to cut your array into many sub arrays on separate processes). However, it has a cost that you have to sync the value in edges after every update of the edge cell. — Bogumił Kamiński, Sep 26 '19 at 11:46

Julia / Cellular Automata: efficient way to get neighborhood

2 Answers2