Example-SSA-delay-telegraph.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing SSA for delayed telegraph model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we show how to use the SSA for the delayed telegraph to produce probability distributions and histograms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"./DSSA-delayed-telegraph.jl\"); using .SSAdelayedtelegraph;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import some necessary packages for plotting and saving data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using PyPlot, DelimitedFiles, LaTeXStrings;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's define:\n",
    "- `par_set`: the set of parameter for which we want to simulate in order $\\{\\rho,\\sigma_{off},\\sigma_{on}\\}$.\n",
    "- `τ`: the deterministic delay time of the transition from nascent to mature RNA.\n",
    "- `total_time`: the total time of the simulation.\n",
    "- `sp`: the sampling period of the SSA (how often the state vector is stored). Total number of data points stored will be $=$`total_time/sp`.\n",
    "- `sims`: the number of independent DSSA trajectories. If equal to 1 only a single DSSA trajectory will be produced.\n",
    "\n",
    "The initial conditions of the SSA are to start in gene state $U^{**}$ with zero $N$ and $M$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "par_set = [2.11,0.609,0.0282]; # params should be float64.\n",
    "τ = 150.0; # must be float64.\n",
    "total_time = 1000000.0; # must be float64.\n",
    "sp = 1.0; # must be float64.\n",
    "sims = 1; # take a single trajectory. Must be integer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.397758 seconds (31.94 M allocations: 800.549 MiB, 7.30% gc time)\n"
     ]
    }
   ],
   "source": [
    "# simulate\n",
    "@time data=SSAdt(sims, par_set, τ, total_time, sp);\n",
    "\n",
    "# take the individual trajectories of each species.\n",
    "n_traj = data[2,1,:]; g_0 = data[1,1,:];\n",
    "g_1 = 1 .- data[1,1,:];\n",
    "# save n and m traj to file\n",
    "writedlm(\"delayTMtest-bimod-N.csv\",n_traj[1:end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAABMvklEQVR4nO3de3xTVbo38N9O0qS3tNzsDQoUBBULDoIyIEdwBFS8zeAFFQUvc1G8dTiKOvi+Vkcp4xkZVEZ8x6PImTkIR0VFjyJ4oeiACC0IwigqtyrUgkCb3nJd7x/p3knatCTNTnZ29u/7+fDR7u4mqztp8uRZz3qWJIQQICIiItIpk9YDICIiIooFgxkiIiLSNQYzREREpGsMZoiIiEjXGMwQERGRrjGYISIiIl1jMENERES6ZtF6APHm8/lw6NAh2O12SJKk9XCIiIgoAkIIOBwOFBUVwWTqOveS8sHMoUOHUFxcrPUwiIiIqBtqamrQr1+/Ls9J+WDGbrcD8F+MnJwcjUdDREREkWhoaEBxcbHyPt6VlA9m5KmlnJwcBjNEREQ6E0mJCAuAiYiISNcYzBAREZGuMZghIiIiXUv5mplIeb1euN1urYcRd2lpaTCbzVoPg4iISDWGD2aEEKitrcWJEye0HkrC9OjRAwUFBey7Q0REKcHwwYwcyOTl5SEzMzOl3+CFEGhubkZdXR0AoLCwUOMRERERxc7QwYzX61UCmd69e2s9nITIyMgAANTV1SEvL49TTkREpHuGLgCWa2QyMzM1Hkliyb+vEWqEiIgo9Rk6mJGl8tRSOEb7fYmIKLUxmCEiIiJdYzBDREREusZghoiIiHSNwYyOnX/++ZAkCa+88krI8eeeew55eXkajYqIiCIhhECLy6v1MFICgxmdEkJg+/btKCwsxOuvvx7yverqapx99tkajYyIiCIx780vccb/XYOvahu0HoruGbrPTDhCCLS4Ex8pZ6SZo1pl9M0338DhcGDBggW4//770dzcrCy5rqqqwtSpU+M1VCIiUsHyzQcBAH/9+Ds8e/1IjUejbwxm2mlxezHs/76f8Pvd/dhFyLRG/nBUVVUhPT0dv/71r/HHP/4R7733Hq666io4nU7s2rULDz/8cBxHS0REavEJofUQdI/TTDpVXV2NESNGwGq14le/+hVee+01AMCOHTvgdrsxatQojUdIRESREAxmYsbMTDsZaWbsfuwiTe43GlVVVUpdzLRp0zBt2jQ4nU5UVVWhV69eGDhwYBxGSUREavP6GMzEKmkyMxUVFZAkCWVlZcoxIQTKy8tRVFSEjIwMTJw4Ebt27YrrOCRJQqbVkvB/0Xbl3bZtm5J9mThxIqxWK95//31UV1dj5EjOvRIR6YXXp/UI9C8pgpktW7bgb3/7G0aMGBFy/Mknn8TChQuxePFibNmyBQUFBZg8eTIcDodGI00Oe/fuxYkTJ5TMjMViweWXX47XX38dVVVVnGIiItIRTjPFTvNgprGxETNmzMALL7yAnj17KseFEFi0aBHmzZuHadOmobS0FMuWLUNzczOWL1+u4Yi1V1VVBavVitLSUuXYVVddhdWrV+PLL7/ksmwiIh1hAXDsNK+ZufPOO3HppZdi0qRJePzxx5Xj+/btQ21tLaZMmaIcs9lsmDBhAjZu3Ijf/e53YW/P6XTC6XQqXzc0pN76/erqapSWlsJqtSrHJk+eDK/XC5fLxWCGiEhHvBrEMh/s/hEPv/klWj0dW5Gckm3D328bg4Lc9MQPrJs0zcysWLEC1dXVqKio6PC92tpaAEB+fn7I8fz8fOV74VRUVCA3N1f5V1xcrO6gk0BFRQWqqqpCjtlsNjQ0NEAIgSFDhmg0MiIiipYW00zv7DiE2oZWnGh2d/j3TV0jNu09mvAxxUKzzExNTQ3uvfderF27FunpnUd/7QtjhRBdFss+9NBDmDNnjvJ1Q0NDSgY0RESUGrSYZmp0egAAcyYPxdThBcrxR9/ejU++OYrGVk/CxxQLzYKZqqoq1NXVhRSrer1ebNiwAYsXL8bXX38NwJ+hKSwsVM6pq6vrkK0JZrPZYLPZ4jdwIiIiFWmxNNvRFqyU9MnCqXl25Xh+jj+54HDqK5jRbJrpwgsvxM6dO7F9+3bl3+jRozFjxgxs374dgwYNQkFBAdatW6f8jMvlQmVlJcaNG6fVsImIiFSlRZsZOTOTnR6a07C3fc3MTITsdnvIahwAyMrKQu/evZXjZWVlmD9/PoYMGYIhQ4Zg/vz5yMzMxA033KDqWIy2LM5ovy8RUTLT4jVZDmZy2gczNkvI9/VC89VMXZk7dy5aWlowe/ZsHD9+HGPGjMHatWtht9tP/sMRSEtLAwA0NzcjIyNDldvUg+bmZgCB35+IiBLLF5SO0WKaSc68ZNtC3weymZmJ3fr160O+liQJ5eXlKC8vj8v9mc1m9OjRA3V1dQCAzMzMqDvx6okQAs3Nzairq0OPHj1gNke3hQIREanD7Qu0/dVimsnRyTSTHNw0MJjRl4ICfxW3HNAYQY8ePZTfm4iIEud/ttRg9ReH0OQKBAudTTMJIfDYO7vxzY+NsJgl3D5hMH4+qHfMY3B6vHB5/MFUtq1dMCNnZpzumO8nkQwfzEiShMLCQuTl5cHt1teD1x1paWnMyBARaWTBmq9wrMkVcszbSTDzTV0jlv5zv/K1BKgSzDQ5A43yOgQzNnOHc/TA8MGMzGw2802eiIjiRgiB+hb/h+a+PTLww4mWtuPhz29xhQYU8s/GSq6HybSaYTaFllZY294H5cyNXmi+NxMREZERtLp9SrHvsKIc5XhnBcCedsfVWmHkaJtCap+VAQCrxR8WuHS2lTeDGSIiogRwtPqDCJME9MwMrCJqH7Qox9sFFGqtMJIb5rUv/gWCghlmZoiIiKg9ZQWRzYJMayCQaB+0KMfbghyr2RTy87GSgyJ7uMxM2305GcwQERFRe0oQkZ4GW1rg7bfTzEzb8R5tWZxGp0eVBnvydJU9vWOvsUBmRl8FwAxmiIiIEqAxKDOTbgksOPF4u55mkoMZIYBmV+xBRnCGqD2bTmtmuJqJiIgIwP9srcF3RxpVvc3BfbJx7TnFAEJrVTKsQcGML3zg4G4LcuzpaTCbJHh9An9a81XIz0Yi3WLGDWP6K5tINqZgzQyDGSIiMrxv6xox97Udcbntswf0wKl59pDMzLDCwGqmzqeZ/AGFxSThlGwbahta8V+bDnRrDA2tbjxy+ZkAgKa2cWSFCYrkmhmf8GeGLGZ9TOAwmCEiIsM74nAC8K8yuursfqrc5qptP+BYkwt1Dqc/mGlbzZSdbsG/DemD2ycMxvOV33U6zSQv2baYJfxl+s/w4b9+jHoMuw83YON3Pym/HxAIntLCBCrBtTwuBjNERET6IS+bHtA7Cw9fNkyV26w6eBzHmlzK9JIjaBWRJEm48ef98Xzld3B3Up8iTzNZTCaMHdwbYwdH3/331a012PjdTyE9arxtGZ/2DfOAQGYG8E81ZVqjvktN6CPkIiIiiqPACh/1PuPLBbZyjUpju8Jbi8n/Ftxp07y2ICfN3P0NkO1hdsGWY6dwwYzFbIJ8WE91MwxmiIjI8NoHGmpQAom2226/U7WlLUjx+ETYJdduXyAz013yLtiO1sgyM0CgCFhPvWYYzBARkeEpK41UDGaUzExbEBPcZwYA0kxd95rxtqVQzDFkZrLbBVTB99VpMGPW3/JsBjNERGR4XTWS6672WRHlPtqCnOAgJdxUk1Ko20nQEdkYLG1jCGxS6WvLApmlzjIz+ttsksEMEREZXle9V7pLnmbadageb39xCAePNYfchyUoSAlXBKwUAMewoih4qkueypJXT3WW8bHpsNcMVzMREZHhBWpmomtI15XcDH9m5pNvjuKTb44qx3PkaaagICXc8mxvUJ+Z7pKDGZ8AWtxeZFot8J40M6O/aSYGM0REZHhuZeWQehMWl44oxJb9x3C82aUc69sjE+eU9AQABMco4WpmApmZ7gczGWlmmCR/MNPY6vEHM5HWzDAzQ0REpB/yYiJTJ9mK7sjPSceSG0d1+n1JkpBmluD2irBbGgQ6AHc/wJIkCdk2CxpaPXA4PcgDTh7M6HCaiTUzRERkeHJRrCmGKZ3ukAOVcNNM8rFY+swAgaJmuS5I6SzMpdlERESpQ36DT3AsowQU4aaZAkuoY3urbr9EXPlduTSbiIgodfjiMM0UCaVxXpjAQY0OwEBg9ZQjyswMp5mIiIh0RJxkhU+8yMuu3WGmmdToAAyEyczIU2onWc3k9Hhjut9EYjBDRESGJ7/BJziWUbIj4ZrmeVVYzQQEdQFua5wXvBt3OMzMEBER6ZDW00zuMKuZ3Cr0mQECHYc71Mx08rvadLg0m8EMEREZnu8ky5XjJSPN36SvxdVxSkc+lmGNrZGf3DhP3ujSc5LpK2ZmiIiIdMin0TSTvGw6eFdrWWC/qNhawrXfIyoQuIU/X48dgBnMEBGR4SmbLyY4M9O+ODdYYCfv2Da/DNTMhGZmOlvyrccOwAxmiIjI8OSSlUTXzLQvzg0W2C8qtsxMZzUzJ8vMsGkeERGRjigdgBM9zdRFZkbOpMQ8zdQuM+M9SWbGZvHX6HCaiYiISEd8J+m9Ei9y1qWrmplYMzPKfbTPzJxs12xmZoiIiPTDq9HS7Ox2K41kPp8IBDNqZWacbX1mTlIfpMdghrtmExGR4Qllo8nE3q+cNVm5pQZvbfshMJ4w53SXPJVVc6wFT639OiV3zWYwQ0REhqfVNNPPinvAbJLg9Qk0hek1c2ZRDmyW2CKsfj0zlf9/vep7ZQuFzoKZNGXzSwYzREREuuHVaDXT6IG9sHXepLA1MwBQ2CMdUoxjyrCa8b/3jMelz3yKRqdH6W3TWTAj76YdZoeFpMVghoiIDE9olJkBgJ5ZVvTMssb1Pk6x2wD4i4rljsKdbZMgXwM5W6UHLAAmIiLD82lUM5Mo9rbGez4RWKLdWeAmxzjhNr9MVin6sBEREUXuZJsv6l16mkmZVpJrczrbNVs+T0eJGQYzREREQqOl2YkiSVKHVVGd/a4Sp5mIiIj0J9B7ReOBxFH7YKbzmhn/fznNREREpCOBXbNTMzMDdNwWobPVTHJnYB0lZhjMEBERabXRZCJltq1ikvXIDL8bN6eZiIiIdEh+4+5sv6JUYA1qvndWv1yl30x7yjQTgxkiIiL9CEwzaTyQOLJaApmZMwpzOj3PrMOmeQxmiIjI8HwpvpoJAKxB1c1dbZFgUmpm9BPNMJghIiLD851k88VUYEsLvOVbuwhmJK5mIiIi0p/ARpMaDySObObIghlOMxEREemQ/MadykuzgwMYq9nc6XnK3kw6ima40SQRERmeEaaZQoKZCKaZ2i/NXv91HTbsORr2Z8YP6Y1fnJ4f+yC7icEMEREZnhGmmayRTjOF6TMjhMDdy7fB4fSE/Zn0NBODGSIiIi0ZYjVThJkZU5iaGafHpwQyvz1/UIetEMYM6q3iSKPHYIaIiAxPbhBnSuHUTHAAY+tiEypTmGkmR6s/kJEk4MGLT0+668QCYCIiMjxhhGmmSDMzYaaZGtuyMtlWS9IFMgCDGSIiIqWnSkpPM0VYMxNYzRQ41tiWmclOT84JHQYzRERkeEaombGFLM2OLjPjcLoBdNx5O1kwmCEiIkMLbtufhDMoqom8ANj/35BpJjkzY2MwQ0RElHSC2/ancmYm2poZb9A0k0OZZgq/07bWGMwQEZGhBS9BTsbiVrVkWgNZlSxr5xmWcBtNygXA9iTNzCTnqIiIiBLEZ5BppvGn9sE1o/oh02rGmUU5nZ5nDjfN5EzuaabkHBUREVGChAYzqRvNZNks+I9rzjrpeZIyzdSxzwxXMxERESWh4GmmVN6bKVKBaabAsca21UzJmplhMENERIYWnJlJ4cRMxMLtzSSvZuLSbCIioiTkM8hqpkjJl8Cro5oZBjNERGRoIdNMDGbCbjTZwJoZIiKi5MVpplDmcEuz2TSPiIgoecnTTJIUWMljZHINdPBqJqXPDDMzREREycfV1uo2zcS3RCAQ0AVPM7W4vQCA9DSzFkM6KT5yRERkaC6PP5jpqsW/kQQvT5enmuRrZEvSa5ScoyIiIkoQOTPDYMYvuNWOPNWkBHxmZmaIiIiSTuCNmm+JQGjdkDzVlOwBn6ajWrJkCUaMGIGcnBzk5ORg7NixeO+995TvCyFQXl6OoqIiZGRkYOLEidi1a5eGIyYiolTDaaZQwdNMPiHg9QklQ5Os10jTUfXr1w8LFizA1q1bsXXrVvziF7/AlVdeqQQsTz75JBYuXIjFixdjy5YtKCgowOTJk+FwOLQcNhERpZBkrwdJtOBpJp8QyvUBGMyEdfnll2Pq1KkYOnQohg4diieeeALZ2dn47LPPIITAokWLMG/ePEybNg2lpaVYtmwZmpubsXz5ci2HTUREKcSZ5FMoiWZqN80UEswk6VRc0ozK6/VixYoVaGpqwtixY7Fv3z7U1tZiypQpyjk2mw0TJkzAxo0bO70dp9OJhoaGkH9ERESd4TRTqOBgxusTcHq9ytdp5uTsw6P5I7dz505kZ2fDZrPh9ttvxxtvvIFhw4ahtrYWAJCfnx9yfn5+vvK9cCoqKpCbm6v8Ky4ujuv4iYhI31gAHCp4mkkETTNZLaakbSqo+SN32mmnYfv27fjss89wxx13YNasWdi9e7fy/fYXTgjR5cV86KGHUF9fr/yrqamJ29iJiEj/mJkJ1dk0ky2Jgz3N+xJbrVaceuqpAIDRo0djy5YtePrpp/HAAw8AAGpra1FYWKicX1dX1yFbE8xms8Fms8V30ERElDLkZccsAPYzmUKnmZJ9WTaQBJmZ9oQQcDqdKCkpQUFBAdatW6d8z+VyobKyEuPGjdNwhERElEqYmelIjmfaTzMlK00zM3/4wx9wySWXoLi4GA6HAytWrMD69euxZs0aSJKEsrIyzJ8/H0OGDMGQIUMwf/58ZGZm4oYbbtBy2ERElEJYM9OR2STB5xUh00wMZjrx448/4qabbsLhw4eRm5uLESNGYM2aNZg8eTIAYO7cuWhpacHs2bNx/PhxjBkzBmvXroXdbtdy2ERElEL0MI2SaP7aVAFvcGYmiYM9TYOZF198scvvS5KE8vJylJeXJ2ZARERkOE4dZB4STZ5m8vmELvrwaF4ATEREpKVk30RRC+a2FU2/+3sVWt3+PjPJHMwk78iIiIgSQA81IYlW1CMDALD7cAP2Hm0CABT3zNRySF1iZoaIiAzN5U3+zEOivfLbn2PbwRPK1xaThHNLemk3oJNgMENERIbGjSY76pNtw+Rhnfd0SzZ85IiIyND0sFqHusZHjoiIDI1Ls/WPjxwRERkaC4D1j48cEREZmpPTTLrHR46IiAyNmRn94yNHRESGxpoZ/eMjR0REhsbMjP7xkSMiIkNT+sywZka3+MgREZGhydNMtjS+JeoVHzkiIjI0bjSpfwxmiIjI0Fgzo3985IiIyNAYzOgfN5okIqKUJ4TA+j1HUNfQ2uF7rR7umq13DGaIiCjlfb7vGG5ZuqXLczLTWDOjVwxmiIgo5X1/vAUA0Cfbip8V9+jw/VEDeqFnljXBoyK1MJghIqKU1+j0AADGlPTGX2ecrfFoSG2cICQiopTnaHUDALJt/AyfihjMEBFRynO0ZWay0xnMpCIGM0RElPIaW9uCGWZmUhKDGSIiSnlyzYydmZmUxGCGiIhSnpyZYTCTmhjMEBFRypNrZrI4zZSSGMwQEVHKa3X7u/xmsDFeSmIwQ0REKY/7L6U2PqpERJTylGDGzLe9VMRHlYiIUp6TmZmUFvOj6vV6sX37dhw/flyN8RAREanO5WUwk8qiflTLysrw4osvAvAHMhMmTMDZZ5+N4uJirF+/Xu3xERERxUyeZrIxmElJUT+qr732Gs466ywAwNtvv419+/bhq6++QllZGebNm6f6AImIiGIVqJnhaqZUFPWC+6NHj6KgoAAA8O677+Kaa67B0KFDcdttt+GZZ55RfYBElNzue/ULvFb1fcy3U5CTjjfvPA8FuekqjIooFKeZUlvUj2p+fj52794Nr9eLNWvWYNKkSQCA5uZmmBnxEhnOezsPq3I7tQ2tqD7I2jtSn9cn4PUJAAxmUlXUmZlbbrkF1157LQoLCyFJEiZPngwA2Lx5M04//XTVB0hEyU3+xPu/94xHfk73siq/X7kdn3xzVGk5T6QmeYoJYDCTqqIOZsrLy1FaWoqamhpcc801sNlsAACz2YwHH3xQ9QESUfLy+QTcXv8n3oKcdPTOtnXrdnpmWgEEWs4TqSkkmGGfmZTUrU0qrr766g7HZs2aFfNgiEhf5KwMENsn3uy2zf+YmaF4cHq9yv+nmSUNR0Lx0q1Xn7vuugvHjh1TeyxEpDNqBTP2ts3/Gp3umMdE1F7wVgaSxGAmFUX86vP994HVCsuXL0djYyMAYPjw4aipqVF/ZESU9NRK32crwQwzM6Q+pccMp5hSVsTTTKeffjp69+6N8847D62traipqUH//v2xf/9+uN38NEVkRMH73cTyiVeeZmrgNBPFgZxBtKUxmElVEQcz9fX1qKqqwieffIJVq1Zh6tSpyM/Ph9PpxPvvv49f/epXSv8ZItK/w/UtmPvaDhxvdnV6zsGfmgHEvkLEnp4GgDUzFB/cZDL1RRzMuN1unHvuuTj33HPx+OOPo6qqCocPH8akSZPw0ksv4b777kO/fv3w9ddfx3O8RJQga3f9iE++ORrRubEGM1lWf4+qZheDGVJfq1vOzLAXWqqKOJjJycnByJEjcd5558HlcqG5uRnnnXceLBYLVq5ciX79+uHzzz+P51iJKIEaWvzTxxNPOwU3jxvY4ft7jzThsXd2A4j9E68cDAXX4BCpRS4st6d3awEv6UDEj+yhQ4ewadMmbNy4ER6PB6NHj8Y555wDl8uF6upqFBcXY/z48fEcKxElkFyMOyQvGxNPy+vw/b49HMr/x5qZkX/eyWCG4sDRNn0pF5pT6on4FahPnz64/PLLUVFRgczMTGzZsgV33303JEnCfffdh5ycHEyYMCGeYyWiBJIb2GXb0sJ+PzvoU27MwUxbZid4qTeRWhqdDGZSXbdfgXJzc3HttdciLS0NH330Efbt24fZs2erOTYi0pDyabaT1HzwG4MQIqb74jQTxdPJnsukf916ZHfs2IG+ffsCAAYMGIC0tDQUFBRg+vTpqg6OiLTT2NpWZ9DJp9ksa+B4rBkVBjMUT/Iquc6ey6R/3Xpki4uLlf//8ssvVRsMESUPJTXfyadZkynQVybWIMRm4TQTxc/Jnsukf1x0T0RhRVM0GWswYzWbVbkdonDk57Lcz4hSD4MZIgqr2eXfnC/L1nlvjv69MgEAZ/fvGdN9cZqJ4qmpLTMj9zOi1MOcGxGF5Wmb8rGYOv/M89LNo7H+6yO4dERhTPclBzMen4DPJ0KmsIhi5fT4A3ObhcFMqmIwQ0RhedtWKJm7CCxOzbPj1Dx7zPcVvLTb5fUh3cQ3HVKPXIsVawsBSl58ZIkoLK/v5MGMWoI7CLNxHqlN2ZuJwUzK6tYj+/e//x3nnXceioqKcODAAQDAokWL8NZbb6k6OCLSTiKDmTSzeiujiNrjRpOpL+pHdsmSJZgzZw6mTp2KEydOwOv1z0X26NEDixYtUnt8RKQRTwKDGUmSAkXAXJ5NKnMyM5Pyon5kn332WbzwwguYN28ezObAvPbo0aOxc+dOVQdHRNrxtQUzlgQV49rMXNFE8cGamdQX9SO7b98+jBw5ssNxm82GpqYmVQZFRNqTMzMmKTHBjPxGs6r6e+w6VJ+Q+yRjkANkG4OZlBX1I1tSUoLt27d3OP7ee+9h2LBhaoyJiJKAr201k8WcmGAmq60537MffYsbXtis1OwQxYoFwKkv6qXZ999/P+688060trZCCIHPP/8cr7zyCioqKvCf//mf8RgjEWlAqZlJUGbm4UvPwBvbfsB7X9aivsWNJpcHOezYSiqQp5mYmUldUQczt9xyCzweD+bOnYvm5mbccMMN6Nu3L55++mlcd9118RgjESWYzycgb4SdiAJgAJhyZgGmnFmAofPeg8vrQ2MrgxlSR2A1E/sXpapuNc37zW9+g9/85jc4evQofD4f8vLy1B4XEWlIbpgHJC6YkWWnW3CsyaVsDkgUK04zpb6YOgD36dNHrXEQURIJrldJeDBj8wcz8uaARLHw+YQyZcpgJnVFFMyMHDkSUoTz5tXV1TENiIi0p3UwA4CZGVJFcN8iBjOpK6Jg5pe//GWch0FEyUTraSYAaGRmhlQQvD0GOwCnroiCmUceeSTe4yCiJOL1BgUzCVrNJMtpC2Ycre6E3i+lpuAmjGkJajNAiccwlYg60DQzw2kmUlFw999IyyVIf6IuADaZun5CyHs1EZF+eZXuv0j4G0C2kplhMEOxc7r970k2TjGltKiDmTfeeCPka7fbjW3btmHZsmV49NFHVRsYEWnHq+zLlPg3gGybv7cMMzOkBu7LZAxRBzNXXnllh2NXX301zjzzTKxcuRK33XabKgMjIu0omRkNXv/tLAAmFbHHjDGo9uiOGTMGH3zwgVo3R0Qa8miamWHNDKmHwYwxqPLotrS04Nlnn0W/fv2i+rmKigqcc845sNvtyMvLwy9/+Ut8/fXXIecIIVBeXo6ioiJkZGRg4sSJ2LVrlxrDJqJOBNfMJJoczDRwNROpILCVAYOZVBb1o9uzZ0/06tVL+dezZ0/Y7Xa89NJL+I//+I+obquyshJ33nknPvvsM6xbtw4ejwdTpkxBU1OTcs6TTz6JhQsXYvHixdiyZQsKCgowefJkOByOaIdORBFSamY0eANQ+swwM0MqcLJmxhCirpn5y1/+ErK6wWQy4ZRTTsGYMWPQs2fPqG5rzZo1IV8vXboUeXl5qKqqwvnnnw8hBBYtWoR58+Zh2rRpAIBly5YhPz8fy5cvx+9+97toh09EEQhkZhKfmrHLmZkWd0iPELWkmSUu0TUQTjMZQ9TBzM033xyHYfjV19cDAHr16gUA2LdvH2prazFlyhTlHJvNhgkTJmDjxo1hgxmn0wmn06l83dDQELfxEqWqwGqmxL/py5mZ7440YejD76l++2MH9cby34xhQGMQnGYyhqgf3aVLl+LVV1/tcPzVV1/FsmXLuj0QIQTmzJmD8ePHo7S0FABQW1sLAMjPzw85Nz8/X/leexUVFcjNzVX+FRcXd3tMREYlN81LdMM8ABh8Sjb69siI2+1v2vsTp7AMhJkZY4g6M7NgwQI8//zzHY7n5eXht7/9LWbNmtWtgdx1113YsWMHPv300w7fa/8JSgjR6aeqhx56CHPmzFG+bmhoYEBDFCWvz/8GoEUwk2WzoPL+iWhyqd+A85zHP4DL60Oj0wN7eprqt0/JR+4zY2Mwk9KiDmYOHDiAkpKSDscHDBiAgwcPdmsQd999N1avXo0NGzaErIgqKCgA4M/QFBYWKsfr6uo6ZGtkNpsNNputW+MgIj95o2EtppkAf+Fxbob6bz7Z6RYca3L5e9jkqn7zlISYmTGGqB/dvLw87Nixo8PxL774Ar17947qtoQQuOuuu7Bq1Sp89NFHHYKkkpISFBQUYN26dcoxl8uFyspKjBs3LtqhE1GEPG2ZGZNGwUy8BJZ9c5rJKFgzYwxRZ2auu+463HPPPbDb7Tj//PMB+JdY33vvvbjuuuuiuq0777wTy5cvx1tvvQW73a7UweTm5iIjIwOSJKGsrAzz58/HkCFDMGTIEMyfPx+ZmZm44YYboh06EUXIp3FmJl7YkM94AtNMZo1HQvEUdTDz+OOP48CBA7jwwgthsfh/3OfzYebMmZg/f35Ut7VkyRIAwMSJE0OOL126VFk1NXfuXLS0tGD27Nk4fvw4xowZg7Vr18Jut0c7dCKKkJKZSbEVP9ncKsFwnJxmMoSogxmr1YqVK1fij3/8I7744gtkZGRg+PDhGDBgQNR3LtpWTHRFkiSUl5ejvLw86tsnou7xCblpXmoFM3YlM8PuwkbBmhljiDqYkQ0cOBBCCAwePFjJ0BBRavB4tWuaF09yZubPa/fgxU/3RfQzA3tnYfENZ/PNUKcYzBhD1FFIc3Mz7r77bqWnzJ49ezBo0CDcc889KCoqwoMPPqj6IIkosZTMTIrVzJx6SjYA4IjDiSMO50nO9tvzYyN2fH8Cowf2iufQKE5aPf4l/lyandqiDmYeeughfPHFF1i/fj0uvvhi5fikSZPwyCOPMJghSgHyrtmptppp9gWnYtypfeB0R9bD5uG3vsTeI01wsMZGt+T6KPYVSm1RBzNvvvkmVq5ciZ///OchjeuGDRuG7777TtXBEZE2tNzOIJ7MJgmjBkS+h1y+Pd0fzHD1k27JK9fkeilKTVHn3Y4cOYK8vLwOx5uamrjXCVGKkIMZLToAJxOuftI/+bGTH0tKTVEHM+eccw7+93//V/laDmBeeOEFjB07Vr2REZFmGMz4cfWT/slZtWxmZlJa1I9uRUUFLr74YuzevRsejwdPP/00du3ahU2bNqGysjIeYySiBFOCGYNnW5mZ0T9Hqz8QZWYmtUWdmRk3bhz++c9/orm5GYMHD8batWuRn5+PTZs2YdSoUfEYIxElmJa7ZicT+dM8a2b0izUzxtCtR3f48OHK0mwiSj2cZvJjZkbfhBBczWQQEQUzDQ0NEd9gTk5OtwdDRMlBbppn9GBGfgPk0mx9cnp8SpsBTjOltoge3R49epx0pZIQApIkweuNrH8DESUvH6eZAAQXADOY0SM5CJUkIDONG02msoiCmY8//jje4yCiJOLhNBMA1szonRyEZlstKdcAkkJFFMxMmDAh3uMgoiTC1Ux+gZoZLs3WI/aYMY5uPcInTpzAiy++iH/961+QJAnDhg3DrbfeitzcXLXHR0QaUDoAp9iu2dHK5jSTrjna+gOxx0zqi3pp9tatWzF48GD85S9/wbFjx3D06FEsXLgQgwcPRnV1dTzGSEQJJgczqbZrdrTsXM2ka8zMGEfUj/Dvf/97XHHFFXjhhRdgsfh/3OPx4Ne//jXKysqwYcMG1QdJRImVqnszRUv+RN/k8sLrE4avIdKbRnb/NYyoH+GtW7eGBDIAYLFYMHfuXIwePVrVwRGRNgJN86JO3qaU4E/0TS4PctirRBWf7f0Juw51bPnRO8uKS0cUIs2szvNOaZjHzEzKi/oRzsnJwcGDB3H66aeHHK+pqYHdbldtYESknUDTPI0HojGbxYw0swS31998jcFM7Opb3Ljpxc1wt/Uyas9qMWHq8EJV7ktemm238XFLdVEHM9OnT8dtt92GP//5zxg3bhwkScKnn36K+++/H9dff308xkhECRYIZgwezQCwmk1we71we31aDyUl1De74fYKWExSSNBSffA4vj/egkMnWlS7LwdrZgwj6kf4z3/+MyRJwsyZM+Hx+J8oaWlpuOOOO7BgwQLVB0hEicfMTIDFbALg7TSTQNFxevyNVbPTLXjm+pHK8T+8sRPLNx9UdeVYI1czGUbUj7DVasXTTz+NiooKfPfddxBC4NRTT0VmZmY8xkdEGmBmJiCtbXm6x8fMjBqcHv91tLaLlOOxciywLxODmVTX7Uc4MzMTw4cPV3MsRJQkPGyap5BXMHmYmVGFq226zmppF8zEoacPVzMZR8SP8K233hrReS+99FK3B0NEycHHpnkKS1t2Sg7wKDYuT/hgJh5bR7BmxjgifoRffvllDBgwACNHjoQQ/KMmSmUeNs1TKNNMLABWhauTaabstpViqk4zMTNjGBE/wrfffjtWrFiBvXv34tZbb8WNN96IXr16xXNsRKQReddsozfNA4KmmZiZUYUczNja7WIdj60j5F42rJlJfRFX9z333HM4fPgwHnjgAbz99tsoLi7Gtddei/fff5+ZGqIUo2RmGMwoDdxYM6MOuWbG1kkBsEOlTT03fnc06LbZZybVRbVUwWaz4frrr8e6deuwe/dunHnmmZg9ezYGDBiAxsbGeI2RiBLMx+0MFHLdkJurmVRxspoZtaaZvvkx8J40+JRsVW6Tkle3111KkgRJkiCEgI9/5EQpRV6GzMxMYHm6l5kZVXQazKSrWwAsT1dNH13MPbUMIKpgxul04pVXXsHkyZNx2mmnYefOnVi8eDEOHjyI7GxGvkSpQq51ZWYGSDOxz4yanN5O+swE1cyoUbrAlUzGEvGjPHv2bKxYsQL9+/fHLbfcghUrVqB3797xHBsRacTb9sbNPjNB00zMzKjiZJkZIYBmlxdZMa5AkmtvuJLJGCJ+lJ9//nn0798fJSUlqKysRGVlZdjzVq1apdrgiEgb8vs20/NBBcDMzKiis2AmI80Ms0mC1yfQ6PTEHMxwx2xjifhRnjlzJiR+SiMyBCUzw2CGHYBV1lkwI0kSsm0W1Le44Wj1ID8ntvuRC4mZmTGGqJrmEZExyG/cDGbYAVht8kaT7WtmACjBjBq9ZuRCYtbMGAN3kSOiDuSmeQxm2AFYbUrTPEvHtx81e80wM2MsDGaIqANlo0kGM+wArDJPF/t+pbd1BW5xeWO+H9bMGAuDGSLqQJ5mSuNGk+wArDIl6xemBlOuo3GpkAULBDPs/msEDGaIqAN325uJXC9iZHKvHXYAVoe3i60y5KkneSoqFpxmMha+UhFRB94upgKMRr4G7ACsjq4yM2oFM06PV8nusADYGBjMEFEHSl0DMzPKNXCzZkYVXWVm1JpmcgTt75RlZTBjBHylIqIOlGkmZmaUa8DVTOqQL6MpXM2MWZ3MjDzFlGU1s4jdIBjMEFEH8qfnNGZmgjoAMzOjBqEs++/4PTkz44w1mGGPGcPhKxURdeBm0zwFOwCry9sWzITNzKhUM+Ng8a/hMJghog7kfYi4NJu7ZqvN20UPI6vZ32cm1pqZQGaGy7KNgsEMEXUgZyEs4eYCDEa+Btw1Wx1ddZdWppncsQYz/g7CdmZmDIOvVETUgZyFsHCaSXnT9TIzowplNVOXTfNi6wAsFwCz+69xMJghog4CmRkGM4G9mZiZUUNXq5nU6jOjbDLJzIxhMJghohBCCPaZCcI+M+rqcjWTykuzuZrJOPhKRUQhvEFv2iwADuoAzGkmVUS0mkmlpnmsmTEOBjNEFCK4nwqXZgdlZjjNpIouVzOpNM3EPjPGw0eaiEK4gz4Vp3E1EzsAq6zL1UzmQNO8PT868P8q98Lpib4YeOv+4wCAbBuXZhsFgxkiChE8zcTVTEEFwKyZUYX8/JLCTDOlp/n7zDjdPvxtw168Xv19TPdV1CM9pp8n/WAwQ0QhgqdTOM0EmNummbiaSR1y6VG4XbPlaaFGpwcnml0AgCt/VoSRxT2ivp9T7Ok4f8gp3R4n6QuDGSIKEdxjJtynZ6NhB2B1+bpYzSQvpW50epQi3snD8nHZiKKEjY/0iRPiRBSCPWZCsQOwurpazWQPysw0slcMRYHBDBGF8HDH7BAWpQMwgxk1+LpYzaRkZloDwQy7+FIk+GpFRCHkVTtmZmYABDJUbq5mUkVXmRm5Zsbl9eFYo79mhiuSKBIMZogohDydwu6/fvJ14GomdSjbGYTJzGRZA1kYB3vFUBT4akVEIeTpFHb/9Qt0AGYwowZlmilMZsZskpBlNYccY80MRYLBDBGFcMurmRjMAAjUzHCaSR3yaqbOEn+9sq3K/1vNJgYzFBE+S4gohIfTTCHkLsjsM6MOuWYmXGYGAB7/5XC888UhCADnDz2FvY4oIgxmiChEcJ8ZCtrOgNNMqpCnmcLVzADAhKGnYMJQNruj6PCjFxGFCPSZ4csDEAjq2DRPHV2tZiLqLr5aEVEI+U2bBcB+Fm5noCplOwNm/khFDGaIKIS8NJtvNn7sM6MubxermYi6i8EMEYXwsgNwCDkzw6XZ6jjZaiai7uDTiYhCyBkILs32Cy4AFoIBTax8rJmhOGAwQ0QhPJxmChGcoeKKpth5u9ibiai7GMwQUYhAB2C+PAChe1Rxqil28jVkZobUxFcrIgrhZp+ZEMGFqgxmYidfQmZmSE0MZogohDzNxMyMX/Cbrpc1MzHznaQDMFF38NWKiEJ4WNMQIiSYYa+ZmMnZLcYypCZNg5kNGzbg8ssvR1FRESRJwptvvhnyfSEEysvLUVRUhIyMDEycOBG7du3SZrBEBuHhaqYQwTEdMzOxUzIzDJZJRZoGM01NTTjrrLOwePHisN9/8sknsXDhQixevBhbtmxBQUEBJk+eDIfDkeCREhmHh31mQkiSpLzxsmYmdlzNRPGg6UaTl1xyCS655JKw3xNCYNGiRZg3bx6mTZsGAFi2bBny8/OxfPly/O53v0vkUIniSgiBZpcXWTbt936V+8yYmZlRmE0SvD7BYCZGQgilAJirmUhNSfvRa9++faitrcWUKVOUYzabDRMmTMDGjRs7/Tmn04mGhoaQf0TJ7t//5wuUlr+Pl/+5T+uhBHUA5puNTC5WZTATm+DLx8wMqSlpg5na2loAQH5+fsjx/Px85XvhVFRUIDc3V/lXXFwc13ESqWHVth8gBLBs0wGth6LszcRdswMsnGZSRfD+VqzJIjUl/auV1C4VKYTocCzYQw89hPr6euVfTU1NvIdIpBp5x2pNx8AC4A5MpsCWBtR9rqBgxspgmVSk/QR9JwoKCgD4MzSFhYXK8bq6ug7ZmmA2mw02my3u4yNKVfIbNpvmBTAzow6Xh8EMxUfSPptKSkpQUFCAdevWKcdcLhcqKysxbtw4DUdGlNo8SgfgpH15SDiuZlKHHMykmSUl20WkBk0zM42Njfj222+Vr/ft24ft27ejV69e6N+/P8rKyjB//nwMGTIEQ4YMwfz585GZmYkbbrhBw1ETqcuXZG+QgQ7AfLORMZhRhxzMMCtDatM0mNm6dSsuuOAC5es5c+YAAGbNmoWXX34Zc+fORUtLC2bPno3jx49jzJgxWLt2Lex2u1ZDJlJdcB1BMgh0AOYbjkwJZtg0Lybyc91q4XOL1KVpMDNx4kSILl4cJElCeXk5ysvLEzcoogQLDmaS4b1SLgBmZiYgkJlJrsBTb5TMDIMZUhmfUUQaCy6KTIZgxs0C4A4CwYzGA9E5J4MZihM+o4g0FhzMuJPg3dLplt9wzBqPJHnITfOSYem8nrFmhuKFzygijQUHM8lQP9PodAMA7OlJ27kh4eTMDGOZ2ARqZhgok7oYzBBpLDiACQ5stNLo9AAAshnMKMwmZmbUwJoZihc+o4g0FpKZSYZgptUfzNiTYNPLZCHXD/mSoahJx+Tnt43TTKQyPqOINOYMCmA8PqF53xlmZjpStjPwMpiJhcvrBcDMDKmPzygijbXPxmhdN+Noy8xkMzOjYGZGHUpmhsEMqYzPKCKNtQ9enBpONbk8PuX+7bY0zcaRbEwSN5pUA2tmKF740YtII98dacSzH36DA8eaQ46XrdjW7Rf7Ptk2PDT1DCWr0ur2ouLdf6G2oTWinw+eRsmyccWJTN5BnNsZxIZ9ZiheGMwQaeQfnx3Am9sPdTj+8ddHYrrdc0t64cqf9QUA/PPbo1i26UDUt1GQkw4LizQVcmaGwUxsnJxmojhhMEOkkRaXvxhy0hn5mDwsD8W9MrHvaFO3b2/llhrs+L4eJ5rdyrHjbf9/Wr4dM8cNiPi2xpT06vY4UpGFG02qQiku5xQmqYzBDJFG3G1TOqMG9MT0c/oDAMYN7tPt29tRU48d39fD0RoIZhrb/v/U/GzMGBN5MEOhuGu2OuRl/1wpR2pjro9II3IDNrU2dJTfIBxtn36BwCdh9oyJDXfNVgefjxQvDGaINOJReUNHuehX/vQLBAIbLrOODTMz6pCzhszMkNoYzBBpxNO2JNusUqGtvJdSY3Bmhml9VZhN/seIwUxs2MOI4oXBDJFG5GXQaSplZpRgpjXMNFM6Cy5jIc8EMpiJDbtLU7zwGUW6VnXgGDZ++5Pqt2sySbhsRCEG9M5S/bZlbnmaSaXMjLxCZE+dA89++A0A4Msf6gGwRiFWzMyogzUzFC98RpFuCSFw27KtIUuR1bRl/zG8fMu5cbltAPCqXADcK8sKAKg51oKn1u0J+z3qHouJHYDVwGlPihc+o0i3nB6fEshcM6qfahmOH060YMOeI0qPlniRl2abVZpmOrekF+6bMhQ/nAjt9pufY8P5Q09R5T6MSt5oUutNQPXOwWlPihMGM6RbwYWuC64aoVpQ8PFXddiw54iSOYkXuQDYYlInCDObJNz1iyGq3BaFYmYmdk6PV9mbiQXApDYWAJNuySnrLKtZtUAGCF6Gq9pNhiXXX6g1zUTxY+au2TFrcnqV/2cwQ2pjMEO6Fa+VEYHW9fGNZtSeZqL4MTMzEzP5w0emyh8+iAAGM6Rjcs8KteffTQlqkBboAMw/w2RnZs1MzBzOtoZ5zMpQHPBVlHSrMU7dbRO1qaDcZ0atDsAUP8zMxI4rmSie+KwiAP437q37j6HZ5T35yRGSJP8mil1lTr6tc6DmWEu3bn/r/mMAAs3i1GJK0D48HpX7zFD8mCVuZxCpVrcXW/cfh7td0dkX358AwB4zFB98VhEA4KVP9+GJd/+l+u2eP/QU/Net4Xu11BxrxuS/bECsMYPawYySmfHGOzMjr2ZiZibZcW+myM1/91/4r00HOv0+l2VTPDCYIQDA98ebAQB5dhsKctNjvr0mpwffHWlCzbHmTs/Z/1MThADS00wYmm/v1v1YzSbc9POB3RxleCYpMZmZQAdgBjPJjtNMkdt7pAkA0L9XJnpkhgYuFpOEW8cP1GBUlOoYzBCAwIv0DWP6o2zS0Jhv74uaE7jyr/9U+kqEI8+hlxbl4rU7xsV8n2qRg4t4fwoPLM3mNFOyYwFw5OTGeP/nsmGYPCxf49GQUfBVlAAE+mfItQGxslr8Ty1nF8FMoBtocsXUiaqPkGsKuEw1+TEzE7nGVv+qpWT7u6bUxmCGAARW1phVmvIIBDOdFxQHVjck1xx6ot64Artm888w2VnYNC9i8VplSNQVvooSgEAWQrXMTNvUSZfTTEn6opeoKQW5zwxrZpKfXEfFzMzJNbYmZ8aVUhuDGQIQKHZVa8rD1paZcXl9EJ18mm1M1mmmRGVmWACsG/JjxJqZrnl9Ak1t7R2S7UMKpTY+2whA4I1VrWBGnmYSwn/b4fYfkjv4JtuL3sn24fF4fTjW7IrpPoSAsiRdrY0mKX4CmZn4bHHh8vhwoiW251QyCNl/Kck+pFBq47ONAAQ+carV80QOZgD/C3W4FTvJPs0ULjPj9vpw0V82YO/RJtXuj5mZ5GeJ4+ajLS4vLnxqPQ7Vt6p/4xqxmk2wWcxaD4MMJLneRUgz8hu3Sa1gxhwazGTZOp7jaFv1kGyf4OS6ISH8QV7wNfmp0aUEMmpcqvNO7cOOqDpgiuPmowePNSuBTKosbLvsrEKth0AGw1dRAqB+ZsZiNsEkAT7hr5sJRykUTLI38+BpH68QMCFwTeSC5kyrGbsfuzjhYyNtKJmZOJTMNLZtwDigdyYq779A/TsgMgBO1hOA4JoZ9Z4S8lRTZyualGmmJMvMBF+C9r1mXF5/TUDwNBqlPnMcMzPJWjtGpCd8RSYAQU3zVHxGyFNNnTXOS9YX8ZDMTLtgRv5drOzaayjx3JspWWvHiPSEr8gEIKhpnqqZGX8B4MkyM8m2NDv4ErQvApZ/F2ZmjMUSz2CGfVmIYsZXZAIQ1GdGpaZ5QGivmfaEEEGfSJOrA3BwZqZ9XxEngxlDMsVxi4tkzVAS6QlfkQlAUAdgFZdTdFUz0+r2KfeZbJ9Igy9Bp5kZTjMZSjw3H3Ukae0YkZ7wr4cAxCmYaXvD/+3ft3Z485ffEyTJvzIomUiSBLNJgtcnOjTOk4MZGzMzhqJkZuKwN5OyR1mSZSiJ9ITBDAEIBDNqLc0GgDMK7fj6RwdONLs7PefMohxIKk5tqUUOZjpkZrycZjIieerRE4e12fLS7GTLUBLpCf96CEAgmFGraR4APHXtz3DHxFO7TM0POiVLtftTk1w71L5mhgXAxiSXUcVj12yuZiKKHf96CEB8MjNmk4TTCuyq3V4iWTrZ0oA1M8akZGbiWADMzAxR9/EVmQAEagFMSTjlowVTJ0txnV65Zia56nwovuTYNR67ZjMzQxQ7BjMEICgzw00PAXTeV4TTTMZkjmNmRikAZmaGqNv4ikwAgmpmmJkB0HlmhsGMMXVWQ6UGpXkkVzMRdRtfkQlAfGpm9IyZGQpm7qSGSg3MzBDFjn89BADwtG2gp2afGT2TM1SPvr0LORmBT8x7jzQCYAGw0ch/Fyea3bj15S1dnnvxmQW49pziiG7X5xNodLFmhihW/OshAIC84wCDGb/8HBt+ONGCrQeOh/1+YW56gkdEWuqdbYXZJMHl9eGjr+q6PHfL/mMRBzMtbi/k1d5ZNhaVE3UXgxkCAHiZmQnx3IxR+OSbIwg3qZBlteDCM/ISPibSTp9sG167fSy+qWvs9BxHqwd/fGc3ml3eiG83eEd5rpAj6j4GMwQgPtsZ6FlBbjquGR3Zp2syhpH9e2Jk/56dfv94kwt/fGc3vD4Br09E9Lck12CZTRL/9ohiwIl/AhAUzHA1E1G3BBeFh9tcNRw2YSRSB/+CCECgaR4/HRJ1T7eCGa+3w88SUfT4F0QA2DSPKFYWkwQ5sen0RlY34+RSfyJV8C+IAHCaiShWkiQp00WcZiJKLBYAq2z913Wo7mQ5bzKTe4Fxmomo+6wWE5weX9TBjI2ZGaKYMJhRUYvLi9/+VxVc3sheyJKN2SQhPY3LQ4m6y2YxwYHQJdddkV8rOM1EFBsGMyo60eKCy+uD2SThxjH9tR5O1EYN7IUsdiEl6rbuTjMxM0MUG75zqUjeYyUn3YJHryzVeDRElGhyhiXS7Cz3+iJSB/+CVORwcsM4IiNTghlOMxElFP+CVKTsfmtLO8mZRJSKog1mnFzNRKQK/gWpqLEtM2Nn3QmRIclBScQFwJxmIlIF/4JUpGRmOM1EZEjR1swEmuZxFSFRLBjMqEipmWFmhsiQ5KCETfOIEot/QSpytLoBMDNDZFS2aAuAOc1EpAr+BalInmayM5ghMqQsqz8z0+zyRHS+fF62jdNMRLFgMKMiFgATGZuclXW0RhbMBKamuQKSKBYMZlTEmhkiY5ODkkiDGS4aIFIHgxkVBV6Y+CmLyIjkKeZGpzui85nNJVKHLoKZ5557DiUlJUhPT8eoUaPwySefaD2ksBqZmSEyNPlvX34tOBlmZojUkfTBzMqVK1FWVoZ58+Zh27Zt+Ld/+zdccsklOHjwoNZD64AFwETGJgcz0dfM8DWDKBZJ/xe0cOFC3Hbbbfj1r38NAFi0aBHef/99LFmyBBUVFZqNq6HVjYaW0FRyfdvXfGEiMib5g8yxJhe+P9580vPrm10hP0dE3ZPUf0EulwtVVVV48MEHQ45PmTIFGzduDPszTqcTTqdT+bqhoSEuY/vHZwfw5Jqvw36PKWMiY5L/9ncdasD4P30c+c/xAxBRTJL6L+jo0aPwer3Iz88POZ6fn4/a2tqwP1NRUYFHH3007mOzmCSlQVaw4X1zMaBXZtzvn4iSz/C+uRian40DP508KyM7sygHA/tkxXFURKkvqYMZmSRJIV8LITockz300EOYM2eO8nVDQwOKi4tVH9Nvzx+M354/WPXbJSL9sqenYe3vJ2g9DCLDSepgpk+fPjCbzR2yMHV1dR2yNTKbzQabzZaI4REREVESSOrVTFarFaNGjcK6detCjq9btw7jxo3TaFRERESUTJI6MwMAc+bMwU033YTRo0dj7Nix+Nvf/oaDBw/i9ttv13poRERElASSPpiZPn06fvrpJzz22GM4fPgwSktL8e6772LAgAFaD42IiIiSgCSEEFoPIp4aGhqQm5uL+vp65OTkaD0cIiIiikA0799JXTNDREREdDIMZoiIiEjXGMwQERGRrjGYISIiIl1jMENERES6xmCGiIiIdI3BDBEREekagxkiIiLSNQYzREREpGtJv51BrOQGxw0NDRqPhIiIiCIlv29HslFBygczDocDAFBcXKzxSIiIiChaDocDubm5XZ6T8nsz+Xw+HDp0CHa7HZIkqXa7DQ0NKC4uRk1NDfd8iiNe58ThtU4MXufE4HVOjHheZyEEHA4HioqKYDJ1XRWT8pkZk8mEfv36xe32c3Jy+IeSALzOicNrnRi8zonB65wY8brOJ8vIyFgATERERLrGYIaIiIh0jcFMN9lsNjzyyCOw2WxaDyWl8TonDq91YvA6Jwavc2Iky3VO+QJgIiIiSm3MzBAREZGuMZghIiIiXWMwQ0RERLrGYIaIiIh0jcFMNz333HMoKSlBeno6Ro0ahU8++UTrIelGRUUFzjnnHNjtduTl5eGXv/wlvv7665BzhBAoLy9HUVERMjIyMHHiROzatSvkHKfTibvvvht9+vRBVlYWrrjiCnz//feJ/FV0paKiApIkoaysTDnG66yeH374ATfeeCN69+6NzMxM/OxnP0NVVZXyfV7r2Hk8Hjz88MMoKSlBRkYGBg0ahMceeww+n085h9c5ehs2bMDll1+OoqIiSJKEN998M+T7al3T48eP46abbkJubi5yc3Nx00034cSJE+r8EoKitmLFCpGWliZeeOEFsXv3bnHvvfeKrKwsceDAAa2HpgsXXXSRWLp0qfjyyy/F9u3bxaWXXir69+8vGhsblXMWLFgg7Ha7eP3118XOnTvF9OnTRWFhoWhoaFDOuf3220Xfvn3FunXrRHV1tbjgggvEWWedJTwejxa/VlL7/PPPxcCBA8WIESPEvffeqxzndVbHsWPHxIABA8TNN98sNm/eLPbt2yc++OAD8e233yrn8FrH7vHHHxe9e/cW77zzjti3b5949dVXRXZ2tli0aJFyDq9z9N59910xb9488frrrwsA4o033gj5vlrX9OKLLxalpaVi48aNYuPGjaK0tFRcdtllqvwODGa64dxzzxW33357yLHTTz9dPPjggxqNSN/q6uoEAFFZWSmEEMLn84mCggKxYMEC5ZzW1laRm5srnn/+eSGEECdOnBBpaWlixYoVyjk//PCDMJlMYs2aNYn9BZKcw+EQQ4YMEevWrRMTJkxQghleZ/U88MADYvz48Z1+n9daHZdeeqm49dZbQ45NmzZN3HjjjUIIXmc1tA9m1Lqmu3fvFgDEZ599ppyzadMmAUB89dVXMY+b00xRcrlcqKqqwpQpU0KOT5kyBRs3btRoVPpWX18PAOjVqxcAYN++faitrQ25xjabDRMmTFCucVVVFdxud8g5RUVFKC0t5ePQzp133olLL70UkyZNCjnO66ye1atXY/To0bjmmmuQl5eHkSNH4oUXXlC+z2utjvHjx+PDDz/Enj17AABffPEFPv30U0ydOhUAr3M8qHVNN23ahNzcXIwZM0Y55+c//zlyc3NVue4pv9Gk2o4ePQqv14v8/PyQ4/n5+aitrdVoVPolhMCcOXMwfvx4lJaWAoByHcNd4wMHDijnWK1W9OzZs8M5fBwCVqxYgerqamzZsqXD93id1bN3714sWbIEc+bMwR/+8Ad8/vnnuOeee2Cz2TBz5kxea5U88MADqK+vx+mnnw6z2Qyv14snnngC119/PQA+p+NBrWtaW1uLvLy8Drefl5enynVnMNNNkiSFfC2E6HCMTu6uu+7Cjh078Omnn3b4XneuMR+HgJqaGtx7771Yu3Yt0tPTOz2P1zl2Pp8Po0ePxvz58wEAI0eOxK5du7BkyRLMnDlTOY/XOjYrV67EP/7xDyxfvhxnnnkmtm/fjrKyMhQVFWHWrFnKebzO6lPjmoY7X63rzmmmKPXp0wdms7lDJFlXV9chcqWu3X333Vi9ejU+/vhj9OvXTzleUFAAAF1e44KCArhcLhw/frzTc4yuqqoKdXV1GDVqFCwWCywWCyorK/HMM8/AYrEo14nXOXaFhYUYNmxYyLEzzjgDBw8eBMDntFruv/9+PPjgg7juuuswfPhw3HTTTfj973+PiooKALzO8aDWNS0oKMCPP/7Y4faPHDmiynVnMBMlq9WKUaNGYd26dSHH161bh3Hjxmk0Kn0RQuCuu+7CqlWr8NFHH6GkpCTk+yUlJSgoKAi5xi6XC5WVlco1HjVqFNLS0kLOOXz4ML788ks+Dm0uvPBC7Ny5E9u3b1f+jR49GjNmzMD27dsxaNAgXmeVnHfeeR3aC+zZswcDBgwAwOe0Wpqbm2Eyhb5tmc1mZWk2r7P61LqmY8eORX19PT7//HPlnM2bN6O+vl6d6x5zCbEByUuzX3zxRbF7925RVlYmsrKyxP79+7Uemi7ccccdIjc3V6xfv14cPnxY+dfc3Kycs2DBApGbmytWrVoldu7cKa6//vqwSwH79esnPvjgA1FdXS1+8YtfGHp5ZSSCVzMJweusls8//1xYLBbxxBNPiG+++Ub893//t8jMzBT/+Mc/lHN4rWM3a9Ys0bdvX2Vp9qpVq0SfPn3E3LlzlXN4naPncDjEtm3bxLZt2wQAsXDhQrFt2zal3Yha1/Tiiy8WI0aMEJs2bRKbNm0Sw4cP59Jsrf31r38VAwYMEFarVZx99tnKsmI6OQBh/y1dulQ5x+fziUceeUQUFBQIm80mzj//fLFz586Q22lpaRF33XWX6NWrl8jIyBCXXXaZOHjwYIJ/G31pH8zwOqvn7bffFqWlpcJms4nTTz9d/O1vfwv5Pq917BoaGsS9994r+vfvL9LT08WgQYPEvHnzhNPpVM7hdY7exx9/HPY1edasWUII9a7pTz/9JGbMmCHsdruw2+1ixowZ4vjx46r8DpIQQsSe3yEiIiLSBmtmiIiISNcYzBAREZGuMZghIiIiXWMwQ0RERLrGYIaIiIh0jcEMERER6RqDGSIiItI1BjNERESkawxmiEhT69evhyRJOHHihNZDISKdYjBDRAk1ceJElJWVKV+PGzcOhw8fRm5ubkLuf/369SgsLASbnxOlDovWAyAiY7NarSgoKEjY/a1evRpXXHEFJElK2H0SUXwxM0NECXPzzTejsrISTz/9NCRJgiRJePnll0OmmV5++WX06NED77zzDk477TRkZmbi6quvRlNTE5YtW4aBAweiZ8+euPvuu+H1epXbdrlcmDt3Lvr27YusrCyMGTMG69ev7zAGOZgBgNdeew3Dhw9HRkYGevfujUmTJqGpqSkRl4KIVMTMDBElzNNPP409e/agtLQUjz32GABg165dHc5rbm7GM888gxUrVsDhcGDatGmYNm0aevTogXfffRd79+7FVVddhfHjx2P69OkAgFtuuQX79+/HihUrUFRUhDfeeAMXX3wxdu7ciSFDhij3VVtbiwsvvBCHDx/G9ddfjyeffBK/+tWv4HA48Mknn3D6iUiHGMwQUcLk5ubCarUiMzNTmVr66quvOpzndruxZMkSDB48GABw9dVX4+9//zt+/PFHZGdnY9iwYbjgggvw8ccfY/r06fjuu+/wyiuv4Pvvv0dRUREA4L777sOaNWuwdOlSzJ8/HwDw1ltv4aKLLkJ6ejoOHz4Mj8eDadOmYcCAAQCA4cOHJ+IyEJHKGMwQUdLJzMxUAhkAyM/Px8CBA5GdnR1yrK6uDgBQXV0NIQSGDh0acjtOpxO9e/dWvn7rrbcwe/ZsAMBZZ52FCy+8EMOHD8dFF12EKVOm4Oqrr0bPnj3j+asRURwwmCGipJOWlhbytSRJYY/5fD4AgM/ng9lsRlVVFcxmc8h5cgBUW1uL6upqXHrppQAAs9mMdevWYePGjVi7di2effZZzJs3D5s3b0ZJSUm8fjUiigMWABNRQlmt1pDCXTWMHDkSXq8XdXV1OPXUU0P+ydNZq1evxtixY9GnTx/l5yRJwnnnnYdHH30U27Ztg9VqxRtvvKHq2Igo/piZIaKEGjhwIDZv3oz9+/cjOztbya7EYujQoZgxYwZmzpyJp556CiNHjsTRo0fx0UcfYfjw4Zg6dSpWr16NK6+8UvmZzZs348MPP8SUKVOQl5eHzZs348iRIzjjjDNiHg8RJRYzM0SUUPfddx/MZjOGDRuGU045BQcPHlTldpcuXYqZM2fi3//933HaaafhiiuuwObNm1FcXIympiZ8+OGHypJsAMjJycGGDRswdepUDB06FA8//DCeeuopXHLJJaqMh4gSRxJch0hEKW7VqlV4+OGHsXv3bq2HQkRxwMwMEaW87Oxs/OlPf9J6GEQUJ8zMEBERka4xM0NERES6xmCGiIiIdI3BDBEREekagxkiIiLSNQYzREREpGsMZoiIiEjXGMwQERGRrjGYISIiIl1jMENERES69v8BqZioHv7BBSkAAAAASUVORK5CYII=",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot first 1000 sample points\n",
    "time = [sp*i for i in 1:1000];\n",
    "PyPlot.plot(time,n_traj[1:1000], label = L\"$N$\");\n",
    "PyPlot.legend(); PyPlot.xlabel(\"time/s\"); PyPlot.ylabel(\"Molecule #s\")\n",
    "PyPlot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGwCAYAAABFFQqPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAmSElEQVR4nO3df3RU9Z3/8ddAkhkVEvllQmpIJsgWaOTXjLIBg9JqMFBXVuqm1iLVlu5UXQxZz+GXXRRXh7Wsy7KQpGDoSukXOD1gl57GhbiFSGWokGaQA1mWHlNCMSkNrQmVmpBwv394mNNxJkhCkpncz/NxzpzjfOZzJ++3qHn5uZ97r8OyLEsAAAAGGhDrAgAAAGKFIAQAAIxFEAIAAMYiCAEAAGMRhAAAgLEIQgAAwFgEIQAAYKyEWBcQzy5fvqwPPvhAgwcPlsPhiHU5AADgGliWpQsXLig9PV0DBlx9zYcgdBUffPCBMjIyYl0GAADohjNnzujWW2+96hyC0FUMHjxY0id/I5OTk2NcDQAAuBYtLS3KyMgI/R6/GoLQVVw5HZacnEwQAgCgn7mWbS1slgYAAMYiCAEAAGMRhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGIsgBAAAjEUQAgAAxiIIAQAAYxGEAACAsQhCAADAWAQhAABgLIIQAAAwFkEIAAAYKyHWBcBsM1+fGXV834J9fVwJAMBErAgBAABjEYQAAICxCEIAAMBYBCEAAGAsNktDEpuWAQBmYkUIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYXD6PHsVl+ACA/oQVIQAAYCyCEAAAMBZBCAAAGIsgBAAAjMVmaVwVm58BAHbGihAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGNx1ZhhOrsKDAAAE7EiBAAAjEUQAgAAxiIIAQAAY8VNECopKZHb7ZbL5ZLH49GBAweuOr+qqkoej0cul0vZ2dkqKyvrdO727dvlcDg0d+7cHq4aAAD0Z3GxWXrHjh0qKipSSUmJpk+fru9///sqKCjQiRMnNGrUqIj5dXV1mj17thYuXKitW7fqnXfe0ZNPPqkRI0Zo3rx5YXNPnz6tZ599Vnl5eX3VjhHYdA0AsIO4CEKvvvqqvvnNb+pb3/qWJGnt2rXas2ePSktL5ff7I+aXlZVp1KhRWrt2rSRp3LhxOnLkiNasWRMWhDo6OvToo4/qhRde0IEDB/Thhx9etY7W1la1traG3re0tFx/czFEWAEA4Opifmqsra1N1dXVys/PDxvPz8/XwYMHox4TCAQi5s+aNUtHjhzRpUuXQmOrVq3SiBEj9M1vfvOaavH7/UpJSQm9MjIyutgNAADoT2IehJqamtTR0aHU1NSw8dTUVDU2NkY9prGxMer89vZ2NTU1SZLeeecdlZeXa9OmTddcy7Jly9Tc3Bx6nTlzpovdAACA/iQuTo1JksPhCHtvWVbE2GfNvzJ+4cIFff3rX9emTZs0fPjwa67B6XTK6XR2oWoAANCfxTwIDR8+XAMHDoxY/Tl37lzEqs8VaWlpUecnJCRo2LBhOn78uH7zm9/ogQceCH1++fJlSVJCQoJOnjyp0aNH93AnAACgv4n5qbGkpCR5PB5VVlaGjVdWVmratGlRj8nNzY2Yv3fvXnm9XiUmJmrs2LE6duyYgsFg6PU3f/M3mjlzpoLBIHt/AACApDhYEZKk4uJizZ8/X16vV7m5udq4caPq6+vl8/kkfbJ35+zZs9qyZYskyefzaf369SouLtbChQsVCARUXl6ubdu2SZJcLpdycnLCfsbNN98sSRHjAADAXHERhAoLC3X+/HmtWrVKDQ0NysnJUUVFhTIzMyVJDQ0Nqq+vD813u92qqKjQ4sWLtWHDBqWnp2vdunUR9xACAAC4God1ZZcxIrS0tCglJUXNzc1KTk6OdTldFk/3Edq3YF/U8c5q7Gw+AACfpSu/v2O+RwgAACBWCEIAAMBYBCEAAGCsuNgsDfuLp/1KAABcwYoQAAAwFkEIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxuHw+hni8BAAAscWKEAAAMBZBCAAAGIsgBAAAjMUeIfQrXd1XxT4sAMDVsCIEAACMRRACAADGIggBAABjsUcIttDZXiAAAK6GIIS4RLABAPQFTo0BAABjEYQAAICxODVmA5xGAgCge1gRAgAAxiIIAQAAYxGEAACAsQhCAADAWAQhAABgLIIQAAAwFkEIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYBCEAAGAsnj4P/IWZr8+MOr5vwb4+rgQA0BdYEQIAAMYiCAEAAGNxagxG6uwUGADALKwIAQAAYxGEAACAsQhCAADAWAQhAABgLIIQAAAwFkEIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYBCEAAGAsghAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGMlxLoAoD+Y+frMqOP7Fuzr40oAAD2JFSEAAGAsghAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGMRhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGCtuglBJSYncbrdcLpc8Ho8OHDhw1flVVVXyeDxyuVzKzs5WWVlZ2Oe7du2S1+vVzTffrJtuukmTJk3SD3/4w95sAQAA9DNxEYR27NihoqIirVixQjU1NcrLy1NBQYHq6+ujzq+rq9Ps2bOVl5enmpoaLV++XIsWLdLOnTtDc4YOHaoVK1YoEAjovffe0+OPP67HH39ce/bs6au2AABAnHNYlmXFuoipU6dqypQpKi0tDY2NGzdOc+fOld/vj5i/ZMkS7d69W7W1taExn8+no0ePKhAIdPpzpkyZojlz5ujFF1+8prpaWlqUkpKi5uZmJScnd6Gja9NTTzTv7HvQ+3j6PADEn678/o75ilBbW5uqq6uVn58fNp6fn6+DBw9GPSYQCETMnzVrlo4cOaJLly5FzLcsS//zP/+jkydPasaMGZ3W0traqpaWlrAXAACwr5gHoaamJnV0dCg1NTVsPDU1VY2NjVGPaWxsjDq/vb1dTU1NobHm5mYNGjRISUlJmjNnjv7jP/5D9913X6e1+P1+paSkhF4ZGRnX0RkAAIh3MQ9CVzgcjrD3lmVFjH3W/E+PDx48WMFgUIcPH9ZLL72k4uJi7d+/v9PvXLZsmZqbm0OvM2fOdKMTAADQXyTEuoDhw4dr4MCBEas/586di1j1uSItLS3q/ISEBA0bNiw0NmDAAN12222SpEmTJqm2tlZ+v1/33HNP1O91Op1yOp3X0Q0AAOhPYr4ilJSUJI/Ho8rKyrDxyspKTZs2Leoxubm5EfP37t0rr9erxMTETn+WZVlqbW29/qIBAIAtxHxFSJKKi4s1f/58eb1e5ebmauPGjaqvr5fP55P0ySmrs2fPasuWLZI+uUJs/fr1Ki4u1sKFCxUIBFReXq5t27aFvtPv98vr9Wr06NFqa2tTRUWFtmzZEnZlGgAAMFtcBKHCwkKdP39eq1atUkNDg3JyclRRUaHMzExJUkNDQ9g9hdxutyoqKrR48WJt2LBB6enpWrdunebNmxea89FHH+nJJ5/Ub3/7W91www0aO3astm7dqsLCwj7vDwAAxKe4uI9QvOI+Qvgs3EcIAOJPv7qPEAAAQKwQhAAAgLEIQgAAwFhxsVka14a9QAAA9CxWhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBabpYEe1lM3ygQA9D5WhAAAgLEIQgAAwFgEIQAAYCz2CAHXgZtcAkD/xooQAAAwFkEIAAAYi1NjcYjTLQAA9A1WhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGIsgBAAAjEUQAgAAxiIIAQAAYxGEAACAsQhCAADAWAQhAABgLIIQAAAwFkEIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYBCEAAGAsghAAADAWQQgAABgr4XoOvnTpkhobG3Xx4kWNGDFCQ4cO7am6AAAAel2XV4T+9Kc/6fvf/77uuecepaSkKCsrS+PHj9eIESOUmZmphQsX6vDhw71RKwAAQI/qUhD6t3/7N2VlZWnTpk364he/qF27dikYDOrkyZMKBAJauXKl2tvbdd999+n+++/XqVOneqtuAACA69alU2MHDx7Uvn37dPvtt0f9/M4779QTTzyhsrIylZeXq6qqSmPGjOmRQgEAAHpal4LQj3/842ua53Q69eSTT3arIAAAgL7CVWMAAMBY3b5q7PDhw1q6dKl+//vf67bbbtOkSZNCr1GjRvVkjQAAAL2i2ytC8+fP18CBA+Xz+ZSdna2qqio9/vjjysrK0rBhw3qyRgAAgF7R7RWhM2fO6Gc/+5lGjx4dNn769GkFg8HrrQsAAKDXdTsI5ebm6re//W1EEMrMzFRmZuZ1FwYAANDbun1qrLi4WC+++KL+8Ic/9GQ9AAAAfabbK0Jf/vKX5XA4NGbMGD3wwAPKzc3V5MmTNXHiRDmdzp6sEbCFma/PjDq+b8G+Pq4EAHBFt4PQqVOndPTo0dBr9erVOn36tBISEjR27Fi99957PVknAABAj+t2EBo9erRGjx6thx56KDTW0tKiYDBICAIAAP1Cl/YI1dfXX/Xz5ORkzZgxQ08//bQk6ezZs92vDAAAoJd1KQjdcccdWrhwod59991O5zQ3N2vTpk3KycnRrl27rrtAAACA3tKlU2O1tbV6+eWXdf/99ysxMVFer1fp6elyuVz64x//qBMnTuj48ePyer363ve+p4KCgt6qGwAA4Lp1aUVo6NChWrNmjT744AOVlZXpr/7qr9TU1KRTp05Jkh599FFVV1frnXfeIQQBAIC4163N0r/73e+UkJCgRx55RHfeeWdP1wQAANAnuhyEtm3bpgULFqi9vV0Oh0OTJ0/Wm2++qREjRvRGfQAAAL2my3eWfuGFF/TYY4/p17/+tX7+859rwIABWrp0aW/UBgAA0Ku6vCL0/vvva8+ePcrMzFR2drZ++MMfasqUKSovL++N+gAAAHpNl1eE2tvbdcMNN4Tef/7zn9fly5fV2NjYo4UBAAD0tm5tln799dc1ffp0TZgwQYMGDVJCQoIuXrzY07UBRuAZZAAQO10OQnfddZf++Z//WRcuXNCAAQPkdrv18ccfq7y8XPfee688Ho+Sk5N7o1YAAIAe1eUg9Pbbb0v65KGr1dXV+tWvfqXq6mqVlpbK7/drwIABGjNmjGpra3u8WAAAgJ7U7YeujhkzRmPGjNFXv/rV0FhdXZ2OHDmimpqaHikOAACgN3U7CEXjdrvldrv18MMP9+TXAgAA9IouXzUGAABgFwQhAABgLIIQAAAwFkEIAAAYK26CUElJidxut1wulzwejw4cOHDV+VVVVfJ4PHK5XMrOzlZZWVnY55s2bVJeXp6GDBmiIUOG6N5779W7777bmy0AAIB+Ji6C0I4dO1RUVKQVK1aopqZGeXl5KigoUH19fdT5dXV1mj17tvLy8lRTU6Ply5dr0aJF2rlzZ2jO/v379cgjj2jfvn0KBAIaNWqU8vPzdfbs2b5qCwAAxDmHZVlWrIuYOnWqpkyZotLS0tDYuHHjNHfuXPn9/oj5S5Ys0e7du8Nu2ujz+XT06FEFAoGoP6Ojo0NDhgzR+vXr9dhjj11TXS0tLUpJSVFzc3Ov3C27s0crABKP2ACA7urK7++Yrwi1tbWpurpa+fn5YeP5+fk6ePBg1GMCgUDE/FmzZunIkSO6dOlS1GMuXryoS5cuaejQoZ3W0traqpaWlrAXAACwr5gHoaamJnV0dCg1NTVsPDU1tdMn2jc2Nkad397erqampqjHLF26VJ/73Od07733dlqL3+9XSkpK6JWRkdHFbgAAQH8S8yB0hcPhCHtvWVbE2GfNjzYuSa+88oq2bdumXbt2yeVydfqdy5YtU3Nzc+h15syZrrQAAAD6mR59xEZ3DB8+XAMHDoxY/Tl37lzEqs8VaWlpUecnJCRo2LBhYeNr1qzRyy+/rLfeeksTJky4ai1Op1NOp7MbXQAAgP4o5itCSUlJ8ng8qqysDBuvrKzUtGnToh6Tm5sbMX/v3r3yer1KTEwMjX3ve9/Tiy++qP/+7/+W1+vt+eIBAEC/FvMgJEnFxcV67bXXtHnzZtXW1mrx4sWqr6+Xz+eT9Mkpq7+80svn8+n06dMqLi5WbW2tNm/erPLycj377LOhOa+88oqee+45bd68WVlZWWpsbFRjY6P+9Kc/9Xl/AAAgPsX81JgkFRYW6vz581q1apUaGhqUk5OjiooKZWZmSpIaGhrC7inkdrtVUVGhxYsXa8OGDUpPT9e6des0b9680JySkhK1tbXpK1/5StjPWrlypZ5//vk+6QsAAMS3uLiPULziPkKIR9xfCACurl/dRwgAACBWCEIAAMBYBCEAAGAsghAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGMRhAAAgLEIQgAAwFhx8awxANeus0ez8OgNAOg6VoQAAICxCEIAAMBYBCEAAGAsghAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGMRhAAAgLEIQgAAwFg8YgOwCR69AQBdx4oQAAAwFkEIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYBCEAAGAsghAAADAWd5YGbI47TgNA51gRAgAAxiIIAQAAYxGEAACAsQhCAADAWAQhAABgLIIQAAAwFkEIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYBCEAAGAsghAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGMRhAAAgLESYl0AgPgy8/WZUcf3LdjXx5UAQO9jRQgAABiLFSHAUJ2t/ACASVgRAgAAxiIIAQAAYxGEAACAsQhCAADAWAQhAABgLIIQAAAwFkEIAAAYiyAEAACMRRACAADG4s7SAK4JzyADYEesCAEAAGMRhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGCtuglBJSYncbrdcLpc8Ho8OHDhw1flVVVXyeDxyuVzKzs5WWVlZ2OfHjx/XvHnzlJWVJYfDobVr1/Zi9QAAoD+KiyC0Y8cOFRUVacWKFaqpqVFeXp4KCgpUX18fdX5dXZ1mz56tvLw81dTUaPny5Vq0aJF27twZmnPx4kVlZ2dr9erVSktL66tWAABAP+KwLMuKdRFTp07VlClTVFpaGhobN26c5s6dK7/fHzF/yZIl2r17t2pra0NjPp9PR48eVSAQiJiflZWloqIiFRUVXbWO1tZWtba2ht63tLQoIyNDzc3NSk5O7kZnV9fZDeqA/oQbKgKINy0tLUpJSbmm398xXxFqa2tTdXW18vPzw8bz8/N18ODBqMcEAoGI+bNmzdKRI0d06dKlbtfi9/uVkpISemVkZHT7uwAAQPyLeRBqampSR0eHUlNTw8ZTU1PV2NgY9ZjGxsao89vb29XU1NTtWpYtW6bm5ubQ68yZM93+LgAAEP/i5lljDocj7L1lWRFjnzU/2nhXOJ1OOZ3Obh8PAAD6l5gHoeHDh2vgwIERqz/nzp2LWPW5Ii0tLer8hIQEDRs2rNdqBXDteEgrgP4g5kEoKSlJHo9HlZWV+tu//dvQeGVlpR588MGox+Tm5uqnP/1p2NjevXvl9XqVmJjYq/UCCMemfwD9Wcz3CElScXGxXnvtNW3evFm1tbVavHix6uvr5fP5JH2yd+exxx4Lzff5fDp9+rSKi4tVW1urzZs3q7y8XM8++2xoTltbm4LBoILBoNra2nT27FkFg0H9+te/7vP+AABAfIr5ipAkFRYW6vz581q1apUaGhqUk5OjiooKZWZmSpIaGhrC7inkdrtVUVGhxYsXa8OGDUpPT9e6des0b9680JwPPvhAkydPDr1fs2aN1qxZo7vvvlv79+/vs94AAED8iov7CMWrrtyHoDs4pQATsUcIQG/rV/cRAgAAiJW4ODUGwBxcTQYgnrAiBAAAjEUQAgAAxiIIAQAAYxGEAACAsQhCAADAWAQhAABgLIIQAAAwFkEIAAAYiyAEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYPH0eQFzgqfQAYoEVIQAAYCyCEAAAMBZBCAAAGIs9QgDiWrS9Q+wbAtBTWBECAADGIggBAABjEYQAAICxCEIAAMBYBCEAAGAsghAAADAWQQgAABiLIAQAAIxFEAIAAMYiCAEAAGPxiA0AthHtcRwSj+QA0DlWhAAAgLFYEQLQ73S28gMAXcWKEAAAMBZBCAAAGIsgBAAAjEUQAgAAxmKzNADb47J6AJ1hRQgAABiLIAQAAIxFEAIAAMYiCAEAAGMRhAAAgLG4agyAsbiaDABBCACuE4EK6L84NQYAAIzFihAAXCOeeg/YDytCAADAWAQhAABgLIIQAAAwFkEIAAAYi83SAPApbIoGzMGKEAAAMBZBCAAAGItTYwAQ57hzNdB7WBECAADGIggBAABjEYQAAICxCEIAAMBYbJYGgF7S1U3O3L8I6HusCAEAAGMRhAAAgLE4NQYA/RT3FwKuH0EIAPoYe4GA+EEQAgCbYaUIuHbsEQIAAMZiRQgADNGVU3KsHsEUrAgBAABjsSIEAIjAPiOYIm5WhEpKSuR2u+VyueTxeHTgwIGrzq+qqpLH45HL5VJ2drbKysoi5uzcuVPjx4+X0+nU+PHj9cYbb/RW+QBghJmvz4z66qn5QF+LiyC0Y8cOFRUVacWKFaqpqVFeXp4KCgpUX18fdX5dXZ1mz56tvLw81dTUaPny5Vq0aJF27twZmhMIBFRYWKj58+fr6NGjmj9/vv7u7/5Ov/zlL/uqLQAAEOcclmVZsS5i6tSpmjJlikpLS0Nj48aN09y5c+X3+yPmL1myRLt371ZtbW1ozOfz6ejRowoEApKkwsJCtbS06M033wzNuf/++zVkyBBt27btmupqaWlRSkqKmpublZyc3N32OsX/FQFA7+AUntm68vs75nuE2traVF1draVLl4aN5+fn6+DBg1GPCQQCys/PDxubNWuWysvLdenSJSUmJioQCGjx4sURc9auXdtpLa2trWptbQ29b25ulvTJ39De0P7n9l75XgAwXV5ZXkx+7s++9rPr/o45/29Oj3x3T31Pf3Tl9/a1rPXEPAg1NTWpo6NDqampYeOpqalqbGyMekxjY2PU+e3t7WpqatLIkSM7ndPZd0qS3+/XCy+8EDGekZFxre0AAAyW8p2UuP/u3qwx3ly4cEEpKVfvN+ZB6AqHwxH23rKsiLHPmv/p8a5+57Jly1RcXBx6f/nyZf3hD3/QsGHDrnpcd7S0tCgjI0NnzpzpldNu8YRe7cukfunVvkzq15ReLcvShQsXlJ6e/plzYx6Ehg8froEDB0as1Jw7dy5iReeKtLS0qPMTEhI0bNiwq87p7Dslyel0yul0ho3dfPPN19pKtyQnJ9v6H8a/RK/2ZVK/9GpfJvVrQq+ftRJ0RcyvGktKSpLH41FlZWXYeGVlpaZNmxb1mNzc3Ij5e/fuldfrVWJi4lXndPadAADAPDFfEZKk4uJizZ8/X16vV7m5udq4caPq6+vl8/kkfXLK6uzZs9qyZYukT64QW79+vYqLi7Vw4UIFAgGVl5eHXQ32zDPPaMaMGfqXf/kXPfjgg/qv//ovvfXWW/rFL34Rkx4BAED8iYsgVFhYqPPnz2vVqlVqaGhQTk6OKioqlJmZKUlqaGgIu6eQ2+1WRUWFFi9erA0bNig9PV3r1q3TvHnzQnOmTZum7du367nnntN3v/tdjR49Wjt27NDUqVP7vL9onE6nVq5cGXEqzo7o1b5M6pde7cukfk3q9VrFxX2EAAAAYiHme4QAAABihSAEAACMRRACAADGIggBAABjEYRioKSkRG63Wy6XSx6PRwcOHIh1ST3i7bff1gMPPKD09HQ5HA795Cc/Cfvcsiw9//zzSk9P1w033KB77rlHx48fj02x18nv9+uOO+7Q4MGDdcstt2ju3Lk6efJk2By79FtaWqoJEyaEbsCWm5sb9jBju/QZjd/vl8PhUFFRUWjMTv0+//zzcjgcYa+0tLTQ53bqVZLOnj2rr3/96xo2bJhuvPFGTZo0SdXV1aHP7dJvVlZWxJ+rw+HQU089Jck+ffYUglAf27Fjh4qKirRixQrV1NQoLy9PBQUFYbcH6K8++ugjTZw4UevXr4/6+SuvvKJXX31V69ev1+HDh5WWlqb77rtPFy5c6ONKr19VVZWeeuopHTp0SJWVlWpvb1d+fr4++uij0By79Hvrrbdq9erVOnLkiI4cOaIvfvGLevDBB0P/4bRLn592+PBhbdy4URMmTAgbt1u/X/jCF9TQ0BB6HTt2LPSZnXr94x//qOnTpysxMVFvvvmmTpw4oX/9138Ne3qAXfo9fPhw2J/plZsLP/zww5Ls02ePsdCn7rzzTsvn84WNjR071lq6dGmMKuodkqw33ngj9P7y5ctWWlqatXr16tDYxx9/bKWkpFhlZWUxqLBnnTt3zpJkVVVVWZZl/36HDBlivfbaa7bt88KFC9aYMWOsyspK6+6777aeeeYZy7Ls9+e6cuVKa+LEiVE/s1uvS5Ysse66665OP7dbv3/pmWeesUaPHm1dvnzZ1n12FytCfaitrU3V1dXKz88PG8/Pz9fBgwdjVFXfqKurU2NjY1jvTqdTd999ty16b25uliQNHTpUkn377ejo0Pbt2/XRRx8pNzfXtn0+9dRTmjNnju69996wcTv2e+rUKaWnp8vtduurX/2q3n//fUn263X37t3yer16+OGHdcstt2jy5MnatGlT6HO79XtFW1ubtm7dqieeeEIOh8O2fV4PglAfampqUkdHR8SDX1NTUyMeEGs3V/qzY++WZam4uFh33XWXcnJyJNmv32PHjmnQoEFyOp3y+Xx64403NH78eNv1KUnbt2/Xr371K/n9/ojP7Nbv1KlTtWXLFu3Zs0ebNm1SY2Ojpk2bpvPnz9uu1/fff1+lpaUaM2aM9uzZI5/Pp0WLFoUe3WS3fq/4yU9+og8//FDf+MY3JNm3z+sRF4/YMI3D4Qh7b1lWxJhd2bH3p59+Wu+9917U59jZpd/Pf/7zCgaD+vDDD7Vz504tWLBAVVVVoc/t0ueZM2f0zDPPaO/evXK5XJ3Os0u/BQUFob++/fbblZubq9GjR+v111/XX//1X0uyT6+XL1+W1+vVyy+/LEmaPHmyjh8/rtLSUj322GOheXbp94ry8nIVFBQoPT09bNxufV4PVoT60PDhwzVw4MCI1H3u3LmIdG43V65EsVvv//AP/6Ddu3dr3759uvXWW0Pjdus3KSlJt912m7xer/x+vyZOnKh///d/t12f1dXVOnfunDwejxISEpSQkKCqqiqtW7dOCQkJoZ7s0u+n3XTTTbr99tt16tQp2/3Zjhw5UuPHjw8bGzduXOhCFbv1K0mnT5/WW2+9pW9961uhMTv2eb0IQn0oKSlJHo8ntIP/isrKSk2bNi1GVfUNt9uttLS0sN7b2tpUVVXVL3u3LEtPP/20du3apZ///Odyu91hn9ut30+zLEutra226/NLX/qSjh07pmAwGHp5vV49+uijCgaDys7OtlW/n9ba2qra2lqNHDnSdn+206dPj7jFxf/93/+FHu5tt34l6Qc/+IFuueUWzZkzJzRmxz6vW4w2aRtr+/btVmJiolVeXm6dOHHCKioqsm666SbrN7/5TaxLu24XLlywampqrJqaGkuS9eqrr1o1NTXW6dOnLcuyrNWrV1spKSnWrl27rGPHjlmPPPKINXLkSKulpSXGlXfdd77zHSslJcXav3+/1dDQEHpdvHgxNMcu/S5btsx6++23rbq6Ouu9996zli9fbg0YMMDau3evZVn26bMzf3nVmGXZq99//Md/tPbv32+9//771qFDh6wvf/nL1uDBg0P/PbJTr++++66VkJBgvfTSS9apU6esH/3oR9aNN95obd26NTTHTv12dHRYo0aNspYsWRLxmZ367AkEoRjYsGGDlZmZaSUlJVlTpkwJXXLd3+3bt8+SFPFasGCBZVmfXJ66cuVKKy0tzXI6ndaMGTOsY8eOxbboborWpyTrBz/4QWiOXfp94oknQv+8jhgxwvrSl74UCkGWZZ8+O/PpIGSnfgsLC62RI0daiYmJVnp6uvXQQw9Zx48fD31up14ty7J++tOfWjk5OZbT6bTGjh1rbdy4MexzO/W7Z88eS5J18uTJiM/s1GdPcFiWZcVkKQoAACDG2CMEAACMRRACAADGIggBAABjEYQAAICxCEIAAMBYBCEAAGAsghAAADAWQQgAABiLIAQAAIxFEAJga7///e+VmJioixcvqr29XTfddFPoieOSlJWVJYfDoUOHDoUdV1RUpHvuuaePqwXQ1whCAGwtEAho0qRJuvHGG1VdXa2hQ4dq1KhRYXNcLpeWLFkSowoBxBJBCICtHTx4UNOnT5ck/eIXvwj99V/6+7//ex06dEgVFRV9XR6AGEuIdQEA0NPq6+s1YcIESdLFixc1cOBA/ed//qf+/Oc/y+Fw6Oabb9bXvvY1lZSUSPrk9JjP59OyZct0//33a8AA/h8RMAX/tgOwnfT0dAWDQb399tuSpEOHDqm6ulpJSUnau3evgsGgVq1aFXbMc889p7q6Ov3oRz+KRckAYoQgBMB2EhISlJWVpf/93//VHXfcoYkTJ6qxsVGpqamaMWOGsrKyNHz48LBjRowYoWeffVb/9E//pLa2thhVDqCvcWoMgO184Qtf0OnTp3Xp0iVdvnxZgwYNUnt7u9rb2zVo0CBlZmbq+PHjEccVFxerpKQkdMoMgP2xIgTAdioqKhQMBpWWlqatW7cqGAwqJydHa9euVTAY7HRT9KBBg/Td735XL730klpaWvq4agCxQBACYDuZmZkaNGiQfve73+nBBx/UqFGjdOLECT300EO67bbblJmZ2emx3/72t5WSkqJt27b1YcUAYoUgBMCW9u/frzvuuEMul0u//OUv9bnPfU7p6emfeVxiYqJefPFFffzxx31QJYBYc1iWZcW6CAAAgFhgRQgAABiLIAQAAIxFEAIAAMYiCAEAAGMRhAAAgLEIQgAAwFgEIQAAYCyCEAAAMBZBCAAAGIsgBAAAjEUQAgAAxvr/XNK0vKvV1GYAAAAASUVORK5CYII=",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot historgram \n",
    "n_bins = floor(Int,maximum(n_traj));\n",
    "PyPlot.xlabel(\"#N\"); PyPlot.ylabel(L\"P(n)\");\n",
    "PyPlot.hist(n_traj[floor(Int,total_time/(2*sp)):end], bins=n_bins, facecolor= \"g\", alpha=0.75, density=true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.2",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}