Merge pull request #31 from MartinuzziFrancesco/fm/st

feat: add stochastic updates

Author: MartinuzziFrancesco

Project.toml

@@ -1,10 +1,11 @@
 name = "CellularAutomata"
 uuid = "878138dc-5b27-11ea-1a71-cb95d38d6b29"
 authors = ["Francesco Martinuzzi"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 Aqua = "0.8"
@@ -20,9 +21,8 @@ julia = "1.10"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "SafeTestsets", "Random", "Aqua", "ForwardDiff", "JET"]
+test = ["Test", "SafeTestsets", "Aqua", "ForwardDiff", "JET"]

README.md

@@ -41,6 +41,23 @@
 CellularAutomata.jl offers lightweight and ready to use implementations of
 cellular automata, one and two dimensional.
 
+## Features
+
+- **Discrete CA** (`DCA`): elementary, general (any radius/state count,
+  symmetric or asymmetric neighborhoods), and totalistic rules, following
+  Wolfram's *A New Kind of Science* numbering.
+- **Continuous CA** (`CCA`): real-valued totalistic rules.
+- **Life-like CA** (`Life`): Conway's Game of Life and any birth/survival
+  variant, in the Golly notation.
+- **Boundary conditions**: `Periodic`, `Reflecting`, `ConstantBoundary`.
+- **Update schemes**: `Synchronous` (default) and `Stochastic(rate)` for
+  per-cell probabilistic updates, given an explicit `rng`.
+- **Functional core** (`next_state`, `rollout`): non-mutating, AD/GPU-friendly
+  building blocks for custom or neural cellular automata, alongside the
+  stateful `CellularAutomaton` wrapper.
+- **Analysis**: `evolution_history`, `lempel_ziv` complexity, and other
+  accessor functions.
+
 ## Installation
 CellularAutomata.jl is registered on the general registry. For the installation use either of:
 

docs/src/ad_design.md

@@ -8,21 +8,29 @@ without changing the existing discrete-rule API.
 ## Functional core API
 
 ```julia
-next_state(rule, state; boundary=Periodic()) -> state
-rollout(rule, state, steps; boundary=Periodic(), save=false) -> result
+next_state(rule, state; boundary=Periodic(), scheme=Synchronous(), rng=nothing) -> state
+rollout(rule, state, steps; boundary=Periodic(), scheme=Synchronous(),
+        rng=nothing, save=false) -> result
 ```
 
-`next_state` is a pure function: it does not mutate `state` or captured parameter
-arrays. `rollout` is the orchestration layer. It returns the final state by default,
-which is the preferred path inside a loss function. With `save=true`, it returns the
-initial state and all subsequent states with time on the last axis.
+`next_state` does not mutate `state` or captured parameter arrays. `rollout` is the
+orchestration layer. It returns the final state by default, which is the preferred
+path inside a loss function. With `save=true`, it returns the initial state and all
+subsequent states with time on the last axis.
 `CellularAutomaton` remains a compatibility wrapper and preserves its original
 history layout.
 
 Boundary behavior is represented by small types — `Periodic`, `Reflecting`, and
 `ConstantBoundary` — rather than being embedded in every rule implementation.
 Neighborhood extraction is likewise shared across rules.
 
+Update schemes follow the same pattern: `Synchronous` (default) applies every
+cell's transition, while `Stochastic(rate)` keeps a cell's previous value unless
+a randomly drawn per-cell mask selects it. Stochastic updates require an explicit
+`rng` argument threaded through `next_state`/`rollout`; drawing a mask advances
+that RNG. The sampled mask is treated as data by the transition, so gradients do
+not differentiate through the random draw or the discrete selection decision.
+
 ## State and rule representation
 
 Discrete automata accept vectors and matrices. Differentiable rules need a documented
@@ -45,9 +53,9 @@ end
 
 spatial_dimensions(::NeuralRule) = 2
 
-function CellularAutomata.next_state(rule::NeuralRule, x; boundary=Periodic())
+function CellularAutomata.__step(rule::NeuralRule, x, boundary)
     features = rule.perceive(x, boundary)
