Cellular Automata

This component implements the cellular automaton as described by [3].

We depend on the box properties being defined. Each Facies should have a viability_range and activation_range defined. Two other parameters are the ca_interval setting how many time steps between every CA advancement, and the ca_random_seed setting the random seed for generating the initial noise.

@kwdef struct Facies <: AbstractFacies
    viability_range::Tuple{Int,Int} = (4, 10)
    activation_range::Tuple{Int,Int} = (6, 10)
end

@kwdef struct Input <: AbstractInput
    ca_interval::Int      = 1
    ca_random_seed::Int   = 0
end

The state of the CA is stored in ca and ca_priority.

@kwdef mutable struct State <: AbstractState
    ca::Matrix{Int}
    ca_priority::Vector{Int}
end

The rules function computes the next value of a cell, given the configured vector of facies, the current facies priority order, and a neighbourhood around the cell.

function rules(facies, ca_priority, neighbourhood)
    cell_facies = neighbourhood[3, 3]
    neighbour_count(f) = sum(neighbourhood .== f)
    if cell_facies == 0
        for f in ca_priority
            n = neighbour_count(f)
            (a, b) = facies[f].activation_range
            if a <= n && n <= b
                return f
            end
        end
        0
    else
        n = neighbour_count(cell_facies) - 1
        (a, b) = facies[cell_facies].viability_range
        (a <= n && n <= b ? cell_facies : 0)
    end
end

The paper talks about cycling the order of preference for occupying an empty cell at each iteration. This means that the rules change slightly every iteration. We need this extra function so that we know the boundary type BT at compile time.

"""
    step_ca(box, facies)

Creates a propagator for the state, updating the celullar automaton in place.

Contract: the `state` should have `ca::Matrix{Int}` and `ca_priority::Vector{Int}`
members.
"""
function step_ca(box::Box{BT}, facies) where {BT<:Boundary{2}}
    tmp = Matrix{Int}(undef, box.grid_size)
    facies_ = facies

    function (state)
        stencil!(BT, Size(5, 5), tmp, state.ca) do nb
            rules(facies_, state.ca_priority, nb)
        end
        state.ca, tmp = tmp, state.ca
        state.ca_priority = circshift(state.ca_priority, 1)
        return state
    end
end

Plot

#| creates: ["docs/src/_fig/ca-long-term.svg"]
#| collect: figures
module Script
using CarboKitten
using CarboKitten.Components: CellularAutomaton as CA
using CairoMakie

function main()
  input = CA.Input(
      box = CarboKitten.Box{Periodic{2}}(
        grid_size=(50, 50), phys_scale=1.0u"m"),
      facies = fill(CA.Facies(), 3)
  )

  state = CA.initial_state(input)
  step! = CA.step!(input)

  for _ in 1:1000
    step!(state)
  end

  fig = Figure(size=(1000, 500))
  axes_indices = Iterators.flatten(eachrow(CartesianIndices((2, 4))))
  xaxis, yaxis = box_axes(input.box)
  i = 1000
  for row in 1:2
    for col in 1:4
      ax = Axis(fig[row, col], aspect=AxisAspect(1), title="step $(i)")

      if row == 2
        ax.xlabel = "x [m]"
      end
      if col == 1
        ax.ylabel = "y [m]"
      end

      heatmap!(ax, xaxis/u"m", yaxis/u"m", state.ca)
      step!(state)
      i += 1
    end
    for _ in 1:996
      step!(state)
      i += 1
    end
  end
  save("docs/src/_fig/ca-long-term.svg", fig)
end

end

Script.main()

Tests

module CellularAutomatonSpec

using Test
using CarboKitten.Components.Common
using CarboKitten.Components: CellularAutomaton as CA

@testset "Components/CellularAutomaton" begin
    let facies = fill(CA.Facies((4, 10), (6, 10)), 3),
        input1 = CA.Input(
            box=Box{Periodic{2}}(grid_size=(50, 50), phys_scale=1.0u"m"),
            facies=facies),
        input2 = CA.Input(
            box=Box{Periodic{2}}(grid_size=(50, 50), phys_scale=1.0u"m"),
            facies=facies,
            ca_random_seed=1)

        state1 = CA.initial_state(input1)
        state2 = CA.initial_state(input2)
        state3 = CA.initial_state(input2)

        @test CA.initial_state(input1).ca == CA.initial_state(input1).ca
        @test state1.ca != state2.ca

        step! = CA.step!(input1)  # inputs have same rules
        for i in 1:20
            step!(state1)
            step!(state2)
            step!(state3)
        end

        @test state1.ca != state2.ca
        @test state2.ca == state3.ca
    end
end

end

Component

@compose module CellularAutomaton
    @mixin Boxes, FaciesBase
    using ..Common
    using ...Stencil: stencil!
    using Random

    <<ca-input>>
    <<ca-state>>
    <<ca-step>>

    function initial_state(input::AbstractInput)
        n_facies = length(input.facies)
        ca = rand(MersenneTwister(input.ca_random_seed), 0:n_facies, input.box.grid_size...)
        return State(ca, 1:n_facies |> collect)
    end

    function step!(input::AbstractInput)
        return step_ca(input.box, input.facies)
    end
end