simulation_scaling_demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "naughty-knock",
   "metadata": {},
   "outputs": [],
   "source": [
    "from simulations import dropletSimulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "several-combination",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "abroad-mission",
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_frequency_vector(length):\n",
    "    unnormalized = np.random.randint(low=0, high=370, size=length)\n",
    "    return unnormalized / sum(unnormalized)\n",
    "\n",
    "def random_copy_numbers(length):\n",
    "    return np.random.randint(low=1, high=16, size=length)\n",
    "\n",
    "A = 2*(np.random.random((100,100)) - 0.5)\n",
    "beta = np.random.random(100)\n",
    "\n",
    "simulation = dropletSimulation(number_species=100, number_droplets=1000, number_batches=5, \n",
    "             copy_numbers=random_copy_numbers(100), frequency_vector=random_frequency_vector(100),\n",
    "             glv_interaction_coefficients=A, glv_baserate_coefficients=beta, noise_scale=8, seed=42,\n",
    "                              timestep=0.01, batch_window=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "continental-teach",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 909 ms, sys: 380 ms, total: 1.29 s\n",
      "Wall time: 1.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "interested-contrary",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation2 = dropletSimulation(number_species=100, number_droplets=1000, number_batches=5, \n",
    "             copy_numbers=random_copy_numbers(100), frequency_vector=random_frequency_vector(100),\n",
    "             glv_interaction_coefficients=A, glv_baserate_coefficients=beta, noise_scale=8,\n",
    "                               timestep = 0.001, batch_window=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "occasional-seventh",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 868 ms, sys: 378 ms, total: 1.25 s\n",
      "Wall time: 3.97 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation2.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "structural-messaging",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation3 = dropletSimulation(number_species=100, number_droplets=1000, number_batches=5, \n",
    "             copy_numbers=random_copy_numbers(100), frequency_vector=random_frequency_vector(100),\n",
    "             glv_interaction_coefficients=A, glv_baserate_coefficients=beta, noise_scale=8,\n",
    "                               timestep = 0.0001, batch_window=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "artificial-digit",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 897 ms, sys: 399 ms, total: 1.3 s\n",
      "Wall time: 33.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation3.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "conceptual-buddy",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "624168203877.475"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation.cells.counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "sunrise-forum",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "23712.272496690533"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation2.cells.counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "comfortable-pilot",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation3.cells.counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "essential-petersburg",
   "metadata": {},
   "source": [
    "Dang OK so nice the carrying capacity limit works a lot better with smaller time steps yay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "revolutionary-travel",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15412.548733897793"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation3.reads.counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "polar-functionality",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation4 = dropletSimulation(number_species=10, number_droplets=1000, number_batches=5, \n",
    "             copy_numbers=random_copy_numbers(10), frequency_vector=random_frequency_vector(10),\n",
    "             glv_interaction_coefficients=A[:,0:10][0:10,:], glv_baserate_coefficients=beta[0:10], noise_scale=8,\n",
    "                               timestep = 0.0001, batch_window=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "neither-construction",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 818 ms, sys: 328 ms, total: 1.15 s\n",
      "Wall time: 32.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation4.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "quantitative-privacy",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000.00000004811"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation4.cells.counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "incorporated-health",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "13234.261582990615"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation4.reads.counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "neither-threshold",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation5 = dropletSimulation(number_species=1000, number_droplets=1000, number_batches=5, \n",
    "             copy_numbers=random_copy_numbers(1000), frequency_vector=random_frequency_vector(1000),\n",
    "             glv_interaction_coefficients=2*(np.random.random((1000,1000))-0.5), glv_baserate_coefficients=np.random.random(1000), \n",
    "                                noise_scale=8,timestep = 0.0001, batch_window=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "embedded-henry",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.88 s, sys: 1.56 s, total: 3.44 s\n",
      "Wall time: 32.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation5.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "continued-thesis",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000.0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation5.cells.counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "running-soccer",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12272.134271466313"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation5.reads.counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "superior-rabbit",
   "metadata": {},
   "source": [
    "HOLY COW! the time is actually nearly independent of the number of strains/species now! yay!\n",
    "\n",
    "I might not have to use any more dumb hacks! yay!\n",
    "\n",
    "oh well I forgot about the fact that the coefficients will still be more expensive to compute though (you know if there are one million of them versus 100 but anyway)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "saved-youth",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation6 = dropletSimulation(number_species=1000, number_droplets=1000, number_batches=5, \n",
    "             copy_numbers=random_copy_numbers(1000), frequency_vector=random_frequency_vector(1000),\n",
    "             glv_interaction_coefficients=2*(np.random.random((1000,1000))-0.5), glv_baserate_coefficients=np.random.random(1000), \n",
    "                                noise_scale=8,timestep = 0.0002, batch_window=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "surprised-quebec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.95 s, sys: 1.55 s, total: 3.5 s\n",
      "Wall time: 17.6 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation6.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "interesting-aircraft",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10425.320558545247"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation6.cells.counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "dirty-google",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "60626.140328808004"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation6.reads.counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "finnish-indian",
   "metadata": {},
   "source": [
    "more inaccurate than the other ones, but not too bad probably, especially since it's about twice as fast or whatever"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "equivalent-fireplace",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation7 = dropletSimulation(number_species=1000, number_droplets=1000, number_batches=5, \n",
    "             copy_numbers=random_copy_numbers(1000), frequency_vector=random_frequency_vector(1000),\n",
    "             glv_interaction_coefficients=2*(np.random.random((1000,1000))-0.5), glv_baserate_coefficients=np.random.random(1000), \n",
    "                                noise_scale=8,timestep = 0.0005, batch_window=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "color-feedback",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.84 s, sys: 1.47 s, total: 3.31 s\n",
      "Wall time: 6.16 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation7.run_simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "affiliated-lecture",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7899.869247461537"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation7.cells.counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "dangerous-arizona",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "29438.231687435455"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation7.reads.counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "altered-talent",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 868 ms, sys: 32.2 ms, total: 900 ms\n",
      "Wall time: 903 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation3.group_droplets()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "knowing-innocent",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 14.3 ms, sys: 1.85 ms, total: 16.1 ms\n",
      "Wall time: 14.8 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation4.group_droplets()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "harmful-underwear",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 19s, sys: 2.96 s, total: 1min 22s\n",
      "Wall time: 1min 22s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation5.group_droplets()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "arranged-intermediate",
   "metadata": {},
   "source": [
    "ok yeah this is about quadratic scaling, and just for the droplet grouping, not even the actual coefficient calculations...\n",
    "\n",
    "So maybe keeping the number of strains small would still be good? Yes maybe but we could still plausibly do 100 probably, or at least 20 or something. Perhaps..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "spare-rebecca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 24s, sys: 22.2 s, total: 1min 46s\n",
      "Wall time: 1min 47s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation6.group_droplets()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "suitable-wrapping",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 19s, sys: 3.61 s, total: 1min 23s\n",
      "Wall time: 1min 23s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation7.group_droplets()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "changed-latitude",
   "metadata": {},
   "source": [
    "keep in mind that even though this doesn't look that bad even though it's a 1000 species, it's also only 1,000 droplets per batch. It will/would get worse with more droplets per batch (not sure about the exact scaling but last time I checked it was super-linear so it's not great). Maybe it wouldn't be 8 hours this time but still"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "billion-relay",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation8 = dropletSimulation(number_species=1000, number_droplets=10000, number_batches=5, \n",
    "             copy_numbers=random_copy_numbers(1000), frequency_vector=random_frequency_vector(1000),\n",
    "             glv_interaction_coefficients=2*(np.random.random((1000,1000))-0.5), glv_baserate_coefficients=np.random.random(1000), \n",
    "                                noise_scale=8,timestep = 0.0001, batch_window=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "plastic-broad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 17.5 s, sys: 29.3 s, total: 46.8 s\n",
      "Wall time: 5min 44s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation8.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "black-condition",
   "metadata": {},
   "source": [
    "notice this is 10,000 droplets per batch this time\n",
    "\n",
    "so will the following only take 10 times longer to run? (which would be 15 minutes) \n",
    "\n",
    "or will it be even worse?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "hollow-celebration",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 17min 9s, sys: 12.3 s, total: 17min 21s\n",
      "Wall time: 17min 24s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation8.group_droplets()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ordinary-season",
   "metadata": {},
   "source": [
    "and then keep in mind too that the cost is not just in/about the droplet grouping -- you then actually have to go and compute the log ratio coefficients (and now also other coefficients) too, and so even though grouping the droplets will amortize you some time, if that's really slow (like in this case), each log ratio variant computation is going to be really slow and expensive too, so..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "liberal-mattress",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 21.6 s, sys: 4.9 s, total: 26.5 s\n",
      "Wall time: 27 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[[ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        ...,\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan]],\n",
       "\n",
       "       [[        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        ...,\n",
       "        [        nan,         nan,         nan, -2.39587673,\n",
       "                 nan],\n",
       "        [ 0.50572676,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan]],\n",
       "\n",
       "       [[        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        ...,\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan]],\n",
       "\n",
       "       ...,\n",
       "\n",
       "       [[        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan, -0.57057329,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        ...,\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan]],\n",
       "\n",
       "       [[        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [-0.52217614,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        ...,\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan]],\n",
       "\n",
       "       [[        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        ...,\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [        nan,         nan,         nan,         nan,\n",
       "                 nan],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ]]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "simulation8.reads.get_fitness_coefficients()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "found-heaven",
   "metadata": {},
   "source": [
    "wow dang, that's actually not bad at all, especially considering that we are talking about _one million_ log ratio coefficients (and then five times) like wow\n",
    "\n",
    "1/32 the compute time of the droplet grouping, impressive"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brilliant-market",
   "metadata": {},
   "source": [
    "Let's also check that the accuracy of these is what we expect from using such a small time step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "noble-denver",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000.000000000005"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation8.cells.counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cutting-product",
   "metadata": {},
   "source": [
    "### 100,000 droplets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "miniature-merit",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation9 = dropletSimulation(number_species=11, number_droplets=100000, number_batches=5, \n",
    "             copy_numbers=random_copy_numbers(11), frequency_vector=random_frequency_vector(11),\n",
    "             glv_interaction_coefficients=2*(np.random.random((11,11))-0.5), glv_baserate_coefficients=np.random.random(11), \n",
    "                                noise_scale=8,timestep = 0.0001, batch_window=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "beautiful-upper",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 4s, sys: 11.2 s, total: 1min 15s\n",
      "Wall time: 51min 11s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation9.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "outside-macintosh",
   "metadata": {},
   "source": [
    "Let's also check that the accuracy of these is what we expect from using such a small time step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "excess-democrat",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000.874344376047"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation9.cells.counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bound-fighter",
   "metadata": {},
   "source": [
    "Increase number of droplets, decrease number of strains..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "ordinary-mauritius",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.08 s, sys: 66.3 ms, total: 1.15 s\n",
      "Wall time: 1.17 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation9.group_droplets()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "authentic-samuel",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 206 ms, sys: 19 ms, total: 225 ms\n",
      "Wall time: 223 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[[ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [-0.01825097, -0.0103619 , -0.03298558, -0.15337361,\n",
       "         -0.17263177],\n",
       "        [-0.03135004,  0.01516   ,  0.03386091, -0.17329343,\n",
       "         -0.12162302],\n",
       "        [ 0.00517588,  0.02271118,  0.02730199,  0.09451613,\n",
       "          0.15043179],\n",
       "        [ 0.0232786 ,  0.04171229,  0.07555189,  0.0844107 ,\n",
       "          0.1625998 ],\n",
       "        [-0.06693257, -0.03944812,  0.16756002,  0.33471494,\n",
       "         -0.03309491],\n",
       "        [-0.00179352, -0.02499574, -0.0277908 , -0.05063811,\n",
       "         -0.00851819],\n",
       "        [-0.00135798, -0.02903899, -0.13542305, -0.01616914,\n",
       "         -0.01249851],\n",
       "        [ 0.0093042 , -0.02840098,  0.0815181 ,  0.15453214,\n",
       "          0.16477393],\n",
       "        [-0.00250099,  0.05505745,  0.0502911 ,  0.0260532 ,\n",
       "          0.10427131],\n",
       "        [ 0.004338  , -0.03608141, -0.00548705,  0.03409703,\n",
       "         -0.11218546]],\n",
       "\n",
       "       [[-0.00172497,  0.0639816 ,  0.01531354,  0.05848579,\n",
       "          0.12380861],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [-0.01070247,  0.01538503,  0.07699368,  0.11317913,\n",
       "          0.27846851],\n",
       "        [ 0.0166332 ,  0.02545433,  0.10365985,  0.18057597,\n",
       "          0.27225476],\n",
       "        [-0.02790372,  0.00512342, -0.02064296, -0.0241811 ,\n",
       "         -0.01644041],\n",
       "        [-0.01591886,  0.0092928 , -0.09470365, -0.1344142 ,\n",
       "          0.06529847],\n",
       "        [-0.00620027, -0.00643632,  0.04379513,  0.05107908,\n",
       "          0.06638549],\n",
       "        [-0.019999  , -0.00875582, -0.03784295,  0.02601414,\n",
       "         -0.07787298],\n",
       "        [-0.00237793, -0.00984054, -0.00936455, -0.02719362,\n",
       "          0.01175842],\n",
       "        [-0.00245831, -0.05024755, -0.06881612, -0.09461539,\n",
       "         -0.10955181],\n",
       "        [ 0.02449929,  0.02336477,  0.10237413,  0.20662202,\n",
       "          0.29017168]],\n",
       "\n",
       "       [[-0.02874556,  0.07660586, -0.00938398, -0.00585859,\n",
       "         -0.01585432],\n",
       "        [-0.0197619 ,  0.04056643, -0.02266376,  0.01616575,\n",
       "         -0.01973922],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [-0.00311762,  0.04273954, -0.00067901, -0.01202645,\n",
       "         -0.09081104],\n",
       "        [-0.01021197, -0.00464385, -0.03109359, -0.02586962,\n",
       "         -0.11967011],\n",
       "        [ 0.00593303,  0.02010219,  0.01335338, -0.07334328,\n",
       "         -0.10962177],\n",
       "        [-0.01746184, -0.02591856, -0.01541491, -0.10729301,\n",
       "         -0.13715425],\n",
       "        [ 0.04761971,  0.03336327,  0.12524195,  0.23068866,\n",
       "          0.26432456],\n",
       "        [-0.00756137,  0.00353805,  0.02435155, -0.00532799,\n",
       "          0.05462445],\n",
       "        [ 0.02263607,  0.03658474,  0.0379495 ,  0.10167241,\n",
       "          0.11439995],\n",
       "        [-0.02926711,  0.03075393,  0.00315539, -0.01050194,\n",
       "         -0.05275215]],\n",
       "\n",
       "       [[-0.03178725, -0.04515939, -0.05201545, -0.06285968,\n",
       "         -0.15970168],\n",
       "        [ 0.01628687,  0.00235046,  0.07187488,  0.10083445,\n",
       "          0.14786609],\n",
       "        [-0.01282703, -0.00168591, -0.02208396,  0.02705657,\n",
       "          0.02685056],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [ 0.00653941, -0.01470541,  0.01911329,  0.05104075,\n",
       "         -0.04483865],\n",
       "        [ 0.00499712,  0.00368569, -0.01205409, -0.03040456,\n",
       "         -0.03738056],\n",
       "        [-0.01347598, -0.04056777, -0.03454016, -0.03553581,\n",
       "         -0.08776208],\n",
       "        [-0.00977455,  0.03284329, -0.0542749 , -0.04689543,\n",
       "         -0.04186985],\n",
       "        [-0.01793357, -0.00814619, -0.01181884, -0.00548424,\n",
       "         -0.06261761],\n",
       "        [ 0.00790587, -0.0020723 , -0.05238689, -0.05353195,\n",
       "         -0.03207072],\n",
       "        [-0.00178487, -0.01233106, -0.01539856, -0.03943893,\n",
       "          0.01493137]],\n",
       "\n",
       "       [[-0.03386464, -0.04299008, -0.06745883, -0.0905388 ,\n",
       "         -0.11616184],\n",
       "        [ 0.0142115 , -0.01975838,  0.01193803, -0.11080343,\n",
       "         -0.04242539],\n",
       "        [ 0.00373523,  0.01877719, -0.00230132, -0.05996908,\n",
       "          0.04686579],\n",
       "        [-0.04819744, -0.05268583, -0.08659912, -0.11778918,\n",
       "         -0.08282316],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [ 0.00395371, -0.02881432,  0.1263989 ,  0.01960831,\n",
       "         -0.03622206],\n",
       "        [-0.00593261, -0.03515896, -0.03480482, -0.10330181,\n",
       "         -0.10153955],\n",
       "        [-0.04005944, -0.03720525, -0.04241362, -0.02250787,\n",
       "          0.02921044],\n",
       "        [ 0.00556044, -0.01985998, -0.04064117, -0.00870518,\n",
       "         -0.06352303],\n",
       "        [-0.00387929, -0.03958916, -0.06886244, -0.04479051,\n",
       "         -0.07267354],\n",
       "        [-0.02040442,  0.00493969, -0.0521997 ,  0.00663409,\n",
       "          0.03472286]],\n",
       "\n",
       "       [[-0.03542633, -0.01840798,  0.00939231, -0.0295982 ,\n",
       "         -0.04755389],\n",
       "        [ 0.02779356, -0.00597425, -0.10143807, -0.09875709,\n",
       "         -0.2697176 ],\n",
       "        [-0.03590956, -0.07671526, -0.18278688, -0.07672435,\n",
       "         -0.24838507],\n",
       "        [-0.01357568, -0.06328916, -0.09403815, -0.08410579,\n",
       "         -0.11704248],\n",
       "        [-0.01873469,  0.05579844,  0.04910796, -0.0654448 ,\n",
       "         -0.08275238],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [ 0.01987792,  0.01624754, -0.01776725,  0.07186232,\n",
       "          0.05383208],\n",
       "        [-0.02166029,  0.0079977 ,  0.09580938,  0.10359828,\n",
       "         -0.03162375],\n",
       "        [-0.02881153, -0.05578067, -0.06178882, -0.10857402,\n",
       "         -0.12935896],\n",
       "        [ 0.00655174, -0.01678589,  0.01018726, -0.13137863,\n",
       "         -0.08893096],\n",
       "        [-0.01746804, -0.06483354,  0.05706119, -0.09550154,\n",
       "         -0.19082639]],\n",
       "\n",
       "       [[ 0.0045219 ,  0.01824575,  0.11587799,  0.11373268,\n",
       "          0.05657277],\n",
       "        [-0.0236208 , -0.03687954, -0.05170003, -0.10242084,\n",
       "         -0.1491951 ],\n",
       "        [ 0.0036473 , -0.00183237, -0.04419255, -0.00374838,\n",
       "          0.01416396],\n",
       "        [ 0.00078155,  0.03527762,  0.06992705,  0.08519993,\n",
       "          0.15345661],\n",
       "        [ 0.02422151,  0.02538854,  0.01579934, -0.02491949,\n",
       "          0.08406237],\n",
       "        [ 0.00530841, -0.02355443,  0.05866898,  0.06983659,\n",
       "          0.2489976 ],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [-0.0231137 , -0.02404875, -0.01623252, -0.07795175,\n",
       "         -0.06088698],\n",
       "        [-0.00314958, -0.00809572, -0.00905638, -0.04048327,\n",
       "         -0.0378336 ],\n",
       "        [-0.01193171, -0.01960488, -0.01726516, -0.00128285,\n",
       "         -0.01699873],\n",
       "        [-0.03568038, -0.04465108, -0.01004902, -0.12364462,\n",
       "         -0.06871892]],\n",
       "\n",
       "       [[-0.02123875, -0.03504875,  0.09939972,  0.04040299,\n",
       "          0.09699384],\n",
       "        [-0.02221508, -0.02911852, -0.02507108, -0.05006609,\n",
       "         -0.08218783],\n",
       "        [-0.01908209, -0.01900972,  0.00167046,  0.00880114,\n",
       "         -0.02069071],\n",
       "        [ 0.00201531,  0.02480695,  0.07903234,  0.06308   ,\n",
       "          0.11590815],\n",
       "        [ 0.00726951, -0.00262689, -0.04231114, -0.06702315,\n",
       "         -0.08670756],\n",
       "        [ 0.03578221,  0.02062698, -0.08835598, -0.10289039,\n",
       "         -0.11250069],\n",
       "        [ 0.00966558, -0.01099901, -0.0281324 , -0.01113943,\n",
       "         -0.08618512],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [-0.01044958,  0.00347148, -0.0334945 , -0.02453698,\n",
       "         -0.08115948],\n",
       "        [-0.00542135, -0.03580308, -0.03408554, -0.07800614,\n",
       "         -0.09427269],\n",
       "        [ 0.02734317, -0.03937146, -0.00661   ,  0.04210474,\n",
       "         -0.03231952]],\n",
       "\n",
       "       [[-0.01217058, -0.08880397, -0.09298823, -0.08372052,\n",
       "         -0.15686606],\n",
       "        [ 0.00977278,  0.02551976,  0.06911932,  0.12660542,\n",
       "          0.15044027],\n",
       "        [ 0.01495126,  0.00370558,  0.0331965 , -0.01767977,\n",
       "         -0.03170234],\n",
       "        [-0.00967639, -0.01044169,  0.0152909 ,  0.03626539,\n",
       "          0.05600491],\n",
       "        [-0.01188443,  0.02005305,  0.04071272,  0.00683079,\n",
       "          0.08987045],\n",
       "        [-0.00268957,  0.04671858, -0.05218644,  0.01473601,\n",
       "          0.04470501],\n",
       "        [-0.00193843,  0.02211943,  0.02278956,  0.03883615,\n",
       "          0.0860794 ],\n",
       "        [ 0.03511145,  0.04728474,  0.01581748,  0.10258971,\n",
       "          0.07088164],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [ 0.02480108, -0.01897226, -0.10091223, -0.06920278,\n",
       "         -0.07931178],\n",
       "        [-0.00309878, -0.0422379 ,  0.00757908, -0.03345665,\n",
       "         -0.06366061]],\n",
       "\n",
       "       [[-0.01081709,  0.08534228, -0.0583901 , -0.0630957 ,\n",
       "         -0.00645193],\n",
       "        [-0.00553945,  0.01199428, -0.02789551,  0.03121219,\n",
       "          0.0034862 ],\n",
       "        [-0.01804381, -0.03929065, -0.08855681, -0.06449319,\n",
       "         -0.14888678],\n",
       "        [ 0.03107587,  0.01341471,  0.06353335,  0.03600807,\n",
       "          0.10210204],\n",
       "        [ 0.02803276,  0.05399963,  0.03779605,  0.06155922,\n",
       "          0.08172269],\n",
       "        [ 0.00390406, -0.05576404, -0.10935606,  0.01782554,\n",
       "         -0.05192889],\n",
       "        [-0.00907845, -0.00143226,  0.00864332, -0.01259945,\n",
       "          0.00087686],\n",
       "        [ 0.01608639, -0.00815361, -0.0273154 ,  0.01821173,\n",
       "          0.06110142],\n",
       "        [ 0.00230748, -0.04512888, -0.02613871, -0.06641315,\n",
       "         -0.09106387],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ],\n",
       "        [-0.01702958, -0.01616833,  0.03094841,  0.08886031,\n",
       "          0.16064726]],\n",
       "\n",
       "       [[-0.07594208, -0.0144625 , -0.15902036, -0.10677819,\n",
       "         -0.09817061],\n",
       "        [ 0.00199943,  0.02092358,  0.00179564, -0.08860236,\n",
       "         -0.10138197],\n",
       "        [-0.01171504, -0.01901548,  0.06436821,  0.00081217,\n",
       "          0.02687616],\n",
       "        [ 0.00828742, -0.03662775, -0.02944251, -0.02599866,\n",
       "         -0.06003979],\n",
       "        [-0.0205805 , -0.01850506, -0.03193337,  0.04310841,\n",
       "         -0.04526514],\n",
       "        [-0.02421393, -0.1373481 ,  0.07006332, -0.11151502,\n",
       "         -0.17427185],\n",
       "        [-0.00195243, -0.03516162, -0.03467714, -0.05570402,\n",
       "         -0.06732243],\n",
       "        [ 0.0444772 , -0.01230211, -0.01670251,  0.02550253,\n",
       "          0.16000009],\n",
       "        [-0.01668879,  0.02286371,  0.0679069 ,  0.12904353,\n",
       "          0.23207559],\n",
       "        [-0.01606663,  0.01539945,  0.10212688,  0.07685599,\n",
       "          0.05631201],\n",
       "        [ 0.        ,  0.        ,  0.        ,  0.        ,\n",
       "          0.        ]]])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "simulation9.reads.get_fitness_coefficients()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "several-atlanta",
   "metadata": {},
   "source": [
    "Allow us to no longer compromise when it comes to the number of strains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "higher-poland",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation10 = dropletSimulation(number_species=1000, number_droplets=100000, number_batches=5, \n",
    "             copy_numbers=random_copy_numbers(1000), frequency_vector=random_frequency_vector(1000),\n",
    "             glv_interaction_coefficients=2*(np.random.random((1000,1000))-0.5), glv_baserate_coefficients=np.random.random(1000), \n",
    "                                noise_scale=8,timestep = 0.0001, batch_window=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecological-documentary",
   "metadata": {},
   "source": [
    "Computer keeps on dying before/instead of running these, which seems like a bad sign...\n",
    "\n",
    "**Update:** it's because my computer keeps running out of memory apparently. I guess this means I wrote bad code that has a major memory leak? uh oh I'm not sure\n",
    "\n",
    "**2nd Update:** Well the code may be bad, but because it's 1000 species with 100,000 droplets with 5 batches we have 500 million entries in the numpy arrays, then two numpy arrays that size, so 1 billion entries, each entry is a float and so probably more than 8 bytes, so works out to more 8GB RAM. I'm guessing/vaguely remember last time I was only doing 10,000 droplets with 1000 species, which would correspond to 1/10 the memory requirement, so that this doesn't run probably isn't too surprising."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "exempt-nursery",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2min, sys: 1min 32s, total: 3min 33s\n",
      "Wall time: 52min 35s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "simulation10.run_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worthy-kentucky",
   "metadata": {},
   "source": [
    "Let's also check that the accuracy of these is what we expect from using such a small time step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "noted-bidder",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000.457629454848"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(simulation10.cells.counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unique-cleaning",
   "metadata": {},
   "source": [
    "is the following only ~10 times as long as previous `group_droplets`? or ~100? \n",
    "\n",
    "(this killed my kernel at least once, so I'm not sure whether it will actually really run or not...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "choice-clerk",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "simulation10.group_droplets()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fourth-blame",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "simulation10.reads.get_fitness_coefficients()"
   ]
  }
 ],
 "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}