EighthSession.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a97ac065",
   "metadata": {},
   "source": [
    "# Eighth Session of Stochastic Modelling\n",
    "\n",
    "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/guglielmopadula/StochasticModelling2022-2023/blob/main/EightSession.ipynb)\n",
    "\n",
    "## A very fast recap on the Gillespie Algorithm\n",
    "\n",
    "Let's suppose we have a CTMC $X_{t}$ with rate matrix $Q$ and initial state $X_{0}$. To simulate it is necessary to build the jump matrix $A$ which has components\n",
    "$$a_{ii}=0$$\n",
    "$$a_{ij}=-\\frac{q_{ij}}{q_{ii}}$$\n",
    "Note that $$\\sum_{j=1}^{n}a_{ij}=1$$\n",
    "With this we can write the Gillespie algoritm:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5fd49644",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nPSEUDOCODE\\ndef Gillespie():\\n    $t=0$\\n    $X=X_{0}$\\n    while:\\n        draw deltat from Exponential(Q[x])\\n        t=d+delta t\\n        draw x from Categorical(A[x])\\n        save (t,x)\\n'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "'''\n",
    "PSEUDOCODE\n",
    "def Gillespie():\n",
    "    $t=0$\n",
    "    $X=X_{0}$\n",
    "    while:\n",
    "        draw deltat from Exponential(Q[x])\n",
    "        t=d+delta t\n",
    "        draw x from Categorical(A[x])\n",
    "        save (t,x)\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54f36ec7",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "- Using scipy.linalg.expm write an algorithm to compute the marginal probability of the CTMC at time t (not the conditional)\n",
    "- Write the python code of the gillespie algorithm \n",
    "- Simulate the marginal probability and the gillespie algorithm for the CTMC associated to the matrix \n",
    "$$\\begin{bmatrix} -1 & 1  \\\\ \n",
    " 2  & -2  \\end{bmatrix}$$\n",
    " \n",
    " \n",
    "### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2dc02a2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def CTMCsteadystate_probability(Q,p0,t):\n",
    "  l=[]\n",
    "  for i in t:\n",
    "    p=p0@scipy.linalg.expm(Q*i) \n",
    "    l.append(p)\n",
    "  return np.array(l)\n",
    "\n",
    "\n",
    "def CtmcSimulation(Q,x0,max_time):\n",
    "  A=Q.copy()\n",
    "  np.fill_diagonal(A,0)\n",
    "  A=(A.T/(np.sum(A,axis=1))).T\n",
    "  times=[]\n",
    "  states=[]\n",
    "  tot_time=0\n",
    "  times.append(tot_time)\n",
    "  x=x0\n",
    "  states.append(x)\n",
    "  while tot_time<max_time:\n",
    "    t=scipy.stats.expon(scale=-1/Q[x0,x0]).rvs(size=1)[0]\n",
    "    tot_time=tot_time+t\n",
    "    times.append(tot_time)\n",
    "    states.append(x)\n",
    "    x=np.random.choice(len(A),p=A[x])\n",
    "    times.append(tot_time)\n",
    "    states.append(x)\n",
    "  return times,states\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d093201f",
   "metadata": {},
   "outputs": [],
   "source": [
    "A=np.array([[-1,1],[2,-2]])\n",
    "p0=np.array([0.5,0.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "32729d84",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f8553674c40>,\n",
       " <matplotlib.lines.Line2D at 0x7f8553674ca0>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuW0lEQVR4nO3df3BU9b3/8dfuJtkQJAk2ZPPDSPihIMUQG2saf9RaVgPTsWjtd2IvLZqrOJemc9HUWnNboKjXWP3Kl3ov11jaFKwzQuu3alstlm6F+7UG0CADtRgIAoHChp/JhpQksHu+f8Rs3EtC9oRkPyF5PmbOJJw95+Szn9E5r3mfz+dzHJZlWQIAABjCnKYbAAAA0BcCCwAAGPIILAAAYMgjsAAAgCGPwAIAAIY8AgsAABjyCCwAAGDII7AAAIAhL850AwZCKBTSoUOHNGbMGDkcDtPNAQAAUbAsSy0tLcrKypLTef4ayrAILIcOHVJOTo7pZgAAgH44cOCALrvssvMeMywCy5gxYyR1fuHk5GTDrQEAANEIBALKyckJ38fPZ1gElq7HQMnJyQQWAAAuMtEM52DQLQAAGPIILAAAYMgjsAAAgCGPwAIAAIY8AgsAABjyCCwAAGDII7AAAIAhj8ACAACGPAILAAAY8ggsAABgyCOwAACAIY/AAgAAhrxh8fLDi0HbmaDqj5xS8+kzEVvH2ZDOhiwFQ5/8DFqyPjnHsqTufwEAYE6c06EffGWaub9v7C+PAEdb2vXnjxq1/m9H9E79UbWdCZluEgAA/ZIQ5ySwDDeHmk6r/FfbtHnvCVmfKpBcOjpBaZckKGVUvFJGxSt5VLwS412KczrkdDgU53TI5XRIDskhhxydvyqKt24DADCoXE6zo0gILAPsbDCkf335A72//6QkKe+yFHmv8sh7lUdXZY6Rg/QBAIBtBJYB9p9v1+v9/Sd1iTtOr377el3hGWO6SQAAXPSYJTSA3t93Qs/5dkuS/v3O6YQVAAAGCIFlgDSfPqOFa7YpZElfuyZbc/KzTTcJAIBhg8AyACzL0g9e3aG/N53W5Zcmaemcz5puEgAAwwqBZQD8361/1++3H5bL6dBP7s7XmMR4000CAGBYIbAMgGfe+kiSVH7rlbrm8rGGWwMAwPBDYLlAJ1o71BholySV3pBrtjEAAAxTBJYLtOfoKUlSduooJSUwSxwAgMFAYLlAe450BpaJ40YbbgkAAMMXgeUCdVVYJqdfYrglAAAMX/0KLCtWrFBubq4SExNVWFioLVu2nPf4pqYmlZWVKTMzU263W1deeaXefPPN8Oc/+tGP5HA4IrapU6f2p2kxV/9JhWXSOAILAACDxfagi7Vr16q8vFxVVVUqLCzU8uXLVVxcrLq6OqWnp59zfEdHh2699Valp6frlVdeUXZ2tvbv36/U1NSI4z772c/qT3/6U3fD4i6O8SB7jrZKosICAMBgsp0Kli1bpvnz56u0tFSSVFVVpTfeeEPV1dV69NFHzzm+urpaJ06c0Lvvvqv4+M71SXJzc89tSFycMjIy7DbHqLYzQR04+Q9JVFgAABhMth4JdXR0qLa2Vl6vt/sCTqe8Xq9qamp6POe3v/2tioqKVFZWJo/Ho+nTp+vJJ59UMBiMOG737t3KysrSxIkTNXfuXDU0NPTajvb2dgUCgYjNhL3HWmVZUsqoeKVdkmCkDQAAjAS2AsuxY8cUDAbl8Xgi9ns8Hvn9/h7P+fjjj/XKK68oGAzqzTff1KJFi/Tss8/qiSeeCB9TWFioVatWad26dXr++ee1d+9e3XTTTWppaenxmpWVlUpJSQlvOTk5dr7GgOkevzJaDofDSBsAABgJBn2gSCgUUnp6un7605/K5XKpoKBAf//73/XMM89oyZIlkqTZs2eHj8/Ly1NhYaHGjx+vX/3qV7rvvvvOuWZFRYXKy8vD/w4EAkZCCzOEAACIDVuBJS0tTS6XS42NjRH7Gxsbex1/kpmZqfj4eLlcrvC+q666Sn6/Xx0dHUpIOPdRSmpqqq688krV19f3eE232y23222n6YOCGUIAAMSGrUdCCQkJKigokM/nC+8LhULy+XwqKirq8ZwbbrhB9fX1CoVC4X27du1SZmZmj2FFkk6dOqU9e/YoMzPTTvNijhlCAADEhu11WMrLy7Vy5UqtXr1aO3fu1IIFC9Ta2hqeNTRv3jxVVFSEj1+wYIFOnDihhQsXateuXXrjjTf05JNPqqysLHzMww8/rI0bN2rfvn169913deedd8rlcukb3/jGAHzFwREMWfr4KBUWAABiwfYYlpKSEh09elSLFy+W3+9Xfn6+1q1bFx6I29DQIKezOwfl5OTorbfe0kMPPaS8vDxlZ2dr4cKF+v73vx8+5uDBg/rGN76h48ePa9y4cbrxxhu1adMmjRs3bgC+4uA41HRa7WdDSnA5lXNpkunmAAAwrDksy7JMN+JCBQIBpaSkqLm5WcnJyTH5m29/dESlq97TFM8YvfXQF2PyNwEAGE7s3L95l1A/MUMIAIDYIbD006fXYAEAAIOLwNJPXRWWSVRYAAAYdASWfmINFgAAYofA0g8nWjt08h9nJBFYAACIBQJLP3RVV7JTR2lUgquPowEAwIUisPQDM4QAAIgtAks/7GH8CgAAMUVg6Yd6KiwAAMQUgaUfwlOaWYMFAICYILDY1HYmqIMnT0uiwgIAQKwQWGz6+GirLEtKTYrXpaMTTDcHAIARgcBiU3j8yrhL5HA4DLcGAICRgcBi08GT/5Akjf8M41cAAIgVAotNpzuCkqRL3CwYBwBArBBYbGo70xlYEuMJLAAAxAqBxaa2MyFJkpvAAgBAzBBYbOqusNB1AADECnddm9rOdlZYEuOosAAAECsEFpu6KixuKiwAAMQMd12b2qmwAAAQcwQWm5glBABA7BFYbGpn0C0AADHHXdemrmnNVFgAAIgdAotNbWepsAAAEGvcdW0KzxJi0C0AADFDYLGJR0IAAMQegcUmVroFACD2uOvaYFlWeB0WHgkBABA7BBYbusKKRIUFAIBY4q5rQ/uZTwcWKiwAAMQKgcWGrinNLqdD8S66DgCAWOGua0N4wG0c3QYAQCxx57WBKc0AAJhBYLGBFx8CAGAGgcWG8Cq3zBACACCmuPPa0PbJtOZE1mABACCmCCw2UGEBAMAM7rw2tFNhAQDACAKLDbxHCAAAM/p1512xYoVyc3OVmJiowsJCbdmy5bzHNzU1qaysTJmZmXK73bryyiv15ptvXtA1TWhnlhAAAEbYDixr165VeXm5lixZoq1bt2rGjBkqLi7WkSNHejy+o6NDt956q/bt26dXXnlFdXV1WrlypbKzs/t9TVNYhwUAADNsB5Zly5Zp/vz5Ki0t1bRp01RVVaWkpCRVV1f3eHx1dbVOnDih1157TTfccINyc3N18803a8aMGf2+pik8EgIAwAxbd96Ojg7V1tbK6/V2X8DplNfrVU1NTY/n/Pa3v1VRUZHKysrk8Xg0ffp0PfnkkwoGg/2+Znt7uwKBQMQWC13vEnIz6BYAgJiyFViOHTumYDAoj8cTsd/j8cjv9/d4zscff6xXXnlFwWBQb775phYtWqRnn31WTzzxRL+vWVlZqZSUlPCWk5Nj52v0G4+EAAAwY9CfbYRCIaWnp+unP/2pCgoKVFJSoh/84Aeqqqrq9zUrKirU3Nwc3g4cODCALe4dj4QAADAjzs7BaWlpcrlcamxsjNjf2NiojIyMHs/JzMxUfHy8XK7uqsRVV10lv9+vjo6Ofl3T7XbL7XbbafqA6Kqw8EgIAIDYslUqSEhIUEFBgXw+X3hfKBSSz+dTUVFRj+fccMMNqq+vVygUCu/btWuXMjMzlZCQ0K9rmtI1hoUKCwAAsWX7zlteXq6VK1dq9erV2rlzpxYsWKDW1laVlpZKkubNm6eKiorw8QsWLNCJEye0cOFC7dq1S2+88YaefPJJlZWVRX3NoaKdMSwAABhh65GQJJWUlOjo0aNavHix/H6/8vPztW7duvCg2YaGBjmd3TkoJydHb731lh566CHl5eUpOztbCxcu1Pe///2orzlUtFNhAQDACIdlWZbpRlyoQCCglJQUNTc3Kzk5edD+zv+qelfv7Tup5+d+TrOvzhy0vwMAwEhg5/5NqcAGpjUDAGAGgcWGrmnNbh4JAQAQU9x5beieJUSFBQCAWCKw2BB+JMQ6LAAAxBSBxQYeCQEAYAZ3XhtYhwUAADMILFEKhSx1BLseCdFtAADEEnfeKLWf7X61ABUWAABii8ASpa7xKxKBBQCAWCOwRKlrSnO8yyGX02G4NQAAjCwEligxpRkAAHMILFHqntJMYAEAINYILFHqCiy8qRkAgNjj7hulrkdCbqY0AwAQc9x9o8R7hAAAMIfAEqX2MwQWAABMIbBEqWvhOMawAAAQe9x9oxQedMu0ZgAAYo7AEqU2XnwIAIAxBJYoda/DQpcBABBr3H2jRIUFAABzCCxRCk9rZgwLAAAxR2CJEo+EAAAwh7tvlHj5IQAA5hBYotTOu4QAADCGu2+UuheOo8ICAECsEViixNuaAQAwh7tvlHj5IQAA5hBYotQ16NbNoFsAAGKOwBIlHgkBAGAOd98odQcWKiwAAMQagSVK3Y+E6DIAAGKNu2+U2hl0CwCAMQSWKPHyQwAAzCGwRKm7wkKXAQAQa9x9oxAMWToTtCTxLiEAAEwgsESha4aQxCMhAABMILBE4dOBhVlCAADEHnffKLR98uLDhDinnE6H4dYAADDy9CuwrFixQrm5uUpMTFRhYaG2bNnS67GrVq2Sw+GI2BITEyOOuffee885ZtasWf1p2qAILxpHdQUAACPi7J6wdu1alZeXq6qqSoWFhVq+fLmKi4tVV1en9PT0Hs9JTk5WXV1d+N8Ox7lVilmzZukXv/hF+N9ut9tu0wZNV2BxM34FAAAjbJcMli1bpvnz56u0tFTTpk1TVVWVkpKSVF1d3es5DodDGRkZ4c3j8ZxzjNvtjjhm7Nixdps2aLrXYKHCAgCACbbuwB0dHaqtrZXX6+2+gNMpr9ermpqaXs87deqUxo8fr5ycHM2ZM0cffvjhOcds2LBB6enpmjJlihYsWKDjx4/3er329nYFAoGIbTC1hx8JUWEBAMAEW4Hl2LFjCgaD51RIPB6P/H5/j+dMmTJF1dXVev311/XSSy8pFArp+uuv18GDB8PHzJo1Sy+++KJ8Pp9+/OMfa+PGjZo9e7aCwWCP16ysrFRKSkp4y8nJsfM1bGtjWX4AAIyyPYbFrqKiIhUVFYX/ff311+uqq67SCy+8oMcff1ySdPfdd4c/v/rqq5WXl6dJkyZpw4YNmjlz5jnXrKioUHl5efjfgUBgUENLO4+EAAAwytYdOC0tTS6XS42NjRH7GxsblZGREdU14uPjdc0116i+vr7XYyZOnKi0tLRej3G73UpOTo7YBhMVFgAAzLIVWBISElRQUCCfzxfeFwqF5PP5Iqoo5xMMBrVjxw5lZmb2eszBgwd1/Pjx8x4TS12Dbt2MYQEAwAjbzzjKy8u1cuVKrV69Wjt37tSCBQvU2tqq0tJSSdK8efNUUVERPv6xxx7TH//4R3388cfaunWrvvnNb2r//v26//77JXUOyP3e976nTZs2ad++ffL5fJozZ44mT56s4uLiAfqaFya8DguPhAAAMML2GJaSkhIdPXpUixcvlt/vV35+vtatWxceiNvQ0CCns/vGfvLkSc2fP19+v19jx45VQUGB3n33XU2bNk2S5HK5tH37dq1evVpNTU3KysrSbbfdpscff3zIrMXSPa2ZCgsAACY4LMuyTDfiQgUCAaWkpKi5uXlQxrP8n/W79BPfbn3zC5friTuuHvDrAwAwEtm5f/OMIwpdg24ZwwIAgBkEligwrRkAALO4A0ehjZVuAQAwisAShfazDLoFAMAkAksUmNYMAIBZ3IGj0BVY3FRYAAAwgsASBdZhAQDALAJLFMLvEoqjuwAAMIE7cBSosAAAYBaBJQrtXWNYqLAAAGAEd+AodM8SosICAIAJBJYotLEOCwAARhFYotDOOiwAABjFHTgKVFgAADCLwNKHM8GQgiFLEu8SAgDAFAJLH7oG3EqSm0dCAAAYwR24D11rsDgcTGsGAMAU7sB9aPvUGiwOh8NwawAAGJkILH1oP9sVWBi/AgCAKQSWPnQvy09XAQBgCnfhPrDKLQAA5hFY+hCusPBICAAAYwgsfegaw8IjIQAAzOEu3IeuCoubR0IAABhDYOkDY1gAADCPwNKHtq5HQiwaBwCAMdyF+9A9rZkKCwAAphBY+vDplW4BAIAZ3IX70M4YFgAAjCOw9KHtLCvdAgBgGnfhPjBLCAAA8wgsfWhn0C0AAMYRWPrQdpZBtwAAmMZduA88EgIAwDwCSx9YhwUAAPMILH3orrDQVQAAmMJduA9d05rdcVRYAAAwhcDSh3YqLAAAGMdduA8MugUAwLx+BZYVK1YoNzdXiYmJKiws1JYtW3o9dtWqVXI4HBFbYmJixDGWZWnx4sXKzMzUqFGj5PV6tXv37v40bcCFB93ySAgAAGNsB5a1a9eqvLxcS5Ys0datWzVjxgwVFxfryJEjvZ6TnJysw4cPh7f9+/dHfP7000/rueeeU1VVlTZv3qzRo0eruLhYbW1t9r/RAOtah4VHQgAAmGP7Lrxs2TLNnz9fpaWlmjZtmqqqqpSUlKTq6upez3E4HMrIyAhvHo8n/JllWVq+fLl++MMfas6cOcrLy9OLL76oQ4cO6bXXXuvXlxpIrHQLAIB5tgJLR0eHamtr5fV6uy/gdMrr9aqmpqbX806dOqXx48crJydHc+bM0Ycffhj+bO/evfL7/RHXTElJUWFhYa/XbG9vVyAQiNgGg3WmTd+1VuuxuF/I7Tw7KH8DAAD0zVZgOXbsmILBYESFRJI8Ho/8fn+P50yZMkXV1dV6/fXX9dJLLykUCun666/XwYMHJSl8np1rVlZWKiUlJbzl5OTY+RpR6wiGdH/cm5oXt16J6hiUvwEAAPo26AMzioqKNG/ePOXn5+vmm2/Wb37zG40bN04vvPBCv69ZUVGh5ubm8HbgwIEBbPGnuBJkySFJGkVgAQDAGFuBJS0tTS6XS42NjRH7GxsblZGREdU14uPjdc0116i+vl6SwufZuabb7VZycnLENhjc8XFyxI/qbLdFYAEAwBRbgSUhIUEFBQXy+XzhfaFQSD6fT0VFRVFdIxgMaseOHcrMzJQkTZgwQRkZGRHXDAQC2rx5c9TXHFRxn0zBPmN+xhIAACNVnN0TysvLdc899+jaa6/Vddddp+XLl6u1tVWlpaWSpHnz5ik7O1uVlZWSpMcee0xf+MIXNHnyZDU1NemZZ57R/v37df/990vqnEH04IMP6oknntAVV1yhCRMmaNGiRcrKytIdd9wxcN+0v+JHSaclnT1tuiUAAIxYtgNLSUmJjh49qsWLF8vv9ys/P1/r1q0LD5ptaGiQ09lduDl58qTmz58vv9+vsWPHqqCgQO+++66mTZsWPuaRRx5Ra2urHnjgATU1NenGG2/UunXrzllgzggqLAAAGOewLMsy3YgLFQgElJKSoubm5oEfz/L8DVLjX6VvvSpN+vLAXhsAgBHMzv2b5Vv7Eufu/EmFBQAAYwgsfYnrnCXEGBYAAMwhsPQlnjEsAACYRmDpS9egWyosAAAYQ2DpS3zXI6F2s+0AAGAEI7D0JTytmQoLAACmEFj6Eq6wMIYFAABTCCx9ocICAIBxBJa+UGEBAMA4AktfWDgOAADjCCx9YeE4AACMI7D0hYXjAAAwjsDSFyosAAAYR2DpCxUWAACMI7D0JY5ZQgAAmEZg6UtXhYXAAgCAMQSWvnRVWFg4DgAAYwgsfaHCAgCAcQSWvlBhAQDAOAJLX7pWuqXCAgCAMQSWvnz6XUKWZbYtAACMUASWvnS9rVmiygIAgCEElr50VVgkxrEAAGAIgaUvrnjJ4er8/Wy72bYAADBCEViiEc/7hAAAMInAEo043icEAIBJBJZoUGEBAMAoAks0qLAAAGAUgSUaXYGFCgsAAEYQWKIRT4UFAACTCCzRiOMFiAAAmERgiUY8L0AEAMAkAks0qLAAAGAUgSUan34BIgAAiDkCSzSY1gwAgFEElmiwcBwAAEYRWKJBhQUAAKMILNGgwgIAgFEElmjEuTt/UmEBAMAIAks04qiwAABgUr8Cy4oVK5Sbm6vExEQVFhZqy5YtUZ23Zs0aORwO3XHHHRH77733Xjkcjoht1qxZ/Wna4GBpfgAAjLIdWNauXavy8nItWbJEW7du1YwZM1RcXKwjR46c97x9+/bp4Ycf1k033dTj57NmzdLhw4fD28svv2y3aYMnjnVYAAAwyXZgWbZsmebPn6/S0lJNmzZNVVVVSkpKUnV1da/nBINBzZ07V0uXLtXEiRN7PMbtdisjIyO8jR071m7TBk88K90CAGCSrcDS0dGh2tpaeb3e7gs4nfJ6vaqpqen1vMcee0zp6em67777ej1mw4YNSk9P15QpU7RgwQIdP36812Pb29sVCAQitkEVx7uEAAAwyVZgOXbsmILBoDweT8R+j8cjv9/f4znvvPOOfv7zn2vlypW9XnfWrFl68cUX5fP59OMf/1gbN27U7NmzFQwGezy+srJSKSkp4S0nJ8fO17CPCgsAAEbFDebFW1pa9K1vfUsrV65UWlpar8fdfffd4d+vvvpq5eXladKkSdqwYYNmzpx5zvEVFRUqLy8P/zsQCAxuaKHCAgCAUbYCS1pamlwulxobGyP2NzY2KiMj45zj9+zZo3379un2228P7wuFQp1/OC5OdXV1mjRp0jnnTZw4UWlpaaqvr+8xsLjdbrndbjtNvzBUWAAAMMrWI6GEhAQVFBTI5/OF94VCIfl8PhUVFZ1z/NSpU7Vjxw5t27YtvH31q1/VLbfcom3btvVaFTl48KCOHz+uzMxMm19nkLA0PwAARtl+JFReXq577rlH1157ra677jotX75cra2tKi0tlSTNmzdP2dnZqqysVGJioqZPnx5xfmpqqiSF9586dUpLly7VXXfdpYyMDO3Zs0ePPPKIJk+erOLi4gv8egOkK7CwcBwAAEbYDiwlJSU6evSoFi9eLL/fr/z8fK1bty48ELehoUFOZ/SFG5fLpe3bt2v16tVqampSVlaWbrvtNj3++OOxfexzPl3vEgp2SKGg5HSZbQ8AACOMw7Isy3QjLlQgEFBKSoqam5uVnJw88H+g/ZRUmd35+78dlhKSBv5vAAAwwti5f/MuoWh0VVgkBt4CAGAAgSUaTpfkjO/8nanNAADEHIElWvG8TwgAAFMILNEKT22mwgIAQKwRWKLF4nEAABhDYIkWFRYAAIwhsEQrjgoLAACmEFiiFc8LEAEAMIXAEq1whaXdbDsAABiBCCzRCk9rpsICAECsEViixRubAQAwhsASLSosAAAYQ2CJFhUWAACMIbBEiwoLAADGEFiiRYUFAABjCCzRCk9rpsICAECsEViiFU+FBQAAUwgs0YrrGsNCYAEAINYILNHibc0AABhDYIlWHO8SAgDAFAJLtKiwAABgDIElWlRYAAAwhsASLSosAAAYQ2CJVrjCQmABACDWCCzRinN3/mThOAAAYo7AEq14KiwAAJhCYIlWHGNYAAAwhcASra4KS+iMFAqabQsAACMMgSVaXRUWianNAADEGIElWp8OLDwWAgAgpggs0XI6JdcnM4WosAAAEFMEFjtYPA4AACMILHawPD8AAEYQWOygwgIAgBEEFju6Bt5SYQEAIKYILHaEF49rN9sOAABGGAKLHV2Lx/E+IQAAYorAYkf4kRBjWAAAiCUCix1UWAAAMKJfgWXFihXKzc1VYmKiCgsLtWXLlqjOW7NmjRwOh+64446I/ZZlafHixcrMzNSoUaPk9Xq1e/fu/jRtcFFhAQDACNuBZe3atSovL9eSJUu0detWzZgxQ8XFxTpy5Mh5z9u3b58efvhh3XTTTed89vTTT+u5555TVVWVNm/erNGjR6u4uFhtbUMsGFBhAQDACNuBZdmyZZo/f75KS0s1bdo0VVVVKSkpSdXV1b2eEwwGNXfuXC1dulQTJ06M+MyyLC1fvlw//OEPNWfOHOXl5enFF1/UoUOH9Nprr9n+QoOKCgsAAEbYCiwdHR2qra2V1+vtvoDTKa/Xq5qaml7Pe+yxx5Senq777rvvnM/27t0rv98fcc2UlBQVFhb2es329nYFAoGILSaosAAAYIStwHLs2DEFg0F5PJ6I/R6PR36/v8dz3nnnHf385z/XypUre/y86zw716ysrFRKSkp4y8nJsfM1+i+u6+WHVFgAAIilQZ0l1NLSom9961tauXKl0tLSBuy6FRUVam5uDm8HDhwYsGufVxwVFgAATIizc3BaWppcLpcaGxsj9jc2NiojI+Oc4/fs2aN9+/bp9ttvD+8LhUKdfzguTnV1deHzGhsblZmZGXHN/Pz8HtvhdrvldrvtNH1gxLPSLQAAJtiqsCQkJKigoEA+ny+8LxQKyefzqaio6Jzjp06dqh07dmjbtm3h7atf/apuueUWbdu2TTk5OZowYYIyMjIirhkIBLR58+Yer2kU7xICAMAIWxUWSSovL9c999yja6+9Vtddd52WL1+u1tZWlZaWSpLmzZun7OxsVVZWKjExUdOnT484PzU1VZIi9j/44IN64okndMUVV2jChAlatGiRsrKyzlmvxbjwoFvGsAAAEEu2A0tJSYmOHj2qxYsXy+/3Kz8/X+vWrQsPmm1oaJDTaW9ozCOPPKLW1lY98MADampq0o033qh169YpMTHRbvMGFxUWAACMcFiWZZluxIUKBAJKSUlRc3OzkpOTB+8PffSGtOafpMs+L93/p8H7OwAAjAB27t+8S8gOFo4DAMAIAosdLBwHAIARBBY7WDgOAAAjCCx2sHAcAABGEFjsYOE4AACMILDY0VVhYVozAAAxRWCxo6vCYgWl4BmzbQEAYAQhsNjRVWGRqLIAABBDBBY74tySHJ2/szw/AAAxQ2Cxw+FgeX4AAAwgsNgVnilEhQUAgFghsNjFTCEAAGKOwGJX12q3VFgAAIgZAotd4fcJEVgAAIgVAotdvLEZAICYI7DYxRubAQCIOQKLXVRYAACIOQKLXVRYAACIOQKLXVRYAACIOQKLXeGF46iwAAAQKwQWu8ILx1FhAQAgVggsdoUXjqPCAgBArBBY7AoPum032w4AAEYQAotd8UmdP9tbzLYDAIARhMBi19jczp/H9xhtBgAAIwmBxa60Kzt/HtslWZbZtgAAMEIQWOz6zGRJDqmtSWo9aro1AACMCAQWu+ITpbHjO38/tstsWwAAGCEILP2RNqXz59E6s+0AAGCEILD0R9oVnT+P7TbbDgAARggCS3+M+6TCcowKCwAAsUBg6Y/wTCEqLAAAxAKBpT+6AkvzAan9lNm2AAAwAhBY+iPpUikprfP34/Vm2wIAwAhAYOmv8DgWpjYDADDYCCz9FZ4pRGABAGCwEVj6i7VYAACIGQJLfzFTCACAmCGw9Ne4TwLL8XopeNZsWwAAGOb6FVhWrFih3NxcJSYmqrCwUFu2bOn12N/85je69tprlZqaqtGjRys/P1+//OUvI46599575XA4IrZZs2b1p2mxk3yZFDdKCp2Rmvabbg0AAMOa7cCydu1alZeXa8mSJdq6datmzJih4uJiHTlypMfjL730Uv3gBz9QTU2Ntm/frtLSUpWWluqtt96KOG7WrFk6fPhweHv55Zf7941ixemU0iZ3/s44FgAABpXtwLJs2TLNnz9fpaWlmjZtmqqqqpSUlKTq6uoej//Sl76kO++8U1dddZUmTZqkhQsXKi8vT++8807EcW63WxkZGeFt7Nix/ftGsZTG1GYAAGLBVmDp6OhQbW2tvF5v9wWcTnm9XtXU1PR5vmVZ8vl8qqur0xe/+MWIzzZs2KD09HRNmTJFCxYs0PHjx+00zYzwwFsCCwAAgynOzsHHjh1TMBiUx+OJ2O/xePTRRx/1el5zc7Oys7PV3t4ul8ul//qv/9Ktt94a/nzWrFn62te+pgkTJmjPnj36t3/7N82ePVs1NTVyuVznXK+9vV3t7e3hfwcCATtfY+CMI7AAABALtgJLf40ZM0bbtm3TqVOn5PP5VF5erokTJ+pLX/qSJOnuu+8OH3v11VcrLy9PkyZN0oYNGzRz5sxzrldZWamlS5fGounn11VhObpLsizJ4TDbHgAAhilbj4TS0tLkcrnU2NgYsb+xsVEZGRm9/xGnU5MnT1Z+fr6++93v6utf/7oqKyt7PX7ixIlKS0tTfX3P7+mpqKhQc3NzeDtw4ICdrzFwLp0kOZxSe7N0qudBxwAA4MLZCiwJCQkqKCiQz+cL7wuFQvL5fCoqKor6OqFQKOKRzv908OBBHT9+XJmZmT1+7na7lZycHLEZEZ8opY7v/P0YM4UAABgstmcJlZeXa+XKlVq9erV27typBQsWqLW1VaWlpZKkefPmqaKiInx8ZWWl1q9fr48//lg7d+7Us88+q1/+8pf65je/KUk6deqUvve972nTpk3at2+ffD6f5syZo8mTJ6u4uHiAvuYg4iWIAAAMOttjWEpKSnT06FEtXrxYfr9f+fn5WrduXXggbkNDg5zO7hzU2tqqb3/72zp48KBGjRqlqVOn6qWXXlJJSYkkyeVyafv27Vq9erWampqUlZWl2267TY8//rjcbvcAfc1BlHaFtGtd5zgWAAAwKByWZVmmG3GhAoGAUlJS1NzcHPvHQ1t/Kf32O9LEW6R5r8X2bwMAcBGzc//mXUIXirVYAAAYdASWC5V2RefPwN+l9hazbQEAYJgisFyopEullJzO3//ynNm2AAAwTBFYBoL3R50//9//lvb9xWhTAAAYjggsA+Hqr0v5cyUrJP3mAen0SdMtAgBgWCGwDJTZP5YunSgFDkq/W9i5VD8AABgQBJaB4h4j3fUzyRkn/e116YOXTLcIAIBhg8AykLILpC8v6vz9D49If99qtj0AAAwTBJaBdv2/ShO+KJ35h7TyFuk/Py+tXyw1bJJCQdOtAwDgosRKt4OhpVH63b9K9X+SQme798cnSUlp0qgUKTFVGpUqxY3qfIzkdHb+dLgkh0OS43/8jEa0xwEAYJPTJRX/+4Be0s792/a7hBCFMR7pn9ZKbc2doaXuD9LuP3b+u7lBajbdQAAAbHK5Bzyw2EFgGUyJKdL0uzq34Bnp5H6prUk63fTJz5NSsKOzChM62/nIKBSUZH0yy8jqZbbRRV8UAwBcbJxmIwOBJVZc8VLaZNOtAADgosSgWwAAMOQRWAAAwJBHYAEAAEMegQUAAAx5BBYAADDkEVgAAMCQR2ABAABDHoEFAAAMeQQWAAAw5BFYAADAkEdgAQAAQx6BBQAADHkEFgAAMOQNi7c1W5YlSQoEAoZbAgAAotV13+66j5/PsAgsLS0tkqScnBzDLQEAAHa1tLQoJSXlvMc4rGhizRAXCoV06NAhjRkzRg6HY0CvHQgElJOTowMHDig5OXlAr41u9HNs0M+xQ1/HBv0cG4PVz5ZlqaWlRVlZWXI6zz9KZVhUWJxOpy677LJB/RvJycn8zxAD9HNs0M+xQ1/HBv0cG4PRz31VVrow6BYAAAx5BBYAADDkEVj64Ha7tWTJErndbtNNGdbo59ign2OHvo4N+jk2hkI/D4tBtwAAYHijwgIAAIY8AgsAABjyCCwAAGDII7AAAIAhj8DShxUrVig3N1eJiYkqLCzUli1bTDfpovbf//3fuv3225WVlSWHw6HXXnst4nPLsrR48WJlZmZq1KhR8nq92r17t5nGXsQqKyv1+c9/XmPGjFF6erruuOMO1dXVRRzT1tamsrIyfeYzn9Ell1yiu+66S42NjYZafHF6/vnnlZeXF15Mq6ioSH/4wx/Cn9PHg+Opp56Sw+HQgw8+GN5HX1+4H/3oR3I4HBHb1KlTw5+b7mMCy3msXbtW5eXlWrJkibZu3aoZM2aouLhYR44cMd20i1Zra6tmzJihFStW9Pj5008/reeee05VVVXavHmzRo8ereLiYrW1tcW4pRe3jRs3qqysTJs2bdL69et15swZ3XbbbWptbQ0f89BDD+l3v/udfv3rX2vjxo06dOiQvva1rxls9cXnsssu01NPPaXa2lq9//77+vKXv6w5c+boww8/lEQfD4b33ntPL7zwgvLy8iL209cD47Of/awOHz4c3t55553wZ8b72EKvrrvuOqusrCz872AwaGVlZVmVlZUGWzV8SLJeffXV8L9DoZCVkZFhPfPMM+F9TU1Nltvttl5++WUDLRw+jhw5YkmyNm7caFlWZ7/Gx8dbv/71r8PH7Ny505Jk1dTUmGrmsDB27FjrZz/7GX08CFpaWqwrrrjCWr9+vXXzzTdbCxcutCyL/54HypIlS6wZM2b0+NlQ6GMqLL3o6OhQbW2tvF5veJ/T6ZTX61VNTY3Blg1fe/fuld/vj+jzlJQUFRYW0ucXqLm5WZJ06aWXSpJqa2t15syZiL6eOnWqLr/8cvq6n4LBoNasWaPW1lYVFRXRx4OgrKxMX/nKVyL6VOK/54G0e/duZWVlaeLEiZo7d64aGhokDY0+HhYvPxwMx44dUzAYlMfjidjv8Xj00UcfGWrV8Ob3+yWpxz7v+gz2hUIhPfjgg7rhhhs0ffp0SZ19nZCQoNTU1Ihj6Wv7duzYoaKiIrW1temSSy7Rq6++qmnTpmnbtm308QBas2aNtm7dqvfee++cz/jveWAUFhZq1apVmjJlig4fPqylS5fqpptu0l//+tch0ccEFmCYKysr01//+teIZ9EYOFOmTNG2bdvU3NysV155Rffcc482btxoulnDyoEDB7Rw4UKtX79eiYmJppszbM2ePTv8e15engoLCzV+/Hj96le/0qhRowy2rBOPhHqRlpYml8t1zgjoxsZGZWRkGGrV8NbVr/T5wPnOd76j3//+93r77bd12WWXhfdnZGSoo6NDTU1NEcfT1/YlJCRo8uTJKigoUGVlpWbMmKGf/OQn9PEAqq2t1ZEjR/S5z31OcXFxiouL08aNG/Xcc88pLi5OHo+Hvh4EqampuvLKK1VfXz8k/nsmsPQiISFBBQUF8vl84X2hUEg+n09FRUUGWzZ8TZgwQRkZGRF9HggEtHnzZvrcJsuy9J3vfEevvvqq/vznP2vChAkRnxcUFCg+Pj6ir+vq6tTQ0EBfX6BQKKT29nb6eADNnDlTO3bs0LZt28Lbtddeq7lz54Z/p68H3qlTp7Rnzx5lZmYOjf+eYzK09yK1Zs0ay+12W6tWrbL+9re/WQ888ICVmppq+f1+0027aLW0tFgffPCB9cEHH1iSrGXLllkffPCBtX//fsuyLOupp56yUlNTrddff93avn27NWfOHGvChAnW6dOnDbf84rJgwQIrJSXF2rBhg3X48OHw9o9//CN8zL/8y79Yl19+ufXnP//Zev/9962ioiKrqKjIYKsvPo8++qi1ceNGa+/evdb27dutRx991HI4HNYf//hHy7Lo48H06VlClkVfD4Tvfve71oYNG6y9e/daf/nLXyyv12ulpaVZR44csSzLfB8TWPrwH//xH9bll19uJSQkWNddd521adMm0026qL399tuWpHO2e+65x7KszqnNixYtsjwej+V2u62ZM2dadXV1Zht9EeqpjyVZv/jFL8LHnD592vr2t79tjR071kpKSrLuvPNO6/Dhw+YafRH653/+Z2v8+PFWQkKCNW7cOGvmzJnhsGJZ9PFg+p+Bhb6+cCUlJVZmZqaVkJBgZWdnWyUlJVZ9fX34c9N97LAsy4pNLQcAAKB/GMMCAACGPAILAAAY8ggsAABgyCOwAACAIY/AAgAAhjwCCwAAGPIILAAAYMgjsAAAgCGPwAIAAIY8AgsAABjyCCwAAGDII7AAAIAh7/8Dj48h7230whYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t=np.linspace(0,50,100)\n",
    "p=CTMCsteadystate_probability(A,p0,t)\n",
    "plt.plot(t,p)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9953b341",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f854f7d29e0>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhL0lEQVR4nO3dfXBU5R328SsJZAMFAjaygRAJ1hekSMCkZFbLqHVrah0qfXkmg44wqeKjJjNA2iqoJFJbQ22haBtNRSn+UQvqFPoCxWIk+FiiSCCjtoqiaBhhExgrwQAJZs/zB8PKhk2yZ5Plx4bvZ2ZHc859n/t37nNyuObs2WyS4ziOAAAAjCRbFwAAAM5thBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYGmBdQDSCwaD27dunoUOHKikpybocAAAQBcdxdPjwYY0ePVrJyV3f/0iIMLJv3z5lZ2dblwEAAGKwd+9ejRkzpsv1CRFGhg4dKunEzgwbNsy4GgAAEI2WlhZlZ2eH/h3vSkKEkZNvzQwbNowwAgBAgunpEQseYAUAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKZch5FXXnlF06dP1+jRo5WUlKR169b12Ke2tlZXXHGFPB6PLrroIq1atSqGUgEAQH/kOoy0trYqNzdXVVVVUbXfs2ePbrzxRl177bVqaGjQvHnzdPvtt+vFF190XSwAAOh/XH83zQ033KAbbrgh6vbV1dUaN26cli5dKkm67LLL9Oqrr+q3v/2tCgsL3Q5/xjiOo6PHOyKuSxuQomNfdGjQwJTT/t5+d/1iFWmcSE4d+2SfaOrpavtd9Y22nng6WduggSmSFNOcu92PzvPb3bjWc+T2PDy13lPn1u15F8t40Wyzu/aRxu/NsZW+/B2Pt3hcQyIdy67Gciuac8NN/V1dp6KpNZY+JwWDjiQpOblv5iMSN8e2c9szee3t63MkVnH/ory6ujr5/f6wZYWFhZo3b16Xfdra2tTW1hb6uaWlJV7lReQ4jn5UXaf6j//Xbbv8sSP0/J2+sF/8aPq51XmcaGrOHztCz/1fn/7PH9zvR6Ttua0nnk6tLW/sCEmKac7d7Efn+ehpXMs5iuU8PFmvpFDfWM47t+N1FYI7b7Or9l2N35tjeybF4xoS6VhGGsutU2tzezy6qzXSdaqnWt2cI50Fg44uvG+DJOnDh78bcyDpaV/dHNtT257pa+/R4x2aUH7inYr//rxQg1Ntvj837g+wBgIBeb3esGVer1ctLS06evRoxD6VlZVKT08PvbKzs+NdZpijxzui+mXa/vH/wtJrtP3c6jxOJJ3H3v7x//TpkfaY9iPS9tzWE0+n1lb/8f9innM3+9F5Pnoa13KOYjkPT9Z7at9Yzju340W7za7adzV+b47tmRSPa0ikYxlpLLeiOTfc1t/VdaqnWt2cI519eqQ94v+71dO+ujm2p7Y9m6+98WQTgXqwcOFClZWVhX5uaWk544HkpO0P+DU49cQt+SPtHcr/xUuu+8XKzXix1BPt9k/27ct64iHaOe/r/Yj1HDkTepqTeJ9jvR3v/91zraY9sjnq8SX1an86j9cXv8eRRDMPbsdOpHOvq1rdHO/e9Olrbq8B0V5TE+Xa2xfiHkYyMzPV1NQUtqypqUnDhg3ToEGDIvbxeDzyeDzxLi0qg1NTYrptFWu/eOlNPWfbvnTFqs6zeX7OZG3xGGuQi3+M+yI0dB7P8tiezedVNGKp383x7k2fvuZ2X6Ntn+jngBtxf5vG5/OppqYmbNmmTZvk8/niPTQAAEgArsPI559/roaGBjU0NEg68dHdhoYGNTY2SjrxFsusWbNC7e+88059+OGHuueee/Tuu+/q8ccf13PPPaf58+f3zR4AAICE5jqMbN++XVOmTNGUKVMkSWVlZZoyZYrKy8slSfv37w8FE0kaN26c1q9fr02bNik3N1dLly7VU089dVZ/rBcAAJw5rt+Muuaaa+Q4TpfrI/111WuuuUY7d+50OxQAADgH8N00AADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAVExhpKqqSjk5OUpLS1NBQYG2bdvWbfvly5fr0ksv1aBBg5Sdna358+fr2LFjMRUMAAD6F9dhZM2aNSorK1NFRYV27Nih3NxcFRYWqrm5OWL7Z599VgsWLFBFRYXeeecdPf3001qzZo3uu+++XhcPAAASn+swsmzZMs2ZM0fFxcWaMGGCqqurNXjwYK1cuTJi+61bt+qqq67SzTffrJycHF1//fWaOXNmj3dTAADAucFVGGlvb1d9fb38fv+XG0hOlt/vV11dXcQ+V155perr60Ph48MPP9SGDRv03e9+t8tx2tra1NLSEvYCAAD90wA3jQ8ePKiOjg55vd6w5V6vV++++27EPjfffLMOHjyob37zm3IcR1988YXuvPPObt+mqays1OLFi92UBgAAElTcP01TW1urhx9+WI8//rh27Nihv/zlL1q/fr0eeuihLvssXLhQhw4dCr327t0b7zIBAIARV3dGMjIylJKSoqamprDlTU1NyszMjNhn0aJFuvXWW3X77bdLki6//HK1trbqjjvu0P3336/k5NPzkMfjkcfjcVMaAABIUK7ujKSmpiovL081NTWhZcFgUDU1NfL5fBH7HDly5LTAkZKSIklyHMdtvQAAoJ9xdWdEksrKyjR79mzl5+dr6tSpWr58uVpbW1VcXCxJmjVrlrKyslRZWSlJmj59upYtW6YpU6aooKBAu3fv1qJFizR9+vRQKAEAAOcu12GkqKhIBw4cUHl5uQKBgCZPnqyNGzeGHmptbGwMuxPywAMPKCkpSQ888IA++eQTnX/++Zo+fbp++ctf9t1eAACAhOU6jEhSaWmpSktLI66rra0NH2DAAFVUVKiioiKWoQAAQD/Hd9MAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYIowAAABThBEAAGCKMAIAAEwRRgAAgCnCCAAAMEUYAQAApggjAADAFGEEAACYiimMVFVVKScnR2lpaSooKNC2bdu6bf/ZZ5+ppKREo0aNksfj0SWXXKINGzbEVDAAAOhfBrjtsGbNGpWVlam6uloFBQVavny5CgsLtWvXLo0cOfK09u3t7fr2t7+tkSNH6oUXXlBWVpY+/vhjDR8+vC/qBwAACc51GFm2bJnmzJmj4uJiSVJ1dbXWr1+vlStXasGCBae1X7lypT799FNt3bpVAwcOlCTl5OT0rmoAANBvuHqbpr29XfX19fL7/V9uIDlZfr9fdXV1Efv87W9/k8/nU0lJibxeryZOnKiHH35YHR0dXY7T1tamlpaWsBcAAOifXIWRgwcPqqOjQ16vN2y51+tVIBCI2OfDDz/UCy+8oI6ODm3YsEGLFi3S0qVL9Ytf/KLLcSorK5Wenh56ZWdnuykTAAAkkLh/miYYDGrkyJF68sknlZeXp6KiIt1///2qrq7uss/ChQt16NCh0Gvv3r3xLhMAABhx9cxIRkaGUlJS1NTUFLa8qalJmZmZEfuMGjVKAwcOVEpKSmjZZZddpkAgoPb2dqWmpp7Wx+PxyOPxuCkNAAAkKFd3RlJTU5WXl6eamprQsmAwqJqaGvl8voh9rrrqKu3evVvBYDC07L333tOoUaMiBhEAAHBucf02TVlZmVasWKFnnnlG77zzju666y61traGPl0za9YsLVy4MNT+rrvu0qeffqq5c+fqvffe0/r16/Xwww+rpKSk7/YCAAAkLNcf7S0qKtKBAwdUXl6uQCCgyZMna+PGjaGHWhsbG5Wc/GXGyc7O1osvvqj58+dr0qRJysrK0ty5c3Xvvff23V4AAICE5TqMSFJpaalKS0sjrqutrT1tmc/n02uvvRbLUAAAoJ/ju2kAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMxRRGqqqqlJOTo7S0NBUUFGjbtm1R9Vu9erWSkpI0Y8aMWIYFAAD9kOswsmbNGpWVlamiokI7duxQbm6uCgsL1dzc3G2/jz76SD/96U81bdq0mIsFAAD9j+swsmzZMs2ZM0fFxcWaMGGCqqurNXjwYK1cubLLPh0dHbrlllu0ePFiXXjhhb0qGAAA9C+uwkh7e7vq6+vl9/u/3EBysvx+v+rq6rrs9/Of/1wjR47UbbfdFtU4bW1tamlpCXsBAID+yVUYOXjwoDo6OuT1esOWe71eBQKBiH1effVVPf3001qxYkXU41RWVio9PT30ys7OdlMmAABIIHH9NM3hw4d16623asWKFcrIyIi638KFC3Xo0KHQa+/evXGsEgAAWBrgpnFGRoZSUlLU1NQUtrypqUmZmZmntf/ggw/00Ucfafr06aFlwWDwxMADBmjXrl362te+dlo/j8cjj8fjpjQAAJCgXN0ZSU1NVV5enmpqakLLgsGgampq5PP5Tms/fvx4vfXWW2poaAi9vve97+naa69VQ0MDb78AAAB3d0YkqaysTLNnz1Z+fr6mTp2q5cuXq7W1VcXFxZKkWbNmKSsrS5WVlUpLS9PEiRPD+g8fPlySTlsOAADOTa7DSFFRkQ4cOKDy8nIFAgFNnjxZGzduDD3U2tjYqORk/rArAACIjuswIkmlpaUqLS2NuK62trbbvqtWrYplSAAA0E9xCwMAAJgijAAAAFOEEQAAYIowAgAATBFGAACAKcIIAAAwRRgBAACmCCMAAMAUYQQAAJgijAAAAFOEEQAAYIowAgAATBFGAACAKcIIAAAwRRgBAACmCCMAAMAUYQQAAJgijAAAAFOEEQAAYIowAgAATBFGAACAKcIIAAAwRRgBAACmCCMAAMAUYQQAAJgijAAAAFOEEQAAYIowAgAATBFGAACAKcIIAAAwRRgBAACmCCMAAMAUYQQAAJgijAAAAFOEEQAAYIowAgAATBFGAACAKcIIAAAwRRgBAACmCCMAAMAUYQQAAJgijAAAAFOEEQAAYIowAgAATBFGAACAKcIIAAAwRRgBAACmCCMAAMAUYQQAAJgijAAAAFOEEQAAYIowAgAATBFGAACAKcIIAAAwRRgBAACmCCMAAMBUTGGkqqpKOTk5SktLU0FBgbZt29Zl2xUrVmjatGkaMWKERowYIb/f3217AABwbnEdRtasWaOysjJVVFRox44dys3NVWFhoZqbmyO2r62t1cyZM7V582bV1dUpOztb119/vT755JNeFw8AABKf6zCybNkyzZkzR8XFxZowYYKqq6s1ePBgrVy5MmL7P/3pT7r77rs1efJkjR8/Xk899ZSCwaBqamp6XTwAAEh8rsJIe3u76uvr5ff7v9xAcrL8fr/q6uqi2saRI0d0/PhxnXfeeV22aWtrU0tLS9gLAAD0T67CyMGDB9XR0SGv1xu23Ov1KhAIRLWNe++9V6NHjw4LNJ1VVlYqPT099MrOznZTJgAASCBn9NM0S5Ys0erVq7V27VqlpaV12W7hwoU6dOhQ6LV3794zWCUAADiTBrhpnJGRoZSUFDU1NYUtb2pqUmZmZrd9f/Ob32jJkiV66aWXNGnSpG7bejweeTweN6UBAIAE5erOSGpqqvLy8sIePj35MKrP5+uy3yOPPKKHHnpIGzduVH5+fuzVAgCAfsfVnRFJKisr0+zZs5Wfn6+pU6dq+fLlam1tVXFxsSRp1qxZysrKUmVlpSTpV7/6lcrLy/Xss88qJycn9GzJkCFDNGTIkD7cFQAAkIhch5GioiIdOHBA5eXlCgQCmjx5sjZu3Bh6qLWxsVHJyV/ecHniiSfU3t6uH/3oR2Hbqaio0IMPPti76gEAQMJzHUYkqbS0VKWlpRHX1dbWhv380UcfxTIEAAA4R/DdNAAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKZiCiNVVVXKyclRWlqaCgoKtG3btm7bP//88xo/frzS0tJ0+eWXa8OGDTEVCwAA+h/XYWTNmjUqKytTRUWFduzYodzcXBUWFqq5uTli+61bt2rmzJm67bbbtHPnTs2YMUMzZszQ22+/3eviAQBA4hvgtsOyZcs0Z84cFRcXS5Kqq6u1fv16rVy5UgsWLDit/aOPPqrvfOc7+tnPfiZJeuihh7Rp0yb9/ve/V3V1dS/Lj53jODp6vCPiuiPtkZf31NZNP7d62nak9Udj3I9YxztT+nLsaLcVy5hWcxTruJH6xfM86Kqvm+XdjR/PY9uX4nENcTu3sWwzlvMlkq6uU7Ec22jGP3W8o+0dOtL+RY99Yh3LzbE9uf5MX3utz/+TXIWR9vZ21dfXa+HChaFlycnJ8vv9qquri9inrq5OZWVlYcsKCwu1bt26Lsdpa2tTW1tb6OeWlhY3ZUbl6PEOTSh/sdfbyf/FS31QTXzGmfbI5rht/0ztd7zFcz8SbY4i1Xs2zU9/PUfjUWdX2+zLsfpqW11dp2LZvts+bq6RsXBTT7RtE+W8dsvV2zQHDx5UR0eHvF5v2HKv16tAIBCxTyAQcNVekiorK5Wenh56ZWdnuymzz+SPHaFBA1NCPw8amKL8sSNc94tVtOOdKm/sCOV16tNVPdFs/9S+sdRzpriZ897sR+f5jfUcOROimZO+qjfa+XczXv7YETpvcKqr37ne7E/n8frq9ziSnuqMZeyuthnpmhBvPdUfqVY3x7s3ffqa22tAtNfUM33tjef5Ho0kx3GcaBvv27dPWVlZ2rp1q3w+X2j5Pffcoy1btuj1118/rU9qaqqeeeYZzZw5M7Ts8ccf1+LFi9XU1BRxnEh3RrKzs3Xo0CENGzYs2nK71d3bNCcNGpiipKSkiP3SBqTo2BeR+0fqF886O48tKaxPd/X0tP3Ofd3WE0+nHgO3cx7rfnSe3+7OEWvRzkmkers7v3szVlfjdbdNt7+rvTm2p47Xl7/HkXRXZ6xjR9pmpGtCb0RzbkRTf+da3Rzv3vQ5MfaJ//bF4XV7DYj2mnqmr73xOt9bWlqUnp7e47/frt6mycjIUEpKymkhoqmpSZmZmRH7ZGZmumovSR6PRx6Px01priUlJWlwqutHZsL6xdK/N+O5EW0ft9uPtZ54ibWW3u5Hd33PtjnqSVf1xmsf4n3O9Xb+z9Txi8c43W2zL8fqi211VWss83K2/c65qSfatmfbPvY1V2/TpKamKi8vTzU1NaFlwWBQNTU1YXdKTuXz+cLaS9KmTZu6bA8AAM4trmNWWVmZZs+erfz8fE2dOlXLly9Xa2tr6NM1s2bNUlZWliorKyVJc+fO1dVXX62lS5fqxhtv1OrVq7V9+3Y9+eSTfbsnAAAgIbkOI0VFRTpw4IDKy8sVCAQ0efJkbdy4MfSQamNjo5KTv7zhcuWVV+rZZ5/VAw88oPvuu08XX3yx1q1bp4kTJ/bdXgAAgITl6gFWK9E+AAMAAM4e0f77zXfTAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAUwnxFYAn/0hsS0uLcSUAACBaJ//d7umPvSdEGDl8+LAkKTs727gSAADg1uHDh5Went7l+oT4bppgMKh9+/Zp6NChSkpK6rPttrS0KDs7W3v37uU7b/oQ8xofzGvfY07jg3mNj0ScV8dxdPjwYY0ePTrsS3Q7S4g7I8nJyRozZkzctj9s2LCEObCJhHmND+a17zGn8cG8xkeizWt3d0RO4gFWAABgijACAABMndNhxOPxqKKiQh6Px7qUfoV5jQ/mte8xp/HBvMZHf57XhHiAFQAA9F/n9J0RAABgjzACAABMEUYAAIApwggAADB1ToeRqqoq5eTkKC0tTQUFBdq2bZt1SQnllVde0fTp0zV69GglJSVp3bp1Yesdx1F5eblGjRqlQYMGye/36/3337cpNkFUVlbqG9/4hoYOHaqRI0dqxowZ2rVrV1ibY8eOqaSkRF/96lc1ZMgQ/fCHP1RTU5NRxYnhiSee0KRJk0J/LMrn8+mf//xnaD1z2ntLlixRUlKS5s2bF1rGvLr34IMPKikpKew1fvz40Pr+OqfnbBhZs2aNysrKVFFRoR07dig3N1eFhYVqbm62Li1htLa2Kjc3V1VVVRHXP/LII3rsscdUXV2t119/XV/5yldUWFioY8eOneFKE8eWLVtUUlKi1157TZs2bdLx48d1/fXXq7W1NdRm/vz5+vvf/67nn39eW7Zs0b59+/SDH/zAsOqz35gxY7RkyRLV19dr+/bt+ta3vqWbbrpJ//nPfyQxp731xhtv6A9/+IMmTZoUtpx5jc3Xv/517d+/P/R69dVXQ+v67Zw656ipU6c6JSUloZ87Ojqc0aNHO5WVlYZVJS5Jztq1a0M/B4NBJzMz0/n1r38dWvbZZ585Ho/H+fOf/2xQYWJqbm52JDlbtmxxHOfEHA4cONB5/vnnQ23eeecdR5JTV1dnVWZCGjFihPPUU08xp710+PBh5+KLL3Y2bdrkXH311c7cuXMdx+FcjVVFRYWTm5sbcV1/ntNz8s5Ie3u76uvr5ff7Q8uSk5Pl9/tVV1dnWFn/sWfPHgUCgbA5Tk9PV0FBAXPswqFDhyRJ5513niSpvr5ex48fD5vX8ePH64ILLmBeo9TR0aHVq1ertbVVPp+POe2lkpIS3XjjjWHzJ3Gu9sb777+v0aNH68ILL9Qtt9yixsZGSf17ThPii/L62sGDB9XR0SGv1xu23Ov16t133zWqqn8JBAKSFHGOT65D94LBoObNm6errrpKEydOlHRiXlNTUzV8+PCwtsxrz9566y35fD4dO3ZMQ4YM0dq1azVhwgQ1NDQwpzFavXq1duzYoTfeeOO0dZyrsSkoKNCqVat06aWXav/+/Vq8eLGmTZumt99+u1/P6TkZRoBEUFJSorfffjvs/WLE7tJLL1VDQ4MOHTqkF154QbNnz9aWLVusy0pYe/fu1dy5c7Vp0yalpaVZl9Nv3HDDDaH/nzRpkgoKCjR27Fg999xzGjRokGFl8XVOvk2TkZGhlJSU055AbmpqUmZmplFV/cvJeWSOY1NaWqp//OMf2rx5s8aMGRNanpmZqfb2dn322Wdh7ZnXnqWmpuqiiy5SXl6eKisrlZubq0cffZQ5jVF9fb2am5t1xRVXaMCAARowYIC2bNmixx57TAMGDJDX62Ve+8Dw4cN1ySWXaPfu3f36XD0nw0hqaqry8vJUU1MTWhYMBlVTUyOfz2dYWf8xbtw4ZWZmhs1xS0uLXn/9dea4G47jqLS0VGvXrtXLL7+scePGha3Py8vTwIEDw+Z1165damxsZF5dCgaDamtrY05jdN111+mtt95SQ0ND6JWfn69bbrkl9P/Ma+99/vnn+uCDDzRq1Kj+fa5aP0FrZfXq1Y7H43FWrVrl/Pe//3XuuOMOZ/jw4U4gELAuLWEcPnzY2blzp7Nz505HkrNs2TJn586dzscff+w4juMsWbLEGT58uPPXv/7VefPNN52bbrrJGTdunHP06FHjys9ed911l5Oenu7U1tY6+/fvD72OHDkSanPnnXc6F1xwgfPyyy8727dvd3w+n+Pz+QyrPvstWLDA2bJli7Nnzx7nzTffdBYsWOAkJSU5//rXvxzHYU77yqmfpnEc5jUWP/nJT5za2lpnz549zr///W/H7/c7GRkZTnNzs+M4/XdOz9kw4jiO87vf/c654IILnNTUVGfq1KnOa6+9Zl1SQtm8ebMj6bTX7NmzHcc58fHeRYsWOV6v1/F4PM51113n7Nq1y7bos1yk+ZTk/PGPfwy1OXr0qHP33Xc7I0aMcAYPHux8//vfd/bv329XdAL48Y9/7IwdO9ZJTU11zj//fOe6664LBRHHYU77Sucwwry6V1RU5IwaNcpJTU11srKynKKiImf37t2h9f11TpMcx3Fs7skAAACco8+MAACAswdhBAAAmCKMAAAAU4QRAABgijACAABMEUYAAIApwggAADBFGAEAAKYIIwAAwBRhBAAAmCKMAAAAU4QRAABg6v8DFQ6d7L/JS8AAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t,s=CtmcSimulation(A,0,50)\n",
    "plt.plot(t,s)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69fd1ec6",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}