mersenneforum.org  

Go Back   mersenneforum.org > Fun Stuff > Puzzles

Reply
 
Thread Tools
Old 2020-06-03, 21:16   #45
uau
 
Jan 2017

2×37 Posts
Default

Here's the first version of the code I used for the challenge (original 10 days before the correction to 19):
Code:
#!/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 27=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:
Code:
#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 228 (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.
uau is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
2020 project goals gd_barnes Conjectures 'R Us 4 2020-09-23 08:59
U.S. Electile Dysentery 2020 ewmayer Soap Box 308 2020-09-14 23:50
March 2020 what Puzzles 1 2020-04-24 05:46
February 2020 what Puzzles 20 2020-03-04 07:55
January 2020 what Puzzles 21 2020-02-02 14:11

All times are UTC. The time now is 16:50.

Wed Sep 30 16:50:22 UTC 2020 up 20 days, 14:01, 0 users, load averages: 1.74, 1.83, 1.80

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.