{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Pytorch Game of Life\n\nConway's Game of Life is an example of 2D cellular automata.\n\nBefore going for the difficult \"reverse\" problem, lets see if we can solve the easier \"forward\" problem first and implement the `life_step()` function as a neural network.\nWe already have a code implemention for this function, so the real challenge is getting a neural network to train to 100% accuracy, which we will indeed see is technically possible. \n\nFor this we will use PyTorch and a pure CNN architure. I have provided a detailed explination of my choices for network sizing, data flows, and state space analysis.\n\nThis network was able to implement the `life_step()` function with 100% accuracy tested over 1 million generated boards. \n\nIn batch mode, the profiler reports that the neural network implemention can be faster than `numpy.roll()`, numba `@njit`, and even `scipy.signal.convolve2d()`.\n\n\n## Update\n\nThis model, with 128 CNN layers, is greatly oversized compared to the theoretical minimum, \nbut it reliably trains from random weight initialization acted as proof of concept that a \nneural network can indeed be trained to 100% accuracy.\n\nIt turns out the Game of Life Forward function can actually be implemented using a minimist network architecture with only 1 or 2 CNN 3x3 channels and only 3 three neurons arranged in 2 layers.\n\nI have written a tutorial on hardcoding neural network weights, and how to implement counting and boolean logic gates using linear algebra.\n- https://www.kaggle.com/jamesmcguigan/pytorch-game-of-life-hardcoding-network-weights/"},{"metadata":{},"cell_type":"markdown","source":"# Classical Implementions of life_step()\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.\n\nHere are three different implementions of the `life_step()` function, using `numpy.roll()`, `scipy.signal.convolve2d()` and an numba @njit optimized version in pure python."},{"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_numpy(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_scipy(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@njit\ndef life_step_njit(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\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@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\n\n# For generating large quantities of data, we can use joblib Parallel() to take advantage of all 4 CPUs available in a Kaggle Notebook\nlife_step = life_step_njit\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 )","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Rather than relying on the public train dataset, we can just reimplement the kaggle dataset generation function"},{"metadata":{"trusted":true},"cell_type":"code","source":"# RULES: https://www.kaggle.com/c/conway-s-reverse-game-of-life/data\ndef generate_random_board() -> np.ndarray:\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(25, 25, density=density, data_rvs=np.ones).toarray()\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: return generate_random_board()  # exclude empty boards and try again\n    return board\n\ndef generate_random_boards(count) -> List[np.ndarray]:\n    generated_boards = Parallel(-1)( delayed(generate_random_board)() for _ in range(count) )\n    return generated_boards","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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# First check if we have a GPU available \n# __file__ is implictly defined when running on localhost, but needs to be manually set when running inside a Kaggle Notebook   \nimport torch\ndevice   = torch.device(\"cuda:0\") if torch.cuda.is_available() else torch.device(\"cpu\")\n__file__ = './notebook.ipynb'","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","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\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\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        \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, discarding saved weights\\n')\n                pass\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\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Pytorch CNN\n\nNow comes the application logic, in the form of a CNN neural network.\n\nWe have one 3x3 convolution, followed by three layers of 1x3 convolutions, we then also remind the final output layer of the original input before making a final prediction.\n\nNotice that there is no nn.Linear layer here that gets the see the entire board state. The effect of this is that we have actually used CNN semantics to created a small 4 layer (128:16:8:1) dense network, with 9 inputs and 1 output, that is reused to predict each pixel individually dependent only on its immediate neighbours. This saves us having to write a looping function over the entire board, and we get `padding_mode='circular'` for free.\n\nIf we think about the core logic implemented in the Game Of Life `life_step()` function, all it is really doing is counting up the number of neighbourhood cells and comparing it to the value of the center cell (hence why we pass in the original board state back into the final layer). T\n- 1 Alive + 2-3 neighbours = 1 Alive\n- 0 Dead  +   3 neighbours = 1 Alive\n- Otheriwse                = 0 Dead \n\nTo correctly size the layers in the neural network, we need to think in terms of state complexity. \n\nThe function only needs to know about its distance 1 nearest neighbours, which requires a single 3x3 convolution. There are 2^9 = 512 possible states given 3x3 binary pixels, however this function only needs to be able to count to either 2 or 3. If we do the maths [(9 choose 2) == 34](https://www.wolframalpha.com/input/?i=9+choose+2) + [(9 choose 3) == 84](https://www.wolframalpha.com/input/?i=9+choose+3) = [120](https://www.wolframalpha.com/input/?i=sum%289+choose+n%29+for+n%3D2..3) and we can round that up to the nearest power of 2 to give the network a little extra capacity. \n\nIn theory it might be possible for the neural network to implement `(n choose 3)`  as `(n choose 2) && (n choose 1) && (n choose 1) not in (n choose 2)` by using a clever choose of state representations, which might require less states, but lets stick to 128 input layers for now.\n\n\nNeural Networks are capable of implementing multi-imput boolean logic gates. NAND, NOR, NOT can all be implemented in a single layer, XNOR requires 2 layers and XOR requires 3.\n- https://medium.com/@stanleydukor/neural-representation-of-and-or-not-xor-and-xnor-logic-gates-perceptron-algorithm-b0275375fea1\n\nThus for a neural network trying to implement a boolean logic function, we need 3 layers in a pyramid design and let the output layer act as a final AND gate to choose between the following rule states, which will have been computed by the previous layers:\n- 1 Alive + <2 neighbous  = 1 Dead\n- 1 Alive + >3 neighbous  = 1 Dead\n- 1 Alive + ==2           = 1 Alive\n- 1 Alive + ==2           = 1 Alive\n- 0 Dead  + <3 neighbous  = 1 Dead\n- 0 Dead  + >3 neighbous  = 1 Dead\n- 0 Dead  + ==3 neighbous = 1 Alive\n\nReLU `max(0,n)` has a very easy to compute derivative function `n if n > 0 else 0`. It essentially implements a greater-than `>` or less-than `<`, with distance given in only one direction (returning 0 for anything on the other side of the line). LeakyReLU `0.01*n if n < 0 else n` however can compute distances in both directions by using a negative slope with a different order of magnitude.\n\nWe also add in dropout, which adds a layer of double-check redundancy to the network, as it needs to achieve 100% accuracy with any arbitrary 10% of the network being removed, thus every logic node must essentially have a backup.\n\nThe `__init__()` function simply defines the layers, and the `forward()` function calls them in sequence, adding `F.leaky_relu()` and `self.dropout()` between them. `torch.cat()` is used to append the original X value to the inputs of the final output layer, which is passed though a `sigmoid()` function to produce an output in the 0-1 domain, which is what we want for predicting `1 = Alive` or `0 = Dead` output cells."},{"metadata":{"trusted":true},"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\n\nclass GameOfLifeForward(GameOfLifeBase):\n    \"\"\"\n    This implements the life_step() function as a Neural Network function\n    Training/tested to 100% accuracy over 10,000 random boards\n    \"\"\"\n    def __init__(self):\n        super().__init__()\n        self.dropout = nn.Dropout(p=0.0)  # disabling dropout as it causes excessive (hour+) runtimes to reach 100% accuracy \n\n        # epoch: 240961 | board_count: 6024000 | loss: 0.0000000000 | accuracy = 1.0000000000 | time: 0.611ms/board | with dropout=0.1\n        # Finished Training: GameOfLifeForward - 240995 epochs in 3569.1s\n        self.conv1   = nn.Conv2d(in_channels=1,     out_channels=128, kernel_size=(3,3), padding=1, padding_mode='circular')\n        self.conv2   = nn.Conv2d(in_channels=128,   out_channels=16,  kernel_size=(1,1))\n        self.conv3   = nn.Conv2d(in_channels=16,    out_channels=8,   kernel_size=(1,1))\n        self.conv4   = nn.Conv2d(in_channels=1+8,   out_channels=1,   kernel_size=(1,1))\n        self.apply(self.weights_init)\n\n\n    def forward(self, x):\n        x = input = self.cast_inputs(x)\n\n        x = self.conv1(x)\n        x = F.leaky_relu(x)\n\n        x = self.conv2(x)\n        x = F.leaky_relu(x)\n        x = self.dropout(x)\n\n        x = self.conv3(x)\n        x = F.leaky_relu(x)\n        x = self.dropout(x)\n\n        x = torch.cat([ x, input ], dim=1)  # remind CNN of the center cell state before making final prediction\n        x = self.conv4(x)\n        x = torch.sigmoid(x)\n\n        return x","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training Loop\n\nTo be honest, this could have been written in a more object oriented and resuable fashion, however this has the advantage of showing the full training loop logic in one place.\n\nFirst we load the model, and call the `.train()` and `.unfreeze()` functions to set the model into training mode. Pytorch by default creates models in this state, but my base class has decided to default to production mode with calls to `.eval()` and `.freeze()` being called as part of the `.load()` function. Note that `.freeze()`, `.unfreeze()`, `.load()` and `.save()` are not native to pytorch and have been custom defined in the `GameOfLifeBase` class. \n\nOptimizer chosen is RMSProp with momentum, using a CyclicLR scheduler with that decays over time. The criterion loss function is MSELoss which is defined in the GameOfLifeBase class.\n\nI have obverseved that when using CyclicLR, it is possible to set a much higher maximum learning rates, because it will start off slow at `min_lr=1e-3`, make sure the gradient is pointing in the right direction before reaching `max_lr=1` after 250 epochs, which helps speed up the training. Combined with dropout, this helps prevent the neural network from getting stuck on solutions such as everything=1 or everything=0.\n\nMost of the code in the loop is just data generation, keeping track of statistics, as well as exit conditions such as timeout or 100% accuracy over the last 10,000 boards.\n\nThe code also contains unused implementions for computing L1 and L2 regularization losses by simply looping over the model parameters and summing up the weights. \n\n\nThe meat of the training loop is simply:\n```\noptimizer.zero_grad()\noutputs = model(inputs)\nloss    = model.loss(outputs, expected, inputs)\nloss.backward()\noptimizer.step()\n```\n\nThe `.loss()` and `.accuracy()` functions defined in the `GameOfLifeBase` class.  \n```\nclass GameOfLifeBase:\n    def loss(self, outputs, expected, input):\n        self.criterion = nn.MSELoss()\n        return self.criterion(outputs, expected)\n        \n    def accuracy(self, outputs, expected, inputs) -> float:\n        return torch.sum(self.cast_int(outputs) == self.cast_int(expected)).cpu().numpy() / np.prod(outputs.shape)\n```\n\nNote that the `accuracy()` metric is checking for an exact match (after int rounding) between the output and expected state, whereas the `loss()` function is providing a Mean Squared Error heuristic for number of incorect pixels and distance from neural network output floats to expected state."},{"metadata":{"trusted":true},"cell_type":"code","source":"import atexit\nimport sys\nimport time\n\nimport numpy as np\nimport torch\nimport torch.optim as optim\nfrom torch import tensor\n\n\ndef train(model, batch_size=25, l1=0, l2=0, timeout=0, reverse_input_output=False):\n    print(f'Training: {model.__class__.__name__}')\n    time_start = time.perf_counter()\n\n    atexit.register(model.save)      # save on exit - BrokenPipeError doesn't trigger finally:\n    model.load().train().unfreeze()  # enable training and dropout\n\n    # NOTE: criterion loss function now defined via model.loss()\n    optimizer = optim.RMSprop(model.parameters(), lr=0.01, momentum=0.9)\n\n    # epoch: 14481 | board_count: 362000 | loss: 0.0000385726 | accuracy = 0.9990336000 | time: 0.965ms/board\n    scheduler = None\n\n    # epoch: 240961 | board_count: 6024000 | loss: 0.0000000000 | accuracy = 1.0000000000 | time: 0.611ms/board\n    # Finished Training: GameOfLifeForward - 240995 epochs in 3569.1s\n    scheduler_config = {\n        'optimizer': optimizer,\n        'max_lr':       1e-2,\n        'base_lr':      1e-4,\n        'step_size_up': 250,\n        'mode':         'exp_range',\n        'gamma':        0.8\n    }\n    scheduler = torch.optim.lr_scheduler.CyclicLR(**scheduler_config)\n\n    success_count = 10_000\n    \n    epoch        = 0\n    board_count  = 0\n    last_loss    = np.inf\n    loop_loss    = 0\n    loop_acc     = 0\n    loop_count   = 0\n    epoch_losses     = [last_loss]\n    epoch_accuracies = [ 0 ]\n    num_params = torch.sum(torch.tensor([\n        torch.prod(torch.tensor(param.shape))\n        for param in model.parameters()\n    ]))\n    try:\n        for epoch in range(1, sys.maxsize):\n            if np.min(epoch_accuracies[-success_count//batch_size:]) == 1.0:   break  # multiple epochs of 100% accuracy to pass\n            if timeout and timeout < time.perf_counter() - time_start:  break\n            epoch_start = time.perf_counter()\n\n            inputs_np   = [ generate_random_board() for _     in range(batch_size) ]\n            expected_np = [ life_step(board)        for board in inputs_np         ]\n            inputs      = model.cast_inputs(inputs_np).to(device)\n            expected    = model.cast_inputs(expected_np).to(device)\n\n            # This is for GameOfLifeReverseOneStep() function, where we are trying to learn the reverse function\n            if reverse_input_output:\n                inputs_np, expected_np = expected_np, inputs_np\n                inputs,    expected    = expected,    inputs\n                assert np.all( life_step(expected_np[0]) == inputs_np[0] )\n\n\n            optimizer.zero_grad()\n            outputs = model(inputs)\n            loss    = model.loss(outputs, expected, inputs)\n            if l1 or l2:\n                l1_loss = torch.sum(tensor([ torch.sum(torch.abs(param)) for param in model.parameters() ])) / num_params\n                l2_loss = torch.sum(tensor([ torch.sum(param**2)         for param in model.parameters() ])) / num_params\n                loss   += ( l1_loss * l1 ) + ( l2_loss * l2 )\n\n            loss.backward()\n            optimizer.step()\n            if scheduler is not None:\n                # scheduler.step(loss)  # only required for\n                scheduler.step()\n\n            # noinspection PyTypeChecker\n            last_accuracy = model.accuracy(outputs, expected, inputs)  # torch.sum( outputs.to(torch.bool) == expected.to(torch.bool) ).cpu().numpy() / np.prod(outputs.shape)\n            last_loss     = loss.item() / batch_size\n\n            epoch_losses.append(last_loss)\n            epoch_accuracies.append( last_accuracy )\n\n            loop_loss   += last_loss\n            loop_acc    += last_accuracy\n            loop_count  += 1\n            board_count += batch_size\n            epoch_time   = time.perf_counter() - epoch_start\n            time_taken   = time.perf_counter() - time_start\n\n            # Print statistics after each epoch\n            if( (epoch <= 10)\n             or (np.log10(board_count) % 1 == 0)\n             or (board_count <  10_000 and board_count %  1_000 == 0)  \n             or (board_count < 100_000 and board_count % 10_000 == 0) \n             or (board_count % 100_000 == 0)\n            ):\n                print(f'epoch: {epoch:6d} | board_count: {board_count:7d} | loss: {loop_loss/loop_count:.10f} | accuracy = {loop_acc/loop_count:.10f} | time: {1000*epoch_time/batch_size:.3f}ms/board | {time_taken//60:2.0f}m {time_taken%60:02.0f}s')\n                loop_loss  = 0\n                loop_acc   = 0\n                loop_count = 0\n        print(f'Successfully trained to 100% accuracy over the last {success_count} boards')\n        print(f'epoch: {epoch:6d} | board_count: {board_count:7d} | loss: {np.mean(epoch_losses[-success_count//batch_size:]):.10f} | accuracy = {np.min(epoch_accuracies[-success_count//batch_size:]):.10f} | time: {1000*epoch_time/batch_size:.3f}ms/board')\n                \n    except (BrokenPipeError, KeyboardInterrupt):\n        pass\n    except Exception as exception:\n        print(exception)\n        raise exception\n    finally:\n        time_taken = time.perf_counter() - time_start\n        print(f'Finished Training: {model.__class__.__name__} - {epoch} epochs in {time_taken:.1f}s')\n        model.save()\n        atexit.unregister(model.save)   # model now saved, so cancel atexit handler\n        # model.eval()                  # disable dropout","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training\n\nHere we run the training loop. \n\nLearnings:\n- There is a long wait time between 99.999% accuracy and 100%, which is due to the network figuring out how to implement robust redundancy through the dropout layer\n- `max_lr=1`, even using CycleLR, is far too high, and was causing weights to get stuck in local minima such as all zeros or all ones. Lowering the learning rate produced more stable training cycles.\n- Conv2D weights should be explictly initialized using `nn.init.kaiming_normal_()`. Before implementing this, I had observed that even with a notebook restart, there would sometimes be very high accuracy levels after epoch 1, leading me to suspect that CUDA memory was being reused.\n- One intresting observation from my localhost experiements, is that restarting training during the era of 99.999% accuracy will instantly result in a 100% accuracy score. I have not be able to explain this effect yet, but please comment below if you understand this effect.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# !rm ./models/GameOfLifeForward.pth\nmodel = GameOfLifeForward()\ntrain(model)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Testing\n\nDid you know that you can write unit tests for Neural Networks?\n\nWe are asserting 100% accuracy for our neural network, so lets see if it can correctly predict a million boards in a row.\n\nNote that the `GameOfLifeBase` class defaults to `.eval()` production mode, which disables the dropout layers (which are only applied during training)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def test_GameOfLifeForward_single():\n    \"\"\"\n    Test GameOfLifeForward().predict() works with single board datastructure semantics\n    \"\"\"\n    model = GameOfLifeForward()\n\n    input    = generate_random_board()\n    expected = life_step(input)\n    output   = model.predict(input)\n    assert np.all( output == expected )  # assert 100% accuracy\n\n\n\ndef test_GameOfLifeForward_batch(count=1000):\n    \"\"\"\n    Test GameOfLifeForward().predict() also works batch-mode datastructure semantics\n    Batch mode is limited to only 1000 boards at a time, which prevents CUDA out of memory exceptions\n    \"\"\"\n    model = GameOfLifeForward()\n\n    # As well as in batch mode\n    for _ in range(max(1,count//1000)):\n        inputs   = generate_random_boards(1000)\n        expected = life_steps(inputs)\n        outputs  = model.predict(inputs)\n        assert np.all( outputs == expected )  # assert 100% accuracy","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# GameOfLifeForward can successfully predict a million boards in a row correctly\ntime_start = time.perf_counter()\ntest_GameOfLifeForward_single()\ntest_GameOfLifeForward_batch(1_000_000)\ntime_taken = time.perf_counter() - time_start\nprint(f'All tests passed! ({time_taken:.0f}s)')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Profiling\n\nSo how fast are GPU neural networks really?\n\nWe can write a profiling loop and speed test them against our range of classical implementions.\n\nIt is also worth noting that neural networks can be used to either predict a single output, or in batch mode to predict an array of input data. Due to the python overhead of function calling and getting the data into and out of CUDA, batch mode is where Neural Network performance really shines, and in this case it can even outperform `numpy.roll()` numba `@njit` and `scipy.signal.convolve2d()`"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Discovery:\n# neural networks in batch mode can be faster than C compiled classical code, but slower when called in a loop\n# also scipy uses numba @njit under the hood\n#\n# number: 1 | batch_size = 1000\n# life_step() - numpy         = 181.7µs\n# life_step() - scipy         = 232.1µs\n# life_step() - njit          = 433.5µs   # includes @njit compile times\n# gameOfLifeForward() - loop  = 567.7µs\n# gameOfLifeForward() - batch = 1807.5µs  # includes .pth file loadtimes\n#\n# number: 100 | batch_size = 1000\n# gameOfLifeForward() - batch =  29.8µs  # faster than even @njit or scipy\n# life_step() - scipy         =  35.6µs\n# life_step() - njit          =  42.8µs\n# life_step() - numpy         = 180.8µs\n# gameOfLifeForward() - loop  = 618.3µs  # much slower, but still fast compared to an expensive function\n\ndef profile_GameOfLifeForward():\n    import timeit\n    import operator\n\n    model  = GameOfLifeForward().load().to(device)\n    for batch_size in [1, 10_000]:\n        boards = generate_random_boards(batch_size)\n        number = 1\n        timings = {\n            'gameOfLifeForward() - batch': timeit.timeit(lambda:   model(boards),                                number=number),\n            'gameOfLifeForward() - loop':  timeit.timeit(lambda: [ model(board)           for board in boards ], number=number),\n            'life_step() - njit':          timeit.timeit(lambda: [ life_step_njit(board)  for board in boards ], number=number),\n            'life_step() - numpy':         timeit.timeit(lambda: [ life_step_numpy(board) for board in boards ], number=number),\n            'life_step() - scipy':         timeit.timeit(lambda: [ life_step_scipy(board) for board in boards ], number=number),\n        }\n        print(f'{device} | batch_size = {len(boards)}')\n        for key, value in sorted(timings.items(), key=operator.itemgetter(1)):\n            print(f'- {key:27s} = {value/number/len(boards) * 1_000:7.3f}ms')\n        print()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"device = torch.device(\"cpu\")\nprofile_GameOfLifeForward()\n\nif torch.cuda.is_available():\n    device = torch.device(\"cuda:0\")\n    profile_GameOfLifeForward()","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 implementations of the Game of Life Forward 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}