{"cells":[{"metadata":{},"cell_type":"markdown","source":"# OuroborosLife - Function Reversal GAN\n\nOuroboros is the ancient symbol of a serpent eating its own tail.\n\n![Ouroboros](https://i.imgur.com/IGOdGbH.jpeg)\n\nThis a novel Ouroboros network architecture inspired by GANs\n\nThe idea is that we have a three headed output, that can predict the Past, Present and Future states of a function sequence.\n\nThe Present function mapping is simply the identity function, which is easy to solve for\n\nThe Future function mapping is a simple boolean function of the 1-distance neighbours for each cell ( Future == Alive + 2-3 neighbours || Dead + 3 neighbours ). There is a many-to-one mapping between inputs and outputs, which can be captured by a static dataset. I was able to solve this function to 100% accuracy with a much simpler neural network:\n- https://www.kaggle.com/jamesmcguigan/game-of-life-forward-in-pytorch-with-100-accuracy \n\nThe Past function mapping is more complicated. There is a one-to-many domain between inputs and outputs. For a Present alive cell, there are (8 choose 2) + (8 choose 3) = 84 perfectly valid Past states. To use a statically defined dataset here, would result in the network being shown only one of these 84 states at random, and recieving a loss error it had predicted a different by still valid past state.\n\nThe novelty here is the Ouroboros loss function. Given that the network will simultaniously output Past, Present and Future function states, it is possible to feed the output for any one of these heads back into the input and generate a second order loss function. Thus we feed the Present into the network to generate a predicted Past, then feed the Past back into the network to predict it's Future. We compare the Past's Future with the Present, which should be the same as the dataset Input. The second-order ouroboros loss is simply the Input compared to itself, using any of the standard MSE / BCE / Focal loss functions.\n\nSpecifically, the training loop generates a sequence dataset of Game of Life timelines. Each timestep in the sequence is used as an individual input. The first order loss is computed by comparing the Ouroboros output to the expected slice from the timeline dataset. Each Past/Present/Future output is then fed back into the Ouroboros network, which is then again compated to the approprate time adjusted slice of the timeline dataset.\n\nThe network thus learns to accept both an external dataset as input, as well as training its own output to be a valid input."},{"metadata":{},"cell_type":"markdown","source":"# Classical Function Implemention\n\n\nUsing the classic ruleset on a 25x25 board with wraparound, the game evolves at each timestep according to the following rules\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."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Functions for implementing Game of Life Forward Play\nfrom typing import List\n\nimport numpy as np\nimport scipy.sparse\nfrom joblib import delayed\nfrom joblib import Parallel\nfrom numba import njit\n\n\n# Source: https://www.kaggle.com/ianmoone0617/reversing-conways-game-of-life-tutorial\ndef life_step_1(X: np.ndarray):\n    \"\"\"Game of life step using generator expressions\"\"\"\n    nbrs_count = sum(np.roll(np.roll(X, i, 0), j, 1)\n                     for i in (-1, 0, 1) for j in (-1, 0, 1)\n                     if (i != 0 or j != 0))\n    return (nbrs_count == 3) | (X & (nbrs_count == 2))\n\n\n# Source: https://www.kaggle.com/ianmoone0617/reversing-conways-game-of-life-tutorial\ndef life_step_2(X: np.ndarray):\n    \"\"\"Game of life step using scipy tools\"\"\"\n    from scipy.signal import convolve2d\n    nbrs_count = convolve2d(X, np.ones((3, 3)), mode='same', boundary='wrap') - X\n    return (nbrs_count == 3) | (X & (nbrs_count == 2))\n\n\n\n# NOTE: @njit doesn't like np.roll(axis=) so reimplement explictly\n@njit\ndef life_neighbours_xy(board: np.ndarray, x, y, max_value=3):\n    size_x = board.shape[0]\n    size_y = board.shape[1]\n    neighbours = 0\n    for i in (-1, 0, 1):\n        for j in (-1, 0, 1):\n            if i == j == 0: continue    # ignore self\n            xi = (x + i) % size_x\n            yj = (y + j) % size_y\n            neighbours += board[xi, yj]\n            if neighbours > max_value:  # shortcircuit return 4 if overpopulated\n                return neighbours\n    return neighbours\n\n\n@njit\ndef life_neighbours(board: np.ndarray, max_value=3):\n    size_x = board.shape[0]\n    size_y = board.shape[1]\n    output = np.zeros(board.shape, dtype=np.int8)\n    for x in range(size_x):\n        for y in range(size_y):\n            output[x,y] = life_neighbours_xy(board, x, y, max_value)\n    return output\n\n\n@njit\ndef life_step(board: np.ndarray) -> np.ndarray:\n    \"\"\"Game of life step using generator expressions\"\"\"\n    size_x = board.shape[0]\n    size_y = board.shape[1]\n    output = np.zeros(board.shape, dtype=np.int8)\n    for x in range(size_x):\n        for y in range(size_y):\n            cell       = board[x,y]\n            neighbours = life_neighbours_xy(board, x, y, max_value=3)\n            if ( (cell == 0 and      neighbours == 3 )\n              or (cell == 1 and 2 <= neighbours <= 3 )\n            ):\n                output[x, y] = 1\n    return output\n\ndef life_steps(boards: List[np.ndarray]) -> List[np.ndarray]:\n    \"\"\" Parallel version of life_step() but for an array of boards \"\"\"\n    return Parallel(-1)( delayed(life_step)(board) for board in boards )\n\n\n@njit\ndef life_step_delta(board: np.ndarray, delta):\n    for t in range(delta): board = life_step(board)\n    return board\n\n\ndef life_step_3d(board: np.ndarray, delta):\n    solution_3d = np.array([ board ], dtype=np.int8)\n    for t in range(delta):\n        board       = life_step(board)\n        solution_3d = np.append( solution_3d, [ board ], axis=0)\n    return solution_3d\n\n\n# RULES: https://www.kaggle.com/c/conway-s-reverse-game-of-life/data\ndef generate_random_board(shape=(25,25)):\n    # An initial board was chosen by filling the board with a random density between 1% full (mostly zeros) and 99% full (mostly ones).\n    # DOCS: https://cmdlinetips.com/2019/02/how-to-create-random-sparse-matrix-of-specific-density/\n    density = np.random.random() * 0.98 + 0.01\n    board   = scipy.sparse.random(*shape, density=density, data_rvs=np.ones).toarray().astype(np.int8)\n\n    # The starting board's state was recorded after the 5 \"warmup steps\". These are the values in the start variables.\n    for t in range(5):\n        board = life_step(board)\n        if np.count_nonzero(board) == 0:\n            return generate_random_board(shape)  # exclude empty boards and try again\n    return board\n\ndef generate_random_boards(count, shape=(25,25)):\n    generated_boards = Parallel(-1)( delayed(generate_random_board)(shape) for _ in range(count) )\n    return generated_boards\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Modified from Source: https://github.com/c0nn3r/RetinaNet/blob/master/focal_loss.py\n# Switched from using F.cross_entropy() to F.binary_cross_entropy()\nimport torch\n\nimport torch.nn as nn\n\n\nclass FocalLoss(nn.Module):\n\n    def __init__(self, focusing_param=2, balance_param=0.25):\n        super(FocalLoss, self).__init__()\n\n        self.focusing_param = focusing_param\n        self.balance_param  = balance_param\n        self.bce            = nn.BCELoss()\n\n    def forward(self, output, target):\n        logpt      = - self.bce(output, target)\n        pt         = torch.exp(logpt)\n        focal_loss = -((1 - pt) ** self.focusing_param) * logpt\n        balanced_focal_loss = self.balance_param * focal_loss\n        return balanced_focal_loss","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Pytorch Base Class\n\nThere is a little bit of core infrasture code that needs to be written, to handle common functionality such as:\n- model save/autoload\n- casting between data formats\n- freezing and unfreezing\n- training loop functions\n\nBy putting this all in a base class, we can seperate out infrasture code from application logic, and makes code reuse easier with the ability to subclass these functions for different usecases."},{"metadata":{"trusted":true},"cell_type":"code","source":"device   = torch.device(\"cuda:0\") if torch.cuda.is_available() else torch.device(\"cpu\")\n__file__ = './notebook.py'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from __future__ import annotations\n\nimport os\nfrom abc import ABCMeta\nfrom typing import List\nfrom typing import TypeVar\nfrom typing import Union\n\nimport humanize\nimport numpy as np\nimport torch\nimport torch.nn as nn\n\n# from neural_networks.device import device\n\n# noinspection PyTypeChecker\nT = TypeVar('T', bound='GameOfLifeBase')\nclass GameOfLifeBase(nn.Module, metaclass=ABCMeta):\n    \"\"\"\n    Base class for GameOfLife based NNs\n    Handles: save/autoload, freeze/unfreeze, casting between data formats, and training loop functions\n    \"\"\"\n    def __init__(self):\n        super().__init__()\n        self.loaded    = False  # can't call sell.load() in constructor, as weights/layers have not been defined yet\n        self.device    = torch.device(\"cuda:0\") if torch.cuda.is_available() else torch.device(\"cpu\")\n        self.criterion = nn.MSELoss()\n\n\n    @staticmethod\n    def weights_init(layer):\n        if isinstance(layer, (nn.Conv2d, nn.ConvTranspose2d)):\n            nn.init.kaiming_normal_(layer.weight)\n            nn.init.constant_(layer.bias, 0.1)\n\n    ### Prediction\n\n    def __call__(self, *args, **kwargs) -> torch.Tensor:\n        if not self.loaded: self.load()  # autoload on first function call\n        return super().__call__(*args, **kwargs)\n\n    def predict(self, inputs: Union[List[np.ndarray], np.ndarray, torch.Tensor]) -> np.ndarray:\n        \"\"\" Wrapper function around __call__() that returns a numpy int8 array for external usage \"\"\"\n        outputs = self(inputs)\n        outputs = self.cast_int(outputs).squeeze().cpu().numpy()\n        return outputs\n\n\n\n    ### Training\n\n    def loss(self, outputs, expected, input):\n        return self.criterion(outputs, expected)\n\n    def accuracy(self, outputs, expected, inputs) -> float:\n        # noinspection PyTypeChecker\n        return torch.sum(self.cast_int(outputs) == self.cast_int(expected)).cpu().numpy() / np.prod(outputs.shape)\n\n\n\n    ### Freee / Unfreeze\n\n    def freeze(self: T) -> T:\n        if not self.loaded: self.load()\n        for name, parameter in self.named_parameters():\n            parameter.requires_grad = False\n        return self\n\n    def unfreeze(self: T) -> T:\n        if not self.loaded: self.load()\n        for name, parameter in self.named_parameters():\n            parameter.requires_grad = True\n        return self\n\n\n\n    ### Load / Save Functionality\n\n    @property\n    def filename(self) -> str:\n        return os.path.join( os.path.dirname(__file__), 'models', f'{self.__class__.__name__}.pth')\n\n\n    # DOCS: https://pytorch.org/tutorials/beginner/saving_loading_models.html\n    def save(self: T) -> T:\n        os.makedirs(os.path.dirname(self.filename), exist_ok=True)\n        torch.save(self.state_dict(), self.filename)\n        print(f'{self.__class__.__name__}.savefile(): {self.filename} = {humanize.naturalsize(os.path.getsize(self.filename))}')\n        return self\n\n\n    def load(self: T) -> T:\n        if os.path.exists(self.filename):\n            try:\n                self.load_state_dict(torch.load(self.filename))\n                print(f'{self.__class__.__name__}.load(): {self.filename} = {humanize.naturalsize(os.path.getsize(self.filename))}')\n            except Exception as exception:\n                # Ignore errors caused by model size mismatch\n                print(f'{self.__class__.__name__}.load(): model has changed dimensions, reinitializing weights\\n')\n                self.apply(self.weights_init)\n        else:\n            print(f'{self.__class__.__name__}.load(): model file not found, reinitializing weights\\n')\n            self.apply(self.weights_init)\n\n        self.loaded = True    # prevent any infinite if self.loaded loops\n        self.to(self.device)  # ensure all weights, either loaded or untrained are moved to GPU\n        self.eval()           # default to production mode - disable dropout\n        self.freeze()         # default to production mode - disable training\n        return self\n\n\n\n    ### Casting\n\n    def cast_bool(self, x: torch.Tensor) -> torch.Tensor:\n        # noinspection PyTypeChecker\n        return (x > 0.5)\n\n    def cast_int(self, x: torch.Tensor) -> torch.Tensor:\n        return self.cast_bool(x).to(torch.int8)\n\n    def cast_int_float(self, x: torch.Tensor) -> torch.Tensor:\n        return self.cast_bool(x).to(torch.float32).requires_grad_(True)\n\n\n    def cast_to_tensor(self, x: Union[np.ndarray, torch.Tensor]) -> torch.Tensor:\n        if torch.is_tensor(x):\n            return x.to(torch.float32).to(device)\n        if isinstance(x, list):\n            x = np.array(x)\n        if isinstance(x, np.ndarray):\n            x = torch.from_numpy(x).to(torch.float32)\n            x = x.to(device)\n            return x  # x.shape = (42,3)\n        raise TypeError(f'{self.__class__.__name__}.cast_to_tensor() invalid type(x) = {type(x)}')\n\n\n    # DOCS: https://towardsdatascience.com/understanding-input-and-output-shapes-in-convolution-network-keras-f143923d56ca\n    # pytorch requires:    contiguous_format = (batch_size, channels, height, width)\n    # tensorflow requires: channels_last     = (batch_size, height, width, channels)\n    def cast_inputs(self, x: Union[List[np.ndarray], np.ndarray, torch.Tensor]) -> torch.Tensor:\n        x = self.cast_to_tensor(x)\n        if x.dim() == 1:             # single row from dataframe\n            x = x.view(1, 1, torch.sqrt(x.shape[0]), torch.sqrt(x.shape[0]))\n        elif x.dim() == 2:\n            if x.shape[0] == x.shape[1]:  # single 2d board\n                x = x.view(1, 1, x.shape[0], x.shape[1])\n            else: # rows of flattened boards\n                x = x.view(-1, 1, torch.sqrt(x.shape[1]), torch.sqrt(x.shape[1]))\n        elif x.dim() == 3:                                 # numpy  == (batch_size, height, width)\n            x = x.view(x.shape[0], 1, x.shape[1], x.shape[2])   # x.shape = (batch_size, channels, height, width)\n        elif x.dim() == 4:\n            pass  # already in (batch_size, channels, height, width) format, so do nothing\n        return x","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# OuroborosLife\n\nThe idea of the Ouroboros Network is that rather than just predicting the next or previous state,\nwe want to past, present and future simultaneously in the same network.\n\nThe dataset is a sequence of 3 consecutive board states generated by life_step().\n\nThe network takes the middle/present board state and attempts to predict all Past, Present and Future states\n\nThe loss function computes the loss against the original training data, but also feeds back in upon itself.\nThe output for Future is fed back in and it's Past is compared with the Present, likewise in reverse with the Past.\n\n---\n\nChangelog:\n- v1 - 3x3 convolution on base layer, was unable to train past 70% accuracy, even after 8 hours\n- v2 - try replacing top layer with 5x5 convolution and increase out_channels in deconvolution layer\n- v3 - refactor loss function to only use the dataset to predict forward play | Past Error is only via second-order ouroboros loss: Past's Future == Present  "},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import atexit\nimport gc\nimport time\nfrom typing import List\nfrom typing import Tuple\nfrom typing import Union\n\nimport numpy as np\nimport torch\nimport torch as pt\nimport torch.nn as nn\n\n# from neural_networks.FocalLoss import FocalLoss\n# from neural_networks.GameOfLifeBase import GameOfLifeBase\n# from utils.game import generate_random_board\n# from utils.game import life_step_3d\n\n\nclass OuroborosLife(GameOfLifeBase):\n    \"\"\"\n    The idea of the Ouroboros Network is that rather than just predicting the next or previous state,\n    we want to past, present and future simultaneously in the same network.\n\n    The dataset is a sequence of 3 consecutive board states generated by life_step().\n\n    The network takes the middle/present board state and attempts to predict all Past, Present and Future states\n\n    The loss function computes the loss against the original training data, but also feeds back in upon itself.\n    The output for Future is fed back in and it's Past is compared with the Present, likewise in reverse with the Paspt.\n    \"\"\"\n    @property\n    def filename(self) -> str:\n        \"\"\" ./models/OuroborosLife3.pth || ./models/OuroborosLife5.pth \"\"\"\n        return super().filename.replace('.pth', f'{self.out_channels}.pth')\n\n\n    def __init__(self, in_channels=1, out_channels=3):\n        \"\"\"\n        TODO:\n        - Create split blocks that return: [identity, avgpool, maxpool ]  # do we need 3x3 convolution with fixed weights\n        - Basically find a way to count neighbours\n        - reduce model size\n        - add in dense layer at end + middle (as opposed to deconvolution???)\n        \"\"\"\n        assert out_channels % 2 == 1, f'{self.__class__.__name__}(out_channels={out_channels}) must be odd'\n\n        super().__init__()\n        self.in_channels  = in_channels\n        self.out_channels = out_channels  # Past, Present and Future\n\n        self.relu    = nn.LeakyReLU()     # combines with nn.init.kaiming_normal_()\n        self.dropout = nn.Dropout(p=0.2)\n\n        # 2**9 = 512 filters and kernel size of 3x3 to allow for full encoding of game rules\n        # Pixels can see distance 5 neighbours, (hopefully) sufficient for delta=2 timesteps or out_channels=5\n        # https://www.youtube.com/watch?v=H3g26EVADgY&feature=youtu.be&t=1h39m410s&ab_channel=JeremyHoward\n        self.cnn_layers = nn.ModuleList([\n            # Previous pixel state requires information from distance 2, so we need two 3x3 convolutions\n            nn.Conv2d(in_channels=in_channels, out_channels=512,  kernel_size=(5,5), padding=2, padding_mode='circular'),\n            nn.Conv2d(in_channels=512,   out_channels=256,  kernel_size=(1,1)),\n            nn.Conv2d(in_channels=256,   out_channels=128,  kernel_size=(1,1)),\n\n            nn.Conv2d(in_channels=1+128, out_channels=128,  kernel_size=(3,3), padding=1, padding_mode='circular'),\n            nn.Conv2d(in_channels=128,   out_channels=512,  kernel_size=(1,1)),\n            nn.Conv2d(in_channels=512,   out_channels=256,  kernel_size=(1,1)),\n            nn.Conv2d(in_channels=256,   out_channels=128,  kernel_size=(1,1)),\n\n            # # Deconvolution + Convolution allows neighbouring pixels to share information to simulate forward play\n            # # This creates a 52x52 grid of interspersed cells that can then be downsampled back down to 25x25\n            nn.ConvTranspose2d(in_channels=1+128, out_channels=512,  kernel_size=(3,3), stride=2, dilation=1),\n            nn.Conv2d(in_channels=512,   out_channels=256,   kernel_size=(1,1)),\n            nn.Conv2d(in_channels=256,   out_channels=64,    kernel_size=(1,1)),\n            nn.Conv2d(in_channels=64,    out_channels=128,   kernel_size=(3,3), stride=2),  # undo deconvolution\n\n            nn.Conv2d(in_channels=1+128, out_channels=64,    kernel_size=(1,1)),\n            nn.Conv2d(in_channels=64,    out_channels=32,    kernel_size=(1,1)),\n            nn.Conv2d(in_channels=32,    out_channels=16,    kernel_size=(1,1)),\n            nn.Conv2d(in_channels=1+16,  out_channels=out_channels, kernel_size=(1,1)),\n        ])\n        self.batchnorm_layers = nn.ModuleList([\n            nn.BatchNorm2d(cnn_layer.out_channels)\n            for cnn_layer in self.cnn_layers\n        ])\n\n\n        # self.criterion = nn.BCELoss()\n        self.criterion = FocalLoss()\n        # self.criterion = nn.MSELoss()\n        self.optimizer = pt.optim.RMSprop(self.parameters(), lr=0.01, momentum=0.9)\n        self.scheduler = torch.optim.lr_scheduler.CyclicLR(\n            self.optimizer,\n            max_lr=1e-3,\n            base_lr=1e-5,\n            step_size_up=10,\n            mode='exp_range',\n            gamma=0.8\n        )\n\n    # def load(self):\n    #     super().load()\n    #     self.apply(self.weights_init)\n\n\n    def forward(self, x):\n        x = input = self.cast_inputs(x)\n        for n, (cnn_layer, batchnorm_layer) in enumerate(zip(self.cnn_layers, self.batchnorm_layers)):\n            if cnn_layer.in_channels > 1 and cnn_layer.in_channels % 2 == 1:   # autodetect 1+in_channels == odd number\n                x = torch.cat([ x, input ], dim=1)                     # passthrough original cell state\n            x = cnn_layer(x)\n            if n != len(self.cnn_layers)-1:\n                x = self.relu(x)\n                if n != 1:               # Don't apply dropout to the first layer\n                    x = self.dropout(x)  # BatchNorm eliminates the need for Dropout in some cases cause BN provides similar regularization benefits as Dropout intuitively\"\n                x = batchnorm_layer(x)   # batchnorm goes after activation\n            else:\n                x = torch.sigmoid(x)  # output requires sigmoid activation\n        return x\n\n\n    # DOCS: https://towardsdatascience.com/understanding-input-and-output-shapes-in-convolution-network-keras-f143923d56ca\n    # pytorch requires:    contiguous_format = (batch_size, channels, height, width)\n    # tensorflow requires: channels_last     = (batch_size, height, width, channels)\n    def cast_inputs(self, x: Union[List[np.ndarray], np.ndarray, torch.Tensor]) -> torch.Tensor:\n        x = self.cast_to_tensor(x)\n        if   x.dim() == 4: pass\n        elif x.dim() == 3 and x.shape[0] == self.out_channels:           # x.shape = (channels, height, width)\n            x = x.view(1, self.in_channels, x.shape[1], x.shape[2])   # x.shape = (batch_size, channels, height, width)\n        else:\n            x = super().cast_inputs(x)\n        return x\n\n\n    def loss_dataset(self, outputs, timeline, inputs, exclude_past=True):\n        # Exclude past losses to avoid the many-pasts to one-future problem\n        if exclude_past:\n            t_present = self.out_channels//2\n            outputs   = outputs[  :, t_present:, :, : ]\n            timeline  = timeline[ :, t_present:, :, : ]\n\n        ### Other ways of computing dataset loss\n        # dataset_loss = torch.mean(torch.mean(( (timeline-outputs)**2 ).flatten(1), dim=1)) # average MSE per timeframe\n        # dataset_loss = torch.sum(torch.tensor([\n        #     self.criterion(outputs[b][t], timeline[b][t])  # NOTE: FocalLoss(outputs, target) needed in correct order\n        #     for b in range(timeline.shape[0])\n        #     for t in range(timeline.shape[1])\n        # ], requires_grad=True))\n\n        dataset_loss = self.criterion(outputs, timeline)\n        return dataset_loss\n\n\n    def loss_accuracy_ouroboros(self, outputs, timeline, inputs) -> Tuple[pt.Tensor, pt.Tensor, pt.Tensor]:\n        \"\"\"\n        Compute simplified losses for each head, only comparing reoutputs with timeline[t_present]\n        reinput    = t1=0 = Past | t1=1 = Present | t1=2 = Future\n        reoutput   = [ Past2, Past, Present ]@t=0, [ Past, Present, Future ]@t=1, [ Present, Future, Future2 ]@t=2\n        \"\"\"\n        losses     = pt.zeros(self.out_channels, dtype=pt.float32, requires_grad=True).to(self.device)\n        acc_boards = pt.zeros(self.out_channels, dtype=pt.float32, requires_grad=False).to(self.device)\n        acc_pixels = pt.zeros(self.out_channels, dtype=pt.float32, requires_grad=False).to(self.device)\n        for t_input in range(self.out_channels):\n            t_present = self.out_channels//2 - t_input\n            reinput   = outputs[:,t_input,:,:].unsqueeze(1)\n            reoutputs = self(reinput)\n\n            # Losses get calculated for present and all future datapoints\n            for d in range(self.out_channels):\n                if reoutputs.shape[1] <= t_present + d: break\n                if timeline.shape[1]  <= t_input   + d: break\n                losses[t_input] += self.criterion(reoutputs[:,t_present+d,:,:], timeline[:,t_input+d,:,:])\n\n            # Accuracy is based only on the self-referential present\n            pixels_correct       = ((reoutputs[:,t_present,:,:] > 0.5) == (timeline[:,t_input,:,:] > 0.5)).to(pt.float).detach()\n            acc_pixels[t_input] +=  pixels_correct.mean()\n            acc_boards[t_input] += (pixels_correct.mean(dim=1) == 1.0).to(pt.float).mean()\n        return losses, acc_pixels, acc_boards\n\n\n\n\n    def fit(self, epochs=100_000, batch_size=25, max_delta=25, timeout=0):\n        gc.collect()\n        torch.cuda.empty_cache()\n        atexit.register(model.save)\n        self.train()\n        self.unfreeze()\n        print(self)\n        try:\n            # timelines_batch = np.array([\n            #     life_step_3d(generate_random_board(), max_delta)\n            #     for _ in range(batch_size)\n            # ])\n            time_start  = time.perf_counter()\n            board_count = 0\n            dataset_accuracies = [0]\n            for epoch in range(1, epochs+1):\n                if np.min(dataset_accuracies[-10:]) == 1.0: break  # we have reached 100% accuracy\n                if timeout and timeout < time.perf_counter() - time_start: break\n\n                epoch_start = time.perf_counter()\n                timelines_batch = np.array([\n                    life_step_3d(generate_random_board(), max_delta)\n                    for _ in range(batch_size)\n                ])\n                epoch_ds_losses  = []\n                epoch_losses     = []\n                epoch_acc_pixels = []\n                epoch_acc_boards = []\n                d = self.out_channels // 2  # In theory this should work for 5 or 7 channels\n                for t in range(d, max_delta - d):\n                    inputs_np   = timelines_batch[:, np.newaxis, t,:,:]  # (batch_size=10, channels=1,  width=25, height=25)\n                    timeline_np = timelines_batch[:, t-d:t+d+1,    :,:]  # (batch_size=10, channels=10, width=25, height=25)\n                    inputs      = pt.tensor(inputs_np).to(self.device).to(pt.float32)\n                    timeline    = pt.tensor(timeline_np).to(self.device).to(pt.float32)\n\n                    self.optimizer.zero_grad()\n                    outputs      = self(inputs)\n                    dataset_loss = self.loss_dataset(outputs, timeline, inputs, exclude_past=True)\n                    orb_losses, acc_pixels, acc_boards = model.loss_accuracy_ouroboros(outputs, timeline, inputs)\n                    loss = pt.mean(orb_losses) + dataset_loss\n                    loss.backward()\n                    self.optimizer.step()\n                    self.scheduler.step()\n\n                    board_count += batch_size\n                    epoch_ds_losses.append(dataset_loss.detach().item())\n                    epoch_losses.append(orb_losses.detach())\n                    epoch_acc_pixels.append(acc_pixels.detach())\n                    epoch_acc_boards.append(acc_boards.detach())\n                    torch.cuda.empty_cache()\n\n                dataset_accuracies.append( pt.stack(epoch_acc_boards).min() )\n                epoch_loss      = f\"{100*np.mean(epoch_ds_losses):.6f} : \" + \" \".join([ f'{100*n:.6f}' for n in pt.stack(epoch_losses).mean(dim=0).tolist()     ])\n                epoch_acc_pixel = \" \".join([ f'{n:.3f}'     for n in pt.stack(epoch_acc_pixels).mean(dim=0).tolist() ])\n                epoch_acc_board = \" \".join([ f'{n:.3f}'     for n in pt.stack(epoch_acc_boards).mean(dim=0).tolist() ])\n\n                epoch_time = time.perf_counter() - epoch_start\n                time_taken = time.perf_counter() - time_start\n                if epoch <= 10 or epoch <= 100 and epoch % 10 == 0 or epoch % 100 == 0:  \n                    print(f'epoch: {epoch:4d} | boards: {board_count:5d} | loss: {epoch_loss} | pixels = {epoch_acc_pixel} | boards = {epoch_acc_board} | time: {time_taken//60:.0f}:{time_taken%60:02.0f} @ {1000*epoch_time//batch_size:3.0f}ms')\n                    # print(f'epoch: {epoch:4d} | boards: {board_count:5d} | loss: {np.mean(epoch_losses):.6f} | ouroboros: {np.mean(ouroboros_losses):.6f} | dataset: {np.mean(dataset_losses):.6f} | accuracy = {np.mean(epoch_accuracies):.6f} | time: {1000*epoch_time//batch_size}ms/board | {time_taken//60:.0f}:{time_taken%60:02.0f}')\n        except KeyboardInterrupt: pass\n        finally:\n            model.save()\n            atexit.unregister(model.save)\n            torch.cuda.empty_cache()\n            gc.collect()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training on Generated Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"!cp -rv ../input/ouroboroslife-function-reversal-gan/models ./\n!cp -v  ../input/ouroboroslife-function-reversal-gan/*.csv  ./","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = OuroborosLife()\nmodel.fit(timeout=1.5*60*60)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Utility Functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"from typing import Dict\n\nimport numpy as np\nimport pandas as pd\nfrom fastcache import clru_cache\n\n\n@clru_cache(None)\ndef csv_column_names(key='start'):\n    return [ f'{key}_{n}' for n in range(25**2) ]\n\n\ndef csv_to_delta(df, idx):\n    return int(df.loc[idx]['delta'])\n\ndef csv_to_delta_list(df):\n    return df['delta'].values\n\n\ndef csv_to_numpy(df, idx, key='start') -> np.ndarray:\n    try:\n        columns = csv_column_names(key)\n        board   = df.loc[idx][columns].values\n    except:\n        board = np.zeros((25, 25))\n    board = board.reshape((25,25)).astype(np.int8)\n    return board\n\n\ndef csv_to_numpy_list(df, key='start') -> np.ndarray:\n    try:\n        columns = csv_column_names(key)\n        output  = df[columns].values.reshape(-1,25,25)\n    except:\n        output  = np.zeros((0,25,25))\n    return output\n\n\n# noinspection PyTypeChecker,PyUnresolvedReferences\ndef numpy_to_dict(board: np.ndarray, key='start') -> Dict:\n    assert len(board.shape) == 2  # we want 2D solutions_3d[0] not 3D solutions_3d\n    assert key in { 'start', 'stop' }\n\n    board  = np.array(board).flatten().tolist()\n    output = { f\"{key}_{n}\": board[n] for n in range(len(board))}\n    return output\n\n\ndef numpy_to_series(board: np.ndarray, key='start') -> pd.Series:\n    solution_dict = numpy_to_dict(board, key)\n    return pd.Series(solution_dict)\n\n\n# Source: https://stackoverflow.com/questions/8290397/how-to-split-an-iterable-in-constant-size-chunks\ndef batch(iterable, n=1):\n    l = len(iterable)\n    for ndx in range(0, l, n):\n        yield iterable[ndx:min(ndx + n, l)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from numba import njit\n\n@njit\ndef is_valid_solution(start: np.ndarray, stop: np.ndarray, delta: int) -> bool:\n    # we are rewriting data, so lets double check our work\n    test_board = start\n    is_valid   = np.count_nonzero(test_board) != 0\n    if not is_valid: return False\n    for t in range(delta):\n        test_board = life_step(test_board)\n        is_valid   = is_valid and np.count_nonzero(test_board) != 0\n        if not is_valid: return False\n    is_valid = is_valid and np.all(test_board == stop)\n    return is_valid\n\n@njit\ndef is_valid_solution_3d(solution_3d: np.ndarray) -> bool:\n    return is_valid_solution(solution_3d[0], solution_3d[-1], delta=len(solution_3d)-1)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\n\ninput_directory  = '../input/conways-reverse-game-of-life-2020/'\noutput_directory = './'\n\ntrain_file             = f'{input_directory}/train.csv'\ntest_file              = f'{input_directory}/test.csv'\nsample_submission_file = f'{input_directory}/sample_submission.csv'\nsubmission_file        = f'{output_directory}/submission.csv'\nexact_submission_file  = f'{output_directory}/submission_exact.csv'\ntimeout_file           = f'{output_directory}/timeouts.csv'\n\ntrain_df              = pd.read_csv(train_file,              index_col='id')\ntest_df               = pd.read_csv(test_file,               index_col='id')\nsample_submission_df  = pd.read_csv(sample_submission_file,  index_col='id')\n\n# try:\n#     submission_df     = pd.read_csv(submission_file,         index_col='id')\n# except:\n#     submission_df     = sample_submission_df.copy()\n\n# try:\n#     exact_submission  = pd.read_csv(exact_submission_file,   index_col='id')\n# except:\n#     exact_submission  = sample_submission_df.copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# import time\n\n# from constraint_satisfaction.fix_submission import is_valid_solution\n# from neural_networks.OuroborosLife import OuroborosLife\n# from utils.datasets import *\n# from utils.util import batch\n# from utils.util import csv_to_delta_list\n# from utils.util import csv_to_numpy_list\n# from utils.util import numpy_to_series\n# import numpy as np\n\ndef ouroborors_dataframe(df: pd.DataFrame):\n    time_start    = time.perf_counter()\n    \n    model         = OuroborosLife()\n    model.load().train().unfreeze()\n    submission_df     = sample_submission_df.copy()\n    exact_submission  = sample_submission_df.copy()\n    \n    stats = {\n        \"boards\":  { \"solved\": 0, \"total\": 0 },\n        \"delta\":   { \"solved\": 0, \"total\": 0 },\n        \"dpixels\": { \"solved\": 0, \"total\": 0 },\n        \"pixels\":  { \"solved\": 0, \"total\": 0 },\n    }\n    for delta in range(1,5+1):\n        df_delta = df[ df.delta == delta ]\n        idxs     = df_delta.index\n        boards   = csv_to_numpy_list(df_delta, key='stop')\n        for idxs, inputs in zip(batch(idxs, 100), batch(boards, 100)):\n            outputs = inputs\n            for t in range(delta):\n                outputs = model.predict(outputs)[:,0,:,:]\n            for idx, output_board, input_board in zip(idxs, outputs, inputs):\n                stats['boards']['total']   += 1\n                stats['delta']['total']    += 1\n                stats['pixels']['total']   += outputs.size\n                stats['pixels']['solved']  += np.count_nonzero( outputs == inputs )\n                stats['dpixels']['total']  += outputs.size\n                stats['dpixels']['solved'] += np.count_nonzero( outputs == inputs )\n                if is_valid_solution(output_board, input_board, delta):\n                    stats['boards']['solved'] += 1\n                    stats['delta']['solved']  += 1\n                    exact_submission.loc[idx] = numpy_to_series(output_board, key='start')\n                submission_df.loc[idx]        = numpy_to_series(output_board, key='start')\n        time_taken = time.perf_counter() - time_start\n        print(f\"delta = {delta} | solved {stats['delta']['solved']:4d}/{stats['delta']['total']:5d} = {100*stats['delta']['solved']/stats['delta']['total']:4.1f}% | {100*stats['dpixels']['solved']/stats['dpixels']['total']:4.1f}% pixels | in {time_taken//60:.0f}:{time_taken%60:02.0f}\")\n        stats['delta']   = { \"solved\": 0, \"total\": 0 }\n        stats['dpixels'] = { \"solved\": 0, \"total\": 0 }\n\n    time_taken = time.perf_counter() - time_start\n    print(f\"ouroborors_dataframe() - solved {stats['boards']['solved']:4d}/{stats['boards']['total']:5d} = {100*stats['boards']['solved']/stats['boards']['total']:4.1f}% | {100*stats['pixels']['solved']/stats['pixels']['total']:4.1f}% pixels | in {time_taken//60:.0f}:{time_taken%60:02.0f}\")\n    submission_df.sort_index().to_csv('submission.csv')\n    exact_submission.sort_index().to_csv('submission_exact.csv')\n    return submission_df\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ouroborors_dataframe(test_df)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Further Reading\n\nI have written an interactive playable demo of the forward version of this game in React Javascript:\n- https://life.jamesmcguigan.com/\n\n\nThis notebook is part of series exploring the Neural Network implementions of the Game of Life Foward Problem\n- [Pytorch Game of Life - First Attempt](https://www.kaggle.com/jamesmcguigan/pytorch-game-of-life-first-attempt)\n- [Pytorch Game of Life - Hardcoding Network Weights](https://www.kaggle.com/jamesmcguigan/pytorch-game-of-life-hardcoding-network-weights)\n- [Its Easy for Neural Networks To Learn Game of Life](https://www.kaggle.com/jamesmcguigan/its-easy-for-neural-networks-to-learn-game-of-life)\n\nThis is preliminary research towards the harder Reverse Game of Life problem, for which I have already designed a novel Ouroboros loss function: \n- [OuroborosLife - Function Reversal GAN](https://www.kaggle.com/jamesmcguigan/ouroboroslife-function-reversal-gan)\n\n\nI also have an extended series of Notebooks exploring different approaches to the Reverse Game of Life problem\n\nMy first attempt was to use the Z3 Constraint Satisfaction SAT solver. This gets 100% accuracy on most boards, but there are a few which it cannot solve. This approach can be slow for boards with large cell counts and large deltas. I managed to figure out how to get cluster compute working inside Kaggle Notebooks, but this solution is estimated to require 10,000+ hours of CPU time to complete.    \n- [Game of Life - Z3 Constraint Satisfaction](https://www.kaggle.com/jamesmcguigan/game-of-life-z3-constraint-satisfaction)\n\nSecond approach was to create a Geometrically Invarient Hash function using Summable Primes, then use forward play and a dictionary lookup table to create a database of known states. For known input/output states at a given delta, the problem is reduced to simply solving the geometric transform between inputs and applying the same function to the outputs. The Hashmap Solver was able to solve about 10% of the test dataset. \n- [Summable Primes](https://www.kaggle.com/jamesmcguigan/summable-primes)\n- [Geometric Invariant Hash Functions](https://www.kaggle.com/jamesmcguigan/geometric-invariant-hash-functions)\n- [Game of Life - Repeating Patterns](https://www.kaggle.com/jamesmcguigan/game-of-life-repeating-patterns)\n- [Game of Life - Hashmap Solver](https://www.kaggle.com/jamesmcguigan/game-of-life-hashmap-solver)\n- [Game of Life - Image Segmentation Solver](https://www.kaggle.com/jamesmcguigan/game-of-life-image-segmentation-solver)"}],"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}