simple_cells.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import re\n",
    "import time\n",
    "import copy\n",
    "import json\n",
    "from collections import defaultdict\n",
    "from tqdm import tqdm\n",
    "from multiprocess import Pool\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.sparse\n",
    "import scipy.spatial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:root:Setting up a new session...\n"
     ]
    }
   ],
   "source": [
    "def fig2data ( fig ):\n",
    "    \"\"\"\n",
    "    @brief Convert a Matplotlib figure to a 4D numpy array with RGBA channels and return it\n",
    "    @param fig a matplotlib figure\n",
    "    @return a numpy 3D array of RGBA values\n",
    "    Ref: http://www.icare.univ-lille1.fr/tutorials/convert_a_matplotlib_figure\n",
    "    \"\"\"\n",
    "    # Get rid of the margins\n",
    "    fig.gca().set_axis_off()\n",
    "    fig.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, \n",
    "                hspace = 0, wspace = 0)\n",
    "    #fig.margins(0,0)\n",
    "    from matplotlib.ticker import NullLocator\n",
    "    fig.gca().xaxis.set_major_locator(NullLocator())\n",
    "    fig.gca().yaxis.set_major_locator(NullLocator())\n",
    "    \n",
    "    # draw the renderer\n",
    "    fig.canvas.draw ( )\n",
    "\n",
    "    # Get the RGBA buffer from the figure\n",
    "    w,h = fig.canvas.get_width_height()\n",
    "    buf = np.frombuffer(fig.canvas.tostring_rgb(), dtype=np.uint8).reshape((h, w, 3))\n",
    "    #buf = buf.transpose((1, 0, 2))\n",
    "    #buf = np.fromstring ( fig.canvas.tostring_argb(), dtype=np.uint8 )\n",
    "    #buf.shape = ( w, h,4 )\n",
    "    \n",
    "    # canvas.tostring_argb give pixmap in ARGB mode. Roll the ALPHA channel to have it in RGBA mode\n",
    "    #buf = np.roll ( buf, 3, axis = 2 )\n",
    "    return buf\n",
    "\n",
    "def frames2web(frames, fps=60):\n",
    "    \"\"\"Converts NxHxWxC NumPu array into Base64-encode video, wrapped in HTML (for Visdom visualization)\"\"\"\n",
    "    import skvideo.io\n",
    "    import tempfile\n",
    "    import mimetypes\n",
    "    import base64\n",
    "    videofile = '/tmp/%s.mp4' % next(tempfile._get_candidate_names())\n",
    "    skvideo.io.vwrite(videofile, frames,\n",
    "                      inputdict={\n",
    "                        \"-r\": \"%d\" % fps\n",
    "                      },\n",
    "                      outputdict={\n",
    "                        \"-b\": \"4000000\"\n",
    "                      })\n",
    "\n",
    "    extension = videofile[-3:].lower()\n",
    "    mimetypes = dict(mp4='mp4', ogv='ogg', avi='avi', webm='webm')\n",
    "    mimetype = mimetypes.get(extension)\n",
    "\n",
    "    with open(videofile, 'rb') as f:\n",
    "        bytestr = f.read()\n",
    "    videodata = \"\"\"\n",
    "        <video width=\"%d\" controls>\n",
    "            <source type=\"video/%s\" src=\"data:video/%s;base64,%s\">\n",
    "            Your browser does not support the video tag.\n",
    "        </video>\n",
    "    \"\"\" % (frames.shape[2]*1.5, mimetype, mimetype, base64.b64encode(bytestr).decode('utf-8'))\n",
    "    return videodata\n",
    "\n",
    "from IPython.display import HTML\n",
    "\n",
    "\n",
    "from visdom import Visdom\n",
    "viz = Visdom()\n",
    "viz_opts = {\n",
    "    'sim_slice': dict()\n",
    "}\n",
    "viz_windows = { k: None for k in viz_opts.keys() }\n",
    "# vizdom .video() modified to use scikit-video instead of OpenCV\n",
    "def viz_video(frames, win=None, env=None, opts=None, fps=60):\n",
    "    return viz.text(text=frames2web(frames, fps), win=win, env=env, opts=opts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch 0, max fitness 16.0\n"
     ]
    }
   ],
   "source": [
    "n_aux_p = 3\n",
    "env_p = [\n",
    "  'aminoacids',\n",
    "#     'atp',\n",
    "]\n",
    "cell_p = [\n",
    "  'genome',\n",
    "  'x',\n",
    "  'y',\n",
    "  'aminoacids',\n",
    "  'atp',\n",
    "  'cycle',\n",
    "]\n",
    "for i in np.arange(n_aux_p):\n",
    "  env_p.append('aux{0}'.format(i))\n",
    "  cell_p.append('aux{0}'.format(i))\n",
    "env_p_map = {k:v for v, k in enumerate(env_p)}\n",
    "cell_p_map = {k:v for v, k in enumerate(cell_p)}\n",
    "\n",
    "# env2cell_p = {cell_p_map[k] for k, v in env_p_map if k in cell_p_map}\n",
    "# cell2env_p = {env_p_map[k] for k, v in cell_p_map if k in env_p_map}\n",
    "# for i in np.arange(n_aux_p):\n",
    "#   env2cell_p[len(env_p) + i] = len(cell_p) + i\n",
    "#   cell2env_p[len(cell_p) + i] = len(env_p) + i\n",
    "\n",
    "resource_cost = {\n",
    "  cell_p_map['cycle']: 1.0,\n",
    "}\n",
    "for i, p in enumerate(cell_p[-n_aux_p:]):\n",
    "  resource_cost[cell_p_map[p]] = 0.01 / (i + 1)\n",
    "\n",
    "cell2cell_diffusion = {\n",
    "  cell_p_map['aminoacids']: 0.2,\n",
    "  cell_p_map['atp']: 0.2,\n",
    "}\n",
    "env2cell_diffusion = {\n",
    "  'aminoacids': 0.05,\n",
    "#     'atp': 0.05,\n",
    "}\n",
    "for i, p in enumerate(env_p[-n_aux_p:]):\n",
    "  env2cell_diffusion[p] = 0.1 / (i + 1)\n",
    "  \n",
    "cell2env_diffusion = {\n",
    "  #\n",
    "}\n",
    "for i, p in enumerate(cell_p[-n_aux_p:]):\n",
    "  cell2env_diffusion[p] = 0.1 / (n_aux_p - i)\n",
    "\n",
    "\n",
    "def simulate(args):\n",
    "  import os\n",
    "  import sys\n",
    "  import re\n",
    "  import time\n",
    "  import copy\n",
    "  import json\n",
    "  from collections import defaultdict\n",
    "  from tqdm import tqdm\n",
    "  from multiprocess import Pool\n",
    "\n",
    "  import numpy as np\n",
    "  import pandas as pd\n",
    "  import scipy.sparse\n",
    "  import scipy.spatial\n",
    "\n",
    "  env_p, env_p_map, cell_p, cell_p_map, resource_cost, cell2cell_diffusion, env2cell_diffusion, cell2env_diffusion, genomes, genes, gene_cell_map, cells, reactions, env, ground_level, steps = args\n",
    "  sim_speed = 0.25\n",
    "  \n",
    "  \n",
    "  def interp_weights(x, y):\n",
    "    xf, xi = np.modf(x)\n",
    "    yf, yi = np.modf(y)\n",
    "    xi = xi.astype(np.int32)\n",
    "    yi = yi.astype(np.int32)\n",
    "    w = np.array([\n",
    "      [xf * yf, xf * (1-yf)],\n",
    "      [(1-xf) * yf, (1-xf) * (1-yf)],\n",
    "    ])\n",
    "    return w, xi, yi\n",
    "\n",
    "  def interp_bilinear(a, x, y):\n",
    "    w, xi, yi = interp_weights(x, y)\n",
    "    w = np.expand_dims(w, -1)\n",
    "    return a[xi, yi] * w[0, 0] \\\n",
    "      + a[xi, yi+1] * w[0, 1] \\\n",
    "      + a[xi+1, yi] * w[1, 0] \\\n",
    "      + a[xi+1, yi+1] * w[1, 1]\n",
    "\n",
    "  def add_interp_bilinear(a, x, y, val):\n",
    "    w, xi, yi = interp_weights(x, y)\n",
    "    w = np.expand_dims(w, -1)\n",
    "  #   print(val.shape)\n",
    "  #   print(a[xi, yi].shape)\n",
    "  #   print(w[0, 0].shape)\n",
    "    a[xi, yi] += w[0, 0] * val\n",
    "    a[xi, yi+1] += w[0, 1] * val\n",
    "    a[xi+1, yi] += w[1, 0] * val\n",
    "    a[xi+1, yi+1] += w[1, 1] * val\n",
    "    return a\n",
    "\n",
    "  def sub_interp_bilinear(a, x, y, val):\n",
    "    w, xi, yi = interp_weights(x, y)\n",
    "    w = np.tile(np.expand_dims(w, -1), (1, val.shape[1]))\n",
    "\n",
    "    # Re-weight by the available resource amounts\n",
    "    w[0, 0] *= a[xi, yi]\n",
    "    w[0, 1] *= a[xi, yi+1]\n",
    "    w[1, 0] *= a[xi+1, yi]\n",
    "    w[1, 1] *= a[xi+1, yi+1]\n",
    "    w /= (np.sum(w, axis=(0,1), keepdims=True) + 1e-6)\n",
    "\n",
    "    a[xi, yi] -= w[0, 0] * val\n",
    "    a[xi, yi+1] -= w[0, 1] * val\n",
    "    a[xi+1, yi] -= w[1, 0] * val\n",
    "    a[xi+1, yi+1] -= w[1, 1] * val\n",
    "    return a\n",
    "\n",
    "\n",
    "  cell_history = []\n",
    "  env_history = []\n",
    "  cell_tree = None\n",
    "  last_neighbour_update = 0\n",
    "  neighbour_update_interval = 4\n",
    "#   cell_steps = [[] for _ in np.arange(cells.shape[0])]\n",
    "#   print('Simulating..')\n",
    "  for step in tqdm(np.arange(steps), disable=True):\n",
    "    # Refill the earth\n",
    "    env[:, :ground_level, env_p_map['aminoacids']] = np.minimum(env[:, :ground_level, env_p_map['aminoacids']] + 0.004, 1.0)\n",
    "\n",
    "    # Sunlight\n",
    "#     env[:, (ground_level-1):, env_p_map['atp']] = np.minimum(env[:, (ground_level-1):, env_p_map['atp']] + 0.1, 1.0)\n",
    "\n",
    "    # Env -> cell diffusion\n",
    "    for pname, v in env2cell_diffusion.items():\n",
    "      diff_amount = np.minimum(\n",
    "        np.minimum(\n",
    "          np.maximum(1.0 - cells[:, cell_p_map[pname]], 0.0), \n",
    "          interp_bilinear(env, cells[:, cell_p_map['x']], cells[:, cell_p_map['y']])[:, env_p_map[pname]]\n",
    "        ),\n",
    "        v\n",
    "      )\n",
    "#       print(diff_amount)\n",
    "#       print(np.sum(np.where(env < 0, env, 0.0)))\n",
    "      cells[:, cell_p_map[pname]] += diff_amount\n",
    "      resource_d = np.zeros((cells.shape[0], env.shape[-1]))\n",
    "      resource_d[:, env_p_map[pname]] = diff_amount\n",
    "      env = sub_interp_bilinear(env, cells[:, cell_p_map['x']], cells[:, cell_p_map['y']], resource_d)\n",
    "        \n",
    "    # Cell -> env diffusion\n",
    "    for pname, v in cell2env_diffusion.items():\n",
    "      diff_amount = np.minimum(\n",
    "        np.minimum(\n",
    "          cells[:, cell_p_map[pname]], \n",
    "          np.maximum(1.0 - interp_bilinear(env, cells[:, cell_p_map['x']], cells[:, cell_p_map['y']])[:, env_p_map[pname]], 0.0)\n",
    "        ),\n",
    "        v\n",
    "      )\n",
    "#       print(diff_amount)\n",
    "#       print(np.sum(np.where(env < 0, env, 0.0)))\n",
    "      cells[:, cell_p_map[pname]] -= diff_amount\n",
    "      resource_d = np.zeros((cells.shape[0], env.shape[-1]))\n",
    "      resource_d[:, env_p_map[pname]] = diff_amount\n",
    "      env = add_interp_bilinear(env, cells[:, cell_p_map['x']], cells[:, cell_p_map['y']], resource_d)\n",
    "\n",
    "    # Express genes\n",
    "    for gene_idx, gene in enumerate(genes):\n",
    "      if gene_idx not in gene_cell_map:\n",
    "        continue\n",
    "        \n",
    "      cell_subset = gene_cell_map[gene_idx]\n",
    "        \n",
    "      def cell_subset_idx(cell_subset, p):\n",
    "        if isinstance(p, list):\n",
    "          return tuple([np.tile(cell_subset, len(p)), np.repeat(p, len(cell_subset))])\n",
    "        return tuple([np.tile(cell_subset, 1), np.repeat([p], len(cell_subset))])\n",
    "        \n",
    "      expression = np.ones((len(cell_subset)))\n",
    "      # Apply activator influence on expression\n",
    "      if len(gene['activators']) > 0:\n",
    "        expression *= 1e-6\n",
    "#         print(cells[cell_subset_idx(cell_subset, gene['activators'])])\n",
    "        expression += np.prod(cells[cell_subset_idx(cell_subset, gene['activators'])], axis=-1)\n",
    "      # Apply repressor influence on expression\n",
    "      if len(gene['repressors']) > 0:\n",
    "#         print(cells[cell_subset_idx(cell_subset, gene['repressors'])])\n",
    "        expression -= np.sum(cells[cell_subset_idx(cell_subset, gene['repressors'])], axis=-1)\n",
    "      expression = np.maximum(expression, 0.0) * sim_speed\n",
    "      # Limit expression by the amount of materials and energy available\n",
    "      expression = np.minimum(expression, cells[cell_subset][:, cell_p_map['aminoacids']] / resource_cost[gene['product']] * sim_speed)\n",
    "      expression = np.minimum(expression, cells[cell_subset][:, cell_p_map['atp']] / resource_cost[gene['product']] * sim_speed)\n",
    "#       print(expression)\n",
    "      cells[cell_subset_idx(cell_subset, gene['product'])] += expression\n",
    "      cells[cell_subset_idx(cell_subset, cell_p_map['aminoacids'])] -= expression * resource_cost[gene['product']]\n",
    "      cells[cell_subset_idx(cell_subset, cell_p_map['atp'])] -= expression * resource_cost[gene['product']]\n",
    "\n",
    "    # Process reactions inside the cells\n",
    "    # TODO: make sure resource amounts do not become negative\n",
    "    for reaction in reactions:\n",
    "#       reaction_power = np.expand_dims(np.minimum(np.min(cells[:, list(reaction['reactants'].keys())] / np.array(list(reaction['reactants'].values())), axis=1), 1.0), -1)\n",
    "      reaction_power = np.expand_dims(\n",
    "        np.minimum(\n",
    "          np.prod(cells[:, list(reaction['reactants'].keys())] / np.array(list(reaction['reactants'].values())), axis=1), \n",
    "          np.min(cells[:, list(reaction['reactants'].keys())] / np.array(list(reaction['reactants'].values())), axis=1)\n",
    "        ), \n",
    "        -1\n",
    "      )\n",
    "      cells[:, list(reaction['reactants'].keys())] = np.maximum(cells[:, list(reaction['reactants'].keys())] - np.array(list(reaction['reactants'].values())) * reaction_power * sim_speed, 0.0)\n",
    "      cells[:, list(reaction['products'].keys())] += np.array(list(reaction['products'].values())) * reaction_power * sim_speed\n",
    "      \n",
    "    # TODO: Variable resource decay?\n",
    "    resource_decay = 0.99\n",
    "    cells[:, cell_p_map['cycle']+1:] *= resource_decay\n",
    "\n",
    "    for ic in np.arange(cells.shape[0]):\n",
    "      # Absorb sunlight\n",
    "      cell_r = 0.5\n",
    "      if cells[ic, cell_p_map['y']] >= ground_level-1 and not np.any((cells[:, cell_p_map['y']] > cells[ic, cell_p_map['y']]) & (np.abs(cells[:, cell_p_map['x']] - cells[ic, cell_p_map['x']]) < cell_r)):\n",
    "        cells[ic, cell_p_map['atp']] += 0.1 * sim_speed\n",
    "    \n",
    "    # Consume resources to stay alive\n",
    "    cells[:, cell_p_map['atp']] -= 0.001 * sim_speed\n",
    "    cells[:, cell_p_map['aminoacids']] -= 0.001 * sim_speed\n",
    "    starving = np.where((cells[:, cell_p_map['atp']] <= 0.0) | (cells[:, cell_p_map['aminoacids']] <= 0.0))\n",
    "#     print(cells.shape)\n",
    "#     print(starving)\n",
    "    had_divisions = False\n",
    "    if len(starving[0]) > 0:\n",
    "      cells = np.delete(cells, starving[0], axis=0)\n",
    "      had_divisions = True\n",
    "      # Rebuild the gene-cell mapping\n",
    "      gene_cell_map = {}\n",
    "      for cell_idx, genome_id in enumerate(cells[:, cell_p_map['genome']]):\n",
    "        for gene_id in genomes[int(genome_id)]['genes']:\n",
    "          if gene_id not in gene_cell_map:\n",
    "            gene_cell_map[gene_id] = []\n",
    "          gene_cell_map[gene_id].append(cell_idx)\n",
    "#     print(cells.shape) \n",
    "        \n",
    "    # Grow up\n",
    "#     cells[np.argmax(cells[:, cell_p_map['y']]), cell_p_map['cycle']] += 0.5 * sim_speed\n",
    "      \n",
    "    next_cells = np.copy(cells)\n",
    "    \n",
    "    for ic in np.arange(cells.shape[0]):\n",
    "      # Divide cells\n",
    "      if cells[ic, cell_p_map['cycle']] >= 1.0:\n",
    "        new_cell = np.zeros((len(cell_p), ))\n",
    "        new_cell[cell_p_map['genome']] = cells[ic, cell_p_map['genome']]\n",
    "        new_cycle = (cells[ic, cell_p_map['cycle']] - 1.0) / 2.0\n",
    "        # Position the new cell randomly around the parent\n",
    "        new_cell[cell_p_map['x']] = np.clip(cells[ic, cell_p_map['x']] + np.clip(np.random.randn(), -3.0, 3.0) / 3.0, 0.5, env.shape[0] - 1.5)\n",
    "        new_cell[cell_p_map['y']] = np.clip(cells[ic, cell_p_map['y']] + np.clip(np.random.randn(), -3.0, 3.0) / 3.0, 0.5, env.shape[1] - 1.5)\n",
    "        # All resources are split semi-equally (up to 10% difference)\n",
    "        resource_diff = 0.05 * np.clip(np.random.randn(*cells[ic, cell_p_map['aminoacids']:].shape) / 3.0, -1.0, 1.0) * cells[ic, cell_p_map['aminoacids']:] / 2.0\n",
    "        new_cell[cell_p_map['aminoacids']:] = cells[ic, cell_p_map['aminoacids']:] / 2.0 + resource_diff\n",
    "        next_cells[ic, cell_p_map['aminoacids']:] = cells[ic, cell_p_map['aminoacids']:] / 2.0 - resource_diff\n",
    "        next_cells[ic, cell_p_map['cycle']] = new_cycle\n",
    "        new_cell[cell_p_map['cycle']] = new_cycle\n",
    "        next_cells = np.concatenate([next_cells, [new_cell]], axis=0)\n",
    "        had_divisions = True\n",
    "#         print('genes')\n",
    "#         print(genomes[int(new_cell[cell_p_map['genome']])]['genes'])\n",
    "        for gene_id in genomes[int(new_cell[cell_p_map['genome']])]['genes']:\n",
    "#           print(gene_id)\n",
    "          gene_cell_map[gene_id].append(next_cells.shape[0] - 1)\n",
    "#         cell_steps.append([])\n",
    "  #       print(cells)\n",
    "\n",
    "#       max_dist = 2.0\n",
    "#       max_neighbours = 6\n",
    "#       neighbour_dist_l1 = np.abs(cells[:, cell_p_map['x']] - cells[ic, cell_p_map['x']]) + np.abs(cells[:, cell_p_map['y']] - cells[ic, cell_p_map['y']])\n",
    "#       neighbour_dist_l1[ic] = sys.float_info.max\n",
    "#       neighbours = list(map(\n",
    "#         lambda x: x[0], \n",
    "#         sorted(\n",
    "#           filter(lambda x: x[1] < max_dist, enumerate(neighbour_dist_l1)), \n",
    "#           key=lambda x: x[1]\n",
    "#         )\n",
    "#       ))[:max_neighbours]\n",
    "# #       neighbour_check = (np.abs(cells[:, cell_p_map['x']] - cells[ic, cell_p_map['x']]) < max_dist) & (np.abs(cells[:, cell_p_map['y']] - cells[ic, cell_p_map['y']]) < max_dist)\n",
    "# #       neighbour_check[ic] = False\n",
    "# #       neighbours = np.where(neighbour_check)[0]\n",
    "\n",
    "#       move_vec = np.zeros((2))\n",
    "#       for ine in neighbours:\n",
    "#         dirv = np.array([cells[ine, cell_p_map['x']] - cells[ic, cell_p_map['x']], cells[ine, cell_p_map['y']] - cells[ic, cell_p_map['y']]])\n",
    "#         dist = np.sqrt(np.power(dirv[0], 2) + np.power(dirv[1], 2))\n",
    "#         diffusion_dist = 1.0\n",
    "#         if dist < diffusion_dist:\n",
    "#           # Cell -> cell diffusion\n",
    "#           for p, v in cell2cell_diffusion.items():\n",
    "#             diff_amount = (cells[ic, p] - cells[ine, p]) / 2.0 * v\n",
    "#   #           diff_amount = min(\n",
    "#   #             np.abs((cells[ic, p] - cells[ine, p]) / 2.0), \n",
    "#   #             v\n",
    "#   #           ) * np.sign(cells[ic, p] - cells[ine, p])\n",
    "#             next_cells[ic, p] = np.maximum(next_cells[ic, p] - diff_amount * sim_speed, 0.0)\n",
    "#             next_cells[ine, p] = np.maximum(next_cells[ine, p] + diff_amount * sim_speed, 0.0)\n",
    "\n",
    "#         # Neighbour attraction/repulsion\n",
    "#         move_vec += dirv * (np.power(max(dist-1.0, 0), 1.5) * 0.2 - 2.0*np.power(1 - min(dist / 1.0, 1), 0.5)) * 0.4\n",
    "\n",
    "#       # Cell movement\n",
    "#       move_vec = np.clip(move_vec * 0.5 * sim_speed, -1.0, 1.0)\n",
    "#       next_cells[ic, cell_p_map['x']] = np.clip(next_cells[ic, cell_p_map['x']] + move_vec[0], 0.5, env.shape[0] - 1.5)\n",
    "#       next_cells[ic, cell_p_map['y']] = np.clip(next_cells[ic, cell_p_map['y']] + move_vec[1], 0.5, env.shape[1] - 1.5)\n",
    "\n",
    "#       cell_steps[ic].append((step, np.copy(next_cells[ic])))\n",
    "    \n",
    "    if cell_tree is None or had_divisions or last_neighbour_update + neighbour_update_interval <= step:\n",
    "      cell_tree = scipy.spatial.KDTree(next_cells[:, cell_p_map['x']:(cell_p_map['y']+1)])\n",
    "      max_dist = 2.0\n",
    "      neighbours = cell_tree.query_pairs(max_dist, p=2)  # cells[:, cell_p_map['x']:cell_p_map['y']], 6\n",
    "      last_neighbour_update = step\n",
    "    \n",
    "    neighbours = np.array(list(neighbours))\n",
    "    \n",
    "    if neighbours.shape[0] > 0:\n",
    "      dirv = np.array([next_cells[neighbours[:, 1], cell_p_map['x']] - next_cells[neighbours[:, 0], cell_p_map['x']], next_cells[neighbours[:, 1], cell_p_map['y']] - next_cells[neighbours[:, 0], cell_p_map['y']]]).T\n",
    "      dist = np.sqrt(np.power(dirv[:, 0], 2) + np.power(dirv[:, 1], 2))\n",
    "      dirv = dirv / (np.linalg.norm(dirv, axis=1, keepdims=True) + 1e-6)\n",
    "      diffusion_dist = 1.0\n",
    "      diffusors = np.argwhere(dist < diffusion_dist)\n",
    "#       for p, v in cell2cell_diffusion.items():\n",
    "#         diff_amount = (cells[ic, p] - cells[ine, p]) / 2.0 * v\n",
    "#       if dist < diffusion_dist:\n",
    "      if diffusors.shape[0] > 0:\n",
    "        # Cell -> cell diffusion\n",
    "        for p, v in cell2cell_diffusion.items():\n",
    "          diff_amount = (next_cells[neighbours[diffusors, 1], p] - next_cells[neighbours[diffusors, 0], p]) / 2.0 * v\n",
    "  #           diff_amount = min(\n",
    "  #             np.abs((cells[ic, p] - cells[ine, p]) / 2.0), \n",
    "  #             v\n",
    "  #           ) * np.sign(cells[ic, p] - cells[ine, p])\n",
    "          np.add.at(next_cells, (neighbours[diffusors, 0], p), diff_amount * sim_speed)\n",
    "          np.add.at(next_cells, (neighbours[diffusors, 1], p), -diff_amount * sim_speed)\n",
    "          next_cells[:, p] = np.maximum(next_cells[:, p], 0.0)\n",
    "\n",
    "#       move_force = 0.1 * (np.power(dist, 0.1) - 0.8 * np.power(1.0 / np.maximum(dist - 0.5, 0.1), 1.5))\n",
    "      move_force = (np.power(np.maximum(dist-1.0, 0), 1.5) * 0.2 - 2.0*np.power(1 - np.minimum(dist / 1.0, 1), 0.5)) * 0.1\n",
    "#       move_force = (np.power(np.maximum(dist-1.0, 0), 1.5) * 0.2 - 2.0*np.power(1 - np.minimum(dist / 1.0, 1), 0.5)) * 0.4\n",
    "    # Neighbour attraction/repulsion\n",
    "#     move_vec += dirv * (np.power(max(dist-1.0, 0), 1.5) * 0.2 - 2.0*np.power(1 - min(dist / 1.0, 1), 0.5)) * 0.4\n",
    "\n",
    "      move_vec = np.zeros((next_cells.shape[0], 2))\n",
    "      np.add.at(move_vec, (neighbours[:, 0], 0), dirv[:, 0] * move_force)\n",
    "      np.add.at(move_vec, (neighbours[:, 0], 1), dirv[:, 1] * move_force)\n",
    "      np.add.at(move_vec, (neighbours[:, 1], 0), -dirv[:, 0] * move_force)\n",
    "      np.add.at(move_vec, (neighbours[:, 1], 1), -dirv[:, 1] * move_force)\n",
    "#       move_vec = move_vec / np.expand_dims(np.maximum(np.bincount(neighbours.ravel(), minlength=move_vec.shape[0]), 1.0).astype(np.float32), -1)\n",
    "      move_vec = np.clip(move_vec, -1.0, 1.0) * sim_speed\n",
    "      move_vec += np.clip(np.random.randn(*move_vec.shape) / 3.0, -1.0, 1.0) * 0.1 * move_vec\n",
    "\n",
    "    # Cell movement\n",
    "#     move_vec = np.clip(move_vec * 0.5 * sim_speed, -1.0, 1.0)\n",
    "      next_cells[:, cell_p_map['x']:(cell_p_map['y']+1)] += move_vec\n",
    "#       next_cells[neighbours[:, 0], cell_p_map['x']] += dirv[:, 0] * move_force * sim_speed\n",
    "#       next_cells[neighbours[:, 0], cell_p_map['y']] += dirv[:, 1] * move_force * sim_speed\n",
    "#       next_cells[neighbours[:, 1], cell_p_map['x']] -= dirv[:, 0] * move_force * sim_speed\n",
    "#       next_cells[neighbours[:, 1], cell_p_map['y']] -= dirv[:, 1] * move_force * sim_speed\n",
    "      next_cells[:, cell_p_map['x']] = np.clip(next_cells[:, cell_p_map['x']], 0.5, env.shape[0] - 1.5)\n",
    "      next_cells[:, cell_p_map['y']] = np.clip(next_cells[:, cell_p_map['y']], 0.5, env.shape[1] - 1.5)\n",
    "\n",
    "    cells = next_cells\n",
    "    if step % 10 == 0:\n",
    "      cell_history.append(np.copy(cells))\n",
    "      env_history.append(np.copy(env))\n",
    "      \n",
    "  return cell_history, env_history\n",
    "\n",
    "  # fix, ax = plt.subplots(len(cell_steps), figsize=(10, 18), sharex=True)\n",
    "  # plot_map = {\n",
    "  #   'cycle': '#44cccc',\n",
    "  #   'aminoacids': '#997744',\n",
    "  #   'atp': '#cccc44',\n",
    "  # }\n",
    "  # for ic in np.arange(len(cell_steps)):\n",
    "  #   for pname, color in plot_map.items():\n",
    "  #     ax[ic].plot(list(map(lambda x: x[0], cell_steps[ic])), list(map(lambda x: x[1][cell_p_map[pname]], cell_steps[ic])), color=color)\n",
    "\n",
    "def render(cell_history, env_history):\n",
    "  def show_frame(cells, env):\n",
    "    fig, ax = plt.subplots(figsize=(6, 2))\n",
    "    ax.fill_between([0, env.shape[0]], 0, ground_level-1, alpha=.4)\n",
    "    ax.plot(0, 0, alpha=0)\n",
    "    ax.plot(env.shape[0], env.shape[1], alpha=0)\n",
    "#     ax.scatter(cells[:, cell_p_map['x']], cells[:, cell_p_map['y']], s=20, c=np.array(list(map(lambda x: '#%02x%02x%02x' % tuple(((np.array([0.2, 0.2, 0.2]) + min(x, 1.0) * np.array([0.0, 0.8, 0.8])) * 250).astype(np.int32)), cells[:, cell_p_map['cycle']]))))\n",
    "    cmap = plt.cm.get_cmap('Set2')\n",
    "    ax.scatter(cells[:, cell_p_map['x']], cells[:, cell_p_map['y']], s=20, c=np.array(list(map(lambda x: '#%02x%02x%02x' % tuple(((np.array([0.0, 0.0, 0.0]) + min(max(x[cell_p_map['cycle']], 0.0), 1.0) * np.array(cmap(x[cell_p_map['genome']] % 20)[:3])) * 250).astype(np.int32)), cells[:]))))\n",
    "#     ax.set_title('Cycle')\n",
    "    f1 = fig2data(fig)\n",
    "    plt.close()\n",
    "\n",
    "    fig2, ax = plt.subplots(figsize=(6, 2))\n",
    "    env_img = np.zeros((env.shape[0], env.shape[1], 3))\n",
    "    env_img[:, ground_level:, :] = 1.0\n",
    "    env_img[:, :ground_level, :] = np.array([0.4, 0.4, 0.0]) + np.minimum(env[:, :ground_level, env_p_map['aminoacids']:env_p_map['aminoacids']+1], 1.0) * np.array([0.6, 0.6, 0.0])\n",
    "    ax.imshow(env_img.transpose((1, 0, 2)), interpolation='bilinear', origin='lower')  # , extent=(50-0.5, (env.shape[0] - 50)-0.5, 50-0.5, (env.shape[1] - 50)-0.5)\n",
    "    ax.plot(0, 0, alpha=0)\n",
    "    ax.plot(env.shape[0], env.shape[1], alpha=0)\n",
    "    ax.scatter(cells[:, cell_p_map['x']], cells[:, cell_p_map['y']], s=20, c=np.array(list(map(lambda x: '#%02x%02x%02x' % tuple(((np.array([0.2, 0.2, 0.2]) + min(x, 1.0) * np.array([0.8, 0.0, 0.8])) * 250).astype(np.int32)), cells[:, cell_p_map['aminoacids']]))))\n",
    "#     ax.set_title('Earth energy')\n",
    "    f2 = fig2data(fig2)\n",
    "    plt.close()\n",
    "\n",
    "    fig3, ax = plt.subplots(figsize=(6, 2))\n",
    "    ax.fill_between([0, env.shape[0]], 0, ground_level-1, alpha=.4)\n",
    "    ax.plot(0, 0, alpha=0)\n",
    "    ax.plot(env.shape[0], env.shape[1], alpha=0)\n",
    "    ax.scatter(cells[:, cell_p_map['x']], cells[:, cell_p_map['y']], s=20, c=np.array(list(map(lambda x: '#%02x%02x%02x' % tuple(((np.array([0.2, 0.2, 0.2]) + min(x, 1.0) * np.array([0.8, 0.8, 0.0])) * 250).astype(np.int32)), cells[:, cell_p_map['atp']]))))\n",
    "#     ax.set_title('Sun energy')\n",
    "    f3 = fig2data(fig3)\n",
    "    plt.close()\n",
    "    \n",
    "    faux = []\n",
    "    for i in np.arange(3):\n",
    "      fig, ax = plt.subplots(figsize=(6, 2))\n",
    "      env_img = np.zeros((env.shape[0], env.shape[1], 3))\n",
    "      env_img = np.array([1.0, 1.0, 1.0]) - np.minimum(env[:, :, (len(env_p) - n_aux_p + i):(len(env_p) - n_aux_p + i)+1], 1.0) * np.array([0.8, 0.8, 0.0])\n",
    "      ax.imshow(env_img.transpose((1, 0, 2)), interpolation='bilinear', origin='lower')\n",
    "#       ax.fill_between([0, env.shape[0]], 0, ground_level-1, alpha=.4)\n",
    "      ax.plot(0, 0, alpha=0)\n",
    "      ax.plot(env.shape[0], env.shape[1], alpha=0)\n",
    "      ax.scatter(cells[:, cell_p_map['x']], cells[:, cell_p_map['y']], s=20, c=np.array(list(map(lambda x: '#%02x%02x%02x' % tuple(((np.array([0.2, 0.2, 0.2]) + min(x, 1.0) * np.array([0.8, 0.8, 0.0])) * 250).astype(np.int32)), cells[:, len(cell_p) - n_aux_p + i]))))\n",
    "  #     ax.set_title('Sun energy')\n",
    "      f = fig2data(fig)\n",
    "      plt.close()\n",
    "      faux.append(f)\n",
    "    \n",
    "    return np.concatenate([\n",
    "      np.concatenate([f1, f2, f3], axis=0),\n",
    "      np.concatenate(faux, axis=0),\n",
    "    ], axis=1)\n",
    "\n",
    "  frames = []\n",
    "  print('Rendering figures..')\n",
    "  for frame_num in tqdm(np.arange(len(cell_history))):\n",
    "      frames.append(show_frame(cell_history[frame_num], env_history[frame_num]))\n",
    "  #     plt.close()\n",
    "\n",
    "  frames = np.array(frames)\n",
    "\n",
    "  # Convert frames to video and send to Visdom\n",
    "  viz_opts = {\n",
    "      'sim_slice': dict(\n",
    "          title='Sim slice',\n",
    "          width=frames.shape[2]*1.5+40,\n",
    "          height=frames.shape[1]*1.5+25\n",
    "      )\n",
    "  }\n",
    "  viz_video(frames, opts=viz_opts['sim_slice'], fps=30)\n",
    "\n",
    "\n",
    "genomes = [\n",
    "  {\n",
    "    'genes': sorted([0]),\n",
    "    'createdTime': 0,\n",
    "    'fitness': [],\n",
    "  }\n",
    "]\n",
    "genome_map = {json.dumps(k['genes']):v for v, k in enumerate(genomes)}\n",
    "\n",
    "genes = [\n",
    "  {\n",
    "    'activators': sorted([\n",
    "      cell_p_map['atp'],\n",
    "      cell_p_map['aminoacids'],\n",
    "    ]),\n",
    "    'repressors': sorted([\n",
    "\n",
    "    ]),\n",
    "    'product': cell_p_map['cycle'],\n",
    "  }\n",
    "]\n",
    "gene_map = {json.dumps(k):v for v, k in enumerate(genes)}\n",
    "\n",
    "reactions = [\n",
    "#     {\n",
    "#       'reactants': {\n",
    "#         cell_p_map['aminoacids']: 1,\n",
    "#         cell_p_map['atp']: 1,\n",
    "#       },\n",
    "#       'products': {\n",
    "#         cell_p_map['cycle']: 1,\n",
    "#       }\n",
    "#     }\n",
    "]\n",
    "\n",
    "\n",
    "def sample_indices(a, p):\n",
    "  if len(a) > 0:\n",
    "    n = np.random.binomial(len(a), p)\n",
    "    if n > 0:\n",
    "      return sorted(np.random.choice(np.arange(len(a)), n, replace=False), reverse=True)\n",
    "  return []\n",
    "\n",
    "\n",
    "def mutate_gene(gene):\n",
    "  new_gene = copy.deepcopy(gene)\n",
    "  add_activator_prob = 0.2\n",
    "  mutate_activator_prob = 0.4\n",
    "  remove_activator_prob = 0.05\n",
    "  add_repressor_prob = 0.1\n",
    "  mutate_repressor_prob = 0.2\n",
    "  remove_repressor_prob = 0.05\n",
    "  mutate_product_prob = 0.3\n",
    "  # Remove/mutate/add activators\n",
    "  for act_idx in sample_indices(new_gene['activators'], remove_activator_prob):\n",
    "    del new_gene['activators'][act_idx]\n",
    "  for act_idx in sample_indices(new_gene['activators'], mutate_activator_prob):\n",
    "    new_gene['activators'][act_idx] = np.random.randint(cell_p_map['cycle'], cells.shape[1])\n",
    "  if np.random.rand() < add_activator_prob:\n",
    "    new_gene['activators'].append(np.random.randint(cell_p_map['cycle'], cells.shape[1]))\n",
    "  new_gene['activators'] = sorted(new_gene['activators'])\n",
    "  # Remove/mutate/add repressors\n",
    "  for rep_idx in sample_indices(new_gene['repressors'], remove_repressor_prob):\n",
    "    del new_gene['repressors'][rep_idx]\n",
    "  for rep_idx in sample_indices(new_gene['repressors'], mutate_repressor_prob):\n",
    "    new_gene['repressors'][rep_idx] = np.random.randint(cell_p_map['cycle'], cells.shape[1])\n",
    "  if np.random.rand() < add_repressor_prob:\n",
    "    new_gene['repressors'].append(np.random.randint(cell_p_map['cycle'], cells.shape[1]))\n",
    "  new_gene['repressors'] = sorted(new_gene['repressors'])\n",
    "  # Mutate the product\n",
    "  if np.random.rand() < mutate_product_prob:\n",
    "    new_gene['product'] = np.random.randint(cell_p_map['cycle'], cells.shape[1])\n",
    "  return new_gene\n",
    "\n",
    "\n",
    "def mutate_genome(genome):\n",
    "  new_genome = copy.deepcopy(genome)\n",
    "  new_genome['fitness'] = []\n",
    "  add_gene_prob = 0.2\n",
    "  mutate_gene_prob = 0.2\n",
    "  remove_gene_prob = 0.05\n",
    "  # Remove some genes\n",
    "  for gene_idx in sample_indices(new_genome['genes'], remove_gene_prob):\n",
    "    del new_genome['genes'][gene_idx]\n",
    "  # Mutate some genes\n",
    "  for gene_idx in sample_indices(new_genome['genes'], mutate_gene_prob):\n",
    "    gene_idx = np.random.randint(0, len(new_genome['genes']))\n",
    "    new_gene = mutate_gene(genes[new_genome['genes'][gene_idx]])\n",
    "    new_gene_id = 0\n",
    "    # Check whether the mutate version already exists\n",
    "    if json.dumps(new_gene) not in gene_map:\n",
    "      genes.append(new_gene)\n",
    "      new_gene_id = len(genes) - 1\n",
    "      gene_map[json.dumps(new_gene)] = new_gene_id\n",
    "    else:\n",
    "      new_gene_id = gene_map[json.dumps(new_gene)]\n",
    "    new_genome['genes'][gene_idx] = new_gene_id\n",
    "    new_genome['genes'] = sorted(new_genome['genes'])\n",
    "  # Add some genes\n",
    "  if np.random.rand() < add_gene_prob:\n",
    "    new_genome['genes'].append(np.random.randint(0, len(genes)))\n",
    "    new_genome['genes'] = sorted(new_genome['genes'])\n",
    "  return new_genome\n",
    "\n",
    "\n",
    "n_threads = 1\n",
    "pool = Pool(n_threads)\n",
    "epochs = 1\n",
    "max_population = 1024\n",
    "for epoch in np.arange(epochs):\n",
    "  # Decide on the parent subpopulation\n",
    "  #sorted_genomes = sorted(genomes, key=lambda x: np.mean(x['fitness'][-3:]) if len(x['fitness']) > 0 else 0, reverse=True)\n",
    "  sorted_genomes = sorted(genomes, key=lambda x: x['createdTime'], reverse=True)[:max_population]\n",
    "  genome_weights = np.arange(len(sorted_genomes))[::-1].astype(np.float32) + 1\n",
    "  genome_weights = genome_weights / np.sum(genome_weights)\n",
    "#   elite_genomes = elite_genomes[:max(int(len(elite_genomes) * 0.2), 5)]\n",
    "  \n",
    "  sim_args = []\n",
    "  start_genome_list = []\n",
    "  for _ in np.arange(n_threads):\n",
    "    # Setup the environment\n",
    "    env_size = (18,6)\n",
    "    env = np.zeros((env_size[0], env_size[1], len(env_p))).astype(np.float32)\n",
    "    ground_level = env_size[1] // 3\n",
    "    # Earth starts full of nutrients\n",
    "    env[:, :ground_level, env_p_map['aminoacids']] = 1.0\n",
    "\n",
    "    # Setup the initial cells\n",
    "    batch_size = 4\n",
    "    cells = np.zeros((batch_size, len(cell_p))).astype(np.float32)\n",
    "    # Sample new genomes\n",
    "    parent_batch = np.random.choice(sorted_genomes, size=batch_size, p=genome_weights)\n",
    "    for i in np.arange(batch_size):\n",
    "      new_genome = parent_batch[i]\n",
    "      # Mutate half of the fresh batch\n",
    "#       if i < batch_size / 2:\n",
    "      new_genome = mutate_genome(new_genome)\n",
    "      new_genome['createdTime'] = int(epoch)\n",
    "  #     print(json.dumps(new_genome))\n",
    "      new_genome_id = 0\n",
    "      # Check whether the mutate version already exists\n",
    "      if json.dumps(new_genome['genes']) not in genome_map:\n",
    "        genomes.append(new_genome)\n",
    "        new_genome_id = len(genomes) - 1\n",
    "        genome_map[json.dumps(new_genome['genes'])] = new_genome_id\n",
    "      else:\n",
    "        new_genome_id = genome_map[json.dumps(new_genome['genes'])]\n",
    "      cells[i, cell_p_map['genome']] = new_genome_id\n",
    "    cells[:, cell_p_map['x']] = np.linspace(0.0, env.shape[0], cells.shape[0], endpoint=False) + float(env.shape[0]) / cells.shape[0] * 0.5\n",
    "    cells[:, cell_p_map['y']] = ground_level - 2\n",
    "    cells[:, cell_p_map['atp']] = 30.0\n",
    "    cells[:, cell_p_map['aminoacids']] = 10.0\n",
    "\n",
    "    gene_cell_map = {}\n",
    "    for cell_idx, genome_id in enumerate(cells[:, cell_p_map['genome']]):\n",
    "      for gene_id in genomes[int(genome_id)]['genes']:\n",
    "        if gene_id not in gene_cell_map:\n",
    "          gene_cell_map[gene_id] = []\n",
    "        gene_cell_map[gene_id].append(cell_idx)\n",
    "\n",
    "    start_genomes = np.bincount(cells[:, cell_p_map['genome']].astype(np.int32))\n",
    "    start_genome_list.append(start_genomes)\n",
    "\n",
    "    n_steps = 150\n",
    "#     simulate((env_p, env_p_map, cell_p, cell_p_map, resource_cost, cell2cell_diffusion, env2cell_diffusion, cell2env_diffusion, genomes, genes, gene_cell_map, cells, reactions, env, ground_level, n_steps))\n",
    "    sim_args.append((env_p, env_p_map, cell_p, cell_p_map, resource_cost, cell2cell_diffusion, env2cell_diffusion, cell2env_diffusion, genomes, genes, gene_cell_map, cells, reactions, env, ground_level, n_steps))\n",
    "  \n",
    "  result = pool.map_async(simulate, sim_args)\n",
    "  all_genome_scores = []\n",
    "  for start_genomes, (cell_history, env_history) in zip(start_genome_list, result.get(timeout=600)):\n",
    "  #   cell_history, env_history = simulate(genomes, genes, gene_cell_map, cells, reactions, env, ground_level, n_steps)\n",
    "    genome_scores = np.bincount(cell_history[-1][:, cell_p_map['genome']].astype(np.int32)) / np.maximum(start_genomes.astype(np.float32), 1.0)\n",
    "    for i, s in enumerate(genome_scores):\n",
    "      if s > 0:\n",
    "        genomes[i]['fitness'].append(float(s))\n",
    "    all_genome_scores.append(genome_scores)\n",
    "  print('epoch {0}, max fitness {1}'.format(epoch, np.max(np.concatenate(all_genome_scores, axis=0))))\n",
    "#   if epoch == epochs-1:\n",
    "#     render(cell_history, env_history)\\\n",
    "#     plt.hist(np.bincount(cell_history[-1][:, cell_p_map['genome']].astype(np.int32)), bins=len(genomes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAAD8CAYAAAB6tolUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF9NJREFUeJzt3W+sZGd9H/DvrzaQBtLarhfk2G6XpFsap1IMWjluqSISGv5GNZFCa6sFixJtX5gWIqpq4Q3pC0tulUCLmlpygoOpEAQBLVZtkbouEkUqkDVQY+MiVuDixa69qQmgIkENv764Z/Hl7v0z9975c2bm85GuZuaZMzPPOec5z3m+Z86cW90dAAAAWLQ/t+gKAAAAQCKgAgAAMBICKgAAAKMgoAIAADAKAioAAACjIKACAAAwCgIqAAAAoyCgAgAAMAoCKgAAAKNw4aIrkCSXXnppHz16dNHVAAAAYAbuu+++P+3uI3tNN4qAevTo0Zw6dWrR1QAAAGAGqup/TTKdU3wBAAAYBQEVAACAURBQAQAAGAUBFQAAgFEQUAEAABgFARUAAIBREFABAAAYBQEVAACAURBQAQAAGAUBFQAApujoybsWXQVYWgIqAAAAoyCgAgAAMAoCKgAAAKMgoAIAADAKAioAAACjIKACAAAwCgIqAAAAoyCgAgAAMAoCKgAAAKMgoAIAADAKAioAAACjsGdAraorq+oTVfVQVT1YVW8eyn+7qr5RVV8Y/l616TVvq6rTVfXlqnr5LGcAAACA1XDhBNM8leSt3f25qvqpJPdV1T3Dc+/q7t/ZPHFVXZXk+iQ/n+Snk/yXqvpr3f2DaVYcAACA1bLnN6jd/Vh3f264/50kDyW5fJeXXJfkg939ve7+WpLTSa6ZRmUBAABYXfv6DWpVHU3ywiSfGYreVFX3V9XtVXXxUHZ5kkc2vexMdg+0AAAAMHlArarnJPlIkrd097eT3JrkZ5NcneSxJL97btJtXt7bvN+JqjpVVafOnj2774oDAACwWiYKqFX1jGyE0/d390eTpLsf7+4fdPcPk/x+nj6N90ySKze9/Iokj259z+6+rbuPd/fxI0eOHGYeAAAAWAGTXMW3krwnyUPd/c5N5ZdtmuzXkzww3L8zyfVV9ayqen6SY0k+O70qAwAAsIomuYrvi5O8LskXq+oLQ9nbk9xQVVdn4/Tdh5P84yTp7ger6kNJvpSNKwDf5Aq+AAAA7GXPgNrdn8r2vyu9e5fX3Jzk5kPUCwAAgDWzr6v4AgAAwKwIqAAAAIyCgAoAAMAoCKgAAACMgoAKAADAKAioAAAAjIKACgAAwCgIqMBKO3ryrkVXAQCACQmoAAAAjIKACgAAwCgIqAAAAIyCgAoAAMAoCKgAAACMgoAKAADAKAioAAAAjIKACgAAa8L/B2fsBFQAAABGQUAFAABgFARUAAAARkFABQAAYBQEVAAAAEZBQAUAAGAUBFQAAABGQUAFAABgFARUlpZ/NA0AAKtFQAUAAGAUBFQAAABGQUAFAABgFARUAAAARmHPgFpVV1bVJ6rqoap6sKrePJRfUlX3VNVXhtuLh/KqqndX1emqur+qXjTrmQAAAGD5TfIN6lNJ3trdP5fk2iQ3VdVVSU4mube7jyW5d3icJK9Mcmz4O5Hk1qnXGgAAgJWzZ0Dt7se6+3PD/e8keSjJ5UmuS3LHMNkdSV4z3L8uyft6w6eTXFRVl0295sDa8C+FAADWw75+g1pVR5O8MMlnkjyvux9LNkJskucOk12e5JFNLzszlAEAAMCOJg6oVfWcJB9J8pbu/vZuk25T1tu834mqOlVVp86ePTtpNQAAAFhREwXUqnpGNsLp+7v7o0Px4+dO3R1unxjKzyS5ctPLr0jy6Nb37O7buvt4dx8/cuTIQesPAADAipjkKr6V5D1JHurud2566s4kNw73b0zysU3lrx+u5nttkm+dOxUYAAAAdnLhBNO8OMnrknyxqr4wlL09yS1JPlRVb0zy9SSvHZ67O8mrkpxO8t0kb5hqjQEAAFhJewbU7v5Utv9daZK8dJvpO8lNh6wXAAAAa2ZfV/EFAACAWRFQAQAAGAUBFQAAgFEQUAEAABgFARUA4ACOnrxr0VUAWDkCKgAAAKMgoAIAADPljAMmJaACAAAwCgIq++YIGAAAMAsCKgCwFBwgBVh9AioAwJwI2QC7E1ABAAAYBQEVAACAURBQAQAAGAUBFQBgjvwOFWBnAioAAACjIKACAKwI384Cy05ABQCALYR9WAwBFQAAgFEQUAFgTfhGCGAc9Mc7E1CBXelAAQCYFwEVAACAURBQAQAAGAUBFYCV4HR0WF22b1gfAiqHYocBAAB7M26ejIAKAMDSMuiH1SKgAjBTBo8AjJn91LgIqAAAAPsk2M6GgAqwAuwkAYBVIKACa0WQAwAYrz0DalXdXlVPVNUDm8p+u6q+UVVfGP5etem5t1XV6ar6clW9fFYVB4B5c4ADAGZrkm9Q35vkFduUv6u7rx7+7k6SqroqyfVJfn54zb+rqgumVVkAYHkse6Bf9vrDsrHNkUwQULv7k0menPD9rkvywe7+Xnd/LcnpJNccon4AADBzwhGMw2F+g/qmqrp/OAX44qHs8iSPbJrmzFB2nqo6UVWnqurU2bNnD1ENAAAAVsFBA+qtSX42ydVJHkvyu0N5bTNtb/cG3X1bdx/v7uNHjhw5YDVWiyN3TNO6t6d1n38AgGV0oIDa3Y939w+6+4dJfj9Pn8Z7JsmVmya9Ismjh6siMA/TCHRCIQAHZR8CJAcMqFV12aaHv57k3BV+70xyfVU9q6qen+RYks8erooAAACsg0n+zcwHkvz3JC+oqjNV9cYk/6qqvlhV9yf55SS/lSTd/WCSDyX5UpKPJ7mpu38ws9oDAADswrfzy+XCvSbo7hu2KX7PLtPfnOTmw1QKAID1IUAA5xzmKr4AAAAwNQLqinIkEjgs/QjA6tLHM1YCKgBzZVAEAOxEQGXUDGThYGw7AKwC+7P1I6ACAAAwCgLqITmqAwDsZJ3GCes0r8DsCKgjo3OfDstxeUxzXVnvG5Z1OUxa72WdP1hVtklgmgRUmDE7bhZJ+wMAlomAug8GesvLumPstFGWmfYLwLQIqLDCDBrHb5XX0SrPG3A4+gdgJwIqAAAAoyCgAsAC+SYJYNz00/MloC6Qxg4ALBvjF2CWBFRgqRgYAQCsLgF1AcY+wD5Xv4PWc+zzBwAsnvECsB0BdUp0spOxnAAAgJ0IqEzksN+qAgAA7EVABQAAYBQEVAAAAEZBQCWJU3e3sjzWw1jX81jrtUwsQwBYTgIqAAArbdoHrRwEg9kRUFk7dirApPQXADBfAiowdQb127NcAAB2J6ACAAAwCgIqAPvm2+DDswwB4HwCKoyIASsAAOtMQAVgaTiIA8ttHtvwvPsJ/dLhWH5sJaDCmlrEDsFOCGB56cOBeRBQAQAAGIU9A2pV3V5VT1TVA5vKLqmqe6rqK8PtxUN5VdW7q+p0Vd1fVS+aZeVhnYz1yPVY6wUAzI/xANMyyTeo703yii1lJ5Pc293Hktw7PE6SVyY5NvydSHLrdKq5vGys57NMWGXaNwDAwe0ZULv7k0me3FJ8XZI7hvt3JHnNpvL39YZPJ7moqi6bVmVZXwb9zIN2xjKbZ/tdlW1lVeYDYJUc9Deoz+vux5JkuH3uUH55kkc2TXdmKDtPVZ2oqlNVders2bMHrMZqsIME2Jk+EpbHtLZX2z2sr2lfJKm2KevtJuzu27r7eHcfP3LkyJSrAQAsgmABjJk+avwOGlAfP3fq7nD7xFB+JsmVm6a7IsmjB6/ecrMBLAfrCQBYtEWNR4yDGJuDBtQ7k9w43L8xycc2lb9+uJrvtUm+de5U4FVkgz6fZQKwvWXtH8de71nXb+zzv0iWDTALF+41QVV9IMlLklxaVWeSvCPJLUk+VFVvTPL1JK8dJr87yauSnE7y3SRvmEGdAQAAWEF7BtTuvmGHp166zbSd5KbDVmpZHT15Vx6+5dWLrgawIvQpAMC6mfZFkhgxp+LA02wP42cdAXBQ9iHLS0AFAABgFARUAADWkm/ZYHwEVABWkoEnq0abBtaBgAprzoAHAGD5rOoYTkAduaMn71rZxgfAcrN/YhG0O1htAiqsADtrABbB/geYNgEVAGDJCYrAqhBQAWBJCCGLZx0AzJaAegD73TmNcWfmt62MlXYJALC+BFTgPELiOFkvrJvNbX4VDg7DLGnzy8u6+3ECKgCwFgwCAcZPQAUAAGAUBFQAAABGQUCFGXAaGatot3atzQMA0yCgAgAAzJmDu9sTUAFWnB0gALAsBNQlZtA5LtYHsA70dQDMkoAKwJ6EkvVifQOwKAIqK80gCwCAVbWKY10BFVgLq9iBAwDns89fbgIqMFV2CgAAHJSACsChOCgBB2PbATifgLrm7Bz3Zz/Ly7KFg5l027GNActGvwV7E1ABAFg6yxb2lq2+sCgCKgBMiQEobM+2AUxKQAUAGAhSAIsloGJnDAB7sK8EmA8BFQAAgFG48DAvrqqHk3wnyQ+SPNXdx6vqkiR/lORokoeT/L3u/ubhqgnAWPgmCQCYlWl8g/rL3X11dx8fHp9Mcm93H0ty7/AYOABBAABgMYzDFmMWp/hel+SO4f4dSV4zg88AAAB2IWCxjA4bUDvJf66q+6rqxFD2vO5+LEmG2+ce8jNgZdlxAADLwJiFeTlsQH1xd78oySuT3FRVvzTpC6vqRFWdqqpTZ8+ePWQ1YD3ZWcBysu0CwPYOFVC7+9Hh9okk/yHJNUker6rLkmS4fWKH197W3ce7+/iRI0cOU421sC6DmXWZTwAA4HwHDqhV9eyq+qlz95O8LMkDSe5McuMw2Y1JPnbYSnJwAh8AcBjGEuthFdbzKswDh/sG9XlJPlVV/yPJZ5Pc1d0fT3JLkl+tqq8k+dXhMSylMXR0m+swhvqweNoBzIZt62mWBYeh/XAYB/4/qN391SS/sE35/0ny0sNUCgAAgPUzi38zA0yBo4+sE+0dAEgEVGbMoBPgafpEFkG7A5aJgAoAAMAoCKj8iCOsLIq2B7BBfwisOwEVNjEwAACAxRFQAWAEJj1ANssDaQ7SAbBoAioAcycI7c7yYSfaBkxukduLbfXgBFRg6c1zJ2CHs3qs09Vm/QIsFwEVAHaxbAFn2eq7aJbXcrLeWGba7+4EVABgVwZT+2N5wTgtYtvUH+yfgMpSsrED03SYPkV/BADTI6AyWgZ9MBu2LYCd+0J9JGO2Du1TQF2QdWhc68K6BGCM7J/Wl3XPMhNQR0inArNh2wI4HP3o7iwfODwBdcnt1hHqJAGYlH0GAGMgoMIOtg7WDN4AWAT7H9iebWM1CajAj9HZj886r5N1nncAWEcCKiyAQTfLSttlLLRFgNUkoM7RtHamLosO2N73x/Jivw7bZhbZ5rT35TX2dTeN+o19Hlk8ARWAtWWgNB7WxeJY9sCYCKjMnB0fbM+2AYth22NZaKusIwGVlaETB9ibvhJWezvYz7xNezms8nKdt3VelgIqDNa5I4CxsB1yjrYAsJ4EVFghBnSwfmz3AKwSARVYOgbkjJF2+bTNy8LVbDfsVpcx1ZP5WoZ1vwx1nBfLYj4EVADWigEGAIyXgAoAMyIMw2qwLcP8CKjA6BgIAGzQHy4H62l5WXfjI6CuGBsZwHzobwFg+mYWUKvqFVX15ao6XVUnZ/U5TGa/AykDL9aNNs+ktBUAmJ2ZBNSquiDJ7yV5ZZKrktxQVVfN4rPYn1UfWC3yn1P73Ol+9jTqu+hlverbG9Mzy7ayju1wGeZ5WnUc+7yOqX5jqstelqmu62brPt4+f/pm9Q3qNUlOd/dXu/v7ST6Y5LoZfRYsxLp3RIed/3VfftN00GW5iHWw9TPH2o7GckBnXgOfSd5/Wutqlbf9VZy3ac7TpO+1SstxbP9eaNqfOY++Y1qWtf0dPXnXnnUaW50P68IZve/lSR7Z9PhMkl+c0WfN1V472KMn78rDt7x634Ow7aZ/+JZX7znddtPOqvPZOmg695nnyrfWd7dlsPn+Xstr8+dPsmwn/dztHk9atttRs0nfc+tzm5fndut+r/faqQ77mefd5nPzst+67nd77W7T7tSmJqnbbmW7Lb/d6rrb+ybnz8t27XGnz9q6fvean2m11UnLt67fret+p9du7Xu2ztuk7zOJSZbPdp+z3XMHXU6TlO+3D9jJXm1/63raOu1++uCdPn+Svmiv/m+Safe7Xnbrb/bqFw6zfvbqK/d63V7P7ba/22+fMMky2m5bn2R97jUfe712mtMeZBlNsq856HrZ6b226yt3e4+95u2cg+yHJvnM/ba/vfr97abfaV8zyfLYrv/bj52WyXZjsYPsa7e+dpI6TrIeJx0rLrvq7um/adVrk7y8u39zePy6JNd09z/ZNM2JJCeGhy9I8uWpV2S6Lk3yp4uuBGtFm2PetDnmTZtj3rQ55kl7+3F/pbuP7DXRrL5BPZPkyk2Pr0jy6OYJuvu2JLfN6POnrqpOdffxRdeD9aHNMW/aHPOmzTFv2hzzpL0dzKx+g/onSY5V1fOr6plJrk9y54w+CwAAgBUwk29Qu/upqnpTkj9OckGS27v7wVl8FgAAAKthVqf4prvvTnL3rN5/AZbmdGRWhjbHvGlzzJs2x7xpc8yT9nYAM7lIEgAAAOzXrH6DCgAAAPsioE6gql5RVV+uqtNVdXLR9WH5VdWVVfWJqnqoqh6sqjcP5ZdU1T1V9ZXh9uKhvKrq3UMbvL+qXrTYOWBZVdUFVfX5qvpPw+PnV9Vnhjb3R8OF7VJVzxoenx6eP7rIerOcquqiqvpwVf3Pob/7m/o5ZqmqfmvYrz5QVR+oqp/QzzFNVXV7VT1RVQ9sKtt3v1ZVNw7Tf6WqblzEvIyVgLqHqrogye8leWWSq5LcUFVXLbZWrICnkry1u38uybVJbhra1ckk93b3sST3Do+TjfZ3bPg7keTW+VeZFfHmJA9tevwvk7xraHPfTPLGofyNSb7Z3X81ybuG6WC//k2Sj3f3X0/yC9loe/o5ZqKqLk/yT5Mc7+6/kY0LdV4f/RzT9d4kr9hStq9+raouSfKOJL+Y5Jok7zgXahFQJ3FNktPd/dXu/n6SDya5bsF1Ysl192Pd/bnh/neyMWi7PBtt645hsjuSvGa4f12S9/WGTye5qKoum3O1WXJVdUWSVyf5g+FxJfmVJB8eJtna5s61xQ8neekwPUykqv5Ckl9K8p4k6e7vd/efRT/HbF2Y5M9X1YVJfjLJY9HPMUXd/ckkT24p3m+/9vIk93T3k939zST35PzQu7YE1L1dnuSRTY/PDGUwFcMpRS9M8pkkz+vux5KNEJvkucNk2iHT8K+T/PMkPxwe/6Ukf9bdTw2PN7erH7W54flvDdPDpH4mydkkfzicVv4HVfXs6OeYke7+RpLfSfL1bATTbyW5L/o5Zm+//Zr+bhcC6t62O5Lm0sdMRVU9J8lHkrylu7+926TblGmHTKyqfi3JE9193+bibSbtCZ6DSVyY5EVJbu3uFyb5v3n6tLftaHMcynCK5HVJnp/kp5M8OxunWG6ln2Nedmpj2t4uBNS9nUly5abHVyR5dEF1YYVU1TOyEU7f390fHYofP3dK23D7xFCuHXJYL07yd6vq4Wz8VOFXsvGN6kXDqXDJj7erH7W54fm/mPNPaYLdnElyprs/Mzz+cDYCq36OWfk7Sb7W3We7+/8l+WiSvxX9HLO3335Nf7cLAXVvf5Lk2HAFuGdm48f2dy64Tiy54Tcu70nyUHe/c9NTdyY5dyW3G5N8bFP564erwV2b5FvnTiWBSXT327r7iu4+mo1+7L929z9I8okkvzFMtrXNnWuLvzFM7+guE+vu/53kkap6wVD00iRfin6O2fl6kmur6ieH/ey5NqefY9b226/9cZKXVdXFwzf/LxvKSFK2w71V1auy8U3DBUlu7+6bF1wlllxV/e0k/y3JF/P07wHfno3foX4oyV/Oxo72td395LCj/bfZ+AH9d5O8obtPzb3irISqekmSf9bdv1ZVP5ONb1QvSfL5JP+wu79XVT+R5N9n4/fRTya5vru/uqg6s5yq6upsXJTrmUm+muQN2Tg4rp9jJqrqXyT5+9m4Wv7nk/xmNn7bp59jKqrqA0lekuTSJI9n42q8/zH77Neq6h9lY+yXJDd39x/Ocz7GTEAFAABgFJziCwAAwCgIqAAAAIyCgAoAAMAoCKgAAACMgoAKAADAKAioAAAAjIKACgAAwCgIqAAAAIzC/wd4hqo7YeqncQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# print(json.dumps(genomes, indent=2))\n",
    "# print(genomes)\n",
    "# print(genes)\n",
    "fig, ax = plt.subplots(figsize=(16, 4))\n",
    "_ = ax.bar(np.arange(len(genomes)), list(map(lambda x: np.mean(x['fitness'][-3:]) if len(x['fitness']) > 0 else 0, genomes)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rendering figures..\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████████████████████████████████████████████████| 150/150 [00:24<00:00,  5.15it/s]\n"
     ]
    }
   ],
   "source": [
    "# Setup the environment\n",
    "env_size = (54,18)\n",
    "env = np.zeros((env_size[0], env_size[1], len(env_p))).astype(np.float32)\n",
    "ground_level = env_size[1] // 3\n",
    "# Earth starts full of nutrients\n",
    "env[:, :ground_level, env_p_map['aminoacids']] = 1.0\n",
    "\n",
    "# Setup the initial cells\n",
    "batch_size = 4\n",
    "cells = np.zeros((batch_size, len(cell_p))).astype(np.float32)\n",
    "best_genomes = np.argsort(list(map(lambda x: np.mean(x['fitness'][-3:]) if len(x['fitness']) > 2 else 0, genomes)))[::-1]\n",
    "best_genomes = [0] * batch_size\n",
    "# Sample new genomes\n",
    "for i in np.arange(batch_size):\n",
    "  cells[i, cell_p_map['genome']] = best_genomes[i]\n",
    "cells[:, cell_p_map['x']] = np.linspace(0.0, env.shape[0], cells.shape[0], endpoint=False) + float(env.shape[0]) / cells.shape[0] * 0.5\n",
    "cells[:, cell_p_map['y']] = ground_level - 2\n",
    "cells[:, cell_p_map['atp']] = 30.0\n",
    "\n",
    "gene_cell_map = {}\n",
    "for cell_idx, genome_id in enumerate(cells[:, cell_p_map['genome']]):\n",
    "  for gene_id in genomes[int(genome_id)]['genes']:\n",
    "    if gene_id not in gene_cell_map:\n",
    "      gene_cell_map[gene_id] = []\n",
    "    gene_cell_map[gene_id].append(cell_idx)\n",
    "\n",
    "n_steps = 1500\n",
    "sim_args = []\n",
    "sim_args.append((env_p, env_p_map, cell_p, cell_p_map, resource_cost, cell2cell_diffusion, env2cell_diffusion, cell2env_diffusion, genomes, genes, gene_cell_map, cells, reactions, env, ground_level, n_steps))\n",
    "\n",
    "result = pool.map_async(simulate, sim_args)\n",
    "cell_history, env_history = result.get(timeout=600)[0]\n",
    "\n",
    "render(cell_history, env_history)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15\n",
      "{'genes': [0, 0, 1, 9, 41], 'createdTime': 8, 'fitness': [153.0, 114.0, 117.0, 191.0, 177.0]}\n",
      "[{'activators': [], 'repressors': [], 'product': 5}, {'activators': [], 'repressors': [], 'product': 5}, {'activators': [], 'repressors': [6], 'product': 6}, {'activators': [8], 'repressors': [], 'product': 5}, {'activators': [8], 'repressors': [], 'product': 6}]\n",
      "\n",
      "{'genes': [0, 1, 11, 16], 'createdTime': 4, 'fitness': [140.0, 94.0, 146.0, 53.0, 148.0, 115.0, 183.0]}\n",
      "[{'activators': [], 'repressors': [], 'product': 5}, {'activators': [], 'repressors': [6], 'product': 6}, {'activators': [], 'repressors': [5], 'product': 8}, {'activators': [6], 'repressors': [], 'product': 5}]\n",
      "\n",
      "{'genes': [0, 1, 35], 'createdTime': 19, 'fitness': [182.0, 170.0, 92.0]}\n",
      "[{'activators': [], 'repressors': [], 'product': 5}, {'activators': [], 'repressors': [6], 'product': 6}, {'activators': [], 'repressors': [5], 'product': 6}]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# print(len(np.unique([plt.cm.get_cmap('Set2')(x) for x in np.arange(21)])))\n",
    "\n",
    "for i in np.arange(3):\n",
    "  print(genomes[best_genomes[i]])\n",
    "  print([genes[g] for g in genomes[best_genomes[i]]['genes']])\n",
    "  print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The line_profiler extension is already loaded. To reload it, use:\n",
      "  %reload_ext line_profiler\n"
     ]
    }
   ],
   "source": [
    "%load_ext line_profiler\n",
    "\n",
    "n_steps = 1000\n",
    "%lprun -f simulate simulate((env_p, env_p_map, cell_p, cell_p_map, resource_cost, cell2cell_diffusion, env2cell_diffusion, cell2env_diffusion, genomes, genes, gene_cell_map, cells, reactions, env, ground_level, n_steps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "66"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "frames.shape\n",
    "len(cells)"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}