notebooks/glaberish_animations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5438b425",
   "metadata": {},
   "source": [
    "# Supporting material for \"Glaberish: Generalizing the Continuously-Valued Lenia Framework to Arbitrary Life-Like Cellular Automata\" -- animations\n",
    "\n",
    "* Read the paper: [https://arxiv.org/abs/2205.10463](https://arxiv.org/abs/2205.10463)\n",
    "* Check out the repository: [https://github.com/rivesunder/yuca](https://github.com/rivesunder/yuca)\n",
    "\n",
    "Some cells take a significant amount of wall time to execute, which may be a problem for free cloud notebooks like mybinder.org. For shorter execution times, experiment with changing the number of simulation steps or CA grid dimensions. In this notebook, these cells can take a bit longer to run: \n",
    "   \n",
    "* Code for animations visualizing mobile patterns and spatial entropy over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12323740",
   "metadata": {},
   "outputs": [],
   "source": [
    "# this section needs to be run for setup on colab\n",
    "# uncomment and run the commands below if on colab\n",
    "# then restart the runtime and run tests in the next cell\n",
    "\n",
    "\n",
    "#! git clone https://github.com/riveSunder/yuca.git my_yuca\n",
    "#%cd my_yuca\n",
    "#! pip install -e ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "029e7777",
   "metadata": {},
   "outputs": [],
   "source": [
    "#! python -m testing.test_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6136a86e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make  cells wider\n",
    "from IPython.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44920194",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import copy\n",
    "import numpy as np\n",
    "import scipy\n",
    "import torch\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.animation \n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "my_cmap = plt.get_cmap(\"magma\")\n",
    "plt.rcParams['figure.dpi'] = 300\n",
    "\n",
    "matplotlib.rcParams[\"pdf.fonttype\"] = 42\n",
    "matplotlib.rcParams[\"ps.fonttype\"] = 42\n",
    "matplotlib.rcParams[\"animation.embed_limit\"] = 256\n",
    "matplotlib.rcParams[\"font.size\"] = 10\n",
    "#matplotlib.rcParams[\"text.usetex\"] = True\n",
    "\n",
    "import IPython\n",
    "\n",
    "import skimage\n",
    "import skimage.io as sio\n",
    "import skimage.transform\n",
    "\n",
    "import notebooks.my_helpers as my_helpers\n",
    "from importlib import reload\n",
    "reload(my_helpers)\n",
    "\n",
    "import yuca\n",
    "from yuca.utils import seed_all\n",
    "from yuca.patterns import get_puffer\n",
    "from yuca.multiverse import CA\n",
    "from yuca.lenia import Lenia\n",
    "from yuca.zoo.librarian import Librarian\n",
    "\n",
    "\n",
    "from yuca.analyze import get_spatial_entropy, \\\n",
    "        get_multiconditional_grid_entropy, \\\n",
    "        get_conditional_grid_entropy,\\\n",
    "        get_grid_entropy\n",
    "\n",
    "my_device = torch.device(\"cpu\")\n",
    "\n",
    "def get_center_of_mass(grid, spacing=20):\n",
    "    \n",
    "    if len(grid.shape) > 2:\n",
    "        my_grid = np.array(grid.cpu().squeeze())\n",
    "    else:\n",
    "        my_grid = np.array(1.0 * grid.cpu())\n",
    "        \n",
    "    edge_l = my_grid[:,:spacing].sum()\n",
    "    edge_r = my_grid[:,-spacing:].sum()\n",
    "    edge_t = my_grid[:spacing,:].sum()\n",
    "    edge_b = my_grid[-spacing:,:].sum()\n",
    "    \n",
    "    max_edge = np.max([edge_l, edge_r, edge_t, edge_b])\n",
    "    \n",
    "    xx_grid, yy_grid = np.meshgrid(np.arange(my_grid.shape[-1]), np.arange(my_grid.shape[-2]))\n",
    "    \n",
    "    if max_edge:\n",
    "        \n",
    "        if edge_r == max_edge:\n",
    "            com_x = xx_grid[:, spacing:] * my_grid[:, spacing:] / my_grid[:, spacing:].sum()\n",
    "            com_y = yy_grid[:, spacing:] * my_grid[:, spacing:] / my_grid[:, spacing:].sum()\n",
    "            \n",
    "        elif edge_l == max_edge:\n",
    "            com_x = xx_grid[:, :-spacing] * my_grid[:, :-spacing] /  my_grid[:, :-spacing].sum()\n",
    "            com_y = yy_grid[:, :-spacing] * my_grid[:, :-spacing] /  my_grid[:, :-spacing].sum()\n",
    "            \n",
    "        elif edge_t == max_edge:\n",
    "            com_x = xx_grid[:-spacing, :] * my_grid[:-spacing, :] / my_grid[:-spacing, :].sum()\n",
    "            com_y = yy_grid[:-spacing, :] * my_grid[:-spacing, :] /my_grid[:-spacing, :].sum()\n",
    "            \n",
    "        else: # edge_b == max_edge:\n",
    "            com_x = xx_grid[spacing:, :] * my_grid[spacing:, :] / my_grid[spacing:, :].sum()\n",
    "            com_y = yy_grid[spacing:, :] * my_grid[spacing:, :] / my_grid[spacing:, :].sum()\n",
    "            \n",
    "    else:\n",
    "        # pattern is not at an edge crossing\n",
    "        \n",
    "        com_x = (xx_grid * my_grid) / my_grid.sum()\n",
    "        com_y = (yy_grid * my_grid) / my_grid.sum()\n",
    "    \n",
    "    return com_x.sum(), com_y.sum()\n",
    "\n",
    "def plot_grid_nbhd_update_next(grid, nbhd, update, next_grid, \\\n",
    "                               my_cmap=plt.get_cmap(\"magma\"), titles=None, vmin=0.0, vmax=1):\n",
    "\n",
    "    global subplot_0\n",
    "    global subplot_1\n",
    "    global subplot_2\n",
    "    global subplot_3\n",
    "\n",
    "    if titles == None:\n",
    "        titles = [\"CA grid time t\", \"Neighborhood\", \"Update\", \"CA grid time t+1\"]\n",
    "\n",
    "    fig, ax = plt.subplots(2,2, figsize=(2,2), facecolor=\"white\")\n",
    "\n",
    "    subplot_0 = ax[0,0].imshow(grid, cmap=my_cmap, vmin=vmin, vmax=vmax, interpolation=\"nearest\") \n",
    "    ax[0,0].set_title(titles[0], fontsize=4)\n",
    "\n",
    "    subplot_1 = ax[0,1].imshow(nbhd, cmap=my_cmap, vmin=vmin, vmax=vmax, interpolation=\"nearest\")\n",
    "    ax[0,1].set_title(titles[1], fontsize=4)\n",
    "\n",
    "    my_vmin = -1.0\n",
    "    my_vmax = 1.0\n",
    "    subplot_2 = ax[1,0].imshow(update, cmap=my_cmap, interpolation=\"nearest\", vmin=my_vmin, vmax=my_vmax)\n",
    "    ax[1,0].set_xlabel(titles[2], fontsize=4)\n",
    "\n",
    "    subplot_3 = ax[1,1].imshow(next_grid, cmap=my_cmap, vmin=vmin, vmax=vmax, interpolation=\"nearest\")\n",
    "    ax[1,1].set_xlabel(titles[3], fontsize=4)\n",
    "\n",
    "    for xxyy in range(4):\n",
    "        xx = xxyy // 2\n",
    "        yy = xxyy % 2\n",
    "        ax[xx,yy].set_xticklabels(\"\")\n",
    "        ax[xx,yy].set_yticklabels(\"\")\n",
    "    #plt.tight_layout()\n",
    "    return fig, ax\n",
    "\n",
    "def update(i):\n",
    "    global grid\n",
    "\n",
    "    n = conv_moore(grid)\n",
    "    m = conv_eye(grid)\n",
    "\n",
    "    if (use_lenia):\n",
    "        # Lenia defines updates solely by a growth function without considering current state\n",
    "        g = combined_growth_fn(n)\n",
    "    else:\n",
    "        # Generally Life-like splits growth into genesis and persistence\n",
    "        g = ((1-m) * growth_fn_b(n) + m * growth_fn_s(n)) * 2 - 1\n",
    "        \n",
    "    subplot_0.set_array(m)\n",
    "    subplot_1.set_array(n)\n",
    "    subplot_2.set_array(g)\n",
    "    subplot_3.set_array(grid)\n",
    "    \n",
    "    grid = np.clip(grid + dt * g, 0, 1.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "daed31a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# execute this cell for gpu support\n",
    "my_device = torch.device(\"cuda\") if torch.cuda.is_available() else torch.device(\"cpu\")\n",
    "print(my_device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc188cdf",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# life works in Lenia \n",
    "\n",
    "# Conway's Game of Life: Implementable in the Lenia Framework\n",
    "\n",
    "gol_b = [3]\n",
    "gol_s = [2,3]\n",
    "\n",
    "b = my_helpers.get_smooth_interval(gol_b) \n",
    "s = my_helpers.get_smooth_interval(gol_s) \n",
    "\n",
    "kernel_args = {\"mu\": 0.99, \"sigma\": 0.3}\n",
    "kernel_moore = my_helpers.get_kernel(my_helpers.gaussian, kernel_args, radius=1)\n",
    "\n",
    "kernel_moore[kernel_moore > 0.01 ] = 1.0 \n",
    "kernel_moore[kernel_moore <= 0.01 ] = 0.0 \n",
    "kernel_moore /= kernel_moore.sum()\n",
    "#kernel_moore *= 0.9\n",
    "\n",
    "kernel_args = {\"mu\": 0.0, \"sigma\": 0.1}\n",
    "kernel_eye = my_helpers.get_kernel(my_helpers.gaussian, kernel_args, radius=1)\n",
    "\n",
    "kernel_eye[kernel_eye < 0.99 ] = 0.0 \n",
    "kernel_eye /= kernel_eye.sum()\n",
    "\n",
    "update_args = {\"a\": b, \"alpha\": 0.0125}\n",
    "growth_fn_b = my_helpers.get_smooth_steps_fn(**update_args)\n",
    "\n",
    "update_args = {\"a\": s, \"alpha\": 0.0125}\n",
    "growth_fn_s = my_helpers.get_smooth_steps_fn(**update_args)\n",
    "\n",
    "combined_growth_fn = lambda x: (growth_fn_b(x) + growth_fn_s(x)) - 1\n",
    "\n",
    "fig = my_helpers.plot_kernel_growth(kernel_moore, combined_growth_fn, vmin=0, vmax=kernel_moore.max(), \\\n",
    "        titles=[\"Moore Kernel\", \"Combined Life \\n Growth Function \\n (Lenia)\"])\n",
    "\n",
    "# uncomment to save figure\n",
    "#plt.savefig(\"assets/life_in_lenia.png\")\n",
    "\n",
    "plt.show()\n",
    "grid = np.zeros((24,24))\n",
    "#glider = my_helpers.get_puffer()\n",
    "\n",
    "glider = np.zeros((3,3))\n",
    "glider[0,1] = 1.\n",
    "glider[1,1:3] = 1.\n",
    "glider[2,0] = 1.\n",
    "glider[2,2] = 1.\n",
    "\n",
    "grid[10:10+glider.shape[1], 4:4+glider.shape[1]] = glider\n",
    "\n",
    "conv_moore = my_helpers.get_conv_fn(kernel_moore)\n",
    "conv_eye = my_helpers.get_conv_fn(kernel_eye)\n",
    "\n",
    "dt = 1.0\n",
    "\n",
    "use_lenia = True\n",
    "num_frames = 64\n",
    "\n",
    "m = grid\n",
    "n = grid\n",
    "g = grid \n",
    "\n",
    "fig, ax = plot_grid_nbhd_update_next(m, n, g, grid, my_cmap=my_cmap)\n",
    "plt.close(\"all\")\n",
    "\n",
    "IPython.display.HTML(matplotlib.animation.FuncAnimation(fig, update, frames=num_frames, interval=100).to_jshtml())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ded8bbe",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Morley: not readily mplementable in the Lenia Framework\n",
    "ca = yuca.lenia.Lenia()\n",
    "\n",
    "ca.restore_config(\"morley.npy\")\n",
    "ca.no_grad()\n",
    "\n",
    "kernel_moore = ca.neighborhood_kernels.squeeze().numpy()\n",
    "\n",
    "growth_fn_b = lambda x: ca.genesis_fns[0](torch.tensor(x).clone()).numpy()\n",
    "growth_fn_s =  lambda x: ca.persistence_fns[0](torch.tensor(x).clone()).numpy()\n",
    "combined_growth_fn = lambda x: (growth_fn_b(x) + growth_fn_s(x)) \n",
    "\n",
    "fig = my_helpers.plot_kernel_growth(kernel_moore, combined_growth_fn, vmin=0, vmax=kernel_moore.max(), \\\n",
    "        titles=[\"Moore Kernel\", \"Combined Morley \\n Growth Function \\n (Lenia)\"])\n",
    "\n",
    "plt.show()\n",
    "\n",
    "grid = np.zeros((24,24))\n",
    "#glider = my_helpers.get_puffer()\n",
    "\n",
    "glider = np.zeros((6,6))\n",
    "glider[2:4,0] = 1.\n",
    "glider[1,1] = 1.\n",
    "glider[4,1] = 1.\n",
    "glider[:,2] = 1.\n",
    "glider[1,3] = 1.\n",
    "glider[4,3] = 1.\n",
    "glider[1,4] = 1.\n",
    "glider[4,4] = 1.\n",
    "glider[2:4,5] = 1.\n",
    "\n",
    "grid[10:10+glider.shape[1], 4:4+glider.shape[1]] = glider\n",
    "\n",
    "conv_moore = my_helpers.get_conv_fn(kernel_moore)\n",
    "conv_eye = my_helpers.get_conv_fn(kernel_eye)\n",
    "\n",
    "dt = 1.0\n",
    "\n",
    "use_lenia = True\n",
    "#rter or longer animation\n",
    "num_frames = 32\n",
    "\n",
    "m = grid\n",
    "n = grid\n",
    "g = grid \n",
    "\n",
    "fig,ax = plot_grid_nbhd_update_next(m, n, g, grid, my_cmap=my_cmap)\n",
    "plt.close(\"all\")\n",
    "\n",
    "\n",
    "IPython.display.HTML(matplotlib.animation.FuncAnimation(fig, update, frames=num_frames, interval=100).to_jshtml())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5d9fc72",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Morley works in Glaberish\n",
    "\n",
    "ca = yuca.multiverse.CA()\n",
    "ca.restore_config(\"morley.npy\")\n",
    "ca.no_grad()\n",
    "\n",
    "kernel_moore = ca.neighborhood_kernels.squeeze().numpy()\n",
    "\n",
    "genesis_fn = lambda x: ca.genesis_fns[0](torch.tensor(x).clone()).numpy()\n",
    "persistence_fn =  lambda x: ca.persistence_fns[0](torch.tensor(x).clone()).numpy()\n",
    "\n",
    "growth_fn_b = genesis_fn\n",
    "growth_fn_s = persistence_fn\n",
    "\n",
    "fig = my_helpers.plot_kernel_genper(kernel_moore, genesis_fn, persistence_fn, vmin=0, vmax=kernel_moore.max(), \\\n",
    "        titles=[\"Moore Kernel\", \"Morley \\n Persistence and genesis \\n (Glaberish)\"])\n",
    "\n",
    "plt.show()\n",
    "\n",
    "grid = np.zeros((36,36))\n",
    "#glider = my_helpers.get_puffer()\n",
    "\n",
    "glider = np.zeros((6,6))\n",
    "glider[2:4,0] = 1.\n",
    "glider[1,1] = 1.\n",
    "glider[4,1] = 1.\n",
    "glider[:,2] = 1.\n",
    "glider[1,3] = 1.\n",
    "glider[4,3] = 1.\n",
    "glider[1,4] = 1.\n",
    "glider[4,4] = 1.\n",
    "glider[2:4,5] = 1.\n",
    "\n",
    "grid[10:10+glider.shape[1], 4:4+glider.shape[1]] = glider\n",
    "\n",
    "conv_moore = my_helpers.get_conv_fn(kernel_moore)\n",
    "conv_eye = my_helpers.get_conv_fn(kernel_eye)\n",
    "\n",
    "dt = 1.0\n",
    "\n",
    "#rter or longer animation\n",
    "num_frames = 256\n",
    "use_lenia = False\n",
    "\n",
    "m = grid\n",
    "n = grid\n",
    "g = grid \n",
    "\n",
    "fig,ax = plot_grid_nbhd_update_next(m, n, g, grid, my_cmap=my_cmap)\n",
    "plt.close(\"all\")\n",
    "\n",
    "\n",
    "IPython.display.HTML(matplotlib.animation.FuncAnimation(fig, update, frames=num_frames, interval=100).to_jshtml())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02cfc881",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Figure showing examples of mobile patterns in s613 and Hydrogeminium\n",
    "\n",
    "num_frames = 512\n",
    "# try fewer frames for faster runtime:\n",
    "# num_frames = 48\n",
    "\n",
    "torch.set_default_dtype(torch.float32)\n",
    "dim_x = 192\n",
    "dim_y = 192\n",
    "\n",
    "lib = Librarian()\n",
    "lib.update_index()\n",
    "\n",
    "grid_s0 = np.zeros((1, 1, dim_x, dim_y))\n",
    "grid_s1 = np.zeros((1, 1, dim_x, dim_y))\n",
    "\n",
    "grid_g0 = np.zeros((1, 1, dim_x, dim_y))\n",
    "grid_g1 = np.zeros((1, 1, dim_x, dim_y))\n",
    "\n",
    "\n",
    "geminium, gem_metadata = lib.load(\"geminium_hydrogeminium_natans000\")\n",
    "gem_config = gem_metadata[\"ca_config\"]\n",
    "gemwobbler, _ = lib.load(\"geminium2_wobble_glider000\")\n",
    "\n",
    "frog, s613_metadata = lib.load(\"s613_s613_frog000\")\n",
    "s613_config = s613_metadata[\"ca_config\"]\n",
    "wobbler, _ = lib.load(\"s613_fast_wobble_glider000\")\n",
    "\n",
    "ca_gem = CA()\n",
    "ca_s613 = CA()\n",
    "\n",
    "ca_gem.restore_config(\"geminium.npy\")\n",
    "ca_s613.restore_config(\"s613.npy\")\n",
    "ca_gem.no_grad()       \n",
    "ca_s613.no_grad()\n",
    "\n",
    "ca_gem.to_device(my_device)\n",
    "ca_s613.to_device(my_device)\n",
    "\n",
    "glider_steps = 64\n",
    "\n",
    "grid_g0[0,0, 60:60+geminium.shape[-2], :geminium.shape[-1]] = geminium\n",
    "grid_g1[0,0, 40:40+gemwobbler.shape[-2], 59:59+gemwobbler.shape[-1]] = gemwobbler\n",
    "                       \n",
    "\n",
    "grid_s0[0,0,:wobbler.shape[-2], :wobbler.shape[-1]] = wobbler\n",
    "grid_s1[0,0, 80:80+frog.shape[-2], 0+wobbler.shape[-1]:0+wobbler.shape[-1]+frog.shape[-1]] = frog\n",
    "                       \n",
    "\n",
    "## update for steps\n",
    "\n",
    "grid_g0 = torch.tensor(grid_g0).to(torch.get_default_dtype()).to(my_device)\n",
    "grid_g1 = torch.tensor(grid_g1).to(torch.get_default_dtype()).to(my_device)\n",
    "\n",
    "grid_s0 = torch.tensor(grid_s0).to(torch.get_default_dtype()).to(my_device)\n",
    "grid_s1 = torch.tensor(grid_s1).to(torch.get_default_dtype()).to(my_device)\n",
    "\n",
    "com_g0 = []\n",
    "com_g1 = []\n",
    "com_s0 = []\n",
    "com_s1 = []\n",
    "\n",
    "def update_gliders(ii):\n",
    "    \n",
    "    global grid_g0\n",
    "    global grid_s0\n",
    "    global grid_g1\n",
    "    global grid_s1\n",
    "    \n",
    "    global subplot_0\n",
    "    global subplot_1\n",
    "    global subplot_2\n",
    "    global subplot_3\n",
    "    \n",
    "    global com_g0\n",
    "    global com_s0\n",
    "    global com_g1\n",
    "    global com_s1\n",
    "    \n",
    "    global ax\n",
    "    \n",
    "    grid_g0 = ca_gem(grid_g0)\n",
    "    grid_g1 = ca_gem(grid_g1)\n",
    "    grid_s0 = ca_s613(grid_s0)\n",
    "    grid_s1 = ca_s613(grid_s1)\n",
    "    \n",
    "    \n",
    "    if (display_center):\n",
    "        \n",
    "        com_g0 = com_g0[-100:]\n",
    "        com_g1 = com_g1[-100:]\n",
    "        com_s0 = com_s0[-100:]\n",
    "        com_s1 = com_s1[-100:]\n",
    "        com_g0.append(get_center_of_mass(grid_g0))   \n",
    "        com_g1.append(get_center_of_mass(grid_g1))\n",
    "        com_s0.append(get_center_of_mass(grid_s0))\n",
    "        com_s1.append(get_center_of_mass(grid_s1))\n",
    "\n",
    "    disp_s0 = 1.0 - my_cmap(grid_s0.cpu().numpy().squeeze())[:,:,:3]\n",
    "    disp_s1 = 1.0 - my_cmap(grid_s1.cpu().numpy().squeeze())[:,:,:3]\n",
    "    disp_g0 = 1.0 - my_cmap(grid_g0.cpu().numpy().squeeze())[:,:,:3]\n",
    "    disp_g1 = 1.0 - my_cmap(grid_g1.cpu().numpy().squeeze())[:,:,:3]\n",
    "\n",
    "    subplot_0.set_array(disp_g0)\n",
    "    subplot_2.set_array(disp_g1)\n",
    "    subplot_1.set_array(disp_s0)\n",
    "    subplot_3.set_array(disp_s1)\n",
    "    \n",
    "    my_green = 1.0 - np.array(my_cmap(grid_g0.cpu().max().numpy()))[:3]\n",
    "    my_blue = 1.0 - np.array(plt.get_cmap(\"afmhot\")(grid_g0.cpu().max().numpy()))[:3]\n",
    "\n",
    "    if (display_center):\n",
    "        for ii in range(len(com_g0)):\n",
    "            progress = ii / len(com_g0)\n",
    "            my_color = progress * my_green + (1.0 - progress) * my_blue\n",
    "            ax[0,0].scatter(com_g0[ii][0],com_g0[ii][1], color=my_color, s=0.5, alpha=0.05) \n",
    "\n",
    "        for ii in range(len(com_g1)):\n",
    "            progress = ii / len(com_g1)\n",
    "            my_color = progress * my_green + (1.0 - progress) * my_blue\n",
    "            ax[1,0].scatter(com_g1[ii][0],com_g1[ii][1], color=my_color, s=0.5, alpha=0.05)    \n",
    "\n",
    "        for ii in range(len(com_s0)):\n",
    "            progress = ii / len(com_s0)\n",
    "            my_color = progress * my_green + (1.0 - progress) * my_blue\n",
    "            ax[0,1].scatter(com_s0[ii][0],com_s0[ii][1], color=my_color, s=0.5, alpha=0.05) \n",
    "\n",
    "        for ii in range(len(com_s1)):\n",
    "            progress = ii / len(com_s1)\n",
    "            my_color = progress * my_green + (1.0 - progress) * my_blue\n",
    "            ax[1,1].scatter(com_s1[ii][0],com_s1[ii][1], color=my_color, s=0.5, alpha=0.05) \n",
    "        \n",
    "    for xx in range(2):\n",
    "        for yy in range(2):\n",
    "            ax[xx, yy].set_xticklabels(\"\")\n",
    "            ax[xx, yy].set_yticklabels(\"\")\n",
    "    \n",
    "\n",
    "disp_s0 = 1.0 - my_cmap(grid_s0.cpu().numpy().squeeze())[:,:,:3]\n",
    "disp_s1 = 1.0 - my_cmap(grid_s1.cpu().numpy().squeeze())[:,:,:3]\n",
    "disp_g0 = 1.0 - my_cmap(grid_g0.cpu().numpy().squeeze())[:,:,:3]\n",
    "disp_g1 = 1.0 - my_cmap(grid_g1.cpu().numpy().squeeze())[:,:,:3]\n",
    "    \n",
    "fig, ax = plot_grid_nbhd_update_next(disp_g0, disp_s0, disp_g1, disp_s1, \\\n",
    "        titles=[\"slow_cucumber\", \"fast, wobbly glider\", \"wobbly glider\", \"frog\"], my_cmap=my_cmap,)\n",
    "\n",
    "plt.close(\"all\")\n",
    "\n",
    "display_center = False\n",
    "IPython.display.HTML(matplotlib.animation.FuncAnimation(fig, update_gliders, frames=num_frames, interval=50).to_jshtml())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a59fdaee",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# visualize spatial entropy over time\n",
    "\n",
    "# static figure used 512 steps, than calculate spatial entropy, 512 more steps and calculate spatial entropy again. \n",
    "num_frames = 384\n",
    "# Use fewer frames for shorter runtime:\n",
    "#num_frames = 128\n",
    "\n",
    "def update_spatial_entropy(ii):\n",
    "    \n",
    "    global grid_g\n",
    "    global grid_s\n",
    "    global grid_hg\n",
    "    global grid_hs\n",
    "    \n",
    "    global subplot_0\n",
    "    global subplot_1\n",
    "    global subplot_2\n",
    "    global subplot_3\n",
    "    \n",
    "    global ax\n",
    "    \n",
    "    grid_g = ca_gem(grid_g)\n",
    "    grid_s = ca_s613(grid_s)\n",
    "    \n",
    "    disp_g = 1.0 - my_cmap(grid_g.cpu().numpy().squeeze())[:,:,:3]\n",
    "    disp_s = 1.0 - my_cmap(grid_s.cpu().numpy().squeeze())[:,:,:3]\n",
    "    \n",
    "    subplot_0.set_array(disp_g)\n",
    "    subplot_2.set_array(disp_s)\n",
    "    \n",
    "    if ii % 32 == 0:\n",
    "        \n",
    "        grid_hg = get_spatial_entropy(grid_g.detach().cpu().numpy().squeeze(), window_size=ca_gem.neighborhood_kernels.shape[-1])\n",
    "        grid_hs = get_spatial_entropy(grid_s.detach().cpu().numpy().squeeze(), window_size=ca_s613.neighborhood_kernels.shape[-1])\n",
    "        disp_hg =  my_cmap(grid_hg.squeeze()/8.)[:,:,:3]\n",
    "        disp_hs =  my_cmap(grid_hs.squeeze()/8.)[:,:,:3]\n",
    "\n",
    "        subplot_1.set_array(disp_hg)\n",
    "        subplot_3.set_array(disp_hs)\n",
    "        print(f\"gem spatial entropy mean {np.mean(grid_hg)} +/- {np.std(grid_hg)} std. dev.\")\n",
    "        print(f\"s613 spatial entropy mean {np.mean(grid_hs)} +/- {np.std(grid_hs)} std. dev.\")\n",
    "        \n",
    "    for xx in range(2):\n",
    "        for yy in range(2):\n",
    "            ax[xx, yy].set_xticklabels(\"\")\n",
    "            ax[xx, yy].set_yticklabels(\"\")\n",
    "            \n",
    "            \n",
    "torch.set_default_dtype(torch.float32)\n",
    "\n",
    "dim_x = 192\n",
    "dim_y = 192\n",
    "half_dim = dim_x // 4\n",
    "\n",
    "grid_s = torch.zeros((1, 1, dim_x, dim_y))\n",
    "grid_g = torch.zeros((1, 1, dim_x, dim_y))\n",
    "\n",
    "my_init = torch.rand(1,1, 2*half_dim, 2*half_dim)\n",
    "\n",
    "grid_s[:,:,half_dim:3 * half_dim, half_dim:3 * half_dim] = my_init \n",
    "grid_g[:,:,half_dim:3 * half_dim, half_dim:3 * half_dim] = my_init \n",
    "\n",
    "grid_s = grid_s.to(my_device)\n",
    "grid_g = grid_g.to(my_device)\n",
    "\n",
    "grid_hs = torch.zeros((1, 1, dim_x, dim_y))\n",
    "grid_hg = torch.zeros((1, 1, dim_x, dim_y))\n",
    "\n",
    "#grid_hg = get_spatial_entropy(grid_g.detach().cpu().numpy().squeeze(), window_size=ca_gem.neighborhood_kernels.shape[-1])\n",
    "#grid_hs = get_spatial_entropy(grid_s.detach().cpu().numpy().squeeze(), window_size=ca_s613.neighborhood_kernels.shape[-1])\n",
    "       \n",
    "ca_gem = CA()\n",
    "ca_s613 = CA()\n",
    "\n",
    "ca_gem.restore_config(\"geminium.npy\")\n",
    "ca_s613.restore_config(\"s613.npy\")\n",
    "ca_gem.no_grad()       \n",
    "ca_s613.no_grad()\n",
    "\n",
    "ca_gem.to_device(my_device)\n",
    "ca_s613.to_device(my_device)\n",
    "\n",
    "grid_hg[0,0,0,0]= 8\n",
    "grid_hs[0,0,0,0]= 8\n",
    "\n",
    "disp_g = 1.0 - my_cmap(grid_g.cpu().numpy().squeeze())[:,:,:3]\n",
    "disp_s = 1.0 - my_cmap(grid_s.cpu().numpy().squeeze())[:,:,:3]\n",
    "\n",
    "disp_hg = my_cmap(grid_hg.squeeze())[:,:,:3]\n",
    "disp_hs = my_cmap(grid_hs.squeeze())[:,:,:3]\n",
    "\n",
    "fig, ax = plot_grid_nbhd_update_next(disp_g, disp_hg, disp_s, disp_hs, \\\n",
    "        titles=[\"${\\it H. natans}$\", \"${\\it H. natans}$ spatial entropy\",\\\n",
    "                \"s613\",\"s613 spatial entropy\"], my_cmap=my_cmap,)\n",
    "plt.show()\n",
    "plt.close(\"all\")\n",
    "\n",
    "\n",
    "IPython.display.HTML(matplotlib.animation.FuncAnimation(fig, update_spatial_entropy, frames=num_frames, interval=50).to_jshtml())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e50007c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}