{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Problem statement\n\nThis Kernel is an attempt to find some clues to generalized problem, that can give solution for this competition. New problem formulation is:\n* Define maximum \"ancestral\" depth (MAD) as maximum number of turns we can go backwards in time in game of life among all possible parent states\n* Find algorithm that for any board state will return one of the immediate parent states with maximum \"ancestral\" depth\n\nSo all states that don't have parent state will have MAD = 0.\nThen all states that have parents only of MAD=0, will have MAD=1. And so on..\n\nFor example if some board state has grandparents but not grand-grandparents (MAD=2), hypotetical algorithm will return any parent state that also has parent.\n\nWe know that all \"stop\" states in the train and test datasets have MAD of at least 6-10 (\"delta(1-5) + warmup(5)\").\nSo algorithm that returns parent state with maximum MAD among all parents and does so for any initial state, will guarantee, that it can be applied recursively at least \"delta\" times, giving us at least one of desired \"start\" states.\nLet's call it Algorithm of God for now ;) \n\nHere is one statement valid for any finite board, that could give a clue for desired algorithm:\n* Any initial state has MAD=\"infinity\", if and only if there is a cycle of states that includes initial state\n(I'll leave proof of this statement to the Reader)\n\nLet's try to find more patterns or properties like the one above. \n\n\n# Reduce board size\n\nIn attempt to find patterns that can help formulate \"Algorithm of God\", I decided to reduce board size to be fully solvable in memory and still meaningful.\n* 1x1 the only cell is own neighbour 8 times\n* 2x2 cells have same other cells as 2 or 4 of their neighbours\n* 3x3 Every cell is a neighbour of any other cell\n* 3x4 Cells in each row share all same neighbours\n\nSo, 4x4 is a first board size that starts to show some properties that possibly can be generalized to larger boards.\nYou can find complete research on it in my previous Kernel:\nhttps://www.kaggle.com/elvenmonk/complete-4x4-board-analysis\n\nNext size to try is 5x5. It is significantly larger and longer to compute, so I had to optimize and cache \"hash\" function before moving on.\n\n**Note on hashing function:**\n\nThere are 2^25 possible states of 5x5 board. We can code each state as a single number in range 0-33554431.\nTo further reduce number of dirrerent states we will consider states equal if we can get one state from another by translation (E.g. using \"numpy.roll\"), rotation or flip around any axis"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"import numpy as np\nfrom scipy.signal import convolve2d\nfrom tqdm.auto import tqdm\nimport matplotlib.pyplot as plt\nimport networkx as nx\n\nSIZE = 5\nSIZE2 = SIZE**2\nM = 2**SIZE\nN = 2**SIZE2\nNM = N//M\n\nempty_board = np.zeros((SIZE,SIZE), dtype=bool)\nm2i = np.array([2**i for i in range(SIZE2)]).reshape(SIZE,SIZE)\nwindow = np.ones((3, 3))\n\ndef plot_3d(solution_3d: np.ndarray, title, size=3, max_cols=2*SIZE, has_mad=False, tile=(1,1), vmin=0, vmax=1):\n    N = len(solution_3d)\n    if N <= 0:\n        return\n    cols = min(N, max_cols)\n    rows = (N - 1) // max_cols + 1\n    plt.figure(figsize=(cols*size, rows*size))\n    if not has_mad:\n        plt.suptitle(title)\n    for t in range(N):\n        board = np.tile(solution_3d[t], tile)\n        plt.subplot(rows, cols, t + 1)\n        plt.imshow(board, cmap='binary', vmin=vmin, vmax=vmax)\n        if has_mad:\n            plt.title(f'MAD {t}')\n    plt.show()\n\ndef life_step(X: np.ndarray):\n    Y = convolve2d(X, window, mode='same', boundary='wrap')\n    return (Y == 3) | (X & (Y == 4)), Y\n\ndef print_base(a):\n    print(np.vectorize(lambda x: x[3:])(np.vectorize(hex)(a+N)))#, base=M)))","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"r = np.arange(SIZE)[:, None]\nc = SIZE*r\nrr = np.cumsum(2**r//2)[:, None] # 0,1,3,7\ncc = np.cumsum(2**c//M)[:, None] # 0, 1, 17, 273\nprint_base(np.hstack([r, c, rr, cc]))\n\nSr = np.sum(2**r)\nSc = np.sum(2**c)\nSd = sum(1 << (SIZE - 1) * (i + 1) for i in range(SIZE))\nprint(hex(Sr), hex(Sc), hex(Sd))\n\nScr = Sc * 2**r\nSrc = Sr * 2**c\n\nSrcc = Sr * cc\nScrr = Sc * rr\nScccSr = (Sc - cc) * Sr\nSrrrSc = (Sr - rr) * Sc\nprint_base(np.hstack([Scr, Src, Srcc, Scrr, ScccSr, SrrrSc]))\n\ndef h_flip(x):\n    return np.sum((x & Scr) << (SIZE - 1) >> 2*r, axis=0)\n\ndef v_flip(x):\n    return np.sum((x & Src) << (SIZE2 - SIZE) >> 2*c, axis=0)\n\ndef d_flip(x):\n    y = x & Sd\n    for i in range(1,SIZE):\n        y |= ((x & (Sd >> (i * SIZE))) << (i * (SIZE + 1))) | ((x & (Sd << (i * SIZE))) >> (i * (SIZE + 1)))\n    return y\n\ndef get_hash(n):\n    h = h_flip(n)\n    v = v_flip(n)\n    hv = h_flip(v)\n    n = np.vstack([n, h, v, hv, d_flip(n), d_flip(h), d_flip(v), d_flip(hv)])[:,None,None,:]\n    n = ((n & Srcc) << (SIZE2 - c)) + ((n & ScccSr) >> c)\n    n = n.T\n    n = ((n & Scrr) << (SIZE - r)) + ((n & SrrrSc) >> r)\n    return np.min(n, axis=(1, 2, 3))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Every state in Hash function has 5x5x8 equivalent states (various translations, rotations and flips). So calculations on entire numpy array don't fit in memory, thus we have to split them into parts:"},{"metadata":{"trusted":true},"cell_type":"code","source":"n = np.arange(N)\nh = np.arange(N)\nfor i in tqdm(range(M), total=M):\n    h[i*N//M:(i+1)*N//M] = get_hash(n[i*N//M:(i+1)*N//M])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# NetworkX DiGraph\n\nTransitions from previous to next board state can be represented as directed graph. Building NetworkX graph will allow us to utilize its various algorithms to analyze and traverse this graph."},{"metadata":{"trusted":true},"cell_type":"code","source":"board = (m2i & n[:,None,None]) != 0\nconv_board = np.zeros(board.shape, dtype=int)\nG = nx.DiGraph()\nfor n0 in tqdm(range(N), total=N, mininterval=1):\n    if G.has_node(h[n0]) and len(list(G.successors(h[n0]))) > 0:\n        continue\n    next_board, conv_board[n0] = life_step(board[n0])\n    n1 = np.sum(next_board * m2i)\n    G.add_edge(h[n0], h[n1])\nn, h = None, None\nprint(nx.info(G))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"First let's extract all cycles. Corresponding board states are the only ones to have infinite MAD. "},{"metadata":{"trusted":true},"cell_type":"code","source":"cyclic = [c for cycle in nx.simple_cycles(G) for c in cycle]\nfor node in cyclic:\n    G.remove_node(node)\nplot_3d(board[cyclic], f'Cyclic boards')\nplot_3d(conv_board[cyclic], f'Convolutions', vmax=8)\nprint(nx.info(G))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Amoung the cyclic boards we can actually see a [Glider](https://en.wikipedia.org/wiki/Glider_(Conway%27s_Life))! This was not possible in 4x4 board, so 5x5 seems to represent more generic case.\n\n##### What is really fascinating about 5x5 cyclic boards is that each of them can be \"tiled up\" to form 25x25 board which will be also cyclic! (with same peroid)"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_3d(board[cyclic], f'Some cyclic 25x25 boards', tile=(5,5))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now our graph is directed and acyclic and we can start extracting longest passes from it!"},{"metadata":{"trusted":true},"cell_type":"code","source":"mad0 = list(nx.isolates(G))\nfor node in mad0:\n    G.remove_node(node)\n\nwith tqdm(total=nx.number_of_nodes(G)) as pbar:\n    path = nx.dag_longest_path(G)\n    while len(path) > 1:\n        for node in path:\n            G.remove_node(node)\n        isolates = list(nx.isolates(G))\n        for node in isolates:\n            G.remove_node(node)\n        mad0 = mad0 + path[:1] + isolates\n        pbar.update(len(path) + len(isolates))\n        plot_3d(board[path], f'{len(path)}-step path', has_mad=True)\n        #plot_3d(conv_board[path], f'Convolutions', vmax=8)\n        path = nx.dag_longest_path(G)\nprint(nx.info(G))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As aresult we can see, which parent state we can choose for every non-cyclic state to be able to get as deep as MAD steps into the past. (Output not always correspond to actual parent, but rather equal state that can be obtained from exact parent by translations and flips)"},{"metadata":{"trusted":true},"cell_type":"code","source":"#L = (len(mad0) + 99) // 100\n#for l in tqdm(range(L), total=L):\n#    plot_3d(board[mad0[l*100:(l+1)*100]], f'MAD=0 boards, part {l}')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}