{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "143a15f7-ac7e-4e8c-a32e-cbbb7139acb5",
   "metadata": {},
   "source": [
    "# Modeling nonhomogeneous Poisson spiking {#sec:variable-spiking}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "631a8658",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| code-fold: true\n",
    "\n",
    "# Colab setup ------------------\n",
    "import os, sys, subprocess\n",
    "if \"google.colab\" in sys.modules:\n",
    "    cmd = \"pip install --upgrade iqplot watermark\"\n",
    "    process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "    stdout, stderr = process.communicate()\n",
    "# ------------------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4f944b2e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "    <style>\n",
       "        .bk-notebook-logo {\n",
       "            display: block;\n",
       "            width: 20px;\n",
       "            height: 20px;\n",
       "            background-image: url(data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABQAAAAUCAYAAACNiR0NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAABx0RVh0U29mdHdhcmUAQWRvYmUgRmlyZXdvcmtzIENTNui8sowAAAOkSURBVDiNjZRtaJVlGMd/1/08zzln5zjP1LWcU9N0NkN8m2CYjpgQYQXqSs0I84OLIC0hkEKoPtiH3gmKoiJDU7QpLgoLjLIQCpEsNJ1vqUOdO7ppbuec5+V+rj4ctwzd8IIbbi6u+8f1539dt3A78eXC7QizUF7gyV1fD1Yqg4JWz84yffhm0qkFqBogB9rM8tZdtwVsPUhWhGcFJngGeWrPzHm5oaMmkfEg1usvLFyc8jLRqDOMru7AyC8saQr7GG7f5fvDeH7Ej8CM66nIF+8yngt6HWaKh7k49Soy9nXurCi1o3qUbS3zWfrYeQDTB/Qj6kX6Ybhw4B+bOYoLKCC9H3Nu/leUTZ1JdRWkkn2ldcCamzrcf47KKXdAJllSlxAOkRgyHsGC/zRday5Qld9DyoM4/q/rUoy/CXh3jzOu3bHUVZeU+DEn8FInkPBFlu3+nW3Nw0mk6vCDiWg8CeJaxEwuHS3+z5RgY+YBR6V1Z1nxSOfoaPa4LASWxxdNp+VWTk7+4vzaou8v8PN+xo+KY2xsw6une2frhw05CTYOmQvsEhjhWjn0bmXPjpE1+kplmmkP3suftwTubK9Vq22qKmrBhpY4jvd5afdRA3wGjFAgcnTK2s4hY0/GPNIb0nErGMCRxWOOX64Z8RAC4oCXdklmEvcL8o0BfkNK4lUg9HTl+oPlQxdNo3Mg4Nv175e/1LDGzZen30MEjRUtmXSfiTVu1kK8W4txyV6BMKlbgk3lMwYCiusNy9fVfvvwMxv8Ynl6vxoByANLTWplvuj/nF9m2+PDtt1eiHPBr1oIfhCChQMBw6Aw0UulqTKZdfVvfG7VcfIqLG9bcldL/+pdWTLxLUy8Qq38heUIjh4XlzZxzQm19lLFlr8vdQ97rjZVOLf8nclzckbcD4wxXMidpX30sFd37Fv/GtwwhzhxGVAprjbg0gCAEeIgwCZyTV2Z1REEW8O4py0wsjeloKoMr6iCY6dP92H6Vw/oTyICIthibxjm/DfN9lVz8IqtqKYLUXfoKVMVQVVJOElGjrnnUt9T9wbgp8AyYKaGlqingHZU/uG2NTZSVqwHQTWkx9hxjkpWDaCg6Ckj5qebgBVbT3V3NNXMSiWSDdGV3hrtzla7J+duwPOToIg42ChPQOQjspnSlp1V+Gjdged7+8UN5CRAV7a5EdFNwCjEaBR27b3W890TE7g24NAP/mMDXRWrGoFPQI9ls/MWO2dWFAar/xcOIImbbpA3zgAAAABJRU5ErkJggg==);\n",
       "        }\n",
       "    </style>\n",
       "    <div>\n",
       "        <a href=\"https://bokeh.org\" target=\"_blank\" class=\"bk-notebook-logo\"></a>\n",
       "        <span id=\"b389e7ec-1a47-41c8-8288-25dab2d4b901\">Loading BokehJS ...</span>\n",
       "    </div>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "'use strict';\n",
       "(function(root) {\n",
       "  function now() {\n",
       "    return new Date();\n",
       "  }\n",
       "\n",
       "  const force = true;\n",
       "\n",
       "  if (typeof root._bokeh_onload_callbacks === \"undefined\" || force === true) {\n",
       "    root._bokeh_onload_callbacks = [];\n",
       "    root._bokeh_is_loading = undefined;\n",
       "  }\n",
       "\n",
       "const JS_MIME_TYPE = 'application/javascript';\n",
       "  const HTML_MIME_TYPE = 'text/html';\n",
       "  const EXEC_MIME_TYPE = 'application/vnd.bokehjs_exec.v0+json';\n",
       "  const CLASS_NAME = 'output_bokeh rendered_html';\n",
       "\n",
       "  /**\n",
       "   * Render data to the DOM node\n",
       "   */\n",
       "  function render(props, node) {\n",
       "    const script = document.createElement(\"script\");\n",
       "    node.appendChild(script);\n",
       "  }\n",
       "\n",
       "  /**\n",
       "   * Handle when an output is cleared or removed\n",
       "   */\n",
       "  function handleClearOutput(event, handle) {\n",
       "    function drop(id) {\n",
       "      const view = Bokeh.index.get_by_id(id)\n",
       "      if (view != null) {\n",
       "        view.model.document.clear()\n",
       "        Bokeh.index.delete(view)\n",
       "      }\n",
       "    }\n",
       "\n",
       "    const cell = handle.cell;\n",
       "\n",
       "    const id = cell.output_area._bokeh_element_id;\n",
       "    const server_id = cell.output_area._bokeh_server_id;\n",
       "\n",
       "    // Clean up Bokeh references\n",
       "    if (id != null) {\n",
       "      drop(id)\n",
       "    }\n",
       "\n",
       "    if (server_id !== undefined) {\n",
       "      // Clean up Bokeh references\n",
       "      const cmd_clean = \"from bokeh.io.state import curstate; print(curstate().uuid_to_server['\" + server_id + \"'].get_sessions()[0].document.roots[0]._id)\";\n",
       "      cell.notebook.kernel.execute(cmd_clean, {\n",
       "        iopub: {\n",
       "          output: function(msg) {\n",
       "            const id = msg.content.text.trim()\n",
       "            drop(id)\n",
       "          }\n",
       "        }\n",
       "      });\n",
       "      // Destroy server and session\n",
       "      const cmd_destroy = \"import bokeh.io.notebook as ion; ion.destroy_server('\" + server_id + \"')\";\n",
       "      cell.notebook.kernel.execute(cmd_destroy);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  /**\n",
       "   * Handle when a new output is added\n",
       "   */\n",
       "  function handleAddOutput(event, handle) {\n",
       "    const output_area = handle.output_area;\n",
       "    const output = handle.output;\n",
       "\n",
       "    // limit handleAddOutput to display_data with EXEC_MIME_TYPE content only\n",
       "    if ((output.output_type != \"display_data\") || (!Object.prototype.hasOwnProperty.call(output.data, EXEC_MIME_TYPE))) {\n",
       "      return\n",
       "    }\n",
       "\n",
       "    const toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n",
       "\n",
       "    if (output.metadata[EXEC_MIME_TYPE][\"id\"] !== undefined) {\n",
       "      toinsert[toinsert.length - 1].firstChild.textContent = output.data[JS_MIME_TYPE];\n",
       "      // store reference to embed id on output_area\n",
       "      output_area._bokeh_element_id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n",
       "    }\n",
       "    if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n",
       "      const bk_div = document.createElement(\"div\");\n",
       "      bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n",
       "      const script_attrs = bk_div.children[0].attributes;\n",
       "      for (let i = 0; i < script_attrs.length; i++) {\n",
       "        toinsert[toinsert.length - 1].firstChild.setAttribute(script_attrs[i].name, script_attrs[i].value);\n",
       "        toinsert[toinsert.length - 1].firstChild.textContent = bk_div.children[0].textContent\n",
       "      }\n",
       "      // store reference to server id on output_area\n",
       "      output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function register_renderer(events, OutputArea) {\n",
       "\n",
       "    function append_mime(data, metadata, element) {\n",
       "      // create a DOM node to render to\n",
       "      const toinsert = this.create_output_subarea(\n",
       "        metadata,\n",
       "        CLASS_NAME,\n",
       "        EXEC_MIME_TYPE\n",
       "      );\n",
       "      this.keyboard_manager.register_events(toinsert);\n",
       "      // Render to node\n",
       "      const props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n",
       "      render(props, toinsert[toinsert.length - 1]);\n",
       "      element.append(toinsert);\n",
       "      return toinsert\n",
       "    }\n",
       "\n",
       "    /* Handle when an output is cleared or removed */\n",
       "    events.on('clear_output.CodeCell', handleClearOutput);\n",
       "    events.on('delete.Cell', handleClearOutput);\n",
       "\n",
       "    /* Handle when a new output is added */\n",
       "    events.on('output_added.OutputArea', handleAddOutput);\n",
       "\n",
       "    /**\n",
       "     * Register the mime type and append_mime function with output_area\n",
       "     */\n",
       "    OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n",
       "      /* Is output safe? */\n",
       "      safe: true,\n",
       "      /* Index of renderer in `output_area.display_order` */\n",
       "      index: 0\n",
       "    });\n",
       "  }\n",
       "\n",
       "  // register the mime type if in Jupyter Notebook environment and previously unregistered\n",
       "  if (root.Jupyter !== undefined) {\n",
       "    const events = require('base/js/events');\n",
       "    const OutputArea = require('notebook/js/outputarea').OutputArea;\n",
       "\n",
       "    if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n",
       "      register_renderer(events, OutputArea);\n",
       "    }\n",
       "  }\n",
       "  if (typeof (root._bokeh_timeout) === \"undefined\" || force === true) {\n",
       "    root._bokeh_timeout = Date.now() + 5000;\n",
       "    root._bokeh_failed_load = false;\n",
       "  }\n",
       "\n",
       "  const NB_LOAD_WARNING = {'data': {'text/html':\n",
       "     \"<div style='background-color: #fdd'>\\n\"+\n",
       "     \"<p>\\n\"+\n",
       "     \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n",
       "     \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n",
       "     \"</p>\\n\"+\n",
       "     \"<ul>\\n\"+\n",
       "     \"<li>re-rerun `output_notebook()` to attempt to load from CDN again, or</li>\\n\"+\n",
       "     \"<li>use INLINE resources instead, as so:</li>\\n\"+\n",
       "     \"</ul>\\n\"+\n",
       "     \"<code>\\n\"+\n",
       "     \"from bokeh.resources import INLINE\\n\"+\n",
       "     \"output_notebook(resources=INLINE)\\n\"+\n",
       "     \"</code>\\n\"+\n",
       "     \"</div>\"}};\n",
       "\n",
       "  function display_loaded(error = null) {\n",
       "    const el = document.getElementById(\"b389e7ec-1a47-41c8-8288-25dab2d4b901\");\n",
       "    if (el != null) {\n",
       "      const html = (() => {\n",
       "        if (typeof root.Bokeh === \"undefined\") {\n",
       "          if (error == null) {\n",
       "            return \"BokehJS is loading ...\";\n",
       "          } else {\n",
       "            return \"BokehJS failed to load.\";\n",
       "          }\n",
       "        } else {\n",
       "          const prefix = `BokehJS ${root.Bokeh.version}`;\n",
       "          if (error == null) {\n",
       "            return `${prefix} successfully loaded.`;\n",
       "          } else {\n",
       "            return `${prefix} <b>encountered errors</b> while loading and may not function as expected.`;\n",
       "          }\n",
       "        }\n",
       "      })();\n",
       "      el.innerHTML = html;\n",
       "\n",
       "      if (error != null) {\n",
       "        const wrapper = document.createElement(\"div\");\n",
       "        wrapper.style.overflow = \"auto\";\n",
       "        wrapper.style.height = \"5em\";\n",
       "        wrapper.style.resize = \"vertical\";\n",
       "        const content = document.createElement(\"div\");\n",
       "        content.style.fontFamily = \"monospace\";\n",
       "        content.style.whiteSpace = \"pre-wrap\";\n",
       "        content.style.backgroundColor = \"rgb(255, 221, 221)\";\n",
       "        content.textContent = error.stack ?? error.toString();\n",
       "        wrapper.append(content);\n",
       "        el.append(wrapper);\n",
       "      }\n",
       "    } else if (Date.now() < root._bokeh_timeout) {\n",
       "      setTimeout(() => display_loaded(error), 100);\n",
       "    }\n",
       "  }\n",
       "\n",
       "  function run_callbacks() {\n",
       "    try {\n",
       "      root._bokeh_onload_callbacks.forEach(function(callback) {\n",
       "        if (callback != null)\n",
       "          callback();\n",
       "      });\n",
       "    } finally {\n",
       "      delete root._bokeh_onload_callbacks\n",
       "    }\n",
       "    console.debug(\"Bokeh: all callbacks have finished\");\n",
       "  }\n",
       "\n",
       "  function load_libs(css_urls, js_urls, callback) {\n",
       "    if (css_urls == null) css_urls = [];\n",
       "    if (js_urls == null) js_urls = [];\n",
       "\n",
       "    root._bokeh_onload_callbacks.push(callback);\n",
       "    if (root._bokeh_is_loading > 0) {\n",
       "      console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n",
       "      return null;\n",
       "    }\n",
       "    if (js_urls == null || js_urls.length === 0) {\n",
       "      run_callbacks();\n",
       "      return null;\n",
       "    }\n",
       "    console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n",
       "    root._bokeh_is_loading = css_urls.length + js_urls.length;\n",
       "\n",
       "    function on_load() {\n",
       "      root._bokeh_is_loading--;\n",
       "      if (root._bokeh_is_loading === 0) {\n",
       "        console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n",
       "        run_callbacks()\n",
       "      }\n",
       "    }\n",
       "\n",
       "    function on_error(url) {\n",
       "      console.error(\"failed to load \" + url);\n",
       "    }\n",
       "\n",
       "    for (let i = 0; i < css_urls.length; i++) {\n",
       "      const url = css_urls[i];\n",
       "      const element = document.createElement(\"link\");\n",
       "      element.onload = on_load;\n",
       "      element.onerror = on_error.bind(null, url);\n",
       "      element.rel = \"stylesheet\";\n",
       "      element.type = \"text/css\";\n",
       "      element.href = url;\n",
       "      console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n",
       "      document.body.appendChild(element);\n",
       "    }\n",
       "\n",
       "    for (let i = 0; i < js_urls.length; i++) {\n",
       "      const url = js_urls[i];\n",
       "      const element = document.createElement('script');\n",
       "      element.onload = on_load;\n",
       "      element.onerror = on_error.bind(null, url);\n",
       "      element.async = false;\n",
       "      element.src = url;\n",
       "      console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n",
       "      document.head.appendChild(element);\n",
       "    }\n",
       "  };\n",
       "\n",
       "  function inject_raw_css(css) {\n",
       "    const element = document.createElement(\"style\");\n",
       "    element.appendChild(document.createTextNode(css));\n",
       "    document.body.appendChild(element);\n",
       "  }\n",
       "\n",
       "  const js_urls = [\"https://cdn.bokeh.org/bokeh/release/bokeh-3.6.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.6.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.6.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.6.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-mathjax-3.6.2.min.js\"];\n",
       "  const css_urls = [];\n",
       "\n",
       "  const inline_js = [    function(Bokeh) {\n",
       "      Bokeh.set_log_level(\"info\");\n",
       "    },\n",
       "function(Bokeh) {\n",
       "    }\n",
       "  ];\n",
       "\n",
       "  function run_inline_js() {\n",
       "    if (root.Bokeh !== undefined || force === true) {\n",
       "      try {\n",
       "            for (let i = 0; i < inline_js.length; i++) {\n",
       "      inline_js[i].call(root, root.Bokeh);\n",
       "    }\n",
       "\n",
       "      } catch (error) {display_loaded(error);throw error;\n",
       "      }if (force === true) {\n",
       "        display_loaded();\n",
       "      }} else if (Date.now() < root._bokeh_timeout) {\n",
       "      setTimeout(run_inline_js, 100);\n",
       "    } else if (!root._bokeh_failed_load) {\n",
       "      console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n",
       "      root._bokeh_failed_load = true;\n",
       "    } else if (force !== true) {\n",
       "      const cell = $(document.getElementById(\"b389e7ec-1a47-41c8-8288-25dab2d4b901\")).parents('.cell').data().cell;\n",
       "      cell.output_area.append_execute_result(NB_LOAD_WARNING)\n",
       "    }\n",
       "  }\n",
       "\n",
       "  if (root._bokeh_is_loading === 0) {\n",
       "    console.debug(\"Bokeh: BokehJS loaded, going straight to plotting\");\n",
       "    run_inline_js();\n",
       "  } else {\n",
       "    load_libs(css_urls, js_urls, function() {\n",
       "      console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n",
       "      run_inline_js();\n",
       "    });\n",
       "  }\n",
       "}(window));"
      ],
      "application/vnd.bokehjs_load.v0+json": ""
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "import iqplot\n",
    "\n",
    "import bokeh.io\n",
    "bokeh.io.output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1661152-e69e-4adf-914c-8d846acdc41b",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "We have previously modeled spiking using a noisy leaky integrate-and-fire model. We saw that in some limits, the noisy LIF model gives Poissonian spiking such that the interspike intervals are Exponentially distributed and the number of spikes in a given time interval are Poisson distributed. Poisson process-based models of spiking are widely used to good effect.\n",
    "\n",
    "However, though well-modelled as Poisson, the rate of spiking can vary over time in what we call a **nonhomogeneous Poisson process**. Here, we derive the PDF for ISIs in this context and show how to sample ISIs.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2aefada",
   "metadata": {},
   "source": [
    "## The PDF for arrival times of a nonhomogenous Poisson process {#sec-pdf-nonhomogeneous-poisson}\n",
    "\n",
    "Consider the waiting time for an arrival of a Poisson process that arrives at rate $\\beta$. We know that the waiting time is Exponentially distributed.\n",
    "\n",
    "$$\\begin{align}\n",
    "f(t\\mid \\beta) = \\beta \\mathrm{e}^{-\\beta t}.\n",
    "\\end{align}\n",
    "$${#eq-exponential-pdf}\n",
    "\n",
    "The probability that a Poisson process has *not* arrived at or before time $t$ is given by the complementary cumulative distribution function (CCDF) of the Exponential distribution.\n",
    "\n",
    "$$\\begin{align}\n",
    "P(\\text{not arrived by time }t) = 1 - P(\\text{arrived by time }t) = 1 - \\int_0^t\\mathrm{d}t'\\,f(t\\mid\\beta) = 1 - (1 - \\mathrm{e}^{-\\beta t}) = \\mathrm{e}^{-\\beta t}.\n",
    "\\end{align}\n",
    "$${#eq-prob-not-arrived-nonhomogeneous-poisson-1}\n",
    "\n",
    "Now, consider the case where the arrival rate varies with time; $\\beta = \\beta(t)$. Consider some small interval, $[0, \\Delta t]$. The probability that no Poisson process has arrived in time $\\Delta t$ is\n",
    "\n",
    "$$\\begin{align}\n",
    "P(\\text{not arrived by time } \\Delta t) \\approx \\mathrm{e}^{-\\beta(\\Delta t) \\Delta t},\n",
    "\\end{align}\n",
    "$${#eq-prob-not-arrived-nonhomogeneous-poisson-2}\n",
    "\n",
    "where $\\Delta t$ is small so that $\\beta(t)$ does not change appreciably over the time interval and $\\beta(t) \\approx \\beta(t+\\Delta t)$. We will soon take the $\\Delta t\\to 0$ limit, so this approximation is justified.\n",
    "\n",
    "Now consider another interval of length $\\Delta t$ right after the one we just considered. The probability that no process arrives in that interval is\n",
    "\n",
    "$$\\begin{align}\n",
    "P(\\text{not arrived in time window } [\\Delta t, 2\\Delta t]) \\approx \\mathrm{e}^{-\\beta(2\\Delta t) \\Delta t}.\n",
    "\\end{align}\n",
    "$${#eq-prob-not-arrived-nonhomogeneous-poisson-3}\n",
    "\n",
    "Therefore, owing to the memorylessness of Poisson processes, the probability that no processes arrived in the interval $[0, 2\\Delta t]$ is the product of the probabilities of not arriving in the respective subintervals,\n",
    "\n",
    "$$\\begin{align}\n",
    "P(\\text{not arrived in time } 2\\Delta t) \\approx \\mathrm{e}^{-\\beta(\\Delta t) \\Delta t}\\,\\mathrm{e}^{-\\beta(2\\Delta t) \\Delta t}.\n",
    "\\end{align}\n",
    "$${#eq-prob-not-arrived-nonhomogeneous-poisson-4}\n",
    "\n",
    "We can consider $m$ such intervals.\n",
    "\n",
    "$$\\begin{align}\n",
    "P(\\text{not arrived between in time } m\\Delta t) \\approx \\prod_{k=1}^m \\mathrm{e}^{-\\beta(k\\Delta t) \\Delta t} = \\exp\\left[-\\sum_{k=1}^m \\Delta t \\,\\beta(k\\Delta t)\\right].\n",
    "\\end{align}\n",
    "$${#eq-prob-not-arrived-nonhomogeneous-poisson-5}\n",
    "\n",
    "In the limit of $\\Delta t\\to 0$, the Riemann sum becomes an integral. Taking $t = m\\Delta t$, we have\n",
    "\n",
    "$$\\begin{align}\n",
    "P(\\text{not arrived by time } t) = \\exp\\left[- \\int_0^t\\mathrm{d}t'\\,\\beta(t')\\right].\n",
    "\\end{align}\n",
    "$${#eq-prob-not-arrived-nonhomogeneous-poisson-6}\n",
    "\n",
    "This is the CCDF of the distribution desribing the arrival of a Poisson process with variable rate $\\beta(t)$. The CDF is then $1 - \\text{CCDF}$,\n",
    "\n",
    "$$\\begin{align}\n",
    "F(t\\mid \\beta(t)) = 1 - \\exp\\left[- \\int_0^t\\mathrm{d}t'\\,\\beta(t')\\right].\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-ccdf}\n",
    "\n",
    "The probability density function is then\n",
    "\n",
    "$$\\begin{align}\n",
    "f(t\\mid \\beta(t)) &= \\frac{\\mathrm{d}F}{\\mathrm{d}t} = -\\frac{\\mathrm{d}}{\\mathrm{d}t}\\,\\exp\\left[- \\int_0^t\\mathrm{d}t'\\,\\beta(t')\\right]\n",
    "= -\\exp\\left[- \\int_0^t\\mathrm{d}t'\\,\\beta(t')\\right]\\,\\frac{\\mathrm{d}}{\\mathrm{d}t}\\left(- \\int_0^t\\mathrm{d}t'\\,\\beta(t')\\right) \\nonumber \\\\[1em]\n",
    "&= \\beta(t)\\,\\exp\\left[- \\int_0^t\\mathrm{d}t'\\,\\beta(t')\\right].\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-pdf-1}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9fe8909-1b8b-4315-9feb-e8893aa8c656",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "### Multiple arrivals\n",
    "\n",
    "Say we observe arrivals of a Poisson process one after another. If we have $n$ arrivals of a Poisson process with a variable rate occuring at times $\\mathbf{t} = \\{t_1, t_2, \\ldots, t_n\\}$, where the times are ordered, the PDF is given by the product of the PDFs for each waiting time.\n",
    "\n",
    "$$\\begin{align}\n",
    "f(\\mathbf{t} \\mid \\beta(t)) = \\prod_{i=1}^n\\beta(t_i)\\exp\\left[- \\int_{t_{i-1}}^{t_i}\\mathrm{d}t'\\,\\beta(t')\\right]\n",
    "= \\left(\\prod_{i=1}^n\\beta(t_i)\\right)\\exp\\left[- \\sum_{i=1}^n\\int_{t_{i-1}}^{t_i}\\mathrm{d}t'\\,\\beta(t')\\right],\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-multiple-arrivals-pdf-1}\n",
    "\n",
    "where $t_0$ is when we started observing (it is not a time of an arrival of a Poisson process), and we take $t_0 = 0$. The integrals add together, giving\n",
    "\n",
    "$$\\begin{align}\n",
    "f(\\mathbf{t} \\mid \\beta(t)) = \\left(\\prod_{i=1}^n\\beta(t_i)\\right)\\exp\\left[-\\int_0^{t_n}\\mathrm{d}t'\\,\\beta(t')\\right].\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-multiple-arrivals-pdf}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a32730b-8d47-4965-b656-1e8050844d1a",
   "metadata": {},
   "source": [
    "### A finite observation time\n",
    "\n",
    "In practice, we start observing at time $t = 0$ and end at time $t = T$. While $t_n$ is the time of the last Poisson process we saw arrive, we should also take into account that we continued watching for time $T - t_n$ and saw no arrivals during that time. We therefore have\n",
    "\n",
    "$$\\begin{align}\n",
    "f(\\mathbf{t}, T;\\beta(t)) &= f(\\mathbf{t} \\mid \\beta(t)) \\times P(\\text{no arrivals between }t_n\\text{ and }T) \\nonumber \\\\[1em]\n",
    "&= \\left(\\prod_{i=1}^n\\beta(t_i)\\right)\\exp\\left[-\\int_0^{t_n}\\mathrm{d}t'\\,\\beta(t')\\right]\\,\\exp\\left[\\int_{t_n}^T\\mathrm{d}t'\\,\\beta(t')\\right] \\nonumber  \\\\[1em]\n",
    "&= \\left(\\prod_{i=1}^n\\beta(t_i)\\right)\\exp\\left[-\\int_0^{T}\\mathrm{d}t'\\,\\beta(t')\\right].\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-pdf}\n",
    "\n",
    "Note that in the special case of a constant rate, $\\beta(t) = \\beta$, we get\n",
    "\n",
    "$$\\begin{align}\n",
    "f(\\mathbf{t}, T;\\beta(t)) &= \\beta^n\\,\\mathrm{e}^{-\\beta T},\n",
    "\\end{align}\n",
    "$${#eq-homogeneous-poisson-pdf}\n",
    "\n",
    "In this case, the PDF is not dependent on the times, but only on the number of arrivals in the time interval. We could have derived by noting that the interrarrival times are Exponentially distributed.\n",
    "\n",
    "$$\\begin{align}\n",
    "f(\\mathbf{t}, T;\\beta(t)) &= \\beta\\,\\mathrm{e}^{-\\beta t_1}\\,\\left(\\prod_{i=2}^{n-1}\\beta\\,\\mathrm{e}^{-\\beta (t_{i+1}-t_i)}\\right)\\mathrm{e}^{-\\beta(T-t_n)} = \\beta^n\\,\\mathrm{e}^{-\\beta T}.\n",
    "\\end{align}\n",
    "$${#eq-homogeneous-poisson-pdf-alt}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed104381",
   "metadata": {},
   "source": [
    "## PMF for number of arrivals in time interval {#sec-nonphomogeneous-poisson-pmf}\n",
    "\n",
    "Consider now a time interval from $0$ to $T$. We break this time interval into $N$ small intervals, all of which are small enough such that, as before, $\\beta(t)$ does not change appreciably across a small interval and furthermore that the probability of more than one spike landing in a given small interval is zero. Let $\\Delta t_i = t_{i+1}-t_i$ be the width of small interval $i$. For reasons that will become clear momentarily, we choose the widths of the small intervals such that the small interval width times the rate is the same for all intervals (which we will define as $\\xi$ for convenience), or\n",
    "\n",
    "$$\\begin{align}\n",
    "\\beta(t_{i+1}) \\,\\Delta t_i = \\xi \\text{ for all }i.\n",
    "\\end{align}\n",
    "$${#eq-nonhomogenous-poisson-interval-width}\n",
    "\n",
    "Since every $\\beta(t_{i+1}) \\,\\Delta t_i$ is the same, we can write\n",
    "\n",
    "$$\\begin{align}\n",
    "\\xi = \\frac{1}{N}\\sum_{i=0}^N\\,\\beta(t_{i+1}) \\,\\Delta t_i = \\frac{1}{N}\\int_0^T\\mathrm{d}t'\\,\\beta(t') \\equiv \\frac{\\Lambda(T)}{N},\n",
    "\\end{align}\n",
    "$${#eq-nonhomogenous-poisson-interval-width}\n",
    "\n",
    "where we have considered small $\\Delta t_i$ and evaluated the Riemann sum, which is valid for unequal intervals provided the intervals are small enough. For convenience, we have defined\n",
    "\n",
    "$$\\begin{align}\n",
    "\\Lambda(T) = \\int_0^T\\mathrm{d}t'\\,\\beta(t').\n",
    "\\end{align}\n",
    "$${#eq-nonhomogenous-poisson-Lambda-1}\n",
    "\n",
    "The probability of a spike not arriving in a small interval spanning from $t_i$ to $t_{i+1}$, as we have already worked out, is\n",
    "\n",
    "$$\\begin{align}\n",
    "P(\\text{not arrived between } t_i \\text{ and } t_{i+1}) \\approx \\mathrm{e}^{-\\beta(t_{i+1})\\,\\Delta t_i} = \\mathrm{e}^{-\\xi},\n",
    "\\end{align}\n",
    "$${#eq-prob-not-arrived-interval-nonhomogeneous-poisson}\n",
    "\n",
    "such that the probability a spike does arrive in that interval is\n",
    "\n",
    "$$\\begin{align}\n",
    "P(\\text{arrived between } t_i \\text{ and } t_{i+1}) \\approx 1 - \\mathrm{e}^{-\\beta(t_{i+1}) \\Delta t_i} = 1-\\mathrm{e}^{-\\xi}.\n",
    "\\end{align}\n",
    "$${#eq-prob-arrived-interval-nonhomogeneous-poisson}\n",
    "\n",
    "Again owing the memorylessness of spiking, the probability of getting a spike in a given small interval constitutes a Bernoulli trial with probability $\\theta$ of success. Then, the probability of getting $n$ spikes in the interval $[0, T]$ is Binomially distributed,\n",
    "\n",
    "$$\\begin{align}\n",
    "f(n;N,\\xi) \\approx \\frac{N!}{n!(N-n)!}\\,(1 - \\mathrm{e}^{-\\xi})^n\\,(\\mathrm{e}^{-\\xi})^{N-n}.\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-binomial}\n",
    "\n",
    "We now consider the limit of large $N$. In this limit,\n",
    "\n",
    "$$\\begin{align}\n",
    "\\mathrm{e}^{-\\xi} \\approx 1 - \\xi,\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-small-limit-1}\n",
    "\n",
    "such that\n",
    "\n",
    "$$\\begin{align}\n",
    "f(n;N,\\xi) \\approx \\frac{N!}{n!(N-n)!}\\,\\xi^n\\,(1-\\xi)^{N-n}.\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-small-limit-2}\n",
    "\n",
    "Now, in the limit of large $N$, $N \\gg n$ such that $N - n \\approx N$. Further, we have that $N = \\Lambda(T) / \\xi$, such that\n",
    "\n",
    "$$\\begin{align}\n",
    "(1-\\xi)^{N-n} \\approx (1-\\xi)^{N} = (1-\\xi)^{\\Lambda(T)/\\xi}.\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-small-limit-3}\n",
    "\n",
    "Because\n",
    "\n",
    "$$\\begin{align}\n",
    "\\lim_{\\xi\\to 0}(1-\\xi)^{1/\\xi} = \\mathrm{e}^{-1},\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-exponential-limit}\n",
    "\n",
    "we have\n",
    "\n",
    "$$\\begin{align}\n",
    "(1-\\xi)^{N-n} \\approx (1-\\xi)^{\\Lambda(T)/\\xi} \\approx \\mathrm{e}^{-\\Lambda(T)}.\n",
    "\\end{align}\n",
    "$${#eq-nonhomogeneous-poisson-small-limit-4}\n",
    "\n",
    "Finally, also have that in the limit of large $N$,\n",
    "\n",
    "$$\\begin{align}\n",
    "N!/(N-n)! \\approx N^n = \\left(\\frac{\\Lambda(T)}{\\xi}\\right)^n.\n",
    "\\end{align}\n",
    "$${#eq-approximate-factorial-ratio}\n",
    "\n",
    "Putting everything together, we have\n",
    "\n",
    "$$\\begin{align}\n",
    "f(n;\\beta(t), T) = \\frac{(\\Lambda(T))^n}{n!}\\,\\mathrm{e}^{-\\Lambda(T)},\n",
    "\\end{align}\n",
    "$${#eq-nonhomogenous-poisson-pmf}\n",
    "\n",
    "where\n",
    "\n",
    "$$\\begin{align}\n",
    "\\Lambda(T) = \\int_0^T\\mathrm{d}t'\\,\\beta(t'),\n",
    "\\end{align}\n",
    "$${#eq-nonhomogenous-poisson-Lambda}\n",
    "\n",
    "giving the PMF for the number of arrivals of a nonhomogenous Poisson process in an interval of length $T$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86c215b9-e1bf-4e55-bc5a-c5947e9d423c",
   "metadata": {},
   "source": [
    "## Sampling arrval times of a nonhomogeneous Poisson process {#sec-nonhomogeneous-poisson-sampling}\n",
    "\n",
    "We can easily sample arrival times of a homogeneous Poisson process by noting that the interarrival times are Exponentially distributed. We keep drawing interarrival times from an Exponential distribution and cumulatively sum them to get the arrival times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d9dc786d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.30127478,  2.12937407,  3.36913533,  3.42966043,  3.75863601,\n",
       "        4.47142917,  4.5084683 ,  4.55892601,  4.83689062,  5.61171643,\n",
       "        5.79705461,  5.94193015,  9.54332889,  9.80945717, 10.56075264,\n",
       "       12.96824828, 14.4488708 , 15.07726874, 16.50660407, 16.88488701,\n",
       "       18.12896431, 18.43671618])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rng = np.random.default_rng()\n",
    "def sample_homogeneous_poisson_process(beta, T):\n",
    "    \"\"\"Draw samples of arrival times of a homogeneous Poisson process.\"\"\"\n",
    "    # Intialize arrival times\n",
    "    arrival_times = []\n",
    "\n",
    "    # Typical arrival time (Numpy likes this)\n",
    "    tau = 1 / beta\n",
    "\n",
    "    # Draw interarrival times and keep adding\n",
    "    t = rng.exponential(tau)\n",
    "    while t < T:\n",
    "        arrival_times.append(t) \n",
    "        t += rng.exponential(tau)\n",
    "\n",
    "    return np.array(arrival_times)\n",
    "\n",
    "# Give it a whirl\n",
    "sample_homogeneous_poisson_process(1, 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c2094b7",
   "metadata": {},
   "source": [
    "However, sampling arrival times of a nonhomogeneous Poisson is more challenging. The arrival time varies, so we cannot directly draw out of an Exponential distribution. To perform the sampling, we can use a **thinning method**. This works as follows. We set a fast rate $\\beta_u$ that is greater than $\\beta(t)$ everywhere on the interval $[0, T]$. We draw a time $t_1$ out of an Exponential distribution with rate $\\beta_u$. We accept $t_1$ as a sample out of the nonhomogeneous distribution with probability $\\beta(t_1)/\\beta_u$. We then draw an interarrival time $\\delta t_2$ from an Exponential distribution parametrized with rate $\\beta_u$. We accept $t_1 + \\delta t_2$ as a draw from the nonhomogeneous distribution with probability $\\beta(t_1 + \\delta t_2) / \\beta_u$. We continue in this manner until the total time of the draws meets or exceeds $T$. This is implemented in the function below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e4c07a60-e7b7-41af-bd79-a826c0cf3b1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_nhpp(beta, beta_u, T, beta_args=()):\n",
    "    \"\"\"Draw arrival times of a nonhomogeneous Poisson process.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    beta : function, call signature beta(t, *beta_args)\n",
    "        The function of the rate of arrivals as a function of time.\n",
    "    beta_u : float\n",
    "        A value of beta that is greater than beta(t) for all time.\n",
    "    T : float\n",
    "        The maximum time of observation.\n",
    "    beta_args : tuple, default ()\n",
    "        Any other arguments passed to the function `beta()`.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    output : Numpy array\n",
    "        Times for arrivals of the nonhomogeneous Poisson process.\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    .. This is an implementation of the algorithm on page 85 of \n",
    "       Simulation, 5th Ed., by Sheldon Ross, Academic Press, 2013.\n",
    "    \"\"\"\n",
    "    samples = []\n",
    "    tau = 1 / beta_u\n",
    "    t = rng.exponential(tau)\n",
    "    while t < T:\n",
    "        r = beta(t, *beta_args) * tau\n",
    "\n",
    "        if r > 1:\n",
    "            raise RuntimeError('beta_u is less than beta; sampling invalid.')\n",
    "\n",
    "        if np.random.uniform() <= r:\n",
    "            samples.append(t)\n",
    "\n",
    "        t += rng.exponential(tau)\n",
    "\n",
    "    return np.array(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa117702",
   "metadata": {},
   "source": [
    "Let's take this for a spin. We will have $\\beta(t)$ given by a sinusoid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7f4eeff2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "  <div id=\"df833f06-6b68-40a8-a32d-eafd0ce8c986\" data-root-id=\"p1002\" style=\"display: contents;\"></div>\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "(function(root) {\n",
       "  function embed_document(root) {\n",
       "  const docs_json = {\"cb21c39a-3a58-4564-9f7d-9636c3e889f0\":{\"version\":\"3.6.2\",\"title\":\"Bokeh Application\",\"roots\":[{\"type\":\"object\",\"name\":\"Figure\",\"id\":\"p1002\",\"attributes\":{\"x_range\":{\"type\":\"object\",\"name\":\"DataRange1d\",\"id\":\"p1004\"},\"y_range\":{\"type\":\"object\",\"name\":\"FactorRange\",\"id\":\"p1001\",\"attributes\":{\"factors\":[\" \"]}},\"x_scale\":{\"type\":\"object\",\"name\":\"LinearScale\",\"id\":\"p1011\"},\"y_scale\":{\"type\":\"object\",\"name\":\"CategoricalScale\",\"id\":\"p1012\"},\"title\":{\"type\":\"object\",\"name\":\"Title\",\"id\":\"p1009\"},\"renderers\":[{\"type\":\"object\",\"name\":\"GlyphRenderer\",\"id\":\"p1047\",\"attributes\":{\"name\":\"hover_glyphs\",\"data_source\":{\"type\":\"object\",\"name\":\"ColumnDataSource\",\"id\":\"p1038\",\"attributes\":{\"selected\":{\"type\":\"object\",\"name\":\"Selection\",\"id\":\"p1039\",\"attributes\":{\"indices\":[],\"line_indices\":[]}},\"selection_policy\":{\"type\":\"object\",\"name\":\"UnionRenderers\",\"id\":\"p1040\"},\"data\":{\"type\":\"map\",\"entries\":[[\"x\",[0.43348730355209997,1.020022968439517,1.4486145933121417,1.5560873962125414,1.6999567522979198,1.8215477644771276,2.2103725698175416,2.350990513300819,3.067354659987503,4.814314034302878,6.318574171873851,6.579402437947184,7.222424953914798,7.942560745494649,8.754048709888593,8.754946273712166,9.839405684399818,10.209614673017477,10.767182888183537,11.52711455427617,11.997558945166414,12.402515371401064,13.566308446872156,13.835333573664261,14.06103454740286,14.84137524718971,14.913338797301085,15.547578578427894,16.535720317820143,17.307540741475552,17.454449108565,18.524132275746794,18.783727844576518,19.30736358707713,20.76584426058035,21.235944963011608,22.474998774239488,25.102960271860542,25.608639892749174,26.405418289469964,28.181792328453465,28.296233544273402,28.4076084540452,28.466381310910318,28.746180385569915,28.782668256600797,28.80884226868342,29.065197286645194,30.56839680938774,31.92541119124079,32.665901048644564,37.86418851347071,45.071411517987684,51.25902739454848,57.749874254759476,58.622799153935304,59.189724350199285,61.23617892310767,63.060426186621115,64.66096242051452,65.20584031742067,65.72392475532668,66.09482911360782,66.23763800341054,68.36177188102792,69.48284746510902,69.49050536081464,70.22122188242912,70.64615120232038,72.08471441222729,73.24858271568166,74.24112084910284,75.64877013189698,75.73225178527417,80.2279577968169,80.52330694989038,81.69691651164307,81.91642957558851,82.36401951545945,82.94586075518573,83.03062558205639,83.0577369623233,83.14685806014131,85.25569099182937,86.06788194653033,86.12890545018982,86.98375998321586,87.15634520007852,88.11314252278233,88.13577840329326,88.45698725115678,88.95617731407408,89.18891262403787,90.93608670333154,91.33408783079685,91.94043522025979,92.04808492838187,92.29852862947945,92.91691154166985,93.14193517925932,95.9782404424489,96.94372090251424,101.71545491672573,102.53344989169736,105.17088503771588,106.2995883338456,106.38041572173107,108.47960978411855,110.10900664064644,110.17972782731226,119.15753557774833,124.06714494216325,126.65079925169668,126.86073007859078,127.89114485118952,128.69018964218503,128.76599433876893,129.3459610205387,129.4307433024083,130.45050636207654,131.5526696352911,131.8638964701547,132.39341562739207,133.58821315522465,133.7564116078935,133.77578743701247,134.10318751067294,134.20319502159734,134.3234079832269,134.69088195518043,135.02381210337026,135.13351988413714,135.3879779392867,136.4319556151832,136.6777026008587,137.23457883387474,138.0798863538993,138.2180863045773,139.21538117536485,139.99197625855507,141.08614192206878,141.1856050879829,141.3309101608946,142.05695660193905,142.0624397903642,142.30510462034098,143.20907580989163,143.45984879815683,143.55203740823748,144.40662738408085,144.54675505019495,145.44925435098952,145.46430563852672,146.08427848362078,146.10443385502623,147.2869845979768,147.54661424506367,148.77111142043918,149.07748772140667,149.784143449754,149.8654475000174,151.33388406221857,152.4134843261774,153.09245140219528,156.01908883372752,156.62563087004682,157.7655142731132,160.3506313041561,160.90705842860973,162.98316388057066,163.92676799908267,166.93017807511595,170.3868461680691,171.56540844437953,181.36819521274148,181.80111849053145,184.33826237356158,185.25970914741373,185.79499563389416,186.8451947656632,187.106409594774,187.64808658305625,188.1006084482192,190.53249988472064,191.1068561021547,192.41079739512176,192.87040314891266,193.38206168947332,194.6775786428782,195.2975768954851,195.94373084423415,197.979770944143,198.14291982779295,199.4172562081756,201.41253356302698,201.8020480584455,201.99981783885,202.53068615378734,202.56741124022318,202.98075257520054,203.3481823766654,203.361604379935,204.68291980868344,204.96907077057088,205.22976520148873,205.24675299143718,205.64472423038265,205.87351525161085,205.97893752614723,206.00910113359788,206.7812460076915,206.9293885921992,207.1368232904111,207.14411223820923,207.8959900559988,207.91874767956452,208.1022257224743,208.18083550337505,208.6335490254542,208.75086601612466,209.68881975014804,210.2937934583255,210.6437781500967,211.1522248715752,211.22831174278937,212.13745431307993,213.68939503192897,213.93636817360638,214.5115422528422,215.56690708397142,215.7654556879544,216.09188646365305,216.5730696747975,216.7223037770032,216.86394222514866,217.37684661138792,218.07806497384067,218.78904103467988,221.86926583598978,223.15440114977727,226.30606991604,226.8643078722171,228.14009212305112,230.10545375897377,232.1359096284597,237.12253209247206,238.80792698678732,241.7960159464189,243.11819601495444,244.91766232208343,245.2931261264572,245.5275318405724,246.0776708128636,246.1631563547306,247.20239064888506,247.7252404677987,247.82114712058674,250.01123244255527,250.6348350368685,252.11792747838783,252.4603258622727,253.80361591316608,253.90173035665765,254.02538071360237,254.20526901486718,254.62938563148236,255.82955977324283,256.474759857216,256.7237540700904,256.78193843106453,257.09999568946864,257.327003901089,257.8850657388703,258.3047079891161,258.3135450408863,258.4311886232356,259.1561894164988,259.2221628412938,259.4852500012693,260.09131343625626,260.0915448013578,260.31687754069435,260.7006900071049,261.7925942817908,261.9726002947012,262.2022444119028,262.28485179445727,262.8729720345592,263.3146458222646,266.09854741017335,266.8071562965696,267.8043858442465,268.61283941508077,268.704057702793,268.76629995586285,269.37626061860624,269.392836064503,269.3933728515582,269.51678656234293,270.57721220915016,271.38499009582665,271.7807888331551,273.456842713199,273.8810691227567,273.8894917222881,275.96681330201045,276.16224106906753,276.4937137984437,277.22617271504186,277.3663082788902,278.7935715082956,279.4402995302093,280.3508908788291,280.5830321449451,281.17360171795326,281.53530952819364,284.1632150624865,284.84874731176023,285.0964950194949,285.3651204393486,285.9079434388548,288.40474081441465,289.0123510007157,289.1397672834361,289.6383969751978,295.3146953861681,305.00206869760365,306.15815179194135,307.64212175548056,307.94015128948195,309.1826562005883,312.6423140485122,314.456390211169,315.3125428635959,315.99870955976917,316.442736995091,317.0248198786499,319.0099653919986,319.22629036580986,319.4520852920668,319.7558342837102,319.99599651359034,320.2220152581898,320.7253862663716,320.99455069998834,322.1550109664916,323.19393455506383,324.35173262346564,324.9447289074548,325.041329972832,325.1623393635105,325.7199617915692,326.4483106288628,327.20152651006805,327.43952633535855,328.62071477316454,328.8989183402467,329.166229399594,329.34284348056997,329.73983390240653,329.9484786550104,330.24239919943346,330.83580910039467,331.58958445653207,331.96529551722233,333.5297055790368,334.07248732163316,334.1299090655457,334.1603947359089,334.52482538066494,335.486334044211,337.3212072108893,337.32616529748753,337.96667910540975,338.2019767312737,338.37947040340856,338.70966876609464,339.2969177808627,339.70839833120624,340.1148326382929,341.6660615224881,341.74721802654364,342.8761523391203,343.3256056273908,343.99669944737497,345.6992189690275,346.3663726535676,347.0769457114124,347.19168396653816,347.650229272251,347.7883540705185,350.9948789136173,351.9284856521489,352.56128889201614,353.3788740364982,355.97063754980127,371.0227437250663,372.4293606557551,374.19333290506574,374.971749263094,375.35488907767984,377.16111136099414,377.1929055908719,377.2667006683954,377.7392950547226,378.2753817974446,379.83567497214233,380.94305754237735,381.4718002862374,382.20270532455834,383.74883833627604,384.38829938071524,384.77834964929076,384.936306330201,386.82354491776294,387.3766039588497,387.7118709745609,387.72216020403266,388.40117272767105,390.02060398911317,390.4879712145096,390.7243970570545,391.35460556658006,392.2968595798365,392.9823203865677,393.4885406383933,394.07217886368244,394.5847012848366,394.748877961789,396.42428623821417,396.5976576272508,396.6328656453799,397.03584438492425,397.40668835360276,397.87109078274256,398.95816195356383,399.2355135485734,399.8720079169583,400.0621790474553,401.27731522636384,401.4249164031407,401.5840311687741,402.48245092396115,403.17855906261855,403.99919631403156,404.55438245640886,407.1015804616998,407.51136655533253,407.9217939651569,408.0071000463723,408.0319834177023,408.0535531813215,408.36130715225715,413.69078015285424,413.9979145435441,415.1774582894793,422.1287980956274,422.33487478227426,429.8365830663819,431.2570955522614,431.7631618371771,435.866293839266,436.94697971814577,437.12523460641506,437.2491607249312,437.53210776603913,438.40996597921935,438.45934173670565,439.12362492462955,439.5716056884266,444.0064277224646,444.389227491881,444.5900060538414,445.4490292181633,446.021055588688,446.3283918068109,446.5679898446221,446.70162895415865,447.41122951911274,447.6404060195616,447.7471118066751,447.80162270667677,448.16371586163683,449.36558816570306,450.01009015212964,450.2880058349035,450.775092634337,451.20979839419635,451.503815850616,453.4307753133725,453.87004642518895,455.0899822989413,455.2147469840774,455.2883501524313,455.31225789179683,455.8025757898682,456.01370473452056,456.7121804018013,456.81376944381935,456.94357206133594,457.2390635774616,458.10155398227596,458.34034157960417,458.359873922517,459.4284260920987,459.74188087701464,460.54158173150563,461.28435241693774,462.6751959193977,463.6820421495,464.4904024514394,464.79173611166743,464.97543261186223,465.86052681626825,465.87542315020727,466.146627636055,467.0150270189673,467.3198903375312,468.0494755291419,469.71649504716555,470.5396629055737,471.6651584797392,472.72238539514944,474.3723142680899,475.0644474832772,476.5478724159839,485.7479325387224,487.34947526499275,498.3310697515958,499.7114837977794]],[\"__dummy_cat\",[\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \"]],[\"cat\",[\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \"]],[\"__label\",[\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \",\" \"]]]}}},\"view\":{\"type\":\"object\",\"name\":\"CDSView\",\"id\":\"p1048\",\"attributes\":{\"filter\":{\"type\":\"object\",\"name\":\"AllIndices\",\"id\":\"p1049\"}}},\"glyph\":{\"type\":\"object\",\"name\":\"Scatter\",\"id\":\"p1044\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"cat\"},\"size\":{\"type\":\"value\",\"value\":68.75},\"angle\":{\"type\":\"value\",\"value\":1.5707963267948966},\"line_color\":{\"type\":\"field\",\"field\":\"cat\",\"transform\":{\"type\":\"object\",\"name\":\"CategoricalColorMapper\",\"id\":\"p1037\",\"attributes\":{\"palette\":[\"#1f77b3\",\"#ff7e0e\",\"#2ba02b\",\"#d62628\",\"#9367bc\",\"#8c564b\",\"#e277c1\",\"#7e7e7e\",\"#bcbc21\",\"#16bdcf\",\"#3a0182\",\"#004201\",\"#0fffa8\",\"#5d003f\",\"#bcbcff\",\"#d8afa1\",\"#b80080\",\"#004d52\",\"#6b6400\",\"#7c0100\",\"#6026ff\",\"#ffff9a\",\"#564964\",\"#8cb893\",\"#93fbff\",\"#018267\",\"#90ff00\",\"#8200a0\",\"#ac8944\",\"#5b3400\",\"#ffbff2\",\"#ff6e75\",\"#798cff\",\"#dd00ff\",\"#505646\",\"#004489\",\"#ffbf60\",\"#ff018c\",\"#bdc8cf\",\"#af97b5\",\"#b65600\",\"#017000\",\"#cd87ff\",\"#1cd646\",\"#bfebc3\",\"#7997b5\",\"#a56089\",\"#6e8956\",\"#bc7c75\",\"#8a2844\",\"#00acff\",\"#8ed4ff\",\"#4b6d77\",\"#00d4b1\",\"#9300f2\",\"#8a9500\",\"#5d5b9e\",\"#fddfba\",\"#00939e\",\"#ffdb00\",\"#00aa79\",\"#520067\",\"#000091\",\"#0a5d3d\",\"#a5e275\",\"#623b41\",\"#c6c689\",\"#ff9eb5\",\"#cd4f6b\",\"#ff07d6\",\"#8a3a05\",\"#7e3d70\",\"#ff4901\",\"#602ba5\",\"#1c00ff\",\"#e6dfff\",\"#aa3baf\",\"#d89c00\",\"#a3a39e\",\"#3f69ff\",\"#46490c\",\"#7b6985\",\"#6b978c\",\"#ff9a75\",\"#835bff\",\"#7c6b46\",\"#80b654\",\"#bc0049\",\"#fd93ff\",\"#5d0018\",\"#89d1d1\",\"#9c8cd3\",\"#da6d42\",\"#8a5700\",\"#3b5069\",\"#4b6b3b\",\"#edcfd8\",\"#cfedff\",\"#aa1500\",\"#dfff4f\",\"#ff2a56\",\"#d1499e\",\"#707cb8\",\"#598000\",\"#00e4fd\",\"#774b95\",\"#67d48c\",\"#3d3a72\",\"#ac413f\",\"#d6a166\",\"#c169cd\",\"#69595d\",\"#87aced\",\"#a0a569\",\"#d1aae6\",\"#870062\",\"#00fddb\",\"#672818\",\"#b342ff\",\"#0e59c4\",\"#168742\",\"#90d300\",\"#cd7900\",\"#f959ff\",\"#5b7466\",\"#8eaeb3\",\"#9c7c8c\",\"#4600c6\",\"#6b4d2d\",\"#a56d46\",\"#9e8972\",\"#a8afca\",\"#cd8ca7\",\"#00fd64\",\"#917900\",\"#ff62a1\",\"#f4ffd8\",\"#018cf0\",\"#13aca0\",\"#5b2d59\",\"#89859e\",\"#cfccba\",\"#d4afc4\",\"#dbdd6d\",\"#cffff4\",\"#006485\",\"#006962\",\"#a84167\",\"#2d97c4\",\"#a874ff\",\"#26ba5d\",\"#57b600\",\"#caffa7\",\"#a379aa\",\"#ffbc93\",\"#89e2c1\",\"#0fc8ff\",\"#d400c4\",\"#626d89\",\"#69858e\",\"#4b4d52\",\"#aa6067\",\"#79b5d4\",\"#2b5916\",\"#9a0024\",\"#bdd1f2\",\"#896e67\",\"#69a56b\",\"#855467\",\"#aecdba\",\"#87997e\",\"#cadb00\",\"#9a0390\",\"#ebbc1a\",\"#eb9cd1\",\"#70006e\",\"#b1a131\",\"#ca6b93\",\"#4146a3\",\"#e48c89\",\"#d44400\",\"#c68aca\",\"#b69597\",\"#d41f75\",\"#724bcc\",\"#674d00\",\"#672138\",\"#38564f\",\"#6ebaaa\",\"#853a31\",\"#a5d397\",\"#b8af8e\",\"#d8e4df\",\"#aa00df\",\"#cac1db\",\"#ffdf8c\",\"#e2524d\",\"#66696e\",\"#ff001c\",\"#522d72\",\"#4d906b\",\"#a86d11\",\"#ff9e26\",\"#5ea3af\",\"#c88556\",\"#915997\",\"#a3a1ff\",\"#fdbaba\",\"#242a87\",\"#dbe6a8\",\"#97f2a7\",\"#6793d6\",\"#ba5b3f\",\"#3a5d91\",\"#364f2f\",\"#267c95\",\"#89959a\",\"#cfb356\",\"#004664\",\"#5e5d2f\",\"#8e8e41\",\"#ac3f13\",\"#69953b\",\"#a13d85\",\"#bfb6ba\",\"#acc667\",\"#6469cf\",\"#91af00\",\"#2be2da\",\"#016e36\",\"#ff7952\",\"#42807e\",\"#4fe800\",\"#995428\",\"#5d0a00\",\"#a30057\",\"#0c8700\",\"#5982a7\",\"#ffebfb\",\"#4b6901\",\"#8775d4\",\"#e6c6ff\",\"#a5ffda\",\"#d86e77\",\"#df014b\",\"#69675b\",\"#776ba1\",\"#7e8067\",\"#594685\",\"#0000ca\",\"#7c002a\",\"#97ff72\",\"#b5e2e1\",\"#db52c8\",\"#777734\",\"#57bd8e\"],\"factors\":[\" \"]}}},\"line_alpha\":{\"type\":\"value\",\"value\":0.3},\"fill_color\":{\"type\":\"field\",\"field\":\"cat\",\"transform\":{\"type\":\"object\",\"name\":\"CategoricalColorMapper\",\"id\":\"p1036\",\"attributes\":{\"palette\":[\"#1f77b3\",\"#ff7e0e\",\"#2ba02b\",\"#d62628\",\"#9367bc\",\"#8c564b\",\"#e277c1\",\"#7e7e7e\",\"#bcbc21\",\"#16bdcf\",\"#3a0182\",\"#004201\",\"#0fffa8\",\"#5d003f\",\"#bcbcff\",\"#d8afa1\",\"#b80080\",\"#004d52\",\"#6b6400\",\"#7c0100\",\"#6026ff\",\"#ffff9a\",\"#564964\",\"#8cb893\",\"#93fbff\",\"#018267\",\"#90ff00\",\"#8200a0\",\"#ac8944\",\"#5b3400\",\"#ffbff2\",\"#ff6e75\",\"#798cff\",\"#dd00ff\",\"#505646\",\"#004489\",\"#ffbf60\",\"#ff018c\",\"#bdc8cf\",\"#af97b5\",\"#b65600\",\"#017000\",\"#cd87ff\",\"#1cd646\",\"#bfebc3\",\"#7997b5\",\"#a56089\",\"#6e8956\",\"#bc7c75\",\"#8a2844\",\"#00acff\",\"#8ed4ff\",\"#4b6d77\",\"#00d4b1\",\"#9300f2\",\"#8a9500\",\"#5d5b9e\",\"#fddfba\",\"#00939e\",\"#ffdb00\",\"#00aa79\",\"#520067\",\"#000091\",\"#0a5d3d\",\"#a5e275\",\"#623b41\",\"#c6c689\",\"#ff9eb5\",\"#cd4f6b\",\"#ff07d6\",\"#8a3a05\",\"#7e3d70\",\"#ff4901\",\"#602ba5\",\"#1c00ff\",\"#e6dfff\",\"#aa3baf\",\"#d89c00\",\"#a3a39e\",\"#3f69ff\",\"#46490c\",\"#7b6985\",\"#6b978c\",\"#ff9a75\",\"#835bff\",\"#7c6b46\",\"#80b654\",\"#bc0049\",\"#fd93ff\",\"#5d0018\",\"#89d1d1\",\"#9c8cd3\",\"#da6d42\",\"#8a5700\",\"#3b5069\",\"#4b6b3b\",\"#edcfd8\",\"#cfedff\",\"#aa1500\",\"#dfff4f\",\"#ff2a56\",\"#d1499e\",\"#707cb8\",\"#598000\",\"#00e4fd\",\"#774b95\",\"#67d48c\",\"#3d3a72\",\"#ac413f\",\"#d6a166\",\"#c169cd\",\"#69595d\",\"#87aced\",\"#a0a569\",\"#d1aae6\",\"#870062\",\"#00fddb\",\"#672818\",\"#b342ff\",\"#0e59c4\",\"#168742\",\"#90d300\",\"#cd7900\",\"#f959ff\",\"#5b7466\",\"#8eaeb3\",\"#9c7c8c\",\"#4600c6\",\"#6b4d2d\",\"#a56d46\",\"#9e8972\",\"#a8afca\",\"#cd8ca7\",\"#00fd64\",\"#917900\",\"#ff62a1\",\"#f4ffd8\",\"#018cf0\",\"#13aca0\",\"#5b2d59\",\"#89859e\",\"#cfccba\",\"#d4afc4\",\"#dbdd6d\",\"#cffff4\",\"#006485\",\"#006962\",\"#a84167\",\"#2d97c4\",\"#a874ff\",\"#26ba5d\",\"#57b600\",\"#caffa7\",\"#a379aa\",\"#ffbc93\",\"#89e2c1\",\"#0fc8ff\",\"#d400c4\",\"#626d89\",\"#69858e\",\"#4b4d52\",\"#aa6067\",\"#79b5d4\",\"#2b5916\",\"#9a0024\",\"#bdd1f2\",\"#896e67\",\"#69a56b\",\"#855467\",\"#aecdba\",\"#87997e\",\"#cadb00\",\"#9a0390\",\"#ebbc1a\",\"#eb9cd1\",\"#70006e\",\"#b1a131\",\"#ca6b93\",\"#4146a3\",\"#e48c89\",\"#d44400\",\"#c68aca\",\"#b69597\",\"#d41f75\",\"#724bcc\",\"#674d00\",\"#672138\",\"#38564f\",\"#6ebaaa\",\"#853a31\",\"#a5d397\",\"#b8af8e\",\"#d8e4df\",\"#aa00df\",\"#cac1db\",\"#ffdf8c\",\"#e2524d\",\"#66696e\",\"#ff001c\",\"#522d72\",\"#4d906b\",\"#a86d11\",\"#ff9e26\",\"#5ea3af\",\"#c88556\",\"#915997\",\"#a3a1ff\",\"#fdbaba\",\"#242a87\",\"#dbe6a8\",\"#97f2a7\",\"#6793d6\",\"#ba5b3f\",\"#3a5d91\",\"#364f2f\",\"#267c95\",\"#89959a\",\"#cfb356\",\"#004664\",\"#5e5d2f\",\"#8e8e41\",\"#ac3f13\",\"#69953b\",\"#a13d85\",\"#bfb6ba\",\"#acc667\",\"#6469cf\",\"#91af00\",\"#2be2da\",\"#016e36\",\"#ff7952\",\"#42807e\",\"#4fe800\",\"#995428\",\"#5d0a00\",\"#a30057\",\"#0c8700\",\"#5982a7\",\"#ffebfb\",\"#4b6901\",\"#8775d4\",\"#e6c6ff\",\"#a5ffda\",\"#d86e77\",\"#df014b\",\"#69675b\",\"#776ba1\",\"#7e8067\",\"#594685\",\"#0000ca\",\"#7c002a\",\"#97ff72\",\"#b5e2e1\",\"#db52c8\",\"#777734\",\"#57bd8e\"],\"factors\":[\" \"]}}},\"fill_alpha\":{\"type\":\"value\",\"value\":0.3},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.3},\"marker\":{\"type\":\"value\",\"value\":\"dash\"}}},\"nonselection_glyph\":{\"type\":\"object\",\"name\":\"Scatter\",\"id\":\"p1045\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"cat\"},\"size\":{\"type\":\"value\",\"value\":68.75},\"angle\":{\"type\":\"value\",\"value\":1.5707963267948966},\"line_color\":{\"type\":\"field\",\"field\":\"cat\",\"transform\":{\"id\":\"p1037\"}},\"line_alpha\":{\"type\":\"value\",\"value\":0.1},\"fill_color\":{\"type\":\"field\",\"field\":\"cat\",\"transform\":{\"id\":\"p1036\"}},\"fill_alpha\":{\"type\":\"value\",\"value\":0.1},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.1},\"marker\":{\"type\":\"value\",\"value\":\"dash\"}}},\"muted_glyph\":{\"type\":\"object\",\"name\":\"Scatter\",\"id\":\"p1046\",\"attributes\":{\"x\":{\"type\":\"field\",\"field\":\"x\"},\"y\":{\"type\":\"field\",\"field\":\"cat\"},\"size\":{\"type\":\"value\",\"value\":68.75},\"angle\":{\"type\":\"value\",\"value\":1.5707963267948966},\"line_color\":{\"type\":\"field\",\"field\":\"cat\",\"transform\":{\"id\":\"p1037\"}},\"line_alpha\":{\"type\":\"value\",\"value\":0.2},\"fill_color\":{\"type\":\"field\",\"field\":\"cat\",\"transform\":{\"id\":\"p1036\"}},\"fill_alpha\":{\"type\":\"value\",\"value\":0.2},\"hatch_alpha\":{\"type\":\"value\",\"value\":0.2},\"marker\":{\"type\":\"value\",\"value\":\"dash\"}}}}}],\"toolbar\":{\"type\":\"object\",\"name\":\"Toolbar\",\"id\":\"p1010\",\"attributes\":{\"tools\":[{\"type\":\"object\",\"name\":\"PanTool\",\"id\":\"p1023\"},{\"type\":\"object\",\"name\":\"WheelZoomTool\",\"id\":\"p1024\",\"attributes\":{\"renderers\":\"auto\"}},{\"type\":\"object\",\"name\":\"BoxZoomTool\",\"id\":\"p1025\",\"attributes\":{\"overlay\":{\"type\":\"object\",\"name\":\"BoxAnnotation\",\"id\":\"p1026\",\"attributes\":{\"syncable\":false,\"line_color\":\"black\",\"line_alpha\":1.0,\"line_width\":2,\"line_dash\":[4,4],\"fill_color\":\"lightgrey\",\"fill_alpha\":0.5,\"level\":\"overlay\",\"visible\":false,\"left\":{\"type\":\"number\",\"value\":\"nan\"},\"right\":{\"type\":\"number\",\"value\":\"nan\"},\"top\":{\"type\":\"number\",\"value\":\"nan\"},\"bottom\":{\"type\":\"number\",\"value\":\"nan\"},\"left_units\":\"canvas\",\"right_units\":\"canvas\",\"top_units\":\"canvas\",\"bottom_units\":\"canvas\",\"handles\":{\"type\":\"object\",\"name\":\"BoxInteractionHandles\",\"id\":\"p1032\",\"attributes\":{\"all\":{\"type\":\"object\",\"name\":\"AreaVisuals\",\"id\":\"p1031\",\"attributes\":{\"fill_color\":\"white\",\"hover_fill_color\":\"lightgray\"}}}}}}}},{\"type\":\"object\",\"name\":\"SaveTool\",\"id\":\"p1033\"},{\"type\":\"object\",\"name\":\"ResetTool\",\"id\":\"p1034\"},{\"type\":\"object\",\"name\":\"HelpTool\",\"id\":\"p1035\"}]}},\"toolbar_location\":\"above\",\"left\":[{\"type\":\"object\",\"name\":\"CategoricalAxis\",\"id\":\"p1018\",\"attributes\":{\"ticker\":{\"type\":\"object\",\"name\":\"CategoricalTicker\",\"id\":\"p1019\"},\"formatter\":{\"type\":\"object\",\"name\":\"CategoricalTickFormatter\",\"id\":\"p1020\"},\"major_label_policy\":{\"type\":\"object\",\"name\":\"AllLabels\",\"id\":\"p1021\"}}}],\"below\":[{\"type\":\"object\",\"name\":\"LinearAxis\",\"id\":\"p1013\",\"attributes\":{\"ticker\":{\"type\":\"object\",\"name\":\"BasicTicker\",\"id\":\"p1014\",\"attributes\":{\"mantissas\":[1,2,5]}},\"formatter\":{\"type\":\"object\",\"name\":\"BasicTickFormatter\",\"id\":\"p1015\"},\"axis_label\":\"x\",\"major_label_policy\":{\"type\":\"object\",\"name\":\"AllLabels\",\"id\":\"p1016\"}}}],\"center\":[{\"type\":\"object\",\"name\":\"Grid\",\"id\":\"p1017\",\"attributes\":{\"axis\":{\"id\":\"p1013\"}}},{\"type\":\"object\",\"name\":\"Grid\",\"id\":\"p1022\",\"attributes\":{\"dimension\":1,\"axis\":{\"id\":\"p1018\"},\"grid_line_color\":null}}],\"frame_width\":375,\"frame_height\":275}}]}};\n",
       "  const render_items = [{\"docid\":\"cb21c39a-3a58-4564-9f7d-9636c3e889f0\",\"roots\":{\"p1002\":\"df833f06-6b68-40a8-a32d-eafd0ce8c986\"},\"root_ids\":[\"p1002\"]}];\n",
       "  void root.Bokeh.embed.embed_items_notebook(docs_json, render_items);\n",
       "  }\n",
       "  if (root.Bokeh !== undefined) {\n",
       "    embed_document(root);\n",
       "  } else {\n",
       "    let attempts = 0;\n",
       "    const timer = setInterval(function(root) {\n",
       "      if (root.Bokeh !== undefined) {\n",
       "        clearInterval(timer);\n",
       "        embed_document(root);\n",
       "      } else {\n",
       "        attempts++;\n",
       "        if (attempts > 100) {\n",
       "          clearInterval(timer);\n",
       "          console.log(\"Bokeh: ERROR: Unable to run BokehJS code because BokehJS library is missing\");\n",
       "        }\n",
       "      }\n",
       "    }, 10, root)\n",
       "  }\n",
       "})(window);"
      ],
      "application/vnd.bokehjs_exec.v0+json": ""
     },
     "metadata": {
      "application/vnd.bokehjs_exec.v0+json": {
       "id": "p1002"
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "beta = lambda t: 1.15 + np.sin(t/10)\n",
    "arrival_times = sample_nhpp(beta, 2.15, 500)\n",
    "\n",
    "# Plot the arrivals\n",
    "p = iqplot.strip(arrival_times, marker='dash', marker_kwargs=dict(alpha=0.3))\n",
    "bokeh.io.show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e0f2e2a",
   "metadata": {},
   "source": [
    "Nice!\n",
    "\n",
    "If we want to speed things up, we can break the interval up into smaller segments and for each interval we make sure that we choose a $\\beta_u$ specific to that interval that is close to the maximum $\\beta(t)$ on that interval. The code (which also includes other speed enhancements) to achieve this is below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7ae07d98-28cc-4b92-afe0-03a83e56b2c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_nhpp(beta, beta_j, t_j, beta_args=()):\n",
    "    \"\"\"Draw arrival times of a nonhomogeneous Poisson process.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    beta : function, call signature beta(t, *beta_args)\n",
    "        The function of the rate of arrivals as a function of time.\n",
    "    beta_j : scalar or array_like\n",
    "        If scalar, a value of beta that is greater than beta(t)\n",
    "        for all time. If array_like, then beta_j[j] > beta(t) for\n",
    "        all times between t_j[j-1] and t_j[j].\n",
    "    t_j : scalar or array_like\n",
    "        If scalar, the maximum time of observation. If array_like, must\n",
    "        be the same length of `beta_j`. beta_j[j] is the value \n",
    "        of the the upper bound of the rate for the interval between\n",
    "        t[j-1] and t[j].\n",
    "\n",
    "    beta_args : tuple, default ()\n",
    "        Any other arguments passed to beta.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    output : Numpy array\n",
    "        Times for arrivals of the nonhomogeneous Poisson process.\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    .. This is an implementation of the algorithm on page 86 of \n",
    "       Simulation, 5th Ed., by Sheldon Ross, Academic Press, 2013.\n",
    "    \"\"\"\n",
    "    # Convert scalar inputs to arrays\n",
    "    if np.isscalar(beta_j):\n",
    "        beta_j = np.array([beta_j])\n",
    "    if np.isscalar(t_j):\n",
    "        t_j = np.array([t_j])\n",
    "\n",
    "    # Make sure dimensions match\n",
    "    if len(beta_j) != len(t_j):\n",
    "        raise RuntimeError(\n",
    "            f'`beta_j` is length {len(beta_j)} '\n",
    "            f'and `t_j` is length {len(t_j)}. '\n",
    "            'They must have the same length.'\n",
    "        )\n",
    "\n",
    "    return _sample_nhpp(beta, beta_j, t_j, beta_args)\n",
    "\n",
    "                                                  \n",
    "def _sample_nhpp(beta, beta_j, t_j, beta_args=()):\n",
    "    \"\"\"\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    beta : function, call signature beta(t, *beta_args)\n",
    "        The function of the rate of arrivals as a function of time.\n",
    "    beta_j : Numpy array\n",
    "        If scalar, a value of beta that is greater than beta(t)\n",
    "        for all time. If array_like, then beta_j[j] > beta(t) for\n",
    "        all times between t_j[j-1] and t_j[j].\n",
    "    t_j : Numpy array\n",
    "        Must be the same length of `beta_j`. beta_j[j] is the value \n",
    "        of the the upper bound of the rate for the interval between\n",
    "        t[j-1] and t[j].\n",
    "\n",
    "    beta_args : tuple, default ()\n",
    "        Any other arguments passed to beta.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    output : Numpy array\n",
    "        Times for arrivals of the nonhomogeneous Poisson process.\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    .. This is an implementation of the algorithm on page 86 of \n",
    "       Simulation, 5th Ed., by Sheldon Ross, Academic Press, 2013.\n",
    "    \"\"\"\n",
    "    # Number of samples to take before concatenating arrays\n",
    "    n_samples = 1000\n",
    "\n",
    "    # Initializations\n",
    "    t = 0.0  # time\n",
    "    j = 0    # Index in beta_j and t_j arrays\n",
    "    i = 0    # index in sample array\n",
    "    n = 0    # total number of samples drawn\n",
    "    x_from_boundary = False  # If we've hit a boundary of\n",
    "\n",
    "    # Array for storing subtrajectory\n",
    "    samples = np.empty(n_samples, dtype=float)\n",
    "\n",
    "    # Output array for all samples\n",
    "    samples_output = np.array([], dtype=float)\n",
    "\n",
    "    # Loop until done (we exceed final time point)\n",
    "    not_done = True\n",
    "    while not_done:\n",
    "        # Take samples until we fill array\n",
    "        # We do it this way for speed to avoid list append operations\n",
    "        while not_done and i < n_samples:\n",
    "            if x_from_boundary:\n",
    "                x = (x - t_j[j] + t) * beta_j[j] / beta_j[j+1]\n",
    "                t = t_j[j]\n",
    "                j += 1\n",
    "            else:\n",
    "                x = np.random.exponential(1 / beta_j[j])\n",
    "\n",
    "            if t + x > t_j[j]:\n",
    "                # If we got here, we went past the edge of this interval\n",
    "                if j == len(t_j) - 1:\n",
    "                    not_done = False\n",
    "                else:\n",
    "                    x_from_boundary = True\n",
    "            else:\n",
    "                t += x\n",
    "                x_from_bounday = False\n",
    "                if np.random.uniform() <= beta(t, *beta_args) / beta_j[j]:\n",
    "                    samples[i] = t\n",
    "                    i += 1\n",
    "                    n += 1\n",
    "        samples_output = np.concatenate((samples_output, samples))\n",
    "        i = 0\n",
    "\n",
    "    return np.array(samples_output[:n])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5a3e682-39bf-4d63-91b8-b450a82cd2bc",
   "metadata": {},
   "source": [
    "## Computing environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "681b1e7d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python implementation: CPython\n",
      "Python version       : 3.12.9\n",
      "IPython version      : 9.1.0\n",
      "\n",
      "numpy     : 2.1.3\n",
      "iqplot    : 0.3.7\n",
      "bokeh     : 3.6.2\n",
      "jupyterlab: 4.4.2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -v -p numpy,iqplot,bokeh,jupyterlab"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}