April 2020

[uau]

Here's the first version of the code I used for the challenge (original 10 days before the correction to 19):

#!/usr/bin/python3

import numpy as np
import itertools
import sys

def calc(adj):
    step = np.identity(128)
    idx = np.arange(128)
    for i in range(8):
        if i == 0:
            rows = np.ones(128, dtype=bool)
        else:
            rows = (idx & (1 << (i-1))).astype(bool)
        for j in range(1, 8):
            if not adj[i, j]:
                continue
            cols = (idx & (1 << (j-1)) == 0)
            x = step[np.ix_(rows, cols)] * 0.1
            step[np.ix_(rows, cols)] -= x
            step[np.ix_(rows, cols==False)] += x
    step = np.linalg.matrix_power(step, 10)
    return step[0, 127]

def marksmall(toosmall, perms, i):
    for j in range(28):
        if not i & (1<<j):
            continue
        toosmall[perms[i & ~(1<<j)]] = 1

def main():
    toosmall = np.zeros(2**28, dtype=np.int8)
    perms = np.fromfile('/tmp/perms.data', dtype=np.int32)
    # performance optimization - "perms[i] < i" slow for Python int 'i'
    canon = perms == np.arange(0, 1<<28, dtype=np.int32)
    for i in np.nonzero(canon)[0][::-1]:
        p = 0
        if not toosmall[i]:
            b = bin(i)[2:].rjust(28, '0')[::-1]
            adj = np.zeros((8, 8), dtype=np.int)
            adj[iu] = [int(x) for x in b]
            adj += adj.T
            p = calc(adj)
            print(i, p)
            sys.stdout.flush()
        if p <= 0.7:
            marksmall(toosmall, perms, i)

# create "/tmp/perms.txt" for C program that generates "/tmp/perms.data"
perms = [[0]+list(x) for x in itertools.permutations(range(1, 8))]
adj = np.zeros((8, 8), dtype=np.int)
iu = np.mask_indices(8, np.triu, 1)
adj[iu] = range(28)
adj += adj.T
perms2 = []
for p in perms:
    a2 = adj[p]
    a2 = a2[:, p]
    perms2.append(a2[iu])
with open('/tmp/perms.txt', 'w') as f:
    for p in perms2:
        f.write(' '.join(str(x) for x in p)+'\n')

main()
# forum editor tends to break last line indentation if it's just main()?

The calc() function calculates the probability for one given adjacency matrix. It creates a 128*128 matrix where each element is the probability to move from state A to state B in one day (there are 2⁷=128 possibilities for the infection states of the 7 nodes other than the initial one). Raising this matrix to power 10 gives the probabilities for 10 days, and the element for 0 to 127 means transition from (all but initial one uninfected) to (all infected). The calculation of the matrix is basically: Start from an identity matrix. For each bit position i: take rows where bit i is set, and for each bit position j where there exists an edge(i, j): take 1/10 of values of columns where bit j=0, and add to where bit j=1.

The program avoids duplicate work by trying all permutations of vertices other than the first one, and only testing one graph from each equivalence class. Since this is kind of slow in Python, I implemented the bit shuffling in a C program:

#include <stdio.h>
#include <stdbool.h>
#include <assert.h>

int a[5040][28];
int seen[1<<28];

int main(void) {
    for (int i = 0; i < 5040; i++)
        for (int j = 0; j < 28; j++) {
            int k;
            scanf("%d", &k);
            a[i][j] = 1 << k;
        }
    for (int i = 1; i < (1<<28); i++) {
        if (seen[i])
            continue;
        for (int j = 0; j < 5040; j++) {
            int s = 0;
            for (int k = 0; k < 28; k++)
                s += (1 << k) * (bool)(i & a[j][k]);
            seen[s] = i;
        }
    }
    int l = fwrite(seen, sizeof(int), 1<<28, stdout);
    assert(l == 1<<28);
    return 0;
 }

All permutations of the vertex order are equivalent to some permutation of adjacency matrix bits. The Python program writes those permutations to file "/tmp/perms.txt"; the C program can then be used to process that to "/tmp/perms.data" (it uses stdin/stdout, so it's called like "./a.out </tmp/perms.txt >/tmp/perms.data"). The data file can then be interpreted as a table of integers, mapping each adjacency matrix to the canonical one of its class (each matrix mapped to an integer by reading the bit values in the upper triangle in order). The Python program thus uses simple table lookup to get the canonical form of any adjacency matrix.

With the original 10 days limit, the problem is solvable within a few seconds. Most adjacency matrices have a probability below the target 0.7, and turning off any bits lowers the probability. Thus the program searches through all possible matrices in descending order, and when a matrix gets a probability <= 0.7, then it knows that all matrices obtained by turning off a bit have an even lower probability. The search space is originally 2²⁸ (28=8*7/2 bits in the adjacency matrix for edges), but most are ruled out as either not canonical after permuting vertex order, or by the <= 0.7 logic (implemented by the marksmall() function).

If you increase the number of days from 10 to 19, the smallness cutoff is no longer as effective, and it'll take minutes instead of seconds.