CODEWITHSUNDEEP
arraysjuliacellular-automata
universal_amateur
universal_amateur

Reputation: 149

Julia / Cellular Automata: efficient way to get neighborhood

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...

Upvotes: 4

Views: 816

Answers (2)

Bogumił Kamiński
Bogumił Kamiński

Reputation: 69949

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

Upvotes: 4

Bill
Bill

Reputation: 6086

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.

Upvotes: 2

Related Questions