{"cells":[{"metadata":{},"cell_type":"markdown","source":"In this kernel I decided to try PicoSAT for few reasons.\n\nFirst I think it's important for this competition to learn as much as possible about states that don't have previous state. That is what is called [Garden of Eden](https://en.wikipedia.org/wiki/Garden_of_Eden_(cellular_automaton)) at least in original version of Game of Life with unbounded board.\n\nOne of the brilliant examples is [Flower of Eden](https://www.conwaylife.com/wiki/Flower_of_Eden) - the smallest known rotationally symmetric Garden of Eden.\n![Flower of Eden](https://www.conwaylife.com/w/images/7/78/Orphan5.png)\n\nMore precisely configuration above is called *orphan* - a subset of cells such as no matter what other cell states are resulting state will be a Garden of Eden. This means that any 25x25 board that contains this or any other orphan configuration inside it, will not have previous state as well.\n\nMore over, for bounded 25x25 board, that we are given, most of the possible board states don't have previous state simply because:\n* there are many board states like empty board with variety of predecessors\n* total number of states stays the same on each step = 2^625\n\nEvery Garden of Eden is a \"dead end\" in solving each task, so we should learn patterns of these states and avoid them as much as possible.\n\nSome research of Gardens of Eden led to solution in 2016 of [Conway's Grandfather problem](https://www.conwaylife.com/wiki/Grandfather_problem) and this is where PicoSAT comes into play as a tool that was used to solve the problem, along with finding\n> \"A father and grandfather, but no great-grandfather\" pattern[7], and a \"father, grandfather and great-grandfather, but no great-great-grandfather\" pattern\n\nWhile link above can give a lot of useful information on its own, one of the take aways is that PicoSAT can be an efficient tool for this competition if used wisely.\n\nAs a python wrapper I'm going to use [pycosat](https://pypi.org/project/pycosat/) - first that I found. It is very easy to use, which is good for a starter Kernel.\n\n\nI took and modified few helper functions (pretty few) from other public notebooks, so don't be surprized, and I apollogize if I didn't mention original authors."},{"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\nimport pandas as pd\nfrom numba import njit, prange\nfrom scipy.signal import convolve2d\nfrom tqdm.auto import tqdm\nimport matplotlib.pyplot as plt\nimport pycosat\nimport random\n\nwindow = np.ones((3, 3))\n\ndef plot_3d(solution_3d: np.ndarray, title, size=3, max_cols=10, has_target=False):\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    plt.suptitle(title)\n    for t in range(N):\n        board = solution_3d[t]\n        plt.subplot(rows, cols, t + 1)\n        plt.imshow(board, cmap='binary')\n        plt.title('target' if t == 0 and has_target else f'state {t}')\n    plt.show()\n\ndef csv_to_numpy_list(df) -> np.ndarray:\n    return df[[ f'stop_{n}' for n in range(25**2) ]].values.reshape(-1,25,25)\n\ndef life_step(X: np.ndarray):\n    C = convolve2d(X, window, mode='same', boundary='wrap')\n    return (C == 3) | (X & (C == 4))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For SAT solvers we only need Test data: there is nothing to train and we always know if solution is correct."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"sample_submission_df = pd.read_csv('../input/conways-reverse-game-of-life-2020/sample_submission.csv', index_col='id')\ntest_df = pd.read_csv('../input/conways-reverse-game-of-life-2020/test.csv', index_col='id')\ndeltas = test_df['delta'].values\nboards = csv_to_numpy_list(test_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, let's try out some of the hepler functions above and look at some sample test board."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"for index in tqdm(range(1), total=1):\n    board = np.tile(boards[index], (2, 1, 1))\n    board[1] = life_step(board[0])\n    plot_3d(board, 'sample board')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, finally we are getting to the solver itself.\n\nLooks like pycosat requires input in [Conjunctive normal form](https://en.wikipedia.org/wiki/Conjunctive_normal_form) and does not support weights, optional clauses or native cardinality constraints. It can present some problems later when we would want to limit the search space or prefer some solutions to another. But for now the simpler the better.\n\nWhile trying to build CNF expression for Game of Life, I found this article, that helped me to reduce number of clauses and variables a lot:\nhttps://sat-smt.codes/SAT_SMT_by_example.pdf (Pages 501-504)"},{"metadata":{"trusted":true},"cell_type":"code","source":"SIZE = 25\n\nempty_board = np.zeros((SIZE,SIZE), dtype=bool)\n\n# Computes number of variables.\n# I negated variable states after some testing, because solver\n# seems to prefer True as default state for solution variables,\n# while I want cells to be empty whenever possible\ndef v(c):\n    return -(SIZE * c[0] + c[1] + 1)\n\ndef dead_clauses(res, c, x):\n    # if cell is dead, there was not exactly 3 alive neighbours (56 clauses)\n    for i1 in range(0, 6):\n        for i2 in range(i1+1, 7):\n            for i3 in range(i2+1, 8):\n                a = [v(x[i]) for i in range(8)]\n                a[i1], a[i2], a[i3] = -a[i1], -a[i2], -a[i3]\n                res.append(a)\n    # if cell is dead and was alive, was not 2-3 alive neighbours (28 clauses)\n    for i1 in range(0, 7):\n        for i2 in range(i1+1, 8):\n            a = [v(x[i]) if i < 8 else -v(c) for i in range(9)]\n            a[i1], a[i2] = -a[i1], -a[i2]\n            res.append(a)\n\ndef live_clauses(res, c, x):\n    # if cell is alive, there was less then 4 alive neighbours (70 clauses)\n    for i1 in range(0, 5):\n        for i2 in range(i1+1, 6):\n            for i3 in range(i2+1, 7):\n                for i4 in range(i3+1, 8):\n                    #from each 4 at least 1 was dead\n                    res.append([-v(x[i1]), -v(x[i2]), -v(x[i3]), -v(x[i4])])\n    # if cell is alive and was dead, there was more than 2 alive (less than 6 dead) neighbours (28 clauses)\n    for i1 in range(0, 7):\n        for i2 in range(i1+1, 8):\n            a = [v(x[i]) if i < 8 else v(c) for i in range(9) if i != i1 and i != i2]\n            res.append(a)\n    # if cell is alive, there was more than 1 alive (less than 7 dead) neighbours (8 clauses)\n    for i1 in range(0, 8):\n        a = [v(x[i]) for i in range(8) if i != i1]\n        res.append(a)\n\ndef board_clauses(board, use_opt = True):\n    res, opt1, opt2 = [], [], []\n    for i in range(SIZE):\n        for j in range(SIZE):\n            x = [((i + k % 3 - 1) % SIZE, (j + k // 3 - 1) % SIZE) for k in range(9) if k != 4]\n            if board[i,j]:\n                live_clauses(res, (i, j), x) # 106 clauses\n            else:\n                dead_clauses(res, (i, j), x) # 84 clauses\n                if use_opt:\n                    y = [((i + k % 5 - 2) % SIZE, (j + k // 5 - 2) % SIZE) for k in range(25)]\n                    if sum(board[ii,jj] for ii,jj in y) < 1: # No alive neighbours\n                        res.append([-v((i, j))]) # Very dead space MUST stay dead! (1 clause)\n                    elif sum(board[ii,jj] for ii,jj in x) < 1: # No alive neighbours\n                        opt1.append([-v((i, j))]) # Dead space should stay dead! (1 clause)\n                    elif sum(board[ii,jj] for ii,jj in x) < 2: # Too few alive neighbours\n                        opt2.append([-v((i, j))]) # Dead space should stay dead! (1 clause)\n\n    return res, opt1, opt2","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Depth(delta)-1 tasks can be solved relatively easy:"},{"metadata":{"trusted":true},"cell_type":"code","source":"N = len(deltas)\nscore = 0\nfor n in tqdm(range(N), total=N):\n    clauses, opt1, opt2 = board_clauses(boards[n], use_opt = False)\n    solution = pycosat.solve(clauses)\n    if isinstance(solution, str):\n        print(f'{n} not solved!')\n        continue\n    board = np.array(solution[:SIZE**2]) < 0\n    sample_submission_df.loc[test_df.index[n]] = 1 * board\n    board = life_step(board.reshape(SIZE,SIZE))\n    d = np.sum(board ^ boards[n])\n    score += d / 625\nprint(score/N)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It takes just a bit more than half a second to solve delta-1 task in average.\n\nNext, I've built step-by-step solver for individual tasks, with ability to gradualy exclude optional clauses and fallback to previous step when there is no solution.\n\nLet's try it out for first tasks."},{"metadata":{"trusted":true},"cell_type":"code","source":"N = 1\nfor n in tqdm(range(N), total=N):\n    T = min(deltas[n], 3)\n    board = np.tile(empty_board, (T+1, 1, 1))\n    board[0] = boards[n]\n    solvers = [None for _ in range(T)]\n    opt = [None for _ in range(T)]\n    os = [0 for _ in range(T)]\n    oe = [0 for _ in range(T)]\n    t = 0\n    while 0 <= t and t < T:\n        if solvers[t] is None:\n            clauses, opt1, opt2 = board_clauses(board[t])\n            solution = pycosat.solve(clauses)\n            if not isinstance(solution, str):\n                if t == T - 1:\n                    print(t, '!!', end=\" \")\n                    t += 1\n                    board[t] = np.array(solution[:SIZE**2]).reshape(SIZE,SIZE) < 0\n                    continue\n                else:\n                    print(t, '??', end=\" \")\n                    random.shuffle(opt1)\n                    random.shuffle(opt2)\n                    opt[t] = opt1 + opt2\n                    os[t] = len(opt[t])+1\n                    oe[t] = len(opt[t])+1\n                    solvers[t] = pycosat.itersolve(clauses + opt[t])\n                    print(len(opt[t]), end=\" \")\n        try:\n            solution = next(solvers[t])\n            if oe[t] - os[t] > 1:\n                os[t] = (os[t]+oe[t])//2\n                solvers[t] = pycosat.itersolve(clauses + opt[t][:(os[t]+oe[t])//2])\n                print((os[t]+oe[t])//2, end=\" \")\n            else:\n                print(t, '++', end=\" \")\n                t += 1\n                board[t] = np.array(solution[:SIZE**2]).reshape(SIZE,SIZE) < 0\n        except Exception as err:\n            if solvers[t] is not None and (os[t]+oe[t])//2 > 0:\n                if os[t] == oe[t]:\n                    os[t] = 0\n                    oe[t] = len(opt[t])+1\n                elif oe[t] - os[t] > 1:\n                    oe[t] = (os[t]+oe[t])//2\n                else:\n                    os[t] -= 1\n                    oe[t] -= 1\n                solvers[t] = pycosat.itersolve(clauses + opt[t][:(os[t]+oe[t])//2])\n                print((os[t]+oe[t])//2, end=\" \")\n            else:\n                print(t, '--', end=\" \")\n                solvers[t] = None\n                opt[t] = None\n                t -= 1\n    if t == T:\n        sample_submission_df.loc[test_df.index[n]] = 1 * board[T].ravel()\n    plot_3d(board, index)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So it does solve some tasks already!\n\nObviously performance is not good enough and a lot of optimization needs to be done for deep backwards search (3+ steps).\n\nBut depth-1 tasks can be solved relatively easy:"},{"metadata":{},"cell_type":"markdown","source":"Finally, let's save resulting submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission_df.to_csv(\"submission.csv\", index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_df = pd.read_csv('./submission.csv', index_col='id')","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}