adventofcode2023/day21/sol.py

import sys
import numpy

def read_map(input):
    map = [x.strip() for x in input]

    if len(map) < 80:
        print(map)

    Y = len(map)
    X = len(map[0])
    assert all(len(row) == X for row in map)
    assert X == Y

    return map

def draw(mask, fill):
    if len(fill) < 50:
        for i, row in enumerate(fill):
            print("".join(".O#!"[f + 2*mask[i,j]] for j,f in enumerate(row)))
    print()


def solve(map):
    Y = len(map)
    X = len(map[0])

    fill = numpy.zeros((Y,X+1), dtype='uint8')
    mask = numpy.zeros((Y,X+1), dtype='uint8')
    for i, row in enumerate(map):
        for j, c in enumerate(row):
            if c == '#':
                mask[i,j] = 1
            if c == 'S':
                start = (i,j)

    # add a border on the side of the map (to prevent wraparound)
    mask[:, -1] = 1

    # start at the start
    fill[start] = 1

    for i in range(64):
        fill = step(mask, fill)
        #draw(mask, fill)
        if i == 6-1:  # sample
            print(fill.sum())

    print("part1 = ", fill.sum())


def solve2(map):
    Y = len(map)
    X = len(map[0])

    fill = numpy.zeros((Y,X), dtype='uint8')
    mask = numpy.zeros((Y,X), dtype='uint8')
    for i, row in enumerate(map):
        for j, c in enumerate(row):
            if c == '#':
                mask[i,j] = 1
            if c == 'S':
                start = (i,j)

    print(fill.shape, mask.shape)

    fill[start] = 1
    values, period, d2 = simulate(fill, mask, max_steps=500)

    for n in 6,10,50,100,500,1000, 5000, 26501365:
        print(n, extrapolate(n, values, period, d2))


def simulate(fill, mask, max_steps):
    Y,X = mask.shape
    assert Y == X, "need a square matrix"
    period = X
    tiles = 1

    original_mask = mask
    prev = [0]
    diffs = []
    for i in range(max_steps):
        # expand the map if necessary
        if fill[0].any() or fill[-1].any() or fill[:,0].any() or fill[:,-1].any():
            draw(mask,fill)
            print("expanding...")
            tiles += 2
            mask = numpy.tile(original_mask, (tiles,tiles))
            fill = numpy.pad(fill, [(Y,Y), (X,X)])

        # take one more step and count the number of reachable squares
        fill = step(mask, fill)
        n = int(fill.sum())

        # look for patterns
        # although the number of squares is somewhat unpreditable from step
        # to step (due to the maze-like structure of the mask), the fact that
        # the mask repeats in tiles (and the fact that there are unobstructed
        # pathways in the orthogonal and diagonal directions) means that,  in
        # the long run, the flood fill will even out and -importantly- it should
        # have some recognizable pattern every N steps (where N is the period
        # of the tiling - 11 in the sample, 131 in the input). this is because
        # we enter new tiles every N steps, and although it's hard to predict
        # exactly how many of the squares in each tile we'll have visited,
        # it should be the same number in every tile (or rather, there should
        # be some small set of repeated tile-states, and we should be able to
        # predict how many of each tile-state there will be).
        #
        # SO the first step is to get the difference between the current
        # number of reachable squares and the number N steps ago.
        #
        #       d1(i) = f(i) - f(i - N)
        #
        # we then have to discover some pattern in that sequence.
        # we know the number of reachable squares will grow roughly as the
        # square of the number of steps (because the map is 2D) so
        # we should be looking for a quadratic relation.
        # we can use second-order differences (see day 9) to do that.
        # d1 is already a first-order difference, so take the difference
        # between d1s to get a second-order difference. if our assumption is
        # correct, then d1 should be a linear sequence and d2 should be a
        # constant.
        #
        #       d2(i) = d1(i) - d1(i-N)
        #             = (f(i) - f(i-N)) - (f(i-N) - f(i-N-N))
        #
        #
        # there may be some unstability at the beginning of the simulation
        # so we need to wait until the d2 values for every step in the period
        # all agree.
        #
        # once it settles down, this gives us a set of N (11, 131, whatever)
        # equations we can use to predict the number of squares after any
        # future number of steps
        #
        d2 = 0
        if len(prev) >= 2*period:
            # find the second differece
            d2 = (n - prev[-period]) - (prev[-period] - prev[-2*period])
            diffs.append(d2)
        prev.append(n)
        print(i, n, d2, sep="\t", flush=True)
        if len(diffs) > period and all_same_value(diffs[-period:]):
            print("gotcha!")
            return prev, period, diffs[-1]
    assert False, "failed to find a stable pattern"
    return prev, period, None

def all_same_value(list):
    if len(list) < 1:
        return False
    x = list[0]
    return all(x == y for y in list)

def extrapolate(n, values, period, d2):
    if n < len(values):
        return values[n]
    quo, rem = divmod(n-len(values)+period, period)
    x = len(values)-period+rem
    y = values[-period+rem]
    d1 = y - values[-2*period+rem]
    while x < n:
        d1 += d2
        y += d1
        x += period
    assert x == n
    return y

def step(mask, old):
    # flood fill
    #draw(fill)
    fill = numpy.zeros(old.shape, dtype='uint8')
    y,x = fill.shape
    for i in range(y):
        f = (old[i] == 1)
        f = numpy.roll(f, 1) | numpy.roll(f, -1)
        if i > 0:    f |= (old[i-1] == 1)
        if i < y-1:  f |= (old[i+1] == 1)
        f &= (mask[i] == 0)
        #print(old, f)
        if f.any():
            fill[i, f] = 1

    assert fill.any()
    return fill


map = read_map(sys.stdin)
solve(map)
solve2(map)