We are simulating a grid of octopi, i.e. of cellular automata. Each is represented by a natural number, and a step of the simulation goes like this:

The energy level of each cell (the number representing it) increases by 1.

If any cell’s energy level is greater than 9, then the energy level of each adjacent (including diagonally adjacent) cell increases by 1. Chain reactions may trigger at this point, but each cell can “flash” at most once each step.

The energy level of each cell that “flashed” becomes 0.

Once again, we have a grid-like type at the center of our puzzle. We require the following operations:

```
neighbors : Grid ℕ → (ℕ × ℕ) → [((ℕ × ℕ) × ℕ)]
inc_all : Grid ℕ → Grid ℕ
```

As usual, we also assume Grid is a functor.

*neighbors (i, j)* gives a list of the neighbors of the cell at
(i, j) paired with their energy levels. *inc_all* increments every
cell in the grid; thus, it describes the first part of a step of the
simulation. [Note: I didn’t actually check what information we were
supposed to be collecting until after writing all of this. We are
looking for the total number of flashes that occur, so some changes
to the following are needed to return this information.]

```
step : Grid ℕ → Grid ℕ
step = inc_all ; ... ?
```

Handling flashing cells and their effects (indeed) is where things get trickier. A flash is a nonlinear phenomenon, so we cannot model it compositionally. Each flashing cell can, in principle, transform the entire grid.

```
resolve_flashes : Grid ℕ → Grid ℕ
reset_flashed : Grid ℕ → Grid ℕ
step = inc_all ; resolve_flashes ; reset_flashed
```

Further complicating matters (though saving us from a truly off-the-rails simulation) is the fact that each cell can flash at most once per step. We therefore need to reference the set of cells that have flashed during the current step.

An important fact for working out *resolve_flashes* is that the set
of cells that will flash during a given step shrinks by only one
cell when that cell’s flash is resolved. A single flash may push
many other cells into the set of flashers, but it can’t pull any
out.

In other words, the set of flashing cells behaves like a stack.

To simplify things, we will consider the propagation of a flash as a computation using the parametric State type. Each computation produces a list of “new” flashing cells.

```
resolve_flashes : Grid ℕ → Grid ℕ
resolve_flashes g = π₁ (execState (resolve (flashers g)) (g, ∅))
flashers : Grid ℕ → [(ℕ × ℕ)]
flashers = keys ∘ filter (> 9)
resolve : [(ℕ × ℕ)] → State ((Grid ℕ) × Set (ℕ × ℕ)) ⊤
```

⊤ here is the unit type, which has a single member, `tt`

.
(Haskell `()`

of type `()`

, not the greatest notation ever, IMO.)

Let’s try to figure out *resolve*. If the list *cs* of flashing
cells is empty, then *resolve* is simply the identity transformer:

```
resolve [] = pure tt
```

Otherwise, at least one flash occurs, so the state must be
transformed. Using monadic `do`

notation,

```
resolve (c ∷ cs) = do new ← resolve1 c
resolve (push_adj new cs)
push_adj : Eq α ⇒ [α] → [α] → [α]
push_adj xs ys = (filter (λ x → not (x ∈ ys)) xs) ⧺ ys
```

*resolve1* (i, j) flashes the cell at (i, j), updates the state,
and returns a list of the neighbors of (i, j) which have become
flashers. We require that the address passed to *resolve1* is
never that of a cell that has already flashed during this step.

```
resolve1 : (ℕ × ℕ) → State ((Grid ℕ) × Set (ℕ × ℕ)) [(ℕ × ℕ)]
resolve1 c =
state (λ (g, s) →
let g' = flash g c
ns = neighbors g' c
cs = list π₁ (filter ((> 9) ∘ π₂) ns)
s′ = adjoin c s
in (filter (not (∈ s′)) cs, (g', s′)))
```

*state* is a library function which embeds a function
σ → (α × σ) (transforming a state σ and returning a value α) into
the state monad; it gives us a pretty compact way to write *resolve1*.
The extensive transformation of the list of neighbors gives us the
list of new flashers; first, we filter for cells with energy above 9,
then we remove those in the “already-flashed” set.

*flash c* increments the neighbors of the cell at *c*:

```
flash : Grid ℕ → (ℕ × ℕ) → Grid ℕ
flash g c =
map_with_key (λ d n →
if in_square c d ∧ c ≠ d
then n + 1
else n)
g
```

*map_with_key* : ((ℕ × ℕ) → α → β) → Grid α → Grid β is a variant of
functorial application in the Grid type that passes the address as
well as the element to the applied function.

The function *in_square* : (ℕ × ℕ) → (ℕ × ℕ) determines whether the
second address is within the 3 × 3 square surrounding and including
the first.

We now have a definition for *step*, with which we can write a
function that runs the simulation for a given number of steps:

```
steps : ℕ → Grid → Grid
steps n = iterateN step n
```

where *iterateN f n x* = f (f … (f x)) (*n* repeated applications
of *f*).

*This is a clean-enough process, but I’d like to be a bit more formal.
We didn’t really come up with any equations describing this
simulation; instead, the above was more a matter of picking a way to
model it that “seemed good”, then elaborating it.*

This is much the same as Part 1, except we will now count the number of steps it takes to reach a certain state; namely, a state where all the cells flash at once. Since a flash resets a cell, to zero, we can express this by the predicate

```
all_flashed : Grid ℕ → Bool
all_flashed = all (== 0)
```

where *all p* is true iff every cell of a grid satisfies *p*.

We can express the step-counting with:

```
steps_until : (Grid → Bool) → Grid → ℕ
steps_until p g = length (takeWhile (¬ p) (iterate step g))
```

where *iterate f x* produces the infinite list of repeated
applications of *f*:
[*x*, *f x*, *f* (*f x*), …].

Then:

```
solution : Grid → ℕ
solution = steps_until all_flashed
```