{"cells":[{"metadata":{},"cell_type":"markdown","source":" # Iterative CNN Approach with optimized thresholds\n \nI use the great [kernel](https://www.kaggle.com/yakuben/crgl2020-iterative-cnn-approach) of Yakubenko Oleksii for model architecture and training.\n\nInstead of setting a common threshold of y_pred = (y_pred_continuous > 0.5) for binarizing the predictions,\nI determine the threshold for each board individually:\n\n- Go through a set of possible thresholds `[0, 0.1, ..., 1]`.\n- Binarize the board according to the threshold, i.e. `y_pred = (y_pred_continuous > threshold)`.\n- Make N steps `y_pred_N = make_move(y_pred, N)` and compute `accuracy(y_pred_N, y_true_N)`. Note that `y_true_N` is the input to the neural network, i.e. the state of the game at time `N`.\n- Use the threshold which maximizes the accuracy.\n \n \n \n \n \n This kernel uses the cythonized make_board method from [here](https://www.kaggle.com/ptyshevs/cnn-for-reversing-game-of-life)\n and [here](https://www.kaggle.com/jpmiller/demo-cython-generator-and-keras-cnn).\n \n ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Task overview/Game Rules\n\n*The game consists of a board of cells that are either on or off. One creates an initial configuration of these on/off states and observes how it evolves. There are four simple rules to determine the next state of the game board, given the current state:*\n\n* #### Overpopulation: if a living cell is surrounded by more than three living cells, it dies.\n* #### Stasis: if a living cell is surrounded by two or three living cells, it survives.\n* #### Underpopulation: if a living cell is surrounded by fewer than two living cells, it dies.\n* #### Reproduction: if a dead cell is surrounded by exactly three cells, it becomes a live cell.\n\n![](https://natureofcode.com/book/imgs/chapter07/ch07_01.png)","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Approach overview\n\n\n### Iterative model\nThe basic idea is to train a model that can predict start state from stop state(with delta=1) and apply this model delta times.\n\nThe model itself can be any (not only neural networks).\n\n\n![](https://i.ibb.co/y0fMSCg/simple-model-diagram.png)\n\n\n### Iterative model with Encoder/Decoder\n\nUsing an encoder and a decoder allows more information to be conveyed to the next step (not just one channel).\n\n![](https://i.ibb.co/MRvZbYj/model-diagram.png)\n\n\n### Iterative CNN model with Encoder/Decoder\n\nThe stop and start state can be viewed as single-channel images(so we can use CV approaches)\n\nInput - 1 channel image\n\nEncoder in - 1 channel image | out - N channel image\n\nReverseOneIterationModel - inp N channel image | out - N channel image\n\nDecoder in - N channel image | out - 1 channel image","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Model","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch\nimport torch.nn as nn\n\n\nclass OneIterationReverseNet(nn.Module):\n    def __init__(self, info_ch, ch):\n        super().__init__()\n        self.relu = nn.ReLU()\n        self.conv1 = nn.Conv2d(info_ch, ch, 5, padding=4, padding_mode='circular')\n        self.conv2 = nn.Conv2d(ch, ch, 3, )\n        self.conv3 = nn.Conv2d(ch, info_ch, 3)\n        \n        \n    def forward(self, inp):\n        x = self.relu(self.conv1(inp))\n        x = self.relu(self.conv2(x))\n        x = self.relu(self.conv3(x))\n        return x\n      \n        \nclass ReverseModel(nn.Module):\n    def __init__(self, info_ch=64, ch=128):\n        super().__init__()\n        self.relu = nn.ReLU()\n        self.encoder = nn.Conv2d(1, info_ch, 7, padding=3, padding_mode='circular')# you can use other model\n        self.reverse_one_iter = OneIterationReverseNet(info_ch, ch)# you can use other model\n        self.decoder = nn.Conv2d(info_ch, 1, 3, padding=1, padding_mode='circular')# you can use other model\n        \n    \n    def forward(self, stop, delta):\n        x = self.relu(self.encoder(stop-0.5))\n        \n        for i in range(delta.max().item()):\n            y = self.reverse_one_iter(x)\n            \n            # this 2 lines allow use samples with different delta in one batch\n            mask = (delta > i).reshape(-1,1,1,1)\n            x = x*(~mask).float() + y*mask.float()\n        \n        x = self.decoder(x)\n        \n        return x ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load data","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\nfrom sklearn.model_selection import train_test_split","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_val = pd.read_csv('/kaggle/input/conways-reverse-game-of-life-2020/train.csv', index_col='id')\ntest = pd.read_csv('/kaggle/input/conways-reverse-game-of-life-2020/test.csv', index_col='id')\n\ntrain, val = train_test_split(train_val, test_size=0.2, shuffle=True, random_state=42, stratify=train_val['delta'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Dataset","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from torch.utils.data import DataLoader, Dataset\nfrom torch import FloatTensor, LongTensor\n\n\ndef line2grid_tensor(data, device='cuda'):\n    grid = data.to_numpy().reshape((data.shape[0], 1, 25, 25))\n    return FloatTensor(grid).to(device)\n\n\nclass TaskDataset(Dataset):\n    def __init__(self, data, device='cuda'):\n        self.delta = LongTensor(data['delta'].to_numpy()).to(device)\n        if data.shape[1] == 1251: \n            self.start = line2grid_tensor(data.iloc[:,1:626], device)\n            self.stop = line2grid_tensor(data.iloc[:,626:], device)\n        else:\n            self.start = None\n            self.stop = line2grid_tensor(data.iloc[:,1:], device)\n        \n    def __len__(self):\n        return len(self.delta)\n\n    def __getitem__(self, idx):\n        if self.start is None:\n            return {'stop': self.stop[idx], 'delta': self.delta[idx]}\n        return {'start': self.start[idx], 'stop': self.stop[idx], 'delta': self.delta[idx]}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dataset_train = TaskDataset(train)\ndataloader_train = DataLoader(dataset_train, batch_size=128, shuffle=True)\n\ndataset_val = TaskDataset(val)\ndataloader_val = DataLoader(dataset_val, batch_size=128, shuffle=False)\n\ndataset_test = TaskDataset(test)\ndataloader_test = DataLoader(dataset_test, batch_size=128, shuffle=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train Loop","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from catalyst.dl import SupervisedRunner\nfrom catalyst.dl.callbacks import CriterionCallback, EarlyStoppingCallback, AccuracyCallback\nfrom catalyst.contrib.nn.optimizers import RAdam, Lookahead\n\nimport collections\n\nrunner = SupervisedRunner(device='cuda', input_key=['stop', 'delta'], )\n\nloaders = {'train': dataloader_train, 'valid': dataloader_val}#collections.OrderedDict({'train': dataloader_train, 'valid': dataloader_val})\n\nmodel = ReverseModel()\n\noptimizer = Lookahead(RAdam(params=model.parameters(), lr=1e-3))\n\ncriterion = {\"bce\": nn.BCEWithLogitsLoss()}\n\nscheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, factor=0.25, patience=2)\n\ncallbacks = [\n        CriterionCallback(input_key='start', prefix=\"loss\", criterion_key=\"bce\"),\n        EarlyStoppingCallback(patience=5),\n    ]\n\nrunner.train(\n    model=model,\n    criterion=criterion,\n    optimizer=optimizer,\n    scheduler=scheduler,\n    loaders=loaders,\n    callbacks=callbacks,\n    logdir=\"./logs\",\n    num_epochs=999,\n    main_metric=\"loss\",\n    minimize_metric=True,\n    verbose=True,\n)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Prediction / LB estimate","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Load best model","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"best_model = ReverseModel().to('cuda')\nbest_model.load_state_dict(torch.load('logs/checkpoints/best.pth')['model_state_dict'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### LB estimate","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"%load_ext Cython","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%cython\ncimport cython\n\nimport numpy as np\n\n@cython.cdivision(True)\n@cython.boundscheck(False)\n@cython.nonecheck(False)\n@cython.wraparound(False)\ncdef int calc_neighs(unsigned char[:, :] field, int i, int j, int n, int k):\n    cdef:\n        int neighs = 0;\n        int i_min = i - 1;\n        int i_pl = i + 1;\n        int j_min = j - 1;\n        int j_pl = j + 1;\n    neighs = 0\n    if i_min >= 0:\n        if j_min >= 0:\n            neighs += field[i_min, j_min]\n        neighs += field[i_min, j]\n        if j_pl < k:\n            neighs += field[i_min, j_pl]\n    if j_min >= 0:\n        neighs += field[i, j_min]\n    if j_pl < k:\n        neighs += field[i, j_pl]\n    if i_pl < n:\n        if j_min >= 0:\n            neighs += field[i_pl, j_min]\n        neighs += field[i_pl, j]\n        if j_pl < k:\n            neighs += field[i_pl, j_pl]\n    return neighs\n\n@cython.cdivision(True)\n@cython.boundscheck(False)\n@cython.nonecheck(False)\n@cython.wraparound(False)\ncpdef make_move_cython(unsigned char[:, :] field, int moves):\n    cdef:\n        int _, i, j, neighs;\n        int n, k;\n        int switch = 0;\n        unsigned char[:, :] cur_field;\n        unsigned char[:, :] next_field;\n    cur_field = np.copy(field)\n    next_field = np.zeros_like(field, 'uint8')\n    n = field.shape[0]\n    k = field.shape[1]\n    for _ in range(moves):\n        if switch == 0:\n            for i in range(n):\n                for j in range(k):\n                    neighs = calc_neighs(cur_field, i, j, n, k)\n                    if cur_field[i, j] and neighs == 2:\n                        next_field[i, j] = 1\n                    elif neighs == 3:\n                        next_field[i, j] = 1\n                    else:\n                        next_field[i, j] = 0\n        else:\n            for i in range(n):\n                for j in range(k):\n                    neighs = calc_neighs(next_field, i, j, n, k)\n                    if next_field[i, j] and neighs == 2:\n                        cur_field[i, j] = 1\n                    elif neighs == 3:\n                        cur_field[i, j] = 1\n                    else:\n                        cur_field[i, j] = 0\n        switch = (switch + 1) % 2\n    return np.array(next_field if switch else cur_field)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from typing import Union, List\n\nimport torch\nimport numpy as np\n\n\n\ndef make_move(board: Union[torch.Tensor, np.ndarray], moves: Union[int, List[int]] = 1):\n    \"\"\"\n    Advance the game by moves steps.\n    Handles:\n     - batch of shape (batch_size, 1, 25, 25)\n     - batch of shape (batch_size, 25, 25)\n     - single board of shape (1, 25, 25)\n     - single board of shape (25, 25)\n\n    Moves can be an int or a list of ints.\n\n    Returns the advanced game with the same shape as the input, with same dtype and device (if torch).\n    \"\"\"\n    if isinstance(board, torch.Tensor):\n        return torch.tensor(_evolve_board(board.detach().cpu().numpy(), moves),\n                            dtype=board.dtype, device=board.device)\n\n    return np.array(_evolve_board(board, moves), dtype=board.dtype)\n\n\ndef _evolve_board(board, moves):\n    # uint8 for make_move_cython\n    board = np.array(board).astype(np.uint8)\n\n    # whole batch with one channel\n    if len(board.shape) == 4:\n        assert board.shape[1:] == (1, 25, 25)\n        board_evolved = _move_board_3dim(board[:, 0], moves)[:, None]\n\n    # whole batch or single board with one channel\n    elif len(board.shape) == 3:\n        board_evolved = _move_board_3dim(board, moves)\n\n    # normal board\n    else:\n        assert board.shape == (25, 25)\n        board_evolved = make_move_cython(board, moves)\n\n    return np.array(board_evolved, dtype=np.float32)\n\n\ndef _move_board_3dim(board, moves):\n    if board.shape == (1, 25, 25):\n        board_evolved = make_move_cython(board[0], moves)[:, None]\n    else:\n        assert board.shape[1:] == (25, 25)\n        if isinstance(moves, int):\n            moves = [moves] * len(board)\n        assert len(moves) == len(board)\n        board_evolved = np.array([make_move_cython(b, move) for b, move in zip(board, moves)])\n    return board_evolved\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch\n\n\ndef batch_accuracy(x, y_pred):\n    return torch.sum(torch.abs(x - y_pred), dim=(1, 2, 3))\n\n\ndef postprocessor(y_pred, x, moves):\n    \"\"\"\n    Find best possible threshold for each prediction.\n    Roughly equivalent to:\n\n    X = x\n    Y = y\n    thresholds = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]\n    best_thresholds = []\n    for y, x, move in zip(Y, X, moves):\n        best_accuracy = 0\n        best_threshold = 0\n        for threshold in thresholds:\n            current_accuracy = accuracy(make_move(y > threshold, move), x)\n            if current_accuracy > best_accuracy:\n                best_accuracy = current_accuracy\n                best_thresholds = threshold\n        best_thresholds.append(best_threshold)\n\n    y_processed = [y > threshold for y, threshold in zip(Y, best_thresholds)]\n\n    :param moves: int or list of int, specifying the number of steps in game of life.\n    :param y_pred: sigmoid predictions of game t in [0, 1], shape (batch_size, 1, 25, 25).\n    :param x: input game (at time t + moves), shape (batch_size, 1, 25, 25)\n    :return: Binarized predictions, optimized for minimal error.\n    \"\"\"\n    thresholds = torch.linspace(0, 1, 11, device=x.device)\n    errors = torch.cat([batch_accuracy(x.float(),\n                                       make_move((y_pred > threshold).float(), moves)\n                                       )[None]\n                        for threshold in thresholds],\n                       0)\n    # errors.shape = (len(threhsolds), len(x))\n    best_thresholds = thresholds[torch.argmin(errors, 0)][:, None, None, None]\n    return (y_pred > best_thresholds).float()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\n\ndef predict_batch(model, batch):\n    model.eval()\n    with torch.no_grad():\n        prediction = model(batch['stop'], batch['delta'])\n        prediction = torch.sigmoid(prediction).detach().cpu().numpy()\n        return prediction\n\n    \ndef predict_loader(model, loader):\n    predict = [predict_batch(model, batch) for batch in loader]\n    predict = np.concatenate(predict)\n    return predict\n\n\ndef validate_loader(model, loader, lb_delta=None, threshold=0.5):\n    prediction_val = predict_loader(best_model, loader)\n    y_val = loader.dataset.start.detach().cpu().numpy()\n    \n\n    score_unoptimized = ((prediction_val > threshold) == y_val).mean(axis=(1,2,3))\n    \n    delta_val = loader.dataset.delta.detach().cpu().numpy()\n    prediction_val = postprocessor(torch.tensor(prediction_val),\n                                   x=loader.dataset.stop.cpu(),\n                                   moves=delta_val).cpu().numpy()\n    score = (prediction_val == y_val).mean(axis=(1,2,3))\n    \n    print(f'All data accuracy (global threshold): {score_unoptimized.mean()}')\n    print(f'All data accuracy (optimized threshold): {score.mean()}')\n        \n    delta_score = {}\n    for i in range(1, 6):\n        delta_score[i] = score[delta_val==i].mean()#print(f'delta={i} accuracy: {score[delta_val==i].mean()}')\n        print(f'delta={i} accuracy: {delta_score[i]}')\n        \n    if lb_delta is not None:\n        lb_delta = lb_delta.value_counts(normalize=True)\n        test_score = sum([lb_delta[i]*delta_score[i] for i in range(1,6)])\n        print(f'VAL score         : {1-score.mean()}')\n        print(f'LB  score estimate: {1-test_score}')\n    \n    \ndef make_submission(prediction, x, moves, sample_submission_path='/kaggle/input/conways-reverse-game-of-life-2020/sample_submission.csv'):\n    prediction = postprocessor(prediction, x, moves).numpy().astype(int).reshape(-1, 625)\n\n    sample_submission = pd.read_csv(sample_submission_path, index_col='id')\n    sample_submission.iloc[:] = prediction\n    return sample_submission","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"validate_loader(best_model, dataloader_val, test['delta'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### submission","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"prediction_test = predict_loader(best_model, dataloader_test)\nsubmission = make_submission(torch.tensor(prediction_test),\n                             x=dataloader_test.dataset.stop.cpu(),\n                             moves=dataloader_test.dataset.delta.detach().cpu().numpy())\nsubmission.to_csv('submission.csv')\nsubmission.head()","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}