{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Reversing game of life problem \n## A NP-Hard problem\nThe problem of reversing game of life (GoL) is NP-hard because it can be framed as a SAT problem. As a consequence, it is very hard to use deep learning to solve it. A lot of public kernels uses various CNNs and ends up at MAE of around .11 at most. \n\n## The approach of this kernel 👀\nThis kernel implements a probabilistic differentiable version of GoL in order to approximate a solution using gradient descent. It is possible to find perfect solutions for simple grids even with multiple deltas. However, most of the time the grids are too complex to converge but the solutions obtained are better than NNs ones. Then for the very complex grids, just submit whatever is better between a blank grid and the grid of origin.\n\n👉 This simple code can easily go in the range 0.09 - 0.10 MAE with 2 hours GPU, and even bellow with hyperparameters tunning and minor modifications."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"\"\"\"\n    Imports\n\"\"\"\n\nimport numpy as np\nimport pandas as pd\nfrom random import randint\nfrom math import sqrt\nimport torch\nimport torch.nn.functional as Func\nimport matplotlib.pyplot as plt\nimport time\nimport gc","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let's keep track of the time to stay in the 2 hour GPU limit\nts = time.time()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Classical bloc of code to use GPU if it's available\nif torch.cuda.is_available(): \n    dev = \"cuda:0\" \nelse:  \n    dev = \"cpu\"\ndevice = torch.device(dev)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n    Hyperparameters\n\"\"\"\n# Max number of forward GoL steps\nmax_delta = 5\n# To debug the kernel, turn it to true and it wont run on the all test set\ntest_mode = False\n# Learning rate\nlr = 1.0\n# Max steps per sample\nmax_steps = 1000\n# Batch size\nbatch_size = 5000\n# Runtime in seconds with a little margin\nt_run = 3600 * 2 - 30 * 60","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n    Game of life functions (With batch version)\n\"\"\"\n\ndef get_padded_version_n(X):\n    \"\"\"\n        Apply circular padding to a batch of grids\n    \"\"\"\n    X_pad = np.zeros((X.shape[0], X.shape[-2] + 2, X.shape[-1] + 2), dtype=X.dtype)\n    X_pad[:, 1:-1,1:-1] += X\n    \n    X_pad[:, 0, 1:-1] = X[:, -1, :]\n    X_pad[:, -1, 1:-1] = X[:, 0, :]\n    \n    X_pad[:, 1:-1, 0] = X[:, :, -1]\n    X_pad[:, 1:-1, -1] = X[:, :, 0]\n    \n    X_pad[:, 0, 0] = X[:, -1, -1]\n    X_pad[:, 0, -1] = X[:, -1, 0]\n    X_pad[:, -1, 0] = X[:, 0, -1]\n    X_pad[:, -1, -1] = X[:, 0, 0]\n    \n    return X_pad\n\ndef nConv2d_sw_3x3(X):\n    \"\"\"\n        Convolve with a 3x3 ones filter a batch of grids\n    \"\"\"\n    X_pad = get_padded_version_n(X)\n    N = np.zeros_like(X_pad)\n    \n    N[:, 1:, 1:] += X_pad[:,:-1,:-1]\n    N[:, 1:, :] += X_pad[:,:-1,:]\n    N[:, 1:, :-1] += X_pad[:,:-1,1:]\n\n    N[:, :, 1:] += X_pad[:,:,:-1]\n    N[:, :, :] += X_pad[:,:,:]\n    N[:, :, :-1] += X_pad[:,:,1:]\n\n    N[:, :-1, 1:] += X_pad[:,1:,:-1]\n    N[:, :-1, :] += X_pad[:,1:,:]\n    N[:, :-1, :-1] += X_pad[:,1:,1:]\n    \n    N = N[:,1:-1,1:-1]\n    \n    return N\n\ndef life_step(X):\n    \"\"\"\n        Forward iteration of game of life\n    \"\"\"\n    N =  nConv2d_sw_3x3(X) - X\n    return np.logical_or(N == 3, np.logical_and(X, N==2)).astype(np.uint8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Probabilistic Forward Differentiable Game of life\nIn order to use gradient descent we need a differentiable and continuous version of the game of life step function. A good candidate for this purpose is the probabilistic version of GoL.\n\n## PFDGol formula 🤓\nPFDGol takes as input a grid of values in the range [0, 1]. Thoses values represent the probability of the given cell to be alive. As the classical version of GoL is computed using the number of neighbors and the current state of the cell, the formula of the probabilistic version is:\n\n$$\n    P(N = 2).I + P(N = 3)\n$$\n\nWith $I$ being the probability of the cell being alive on the initial grid and $N$ being the number of alive neighbors of the cell. As the first one is given, we just have to find a way to compute $p(N=n)$.\n\n## Probability of a cell having n alive neighbors\nTo compute $p(N=n)$, we just consider the likelyhood of all possible scenariis for the neighborhood of the cell. By scenario, I mean all possible configuration of neighborhood to obtain n alive neighbors.\nUnder independancy of cells in initial state assumption, we have:\n\n$$\n    P(S^{x,y}) = \\prod_{(i, j) \\in \\{-1, 0, 1\\}^2 , (i, j) \\neq (0, 0)} S_{i,j} * I_{x + i, y + j}\n$$\n\nWith $S^{x,y}$ being the occurence of a scenario at cell $x, y$ and $S_{i,j}$ being the state of the cell in the considered scenario.\nThen to obtain $p(N=n)$, just sum up among all possible scenariis (again under assumptions of independancy)\n\n## Little trick to go faster 🚀\nIt's possible to unlock a considerable speedup of the function by using to our advantage the highly optimized pytorch's conv2d function. This quite simple trick consist of convolving a 3x3 ones filter to the $log$ of the probabilities. To do so, we need to separate the calulation in three phases:\n1. calculate $P(N \\geq n)$:\n\n$$\n    P(N \\geq n) = Conv2d(log(I), f)\n$$\n\nWith $f$ being a 3x3 filter representing the considered scenario. ($f_{i, j} = 1$ when $cell_{i, j}$ is alive in the scenario, and 0 otherwise). Here, you multiply (add in log space) the probabilities of the cells at the locations of the alive cells in the scenario being actually alive.\n2. calculate $P(N \\leq n)$:\n\n$$\n    P(N \\leq n) = Conv2d(log(1 - I), 1 - f)\n$$\n\nHere it is the opposit. You multiply the probabilities of the cells at the locations of the dead cells in the scenario being actually dead.\n3. Merge\n$$\n    p(N = n) = P(N \\leq n) \\cap P(N \\geq n)\n$$\n\n## Limitations 😕\nIf this model is correct for one forward iteration, it is not (I think so, but didn't try to prove it) when you apply it recursively multiple times (as our independancy assomptions become wrong).\n\nAlso as gradient descent is based on partial derivatives there are a lot of non-optimum local minimas.\nFor example in the case of:\n```python\nI = [\n    [0, 0, 0],\n    [0, 0, 0],\n    [0, 0, 0]\n]\n```\nyou will obtain $\\frac{\\partial MAE(PFDGol(I), F)}{\\partial I_{x, y}} = 0, \\forall F, x, y$ as you need at least 3 alive neighbors to turn a dead cell alive in the next step."},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n    Probabilistic (differentiable) Game of life functions with pytorch\n\"\"\"\n\n# All possible configurations of neighborhood (flattened)\nbits=np.unpackbits(np.expand_dims(np.arange(2**8, dtype=np.uint8), axis=1), axis=1)\n\nfilters, filters_i = {}, {}\nfor k in [2, 3]:\n    # All possible configurations of a neighborhood containing k cells alive\n    poss_n = bits[bits.sum(axis=-1) == k]\n    # All possible neighborhood for a cell with k neighbors (reshaped on 3x3)\n    poss_neigh = np.insert(poss_n, 4, 0, axis=1).reshape(len(poss_n), 3, 3)\n    # Save it as a torch tensor\n    filters[k] = torch.tensor(poss_neigh, dtype=torch.float32, device=device).view(len(poss_n), 1, 3, 3)\n    # Negation of the previous\n    poss_neigh_comp = np.insert(1 - poss_n, 4, 0, axis=1).reshape(len(poss_n), 3, 3)\n    # Save it as a torch tensor\n    filters_i[k] = torch.tensor(poss_neigh_comp, dtype=torch.float32, device=device).view(len(poss_n), 1, 3, 3)\n\ndef p_n_f(n, I):\n    \"\"\"\n        Return the probability of each cell having n neighbors\n    \"\"\"\n    ## Circular padding + log (Trick: addition becomes multiplication)\n    # Complementary of alive probs in initial state\n    p_F_0_ln = Func.pad(torch.log(1 - I), pad=(1, 1, 1, 1), mode='circular')\n    # Probability of cells to be alive on initial state\n    p_F_1_ln = Func.pad(torch.log(I), pad=(1, 1, 1, 1), mode='circular')\n    \n    ## As we are in log space, convolution now become a likelyhood calculation among neighbors\n    ## Each filter represent a configuration of neigborhood\n    ## Filters are NOT in log space so we are just calulcating the likelyhood of partial neighborhood\n    # Probability of a cell having less or n neighbors (convolve complement of I to complement of neighborhood filters)\n    p_Nf_lte_n_ln = Func.conv2d(p_F_0_ln, filters_i[n])\n    ## As we are convolving with the complement\n    # Probability of a cell having n or more neighbors (convolve I to neighborhood fliters)\n    p_Nf_se_n_ln = Func.conv2d(p_F_1_ln, filters[n])\n    # Now merge the twos : p(N <= n) * p(N >= n) -> p(N == n)\n    p_Nf_n_ln = p_Nf_lte_n_ln + p_Nf_se_n_ln\n    # Now we sum up the likelyhood of all possible n neighbors configurations, and obtain the good result\n    return torch.sum(torch.exp(p_Nf_n_ln), dim=-3)\n\ndef p_life_step(I):\n    \"\"\"\n        Probabilistic differentiable forward iteration of game of life\n    \"\"\"\n    # I use logs in here so I tried to handle 0s and 1s by replacing them with very little values, but it doesn't seems to improve the performances\n    n = len(I)\n    I = I.view((n, 1, 25, 25))\n    p_n_2 = p_n_f(2, I)\n    p_n_3 = p_n_f(3, I)\n    # The well known formula\n    return p_n_2 * I.view(n, 25, 25) + p_n_3\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n    Test set\n\"\"\"\n# Let's read the dataset\ntest = pd.read_csv(\"../input/conways-reverse-game-of-life-2020/test.csv\", index_col='id').values\n# Those are the deltas\ntest_delta = test[:,0]\n# Thoses are the states to reach\nX_test = test[:, 1:626].reshape((-1, 25, 25))\n\nprint(test_delta.shape, X_test.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n    Choose the input of the function based on mode of the kernel\n\"\"\"\n\nif test_mode:\n    n = 1000\n    sample_idx = np.random.permutation(np.arange(len(X_test)))[:n]\n    F = X_test[sample_idx]\n    d = test_delta[sample_idx]\nelse:\n    F = X_test\n    d = test_delta","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n    Core of the code (Gradient descent)\n\"\"\"\n\n# Initial states\nIf = np.zeros_like(F)\n# Selector\nselector = np.arange(len(If))\n# Score of each saved initial state to keep only if we found a better one\nS = np.ones_like(selector, dtype=np.float32)\n\nwhile 1:\n    # As memory is limited by Autograd, we have to execute batch per batch\n    # This is a boolean mask to keep track of wich grid is in the current batch\n    runnings = np.zeros(len(S)).astype(np.bool_)\n    # Choose the boards that have the worst score\n    runnings[np.arange(len(S))[S >= np.quantile(S, 1 - (batch_size / len(S)))][:batch_size]] = True\n    # Targets of the current batch\n    F_step = F[runnings]\n    # Initial states (autograd variable) (initialized randomly as zero grid is a local optimum)\n    I = torch.rand(F_step.shape, dtype=torch.float32, requires_grad=True, device=device)\n    # Targets of the current batch as tensors\n    Ft = torch.tensor(F_step, dtype=torch.float32, device=device)\n    # Optimizer is classic gradient descent with momentum (but it's not stochastic here)\n    optimizer = torch.optim.SGD([I], lr=lr)\n    \n    for s in range(max_steps):\n        Fh = torch.zeros(I.size(), dtype=torch.float32, device=device)\n        # Uses sigmoid to squeeze the values in ]0, 1[ as we use log\n        current = torch.sigmoid(I)\n        # Forward passes\n        for i in range(max_delta):\n            current = p_life_step(current)\n            current_m = d[runnings] == i + 1\n            # Save only those that needs i + 1 forward steps\n            Fh[current_m] = current[current_m]\n        ## There are numerous NaNs because of log so we want to backpropagate only the non-NaNs\n        calculable_m = (torch.sum(torch.sum(torch.isnan(Fh), dim=-1), dim=-1) == 0).cpu().numpy()\n        loss = torch.sum(torch.square(Fh[calculable_m] - Ft[calculable_m]))\n        loss.backward()\n\n        with torch.no_grad():\n            optimizer.step()\n            optimizer.zero_grad()\n            #not_calculable = selector[runnings][np.logical_not(calculable_m)]\n            #I[not_calculable] = torch.rand((len(not_calculable), F.shape[-2], F.shape[-1]), dtype=torch.float32, device=device)\n            # Evaluate our grids at a given frequecy\n            if s % 100 == 0:\n                # Compute absolute error\n                AE = torch.abs(Fh[calculable_m]  - Ft[calculable_m] )\n                print(\"Step {}/{} MAE of step: {}\".format(s, max_steps, torch.mean(AE).item()))\n                # Apply deterministic forward GoL in the same way as before\n                Is = torch.sigmoid(I[calculable_m]).cpu().numpy()\n                Iv = (Is > 0.5).astype(np.uint8)\n                Fv = np.zeros_like(Iv)\n                current = Iv\n                for i in range(max_delta):\n                    current = life_step(current)\n                    current_m = d[runnings][calculable_m] == i + 1\n                    Fv[current_m] = current[current_m]\n                # Compute MAE with deterministic version of GoL\n                mae = np.mean(np.abs((Fv - F[runnings][calculable_m])), axis=(-2, -1))\n                # Mask of the better ones\n                local_m = mae < S[runnings][calculable_m]\n                # Uses the selector to compute a mask of a mask\n                If[selector[runnings][calculable_m][local_m]] = Iv[local_m]\n                S[selector[runnings][calculable_m][local_m]] = mae[local_m]\n                print(\"MAE w/ round: {}\".format(np.mean(mae)))\n                print(\"Overall MAE so far: {}\".format(S[S < 1].mean()))\n        # Some memory management   \n        del Fh\n    # Stop if it's time\n    if time.time() - ts >= t_run:\n        break\n        \n    # Some other memory management   \n    del I\n    del Ft\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n    Compute various naive submissions\n\"\"\"\n\n# Score of a blank submission\nblank_score = X_test.mean(axis=(-2, -1))\n\n# Stop submission\nstop_scores = np.zeros_like(blank_score)\n# Forward passes on the final states \ncurrent_step = np.copy(X_test)\nfor i in range(max_delta):\n    current_step = life_step(current_step)\n    stop_scores[test_delta==i+1] = 1 - (current_step[test_delta==i+1] == X_test[test_delta==i+1]).mean(axis=(-2, -1))\n\nprint(blank_score.mean())\nprint(stop_scores.mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n    Merge the naive submissions\n\"\"\"\n\nnaive_submission = np.zeros_like(X_test)\nstop_best_mask = stop_scores < blank_score\nnaive_submission[stop_best_mask] = X_test[stop_best_mask]\nnaive_scores = np.copy(blank_score)\nnaive_scores[stop_best_mask] = stop_scores[stop_best_mask]\n\nprint(naive_scores.mean())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n    Compute the final submission using main and naive ones\n\"\"\"\n\ndef get_optimized_solution_and_score(y, delta, labels, naive_submission, naive_scores):\n    \"\"\"\n        Return best submission using ys and naive submission\n    \"\"\"\n    print(\"Scores from prediction\")\n    res = np.zeros_like(y)\n    scores = []\n    current_step = np.copy(y)\n    for i in range(5):\n        current_step = life_step(current_step)\n        scores.append(current_step[delta==i+1] == labels[delta==i+1])\n        print(\"Actual score for delta={}: {}\".format(i+1, scores[-1].mean()))\n        res[delta == i+1] = current_step[delta == i+1]\n    print(\"Actual LB score = \", 1 - (res == labels).mean())\n    y_scores = 1 - (labels == res).mean(axis=(-2, -1))\n    naive_mask = naive_scores < y_scores\n    y_final = np.copy(y)\n    y_final[naive_mask] = naive_submission[naive_mask]\n    f_score = np.copy(y_scores)\n    f_score[naive_mask] = naive_scores[naive_mask]\n    print(\"LB score estimation = \", f_score.mean())\n    return y_final, f_score.mean()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n    Make the submission\n\"\"\"\ny_final, score = get_optimized_solution_and_score(If, test_delta, X_test, naive_submission, naive_scores)\nprint(score, naive_scores.mean() - score)\n\nsubmission = pd.read_csv(\"../input/conways-reverse-game-of-life-2020/sample_submission.csv\", index_col='id')\nsubmission_values = y_final.reshape((len(submission), 625))\nfor i, col in enumerate(submission.columns):\n    submission[col].values[:] = submission_values[:, i]\nsubmission.to_csv(\"submission.csv\")","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}