# Adjacency Matrices in Advent Of Code

## The iterate-step pattern
One common pattern for solvers on advent of code to use is what I call `iterate step`. In this pattern you write some function $f : S \rightarrow S$, and then apply $f$ to your initial state $s_0$ repeatedly. For instance, in [2021 day 14](https://adventofcode.com/2021/day/14) you might write

type S = String
type Rules = Map (Char,Char) Char
step :: Rules -> S -> S
step rs = \s -> (head s :) $ concat $ zipWith go s (tail s)
    where
    go :: Char -> Char -> String
    go a b = case Map.lookup (a,b) rs of 
        Nothing -> [b]
        Just c -> [c,b]
        
part1 s rs = fin $ count $ iterate (step rs) s !! 10

Where `fin` above is just the function which goes from a count per letter to whatever specific thing they ask for (in this case the max value minus the min value).

## Improving iterate-step for Part 2
There are a few typical failure modes for the `iterate step` pattern. By far the most common failure mode is that your state type $S$ grows too large. For instance, consider this naive solution to [2021 day 6](https://adventofcode.com/2021/day/6):

def step(fish):
    return [
        f 
        for old in fish 
        for f in ([6,8] if old==0 else [old-1])
        ]

def part1(fish):
    for _ in range(80): 
        fish = step(fish) # This is the iterate-step pattern in python
    return len(fish)

Here we can go to part 2 by changing from storing a list of all fish to a list of the counts of all fish:

def step(fish_counts):
    return {
        age: fish_counts.get(age+1,0) + (fish_counts[0] if age in [6,8] else 0)
        for age in range(9)
        }

def part2(fish):
    fish_counts = {
        age:len([f for f in fish if f==age])
        for age in range(9)
    }
    for _ in range(256): 
        fish_counts = step(fish_counts) # This is the iterate-step pattern in python
    return sum(fish_counts.values())

## But what if there was a Part 3?
This is all well in good if you want to do 256 steps. But what happens if you want to do $10^{12}$ steps? Obviously, you cannot actually run the step function that often.

This is where adjacency matrices can help. An _adjacency matrix_ is a way to encode state changes that is particularly useful if your problem is _linearly reducible_. For instance, in day 6 2021 (the fish problem), it is clear that if you knew `part2([0])`, `part2([1])`, up to `part2([8])`, you could compute `part2` of arbitrary inputs.

To construct an adjacency matrix you must first encode your state into a state _vector_. I will call the function from a state to a vector $V$. For some cases this is easy. For instance, in the fish example, an obvious vector is $V(s)=(a_i)$ where $a_i=\text{# of fish aged }i$. For the example input `3,4,3,1,2`, this would create

$$
V(s_0)=
\begin{pmatrix}
0&1&1&2&1&0&0&0&0
\end{pmatrix}
^T
$$

After one day you end up `2,3,2,0,1`

$$
V(s_1)=
\begin{pmatrix}
1&1&2&1&0&0&0&0&0
\end{pmatrix}
^T
$$

And after two days you end up with `1,2,1,6,0,8`

$$
V(s_2)=
\begin{pmatrix}
1&2&1&0&0&0&1&0&1
\end{pmatrix}
^T
$$

I claim that there exists a matrix $A$ such that $V(s_{i+1})=AV(s_i)$, and name that matrix the adjacency matrix.

### Why Linear Algebra?
Before constructing this matrix, note why this would immediately solve our problem for even extremely large numbers of steps. The proof is somewhat simple: recursively applying the above formula, we see that $V(s_4)=A(A(A(AV(s))))$. But matrix multiplication is associative, so this is equivalent to $V(s_4)=(A*A*A*A)V(s)=A^4V(s)$. We can compute $A^4$ by repeated squaring in 3 moves; while each move may take longer than the `step` function above, you will in general compute $V(s_i)$ in $O(\log i)$ moves. What's more interesting, since the entire work is in computing $A^i$, if you had two different initial states $s_0$ and $s_0'$, computing both answers takes as much time as computing either answer.

### Computing the adjacency matrix
Computing the adjacency matrix is simple. To encode our earlier intuition that if we knew `part2([0])`, `part2([1])`, up to `part2([8])` we would know the entire behavior, we look at the result of `step([0])` etc. We have the following (note the transposition):

$$
V(step([0])) =
\begin{pmatrix}
0&0&0&0&0&0&1&0&1
\end{pmatrix}
^T \\
V(step([1])) =
\begin{pmatrix}
1&0&0&0&0&0&0&0&0
\end{pmatrix}
^T \\
V(step([2])) =
\begin{pmatrix}
0&1&0&0&0&0&0&0&0
\end{pmatrix}
^T \\
\dots \\
V(step([8])) =
\begin{pmatrix}
0&0&0&0&0&0&0&1&0
\end{pmatrix}
^T
$$

Then, the matrix is simply

$$
A =
\begin{pmatrix}
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{pmatrix}
$$

## Final Code
Here is code that would compute `part2` for the 2021 day 6:

adjacency = np.array([
    [(1 if j == i+1 or (j==0 and i in [6,8]) else 0) for j in range(9)]
    for i in range(9)
], dtype = np.longlong)
p2m = np.linalg.matrix_power(adjacency, 256)

def part2(fish):
    fish_counts = np.array([
        len([f for f in fish if f==age])
        for age in range(9)
    ], dtype = np.longlong)
    return sum(np.matmul(p2m, fish_counts))

part2([3,4,3,1,2]) # 26984457539

And here is code that will compute the hypothetical part 3, which wanted you to find the value at $10^{12}$. Note that we overflow the `np.ulonglong` data type: we leave it as an exercise to the reader to implement finding ie the precise last 15 digits of this value.

N = 10 ** 12

adjacency = np.array([
    [(1 if j == i+1 or (j==0 and i in [6,8]) else 0) for j in range(9)]
    for i in range(9)
], dtype = np.ulonglong)
p3m = np.linalg.matrix_power(adjacency, N)

def part3(fish):
    fish_counts = np.array([
        len([f for f in fish if f==age])
        for age in range(9)
    ], dtype = np.ulonglong)
    return sum(np.matmul(p3m, fish_counts)) % np.iinfo(np.ulonglong).max

part3([3,4,3,1,2]) # 1.178868630105016e+19 % np.iinfo(np.ulonglong).max