-    return x + rule.update(features)
+    return x + rule.update(features) # proposed next state
 end
 ```
 
@@ -59,10 +67,10 @@ that narrow boundary.
 
 ## Migration sequence
 
- 1. **Implemented:** pure `next_state`, functional `rollout`, explicit boundary types,
-    shared neighborhood access, compatibility through `CellularAutomaton`, and
-    multi-step ForwardDiff coverage. Lookup-based DCA/TCA and thresholded Life rules
-    remain intentionally nondifferentiable.
+ 1. **Implemented:** non-mutating `next_state`, functional `rollout`, explicit
+    boundary types, shared neighborhood access, compatibility through
+    `CellularAutomaton`, and multi-step ForwardDiff coverage. Lookup-based DCA/TCA
+    and thresholded Life rules remain intentionally nondifferentiable.
  2. Test reverse-mode AD and GPU arrays after selecting the package's supported ML
     stack.
  3. Add a generic continuous multidimensional local rule with a documented

docs/src/api/general.md

@@ -4,6 +4,9 @@
     CellularAutomaton
     next_state
     rollout
+    AbstractUpdateScheme
+    Synchronous
+    Stochastic
     cellular_automaton_rule
     evolution_history
     generation_count

src/CellularAutomata.jl

@@ -1,6 +1,7 @@
 module CellularAutomata
 
 using ConcreteStructs: @concrete
+using Random: AbstractRNG
 
 include("generics.jl")
 include("cellular_automaton.jl")
@@ -12,6 +13,7 @@ include("measures.jl")
 
 export AbstractCellularAutomaton, AbstractCellularAutomatonRule
 export AbstractBoundaryCondition, Periodic, Reflecting, ConstantBoundary
+export AbstractUpdateScheme, Synchronous, Stochastic
 export AbstractDiscreteCellularAutomatonRule
 export AbstractContinuousCellularAutomatonRule
 export AbstractTotalisticCellularAutomatonRule

src/cellular_automaton.jl

@@ -6,10 +6,12 @@ Supertype for cellular automata that retain an evolution history.
 abstract type AbstractCellularAutomaton end
 
 """
-    next_state(rule, state; boundary=Periodic())
+    next_state(rule, state; boundary=Periodic(), scheme=Synchronous(), rng=nothing)
 
-Apply one synchronous cellular-automaton transition without mutating `state`.
-Rules may specialize this function to define differentiable transitions.
+Apply one cellular-automaton transition without mutating `state`. Rule subtypes
+define their transition by implementing
+`CellularAutomata.__step(rule, state, boundary)`; this wrapper then applies the
+selected update scheme.
 
 # Arguments
 
@@ -20,6 +22,11 @@ Rules may specialize this function to define differentiable transitions.
 
   - `boundary`: Boundary condition used outside the state array.
     Defaults to `Periodic()`.
+  - `scheme`: Update scheme controlling which cells apply the transition.
+    Defaults to `Synchronous()`; use `Stochastic(rate)` for a per-cell random mask.
+  - `rng`: Random number generator used by stochastic schemes. It must be supplied
+    explicitly with `Stochastic` and its state advances when the mask is drawn.
+    It is not consulted by `Synchronous`.
 
 # Examples
 
@@ -38,23 +45,33 @@ julia> next_state(DCA(30), [0, 0, 1, 0, 0])
 function next_state(
         rule::AbstractCellularAutomatonRule,
         state;
-        boundary::AbstractBoundaryCondition = Periodic()
+        boundary::AbstractBoundaryCondition = Periodic(),
+        scheme::AbstractUpdateScheme = Synchronous(),
+        rng::Union{Nothing, AbstractRNG} = nothing
     )
     __validate_state(rule, state)
-    return __step(rule, state, boundary)
+    __validate_scheme_rng(scheme, rng)
+    new_state = __step(rule, state, boundary)
+    return __apply_scheme(scheme, rng, state, new_state)
 end
 
 @concrete struct __Transition
     rule
     boundary
+    scheme
+    rng
 end
 
 function (transition::__Transition)(state, _)
-    return next_state(transition.rule, state; boundary = transition.boundary)
+    return next_state(
+        transition.rule, state;
+        boundary = transition.boundary, scheme = transition.scheme, rng = transition.rng
+    )
 end
 
 """
-    rollout(rule, initial_state, steps; boundary=Periodic(), save=false)
+    rollout(rule, initial_state, steps; boundary=Periodic(), scheme=Synchronous(),
+        rng=nothing, save=false)
 
 Apply `steps` transitions of `rule` to `initial_state`. By default only the final
 state is returned, which is the preferred path inside a loss function. With
@@ -70,6 +87,10 @@ on the last axis.
 # Keyword arguments
 
   - `boundary`: Boundary condition used by each transition. Defaults to `Periodic()`.
+  - `scheme`: Update scheme applied at every step. Defaults to `Synchronous()`.
+  - `rng`: Random number generator passed to every step. It must be supplied
+    explicitly with `Stochastic`; the same object is reused and its state advances
+    across steps. It is not consulted by `Synchronous`.
   - `save`: Retain the initial state and every subsequent state. Defaults to `false`.
 
 # Examples
@@ -96,11 +117,14 @@ function rollout(
         initial_state,
         steps::Integer;
         boundary::AbstractBoundaryCondition = Periodic(),
+        scheme::AbstractUpdateScheme = Synchronous(),
+        rng::Union{Nothing, AbstractRNG} = nothing,
         save::Bool = false
     )
     steps >= 0 || throw(ArgumentError("steps must be nonnegative"))
     isempty(initial_state) && throw(ArgumentError("initial_state cannot be empty"))
-    transition = __Transition(rule, boundary)
+    __validate_scheme_rng(scheme, rng)
+    transition = __Transition(rule, boundary, scheme, rng)
     if save
         states = accumulate(transition, 1:steps; init = initial_state)
         return stack(Iterators.flatten(((initial_state,), states)))

src/generics.jl

@@ -2,8 +2,9 @@
     AbstractCellularAutomatonRule
 
 Supertype for cellular-automaton transition rules. Subtypes must implement
-`__step(rule, state, boundary)` returning the next state without mutating `state`,
-or specialize [`next_state`](@ref) directly.
+`CellularAutomata.__step(rule, state, boundary)` returning the proposed next state
+without mutating `state`. [`next_state`](@ref) applies validation and the selected
+update scheme around that transition.
 """
 abstract type AbstractCellularAutomatonRule end
 
@@ -79,6 +80,71 @@ function Base.show(io::IO, boundary::ConstantBoundary)
     return print(io, "ConstantBoundary(", boundary.value, ")")
 end
 
+"""
+    AbstractUpdateScheme
+
+Supertype for cellular-automaton update schemes, controlling which cells apply a
+transition each generation.
+"""
+abstract type AbstractUpdateScheme end
+
+"""
+    Synchronous()
+
+Apply the transition to every cell each generation. This is the default update
+scheme and matches classical (non-stochastic) cellular automata.
+"""
+struct Synchronous <: AbstractUpdateScheme end
+
+"""
+    Stochastic(rate)
+
+Apply the transition to each cell independently with probability `rate`; cells
+that are not selected keep their previous value. An explicit `rng` must be passed
+to [`next_state`](@ref) or [`rollout`](@ref). Drawing the per-cell mask advances
+the state of that random number generator.
+
+# Examples
+
+```jldoctest
+julia> using CellularAutomata, Random
+
+julia> next_state(DCA(30), [0, 0, 1, 0, 0]; scheme=Stochastic(0.5), rng=Xoshiro(1))
+5-element Vector{Int64}:
+ 0
+ 1
+ 1
+ 0
+ 0
+```
+"""
+struct Stochastic{T <: Real} <: AbstractUpdateScheme
+    rate::T
+    function Stochastic{T}(rate::T) where {T <: Real}
+        0 <= rate <= 1 || throw(ArgumentError("rate must be between 0 and 1"))
+        return new{T}(rate)
+    end
+end
+
+Stochastic(rate::T) where {T <: Real} = Stochastic{T}(rate)
+
+__validate_scheme_rng(::Synchronous, rng) = nothing
+__validate_scheme_rng(::Stochastic, ::AbstractRNG) = nothing
+function __validate_scheme_rng(::Stochastic, ::Nothing)
+    throw(ArgumentError("an explicit rng is required when using Stochastic"))
+end
+
+__apply_scheme(::Synchronous, rng, state, new_state) = new_state
+
+function __apply_scheme(::Stochastic, ::Nothing, state, new_state)
+    throw(ArgumentError("an explicit rng is required when using Stochastic"))
+end
+
+function __apply_scheme(scheme::Stochastic, rng::AbstractRNG, state, new_state)
+    mask = rand(rng, size(state)...) .< scheme.rate
+    return ifelse.(mask, new_state, state)
+end
+
 @inline __boundary_index(i, n, ::Periodic) = mod1(i, n)
 @inline function __boundary_index(i, n, ::Reflecting)
     n == 1 && return 1

test/regression_test.jl

@@ -34,6 +34,8 @@ end
         CCA(0.1),
         Life(((3,), (2, 3))),
         ConstantBoundary(0.0),
+        Synchronous(),
+        Stochastic(0.5),
     )
     for object in objects
         @test all(isconcretetype, fieldtypes(typeof(object)))

test/runtests.jl

@@ -35,6 +35,10 @@ if GROUP in ("All", "Core")
         @safetestset "Edge cases and hierarchy" include("regression_test.jl")
     end
 
+    @testset "Update schemes" begin
+        @safetestset "Synchronous and Stochastic" include("update_scheme_test.jl")
+    end
+
     @testset "Life-like" begin
         @safetestset "Life glider" include("glider_test.jl")
         @safetestset "Life blinker" include("blinker_test.jl")

test/update_scheme_test.jl

@@ -0,0 +1,47 @@
+using CellularAutomata
+using Random
+using Test
+
+@testset "Stochastic rate validation" begin
+    @test_throws ArgumentError Stochastic(-0.1)
+    @test_throws ArgumentError Stochastic(1.1)
+    @test_throws ArgumentError next_state(DCA(30), [0, 0, 1]; scheme = Stochastic(0.5))
+    @test_throws ArgumentError rollout(DCA(30), [0, 0, 1], 0; scheme = Stochastic(0.5))
+end
+
+@testset "Stochastic edge rates match Synchronous" begin
+    rule = DCA(30)
+    state = [0, 0, 1, 0, 0]
+    @test next_state(rule, state; scheme = Stochastic(0.0), rng = Xoshiro(1)) == state
+    @test next_state(rule, state; scheme = Stochastic(1.0), rng = Xoshiro(1)) ==
+        next_state(rule, state)
+end
+
+@testset "custom rule transitions still apply schemes" begin
+    struct TestRule <: AbstractCellularAutomatonRule end
+    CellularAutomata.__step(::TestRule, state, boundary) = one.(state)
+
+    state = zeros(Int, 4)
+    @test next_state(TestRule(), state) == ones(Int, 4)
+    @test next_state(
+        TestRule(), state; scheme = Stochastic(0.0), rng = Xoshiro(1)
+    ) == state
+end
+
+@testset "Stochastic is reproducible given the same rng" begin
+    rule = DCA(30)
+    state = [0, 0, 1, 0, 0]
+    result = next_state(rule, state; scheme = Stochastic(0.5), rng = Xoshiro(1))
+    @test result == next_state(rule, state; scheme = Stochastic(0.5), rng = Xoshiro(1))
+end
+
+@testset "rollout threads a single rng across steps" begin
+    rule = DCA(30)
+    state = [0, 0, 1, 0, 0]
+
+    rng = Xoshiro(1)
+    first_step = next_state(rule, state; scheme = Stochastic(0.5), rng = rng)
+    second_step = next_state(rule, first_step; scheme = Stochastic(0.5), rng = rng)
+
+    @test rollout(rule, state, 2; scheme = Stochastic(0.5), rng = Xoshiro(1)) == second_step
+end