scripts/analysis_1d.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.pyplot\n",
    "import matplotlib.style\n",
    "import numpy\n",
    "import sympy\n",
    "import sympy.printing\n",
    "import sympy.utilities\n",
    "import scipy.io as sio\n",
    "\n",
    "matplotlib.rcParams[\"figure.figsize\"] = [12, 6]\n",
    "# matplotlib.style.use(\"dark_background\")\n",
    "sympy.init_printing(use_latex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We introduce normalised variables to simplify the equations:\n",
    "- normalised pixel position `p = x / w`\n",
    "- normalised shear `s = delta * v / w`\n",
    "- count `c = delta * rho`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# c0 := s < 1\n",
    "\n",
    "# c, p in ]0, +oo[\n",
    "c, p = sympy.symbols(\"c p\", real=True, positive=True)\n",
    "\n",
    "# s in ]0, 1[\n",
    "s0 = sympy.symbols(\"s0\", real=True, positive=True)\n",
    "s1 = s0 + 1\n",
    "s  = 1 / s1\n",
    "\n",
    "f = sympy.Piecewise(\n",
    "    (0, p < 0),\n",
    "    (p * c / s, p < s),\n",
    "    (c, p < 1),\n",
    "    (c * (1 - (p - 1) / s), p < 1 + s),\n",
    "    (0, True),\n",
    ")\n",
    "\n",
    "mu_with_assumptions = sympy.integrate((f / (1 + s)).simplify(), (p, 0, 1 + s)).simplify()\n",
    "variance_with_assumptions = sympy.integrate(\n",
    "    (((f - mu_with_assumptions) ** 2) / (1 + s)).simplify(), (p, 0, 1 + s)\n",
    ").simplify()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# re-introduce variables to clear bindings\n",
    "s0, s1, s = sympy.symbols(\"s0 s1 s\", real=True, positive=True)\n",
    "\n",
    "mu_c0 = mu_with_assumptions.subs(s0, 1 / s - 1).simplify()\n",
    "variance_c0 = variance_with_assumptions.subs(s0, 1 / s - 1).simplify().factor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\frac{c}{s + 1}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAC4AAAAnCAYAAABwtnr/AAAACXBIWXMAAA7EAAAOxAGVKw4bAAACNElEQVRYCe2Z2VECQRCGwTIAJAMJQVPQDFAjUDOQR3jFDMQMlAzEDMAMJAOPDPD7xx1ql+XYnQPXqu2qdu6ez7Z7YMbmYrFoxJLBYDDE9kfK/oS+t1TbuXrovHLLQuA6DL+gPepjTaVU+wo9VdtXooADJcixhU4gvykfkrp30QwdKsAK7gY9oi7YKHIQweolNhXL0aDFHBQc2BY2pUESUICbJDS49fL7ug35xbrr+l36YiTnCJBzVKWR5C+ho7G6ySlSQAUpsZ5v0Xf/2xXmZ/BTJQzWbitBY3z3duFm1ODhfFnMUu3xYn4KN6v2eDhfFrPU7Pf78W4SxRicZtUfQE5u81hUJ6eH85yW1h53cpvHohgXCQ+c/FK+x+vW9IgeU7c3rEY08GTDjssFgjW6twp2jnZQtTMSDTyzS8lG4tkLLaOup47cXfXfJmfO4/yG6fc+xdSUvujPDfJuGcmAAzhj8bUFpdRTmt76bssY3cfcJTiQSoITC51s3qNUglROluCQfYoO8C+KJ/SZ+kR924Q5Cq2zNXN0ErQZ1wvtqszpN8m3OlC0nfl2iDEByMsWRG+AetwpLazTSeB0HKY3w45OFT0kZR5RjccZtGGit2zjZfruqA81hlYuXOxx+CxIdCnAmpenKkIL0oIrvjMhAbA8rrD5a8l9agrIhAqljrsusBQNnd2arAQy/wZR576FvRXXbdTm2yt9U9ozylEmOUPCYTxIcm5ish7fNO7Tv/Mo9TH+A+aFrTmt6/m7AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{c}{s + 1}$"
      ],
      "text/plain": [
       "  c  \n",
       "─────\n",
       "s + 1"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(sympy.printing.latex(mu_c0))\n",
    "mu_c0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- \\frac{c^{2} s \\left(s - 2\\right)}{3 \\left(s + 1\\right)^{2}}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHAAAAA3CAYAAADUvc6fAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAICElEQVR4Ae2c7XXUOBSGJzlbQJYOFjoAtoPQASwVJHSwHH4l/zjQAdAB0MFCBYR0AFSwbDrIvo9GErJGGn9JtieZe44jWSNd3W9dy1YOrq+vV1PD+fn5K815pOuuru+6nqvtSuUeekrgt579R3e3ynujEsWtVH5Q8VnXA+730E8Ch/26F+n9OMLyUvf3pUi8cQ89JTCHAiExVJYLnWFbTzZub/c5Qui9SNz37f1F1D7q1nr0scq3oxDNOFi0E60uVZrlJkXKwRxJTEiIiPume9bE12H7mLpwkSB9UPloDJ6aYy2NJHOAiz4kc5frpvVf3ZMj0J5U4uQeGBEHA59EXDHlWfwkRU/CuZZUF78Y2DuVnkbVkcVXlY90fQroPVE9m+TNtQauROSpCDtS+SwgdnTV4r1QmbTY0ROUQYCyUIwH0ftcN+QDeJwHtdOGkf/tG4PKLAoUMcei4Z5Kozzudbm1MCBvUBVBuNA0CMEEg/7SHF8T8+B5GLULqa4LmfoLdxOWkyvQKgoBf1H9MZfqCH20xzjGVY7GFQqpQh367opOQmkKGu3qhxd+V4nhN6DaGqjJUNK/wWyEARZo4jkEpkJF0H1djfDACOGxsdAHgzCGcP0IfvpV7Ynz18BCNc2f27RwUShlgPDFmtngr7gCRRzu/48uMqePKlcquX+q64Hqv9PWBdSXMHOi0ijM4oH53LpJ1tkwjHieAThjFFXuRRfehexeq46hxkC2vsF3cQVqEpT1UUQY5VkqIOiNrXcqNB5m2KEJva0t1DLmZ26CgThz6Eq3Ix/kBo8pgC8iVwOKKlCTQwRCZNH1oHafLvvG9opRhMb+p67vdfFc1wgfCRR31JayXtd1CM6V5kVwhP4+4CNH2yDhR24sMRseFow162Zwb6pFH+RFAMJmjSryAC08hBUskhKAySxuO/8TlVlF98W5nrbeX9HD44HPyHMzqR/rI8+JB2GfYh4oxFgp12U4wZC6cLnwSRg2ylAbjL7iN12pRZ6pkmGGHxingpDcFyfDq4BoIelqKM/SuVIZ80h0idtWh6Uo04QudLHYboB+h9iuQCLSeJbTeLNbk2AsxAmDMJqCoThTuEa3iQ886k+VcdjMyQkDdDL28xfzQIuRjWNCnN9AFoF4Jcrok8TgSY11U3jwwNwCr58M4P25FH0oTou6XCFeUAYGRdLiDBU5YXxEidTWInxtLA1FFaiJn0GQJcp54pDtMqySh3wVxupgjgfZMLPltxjIgHOGMhRnPEeJe5SHEjHKGDbCpO3wUOVJ3LloEhMjn+NeSsZw2BDOCWIOskbNKV4wYBKY+FVcuTVwFIVlBxOS2kJt2RnrY3uhKZKR5cZ5ILKUpRJKCec774XiAe/7rDK5thfLQhHcgoAEKGmxC6KxKymsl42ELhx4Iz0QBmWxJAnHKn1GHDK+C3XRzjtTNi+ykeTGKnAXFFSCxpsaQkvIZidw7BW4E2rKE3lwdnZmspx8l+QvGzvtitPTf6OfJO12Ne7XwB3Xd9GttB2XxezkK4qxCUFEJIMm82w99LNXoKS0BLDK633oZ5/ELEF7axri10gv1dx66GevwOUoEEoInQ6ubCVsc7/5ch9CvSjmrSiExm8aeOELXKyL9N+9B6blsoRW81ZFinWemKRpJxQoJvgOhn3BnQXRzwvqreHQMad+KI890NSbedfNlNWfAy3R4XcfhAayrba364ZA9SOtXvRRMUNohz/ihTcL2aNioFAfDJUPoEOZ8VMSqq6BIgKLg2BPjOpkWyiE93Vd3hTwPWb2dUqSq5kaLW/vNP0fqqdCH59EwE/y3Z7GHOs3/5Wavf+pkm99klA7hKKsUxHgU2TVnee1vjVXX6yR70yzr1OSXPVshD5dqe9TWjFpHN/8YJCEvae6iBhJUB+USmjcmEttRCZwfFG986Gfqh4oYrAciE5ZY6pNXRuAkrMf8jZ6znRjlWIihOoYnDfWDEk83/3QFa9veOaRLsKsB4vf38eVqgrU5HwG1zjMojbH4NY35upnFnyVVb0vFkjte/FzBU+6eNnsPxNUvSGnrnRUVWBMhIjEwggTnMBpW/9QtGcwxuXuhQd87hgbXk3Iza4ZbtzMJXzhta38tdE5iQIlULyJRMYphTDSBoTORjiJBwhv3+NnMYq57vn00Sd2Y4ioncQY2iRoQgbZKLsNEP9D9eMWwlE6X1MnQeP5PXX8rDU5SiKcthG+iEajYRIPDKmU4AmffOdI5pZLtxnCZ+aExBwY5QpH5+Nn6ku4TRkOwryj38kiY8D4Sj/GsK5jgKOhqgLFuCFSZZyIsL+HILncY0UvZoSTZIAwi8eR/Z3qvu34WdI7NY7Qzm5PnBmqedlwWJk8wuU3CWZIuMiGGeFD2DwrGYWp5MwcyiGzK2LZleVCdImNetCUtRVICETIcSh0Qt6WhcEgjKaA5IZw6EFzGO9RWUQwHnGdCvzHMhk0U9UQKoo2QpYEbMKVfmMrbRsTl+qT3HJSO97ZWJeEq8vxMw2rCkcdsSePinUc2+hWVYES6ltdhLXwoR3r4/TQNu+DyF05Kray/BEtWNMBzjKwznOiKPW8+1C/sS86Gqq/jRhDoZhnDa1+VEzzTJbEaC68NHlUbIisDocMmnAM69xGGK4wP9Eg5SkVpjL/MiuMSKPmWLQHwpks9tYcFRuiyaV7IDyRrBSz2CFCKjiG7LmRfI3FvXgPhEF5IYkPydBUYW6sXDfGi3Y2G3ikKvqY8z9AoAonbZ+2lwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle - \\frac{c^{2} s \\left(s - 2\\right)}{3 \\left(s + 1\\right)^{2}}$"
      ],
      "text/plain": [
       "  2           \n",
       "-c ⋅s⋅(s - 2) \n",
       "──────────────\n",
       "           2  \n",
       "  3⋅(s + 1)   "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(sympy.printing.latex(variance_c0))\n",
    "variance_c0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# c1 := s > 1\n",
    "\n",
    "# c, p in ]0, +oo[\n",
    "c, p = sympy.symbols(\"c p\", real=True, positive=True)\n",
    "\n",
    "# s in ]1, +oo[\n",
    "s0 = sympy.symbols(\"s0\", real=True, positive=True)\n",
    "s = s0 + 1\n",
    "\n",
    "f = sympy.Piecewise(\n",
    "    (0, p < 0),\n",
    "    (p * c / s, p < 1),\n",
    "    (c / s, p < s),\n",
    "    (c * (1 + s - p) / s, p < 1 + s),\n",
    "    (0, True),\n",
    ")\n",
    "\n",
    "mu_with_assumptions = sympy.integrate((f / (1 + s)).simplify(), (p, 0, 1 + s)).simplify()\n",
    "variance_with_assumptions = sympy.integrate(\n",
    "    (((f - mu_with_assumptions) ** 2) / (1 + s)).simplify(), (p, 0, 1 + s)\n",
    ").simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# re-introduce variables to clear bindings\n",
    "s0, s = sympy.symbols(\"s0 s\", real=True, positive=True)\n",
    "\n",
    "mu_c1 = mu_with_assumptions.subs(s0, s - 1).simplify()\n",
    "variance_c1 = variance_with_assumptions.subs(s0, s - 1).simplify().factor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\frac{c}{s + 1}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAC4AAAAnCAYAAABwtnr/AAAACXBIWXMAAA7EAAAOxAGVKw4bAAACNElEQVRYCe2Z2VECQRCGwTIAJAMJQVPQDFAjUDOQR3jFDMQMlAzEDMAMJAOPDPD7xx1ql+XYnQPXqu2qdu6ez7Z7YMbmYrFoxJLBYDDE9kfK/oS+t1TbuXrovHLLQuA6DL+gPepjTaVU+wo9VdtXooADJcixhU4gvykfkrp30QwdKsAK7gY9oi7YKHIQweolNhXL0aDFHBQc2BY2pUESUICbJDS49fL7ug35xbrr+l36YiTnCJBzVKWR5C+ho7G6ySlSQAUpsZ5v0Xf/2xXmZ/BTJQzWbitBY3z3duFm1ODhfFnMUu3xYn4KN6v2eDhfFrPU7Pf78W4SxRicZtUfQE5u81hUJ6eH85yW1h53cpvHohgXCQ+c/FK+x+vW9IgeU7c3rEY08GTDjssFgjW6twp2jnZQtTMSDTyzS8lG4tkLLaOup47cXfXfJmfO4/yG6fc+xdSUvujPDfJuGcmAAzhj8bUFpdRTmt76bssY3cfcJTiQSoITC51s3qNUglROluCQfYoO8C+KJ/SZ+kR924Q5Cq2zNXN0ErQZ1wvtqszpN8m3OlC0nfl2iDEByMsWRG+AetwpLazTSeB0HKY3w45OFT0kZR5RjccZtGGit2zjZfruqA81hlYuXOxx+CxIdCnAmpenKkIL0oIrvjMhAbA8rrD5a8l9agrIhAqljrsusBQNnd2arAQy/wZR576FvRXXbdTm2yt9U9ozylEmOUPCYTxIcm5ish7fNO7Tv/Mo9TH+A+aFrTmt6/m7AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{c}{s + 1}$"
      ],
      "text/plain": [
       "  c  \n",
       "─────\n",
       "s + 1"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(sympy.printing.latex(mu_c1))\n",
    "mu_c1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\frac{c^{2} \\cdot \\left(2 s - 1\\right)}{3 s^{2} \\left(s + 1\\right)^{2}}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGoAAAA3CAYAAADkMd++AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAIZUlEQVR4Ae2c63HUSBDH1y4CMM7gIAO4y8BkAEcENhkcxSf7mwsyAGdwkMFBBBhncBDBgTPw/X+zM1PzkjTalbTa9XbVeB7qmX5N97S0kg/u7u4WtXBxcfFWuEcqj1S+q7zW2K3qPYysgQe161sjvVeNgRaqP6r6ovKU/h7G1cBhj+WfJ7iX6j+RwfCuPYysgT6GgpXQKC7khWMjs3t/l+8T+h4nanpi+9fJ+GBd660nqj8MtuhMF5KMRKwb1eZoSdk86JNMhJO14L/qc2a9C8eHamtdkpaPqp+Fa9pxkhrAeTNJzc1yaHv/SgbOfWTJjFXtUaH4WghFfVY9ipEsLRKVFwldjHclun7c8vJN9TOVzyH+3NriD6+5UvlNbXd0hGyeqlNM0PqeUQsRONNiR6pfhRSGbFsa16rTncUGQRgPwnmtDkKzG2cH4g9dERng/aUKm60IwkEOHOCvFKGXobTAiRZ4rNoYib6KO6vStdfpo3wX3sJ1/lTnWzhg23gSCnGhsICymSHxdKvyQgWZ/qnggmz6TYpXbSgRwiAo76vazylqQzzd9RpaHbSuUbbq0rqMPdK1pl3ZNL46QxPPlGx41XfVOIWH7IwSAsb4z2MsXZGDmtiJIqIQYxcO0NdusgGKZ41oNd1cO6/OjJvIgxIIqXNPPJCfc9jrwRtKzLOTcU2yjk+qF6rpE1efqv2QsQmALC/aDG00xRc7D97fqY0hPKhPmDxVbQyjGnkw9mjnqye+XoOMOuIxDH0I8UnCGCNZOghOzFwLtCaeWAso/WctsvDeq8A3YdiD+qzDk5PQe8CJ8PyEeTWQP9KZ8SgJg7AIFhlF4z4NXlUOrQHBXyhMpSl0hcsfq8MG6QStB99kSdHusxONsXXtl/p/q5B5+VBicbJKOPBLmO8D3mv7TGrBNWdxeN2FPrIpBK5SULhAV5s1VSCcnR9dc9uua02TwqouGWmhcegSRvEgbinO1EdGxhqBebpYs6Ea1xjjwgMxxg6ihCFiUFqikT5+als/c/sUWeuRcPjbBK5rjIhATcbkwh5h3HiRxjDsW66BA+6MgagS8Xgopp0XcYBloOsoZUqAQRgtgvghw/tDdepJIZ8kI9F9mPDNUxTVkQKKRDY/yEZzdjHcuNDHQ09Cgn/4KYHwMoTlHJgS8Oxi6BFPCIARSB6cIeATw5I4uEdaeGV0vuoaHrXpRAJeawD5o/PUJROvENwK7zyLO/1019YQWReH7LNpc2AkjGXOp4RQ6CnwzU05KOxMFERIDDNark0Coos8bKYTS/CLxq7V5hmldw57jep3ldOgv1j56Xm4yNBtMc9m4SFrqPyhycxyPcnMpsKA0bl+OEtulyF302FqU6p5I8JZRJmlR6Eh7ShCICH53niVZMWbCIvZGT1Xj8JWJAPZzuLCDgNncJQEOVln61EwqJ1F4nCiunTgOhl2opaM3JRzQ16MILM21E5YYCAh5hz6BhJxN5bZG2pL7Hhwfn5e/07zlgi1i2zuz6gtsap71rcl7O4um8r2eHbJfRSZLplf9AHG3lDSyKbBGqn1A4x9MrFpKy3phz/RMHKpEn2AsTfUUlFz+EvIc3BrG35s1NDXFXcdV/e9lp6iJ+XSBz+OAvwUYmA0Q9XEXcfEvs40QGIRJRNjhr7OuJuxlwzI2LzfwDOwrQXxzw+YPoR1CSJcjMQzP/drtZkSeZRdMPxVFxckG1n1l1EYdA8Zw7jrxgwTpT+iSaoK7da3hkpz5zSG7lR4Va34OU3Iq3AaP8DwhhISSmUxbyi18QqI8LtQryfYwu+MuyGThTbv1hUf+RdwNzpk9XQlJnp/TuMY1xonavs3q2z/p2rzdlgY+jAK7775kKW286Qhfm3N4q5jMq1Fl511rbrT89K5ffvIq1J6B6N1Kc3hnRI2MXK9VCECFEE4RBPCWZGOxolcrPNV7eIHGN6jhITlWJCSQmksxWnsi3gx7jZOWL4tNOuQJ5nQifF4tdlYfoM3yMW90Q+V6OyxuESPIxV+OPRgaZi+N5QGeT0p+hBAY474yr+0WiGq32gSvjl4VY/uTV4jEzQkD2/ums9pVEevgqkf6b3EjjdUelGTsTCewFcSxfNJ41x3n+iwwwhXJqaqvVC7Ne6CUwA2RyRIAYe1W2mX5sxgDLnwwk75Ul4zQ0kB7GgSCqcwXDYD4bV+0qLrLu5equ08k3W7EgRCXhQCUuJdtFP8GfV5DQ4d9IaSoQg5JA9kgBx+P1TzaaPfBWpjzNInLWG46oy7WqMErM2brkWopF2cO4NB5CJS9YbMUOEKUgphj/fMyG7C1NMoUmONn7ToWmfcDWkF7WO1CaNN0Em7NFH8ECoJxSmguGNdJ3NLgTOlKwKkc9r6bGRzBrchla55Q4khs4Dq0CuYw/MmBKSYdF04HIyEKDyPjOdM/c5PWoS3NqxKW/PgNQONE5Z5AlLKxjL8TQ14Q4kB4udCDD9UadzRuubCHkYz4VBjhMihPmlpDA8T0EYFYwLRInWEKnqHARbGwStSIzlXdWcUBz1hxIPmmN2oeiUm/ELLBmsgUAnGpl2iOeQYukz1W7V+6FFZaJDiTVjQSjxCcgTY8VHc1jU8KptfxUGOdKOh7JVeizY27ZybupGjOjQjl9vwlVOWaN5QUvYHFd5KDW9u2QHpv64xqbvwWAHjweSQn7S0fXYzNm1kqgarK7z/xE7q/TlNLbFZvoUkBXBeTvLZjWhNkkyIDhs6+5ym1lCHtYgT43EGDhVKu1gnFBWfvHRN7Hmd25wwWvWaPkuPQgLtwJ357MZ6U/FzmlprzdWj4J+EZeUdWKuAifDIVqMErC/d2XoUgmgnksyQ4EwRmvrqrgpfvPNAgNuetW5d/geNPWXEuEXe0wAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\frac{c^{2} \\cdot \\left(2 s - 1\\right)}{3 s^{2} \\left(s + 1\\right)^{2}}$"
      ],
      "text/plain": [
       "  2          \n",
       " c ⋅(2⋅s - 1)\n",
       "─────────────\n",
       "   2        2\n",
       "3⋅s ⋅(s + 1) "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(sympy.printing.latex(variance_c1))\n",
    "variance_c1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAFlCAYAAAAZA3XlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9uElEQVR4nO3dd3hUZcKG8TuNkoQUCKE3AREURUDsa1dsoKvYFQtWbGvFrlhQVNS1F+xrd12xLboqKuqHNBWl995CQkhCQtr3xwwYI2iAJGeS3L/rmmtmzjmZeeJcIY9v3vOeqNLSUiRJkiT9JjroAJIkSVKksSRLkiRJ5ViSJUmSpHIsyZIkSVI5lmRJkiSpHEuyJEmSVE5s0AHKW7VqVemCBQuCjiFJkqRarnfv3quBppvbF3ElecGCBeyxxx5Bx5AkSVItV1pausWRWadbSJIkSeVYkiVJkqRyLMmSJElSOZZkSZIkqRxLsiRJklSOJVmSJEkqx5IsSZIklWNJliRJksqxJEuSJEnlWJIlSZKkcizJkiRJUjmxQQdQ9WuU1oTGrVpQr2FDigsLyVmTScbipRQXFgYdTZIkKSJYkuuA6JgYuh2wH7sediBd9tmTxMapfzimqLCQZbPmMP2b75nyvzEsmT4zgKSSJEmRwZJci0VFR9PnuKM55PyzadK6JTlrMpk+9v9Y9OtUVi9aQkFuHrH16tGoSSotOnek3W7dOWTQWRx24TksnDKVr199k59Gf05JcXHQ34okSVK1iiotLQ06w+9MmDChdI899gg6Ro3XrGMHTrnzZtp278b8n6bw5fOv8uuYsZSWlPzp1yWkptDjiEPY77QBpHdox/I58/hwxONM+/rbakouSZJUPUpLSycCvTe3z5JcC/U69khOuPlaNqxfz/vDH2Hyx59u9WtERUWx80F/4+grLya9QzumfvUtbw+9j+yVq6ogsSRJUvWzJNchh110Ln0Hn8/s8ZN49bpbWbc6Y7teLzo2hv1OHcCRl11IcWEh/x72IJM+HF1JaSVJkoLzZyXZJeBqkWOvvoy+g89n/Psf8/T5l293QQYoKSrm61fe4IETzmTZrDmcPux2Trz1emLi4iohsSRJUmSyJNcSB517BgeefRpjX3+HN2+5q9JPtstYtJgnz7uUz597mb0HHMelLz1FUnrTSn0PSZKkSGFJrgV6Hn04x/xjMJM+Gs1/ho2gqqbQlBQX8/EjT/LCFdeTvkM7LnvladI7tKuS95IkSQqSJbmGa965IwNuu4HZP0zkjZvvqrKCXNYvX3zNE2dfQmxcHJe98gwddt+1yt9TkiSpOlmSa7D6CfEMfPBu1q9bxyvX3UJxUVG1vfeS6TN59MwLyFmTyYXP/JPOe252zrskSVKNZEmuwfpdczlpbVvz6nW3kpORWe3vv2bJMh4beBGrFi7ivMcesChLkqRaw5JcQ3Xeaw/2OrE/Y156jbkTfwwsR25mFk8NusyiLEmSahVLcg1Ur2EDBtw2hJXzFjD6iZFBx/ldUT7nn/fRZueuQUeSJEnaLpbkGujg886iSeuWvHX7MIoKCoKOA4SK8jMXXEHOmkwGPfEgae3aBB1JkiRpm1mSa5iU5s04cOBpTPr4U+ZN+inoOL+zLmMNz1x4JaWlpVz49CMkNU0LOpIkSdI2sSTXMEdfeTEAHz30RMBJNm/1wsU8d8nVJKQmM+iJB6nXsGHQkSRJkraaJbkGad2tCz2PPoIxL71G1vIVQcfZosVTp/Py1TfRonNHTr37FqKiooKOJEmStFUsyTXIYRedS97abL584dWgo/yl6WP/jw8efIxdDzuIwy8ZFHQcSZKkrWJJriFa7bQjuxz0N756+XUKcvOCjlMhX7/yBj+89yGHX3QuPY44JOg4kiRJFWZJriEOvfAc1mevY+xrbwcdZau8c+dw5k36iVPuuoUWO3YKOo4kSVKFWJJrgGYdO7DroQfyzb/eIj8nN+g4W6W4sJAXr7qBvOxsBo64hwaJCUFHkiRJ+kuW5Bpg/9NPojC/oMaNIm+Uk5HJK9feQuNWLTjpjhuDjiNJkvSXLMkRrmFSEr2O6cvED/9LbtbaoONss3mTfuLjh59kt8MPZv8zTg46jiRJ0p+yJEe4vU7sR72GDfimho4ilzXmpdeY8vlXHHvVpbTbbZeg40iSJG2RJTmCRcfEsO8pJzDr/yawfNacoONUijduuYvM5cs5Y/hQGjRKDDqOJEnSZlmSI1iXffcitUVzxr7+TtBRKk3+uhz+df1tJKc35YSbrw06jiRJ0mZZkiNYn+OPYV3GGqZ+PTboKJVq4ZSpfPrkSHoedTg9jzki6DiSJEl/UNGS3BeYAcwGhmxm/9+ASUARcGK5fQOBWeHbwG2LWfckNk5l5wP2Y8KoTygpKg46TqX7/LmXmTvxR0646Voat24ZdBxJkqTfqUhJjgEeB44EugGnhu/LWgicDbxWbntj4DZgT6BP+HHqtsetO3oecwQxcbGM/8+HQUepEqUlJbx2wx2UlpZy+rDbiY6JCTqSJEnSJhUpyX0IjSDPBTYAbwD9yx0zH/gZKCm3/QjgM2ANkBl+3Hfb49Ydex5/LPN/msKKufODjlJlMpct5907h9O+R3cOOd8/MkiSpMhRkZLcClhU5vni8LaKqOjXXgBMACakpaVV8KVrr9bdutC80w6M/89HQUepcpM/+YxJH3/KYRecQ8sunYOOI0mSBETOiXvPAL2B3qtXrw46S+B69D2MosJCfvr0i6CjVIv37nmQ3LVrOeWum4mJjQ06jiRJUoVK8hKgTZnnrcPbKmJ7vrbO2u2Ig5n53Q+sz14XdJRqkbc2m3fvHE6rnXZ02oUkSYoIFSnJ44HOQAegHnAKMKqCrz8aOJzQyXqp4cejtz5m3dFut11o3LIFP47+POgo1eqXL75m4of/5dDzz3bahSRJClxFSnIRcCmhcjsNeAv4FRgK9Asfsweh+cYDgKfD+yF0wt6dhIr2+PDXrKmk7LXSbkccQtGGDfz65ddBR6l27w17iNysLE69+xanXUiSpEBVdE7yx8COQEfg7vC2W/ltRHk8oakUCUATYOcyX/s80Cl8e2E789ZqUVFR7Hb4wUwf+z35OblBx6l267OzeWfofbTs0plDLzg76DiSJKkOi5QT9wS02607Kc3S+XF03Thhb3N+HTOWCR98wiGDBtK80w5Bx5EkSXWUJTmC7HzQfhQXFjHt62+DjhKoUcMfIT8nhwG3DSEqKiroOJIkqQ6yJEeQnQ/cnzkTJtXJqRZl5Wat5f37/0n7Ht3Za8BxQceRJEl1kCU5QjRp05pmO7Tn1zFjg44SESZ+8Akz/288R195CUlNvcCMJEmqXpbkCLHzgfsBMPUrS/JG7wwdTmxcHMcN+UfQUSRJUh1jSY4QOx+4H8tmzWHNkmVBR4kYGYsW8+lTI9nt8IPZ+aD9g44jSZLqEEtyBGiY1IgOPXdzqsVmjHnpNZbOnM3fb7qG+gnxQceRJEl1hCU5AnTZuw8xsbFOtdiMkqJi3r59GElN0+h76QVBx5EkSXWEJTkC7LjPnuRlZ7NwytSgo0SkhVOm8v1b77HfqSfSYsdOQceRJEl1gCU5Auy49x7MHjeR0pKSoKNErE8efZq8tdmccNM1rp0sSZKqnCU5YE3btyW1RXNmfj8+6CgRbX32Oj586HE69NyNXsceGXQcSZJUy1mSA7bj3n0AmPn9DwEniXwT3v+Y+T9O4ZirBtOgUWLQcSRJUi1mSQ7YjnvvQcbiJWQsXhJ0lIhXWlrKv+9+gISUZI70JD5JklSFLMkBio6NodMevZxqsRWWTJ/Jd2/+m31O/jutdtox6DiSJKmWsiQHqF33nWmQmOBUi630yWPPkJu1lr97Ep8kSaoiluQAdezTi5KSEmaNmxh0lBolf10OH454nPY9utO7/1FBx5EkSbWQJTlAHXv1YPmsOazPzg46So0z8YNPmDfpJ475hyfxSZKkymdJDkh0bAztdtuFuRN/DDpKjVRaWsp7w0YQn5LMYReeE3QcSZJUy1iSA9Jqpy7Uj49njiV5my2ZPpMf/v0B+592Ek3btw06jiRJqkUsyQHp2KsHAPMsydvlk8eeprCggH7XXB50FEmSVItYkgOyQ68erJy3gHUZa4KOUqPlZGTy2VPP0+2Afemy715Bx5EkSbWEJTkAUVFRdOi5m/ORK8k3/3qLVQsW0f+6K4iOjQk6jiRJqgUsyQFo1mkH4pOTLMmVpLioiA8e+CfNdmjPPif9Peg4kiSpFrAkB2DjfGRLcuX5dcxYZnw3jiMGDyIhJTnoOJIkqYazJAegw+67krV8BZnLlgcdpVZ5f/gj1I+P54jB5wcdRZIk1XCW5AC03XVnFvz8a9Axap0Vc+bx/VvvsfeA42jeuWPQcSRJUg1mSa5miY1TadK6FQstyVVi9BPPkZ+TS//rrgg6iiRJqsEsydWsbfedAVjw8y8BJ6md8tZmM/qJZ9lxrz3o+rd9g44jSZJqKEtyNWu7azeKi4pYPG1G0FFqre/eeo+V8xZw7NWXEh3jknCSJGnrWZKrWbtdd2HpzNkU5hcEHaXWKikq5sOHHqfZDu3Z84R+QceRJEk1kCW5GkVFR9N2l27OR64Gv375DXMmTOaISwZRPyE+6DiSJKmGsSRXo/QO7WiQmODKFtVk1P3/pFGTxhx87plBR5EkSTWMJbkatQuftLdwiiW5OiyeOp2JH/6XA846lZRm6UHHkSRJNYgluRq13W1n8rKzWb1gUdBR6oxP/vk0RMGRl18UdBRJklSDWJKrUdudu7FoylRKS0uDjlJnZC5bzjevvknvfkfSquuOQceRJEk1hCW5msTExdG80w4smurSb9Xt8+deJmdNJv2uuTzoKJIkqYawJFeTFp13ICYuliWuj1zt8nNy+fTJkXTq04tuB+wXdBxJklQDWJKrSauuXQC8iEhAvn/nP6yct4ALh15OvQb1go4jSZIinCW5mrTu2oX12etYs3hp0FHqpJKiYuZ++A6X/a0hN9w1KOg4kiQpwlmSq0mrrl0cRQ7Y28+8w/c/zOaiM/YiMbFh0HEkSVIEsyRXg+jYGFp26cSSaTODjlLnXTH4CZo1S+W66/4edBRJkhTBLMnVIL1De+Lq12fJdEeSgzZhwixef/0rrrr6eFq2bBx0HEmSFKEsydWg9caT9lz+LSLcdOMrxMREM3To6UFHkSRJEcqSXA1add2Rgrz1rPJKexFh/vwVPP7Yh5x9zqHssku7oONIkqQIVNGS3BeYAcwGhmxmf33gzfD+cUD78PY44CVgCjANuGE7stZYrbt2YemMWZSWlAQdRWF33/0Wa9fmcd/wc4KOIkmSIlBFSnIM8DhwJNANODV8X9Z5QCbQCXgIuC+8fQChAt0d6AVcyG8Fuk6Iioqi5U6dWTLdk/YiSWZmDnff9SZHHtmLQw7ZLeg4kiQpwlSkJPchNEI8F9gAvAH0L3dMf0IjxgDvAIcAUUApkADEAg3DX5+93alrkNRWLWiQkMDSGbOCjqJyHnvsQ+bNW8Hw+88hKioq6DiSJCmCVKQktwLKTqZdHN62pWOKgLVAE0KFORdYBiwEHgDWbOY9LgAmABPS0tIqmr1GaNG5IwDLZs0JOInK27ChiJtufJndd+/IGWccGHQcSZIUQar6xL0+QDHQEugAXA3ssJnjngF6A71Xr15dxZGq18aSvGL2vICTaHPefPMbxo+fxZ13nUkDL1ctSZLCKlKSlwBtyjxvHd62pWNigWQgAzgN+C9QCKwEviVUhuuMFp07krF4CQV5eUFH0WaUlpZy7TXP07ZtU664ol/QcSRJUoSoSEkeD3QmNBJcDzgFGFXumFHAwPDjE4EvCM1HXggcHN6eAOwFTN++yDVL884dnWoR4b7++hdGjRrHkBtOJC0tKeg4kiQpAlSkJBcBlwKjCS3j9hbwKzAU2Dj0NpLQHOTZwFX8tkzc40Bi+PjxwAvAz5WUPeLFxMXRtF0bS3INMOT6F0lIaMAtt5wSdBRJkhQBYit43MfhW1m3lnmcT2i5t/JytrC9Tmi2Q3tiYmNZPtOSHOmmT1/MyOc+5aKLj+TRRz9g9uxlQUeSJEkB8op7VciVLWqW229/jYKCQu4ZNvCvD5YkSbWaJbkKtejckaING1i10MtR1wQrVmRx//B/c+KJ+7LPPl2DjiNJkgJkSa5CzXfsyIq58ykpKg46iirowQffY+nSDIbf7+WqJUmqyyzJVaiFK1vUOHl5Bdx262vss09XTjhhn6DjSJKkgFiSq0jDpEakNEtnuSW5xnnhhf/xyy8LGHbvQOLiKnpuqyRJqk0syVWk2Q4dAFjulfZqnJKSEq6/7gU6dWrJRRcdGXQcSZIUAEtyFWm2QzsAVsybH2wQbZNPPpnI55//xC23nkJyckLQcSRJUjWzJFeR9B3aU5hfQObS5UFH0Ta67toXSEtLYsiQE4OOIkmSqpkluYqkd2jHyvkLKC0pCTqKttHkyXN4+eUvuOLKfrRp0zToOJIkqRpZkqtIsx3as3LegqBjaDvdcvOrANx51xkBJ5EkSdXJklwFYuvXJ7VlC1bOnR90FG2nRYtW8cjDozjjjAPp0WOHoONIkqRqYkmuAk3btSE6OtqR5Fpi2LC3WbMmxwuMSJJUh1iSq0CzDq5sUZtkZ+dx59A3OPTQHvTt2yvoOJIkqRpYkqtA+g7tKSkpYdWCxUFHUSV56qlPmDVrKfcNP5voaH9sJEmq7fxtXwWa7dCeNUuWUlRQEHQUVZLCwiJuvOElundvz9lnHxJ0HEmSVMUsyVUgvUM75yPXQu+++x3ffTeNoXeeTnx8/aDjSJKkKmRJrmRR0dE0bd+WlXMtybXRtdc8T8uWTbjqquOCjiJJkqqQJbmSpbZsTlz9+qxw+bda6fvvp/Puu99x3fUn0KxZStBxJElSFbEkV7JmHdoDON2iFrthyIvUrx/H7befFnQUSZJURSzJlSw9vPzbSpd/q7Vmz17G0099wnmDDqdr1zZBx5EkSVXAklzJ0tq1ITdrLXlrs4OOoio0dOgb5OSsZ8RDg4KOIkmSqoAluZI1bduG1QsWBR1DVWz16mzuuP11jjiiJ8ce2yfoOJIkqZJZkitZk7atWLXQklwXPP74R0ydupARDw2ifv24oONIkqRKZEmuRLH16pHSvBkZC73SXl1QVFTMlVc8S8eOLfjHP/oHHUeSJFUiS3IlatK6JdHR0ayyJNcZ//vfj7z33vfceNNJtGzZOOg4kiSpkliSK1Fau9BKB6styXXKNVePJDY2hnvvOzvoKJIkqZJYkitR07YbS7JzkuuSefNW8MD9/+aMMw5in326Bh1HkiRVAktyJWrStjW5WWtZn70u6CiqZvfe+w6LFq3ikX9eQHS0P1aSJNV0/javRC7/Vnfl5RVw3bUv0KtXJ84999Cg40iSpO1kSa5ETdq2YvUi5yPXVW+++Q1ff/0Ld99zFikpCUHHkSRJ28GSXEk2Lv/mSHLddsXlz9C4cSK3335a0FEkSdJ2sCRXko3LvzmSXLf99NM8nn7qvwy+9Gh69Ngh6DiSJGkbWZIrycbl31YtsCTXdTfd9AqrV2fz5FOXeBKfJEk1lL/BK0lam9aAy78J1q7N5Zqrn2fPPbtw/vmHBx1HkiRtA0tyJUlr18bl37TJv/41hi+++Il7hg0kPT0l6DiSJGkrWZIrSVrb1p60p98ZfMlTJCTUZ/j95wQdRZIkbSVLciVp0roVGUuWBh1DEWTGjMXcP/zfnHXWwRxwwC5Bx5EkSVvBklwJomNiSGmezprFlmT93t13v8Xcuct5/ImLiYuLDTqOJEmqIEtyJUhpnk5MbCxrHElWOfn5G7js0qfo1q0tV199XNBxJElSBVmSK0HjVi0BWLNkWcBJFIk++WQi7777HTffcgrt2zcLOo4kSaoAS3Il2FiSnZOsLfnHlc9SUlLC409cHHQUSZJUAZbkStC4dQuKi4rIWr4i6CiKUIsXr+amG1/hyCN7cfrpBwYdR5Ik/QVLciVo0qolWctXUlJUHHQURbDHH/+I77+fzkMPn09aWlLQcSRJ0p+wJFeCxq1aetKe/lJJSQnnD3qUpKSGPPzIBUHHkSRJf6KiJbkvMAOYDQzZzP76wJvh/eOA9mX27Qp8D/wKTAEabGPWiNW4VQtP2lOFTJ26kHvufovTTjuAo47qHXQcSZK0BRUpyTHA48CRQDfg1PB9WecBmUAn4CHgvvD2WOBV4CJgZ+BAoHB7Q0eS2Pr1SWqa5kiyKmzYsHf45ZcFPPnUJTRq1DDoOJIkaTMqUpL7EBohngtsAN4A+pc7pj/wUvjxO8AhQBRwOPAz8FN4XwZQqybuNm7ZHMCSrAorLCzi/EGP0qpVE4YNGxh0HEmStBkVKcmtgEVlni8Ob9vSMUXAWqAJsCNQCowGJgHXbU/YSNS4dXiN5MVOt1DFjRs3g0f/+QGXDD6affct/4cZSZIUtKo+cS8W2A84PXx/PKFR5vIuACYAE9LS0qo4UuVq4hrJ2kY33/wq8+ev4NnnLqNBg3pBx5EkSWVUpCQvAdqUed46vG1Lx8QCyYSmViwGvgZWA3nAx0DPzbzHM0BvoPfq1asrmj0iNG7VksL8Atatzgg6imqY3Nx8zh/0KDvt1Jq77joj6DiSJKmMipTk8UBnoANQDzgFGFXumFHAxsmVJwJf8Ns0i+5APKHyfAAwdbtTR5DGrVqwZqlTLbRtPv/8J554/COu/Ed/9tvPaReSJEWKipTkIuBSQoV3GvAWoeXchgL9wseMJDQHeTZwFb8tE5cJjCBUtH8kNC/5o8qJHhlcI1nb6/rrX2TevBW88OKVJCTUuhUSJUmqkSo6J/ljQifhdQTuDm+7ld9GlPOBAYSWgOtDaCWMjV4ltPzbLtTKE/dcI1nbJzc3n3POfpgOHZoxfPg5QceRJEl4xb3t0qBRIvFJSaxZ7Eiyts/YsVN5aMT7XHzJURx6aI+g40iSVOdZkrdDaotmAM5JVqW45ZZXmTZtESOfv5ykpPig40iSVKdZkrdDSvPQhUQyl60IOIlqg/z8DZw98CFatGjMQw+fH3QcSZLqNEvydtg4kpy13JKsyjF+/CzuHfY255xzKP377xV0HEmS6ixL8nZIbdGMosJCcjLWBB1Ftcidd77JhAmzePa5y2jRonHQcSRJqpMsydshpUVzspavoLS0NOgoqkUKC4s44/QHadiwHi+9/A+ioqKCjiRJUp1jSd4Oqc2bkeV8ZFWBmTOXcOUVz3LooT24+urjgo4jSVKdY0neDsnN0z1pT1Vm5MhPeffd77jr7jPZffeOQceRJKlOsSRvo+iYGJLTm3rSnqrUBec/ysqVa3nt9WuIj68fdBxJkuoMS/I2SkprQkxsrCVZVSozM4ezzhxB584teeihQUHHkSSpzrAkb6OUFq6RrOoxZswU7h/+b86/oC9///s+QceRJKlOsCRvo01rJC9bHnAS1QW33vovxo+fxXMjL6NDh2ZBx5EkqdazJG+jlE0XElkZcBLVBYWFRZx80n0AvPnW9dSrFxtwIkmSajdL8jZKbdGcvLXZFOTlBR1FdcT8+Ss45+xH6N27Mw88cG7QcSRJqtUsydsopXkzMp1qoWr2/vv/x4gH3+PSy45lwID9go4jSVKtZUneRinN051qoUAMGfIS338/nWefu4xOnVoEHUeSpFrJkryNUsOXpJaqW1FRMaecPJzCwiLeensIDRrUCzqSJEm1jiV5G9SPjyc+OcnpFgrMokWrOOvMEfTosQOPPHJ+0HEkSap1LMnbIKV5OgBZrpGsAH3yyUSG3fMW51/Ql3PPPSzoOJIk1SqW5G3ghUQUKW655V+MHj2Jx5+4mD59dgw6jiRJtYYleRtsupCIc5IVsJKSEk479X6WLMng3X/fSLNmKUFHkiSpVrAkb4PkZumUlJSQvXp10FEkMjNzOP64u0lJSeDtd24gLs4LjUiStL0sydsgOb0pORlrKCkqDjqKBMCUKfM595xH2G+/bjz8sCfySZK0vSzJ2yA5vSlrV64KOob0O2+/PZbh973DxZcc5Yl8kiRtJ0vyNkhuZklWZLrxxlc2nci3115dgo4jSVKNZUneBsnpTVm7wpKsyLPxRL5Fi1bxn/dvpl279KAjSZJUI1mSt1Js/frEJyc5kqyIlZmZwzFHDyUuLpYPPryVpKT4oCNJklTjWJK3UnJ6UwCyLcmKYDNnLuHEE4bRpUsr3njzOmJi/FGXJGlr+JtzKyU3C5VkR5IV6b788mcuufhJ+vbt5YoXkiRtJRdU3UobR5Kdk6yaYOTIT+nSpRXXXPt3ZsxYwmOPfRh0JEmSagRL8lbaVJIdSVYNMWTIS3TesSUPPTyIOXOW8cknE4OOJElSxHO6xVZKTm9KQV4e+Tm5QUeRKqSkpIQzTn+Qn36az5tvXU+vXp2CjiRJUsSzJG+l5GYu/6aaJzc3n2OOvoNVq7L56OPb6NixRdCRJEmKaJbkreTV9lRTLV+eSd8jbiUqKor/jr6D9PSUoCNJkhSxLMlbKSk9zZKsGmvWrKUce8xQmjdP5aOPbyMxsWHQkSRJikiW5K0QFRVFctOmrpGsGu2HH2Zy0oB72W23Drzz7g3ExXn+riRJ5VmSt0JC4xRi4mIdSVaN98knEzl/0KMcfvjuPP/CFURFRQUdSZKkiOIQ0lZwjWTVJi+99DktWqRyz7CBrM3K5dJLnwo6kiRJEcOSvBWS09MB10hW7XHvve+QkpLAddefSG5uPtdf/2LQkSRJigiW5K3ghURUGw0Z8hKJiQ259roTWLduPXfd9WbQkSRJCpwleSskN2tKSXEx61avCTqKVKkuu+xpEhIbMPTOM8jJyefhh98POpIkSYGyJG+F5PSmrMtYQ0lxcdBRpEpVWlrKoPP+SXx8fUY8NIjc3HyefXZ00LEkSQqMJXkrJLtGsmqx4uLQ5avj4+vz5FOXkJ+/gVde+TLoWJIkBcIl4LZCo6ZpZK9aHXQMqcoUFhYx4MR7+eKLn3nhxSsZOPCQoCNJkhQIS/JWaNSksfORVevl52+g37F38r///cTI5y/nnHMODTqSJEnVrqIluS8wA5gNDNnM/vrAm+H944D25fa3BXKAa7YpZQSIjo2hUZPGjiSrTsjP38Bx/e/i008nM/L5Kxg06PCgI0mSVK0qUpJjgMeBI4FuwKnh+7LOAzKBTsBDwH3l9o8APtmupAFLbNwYwJFk1Rn5+Rs4/ri7+eij8Tzz7GVceGHfoCNJklRtKlKS+xAaIZ4LbADeAPqXO6Y/8FL48TvAIcDG69weB8wDft3OrIFKSmsCQPZqR5JVdxQUFHLC3+/hww/H8+RTg7nkkqOCjiRJUrWoSEluBSwq83xxeNuWjikC1gJNgETgeuCOv3iPC4AJwIS0tLQKRKp+jcIled3qjICTSNVrw4YiTjzhHt5///947PGLGTLkxKAjSZJU5ar6xL3bCU2/yPmL454BegO9V0foSG1S0/BI8ipLsuqeDRtCq178619juGfYQO677+ygI0mSVKUqsk7yEqBNmeetw9s2d8zi8GsmAxnAnsCJwHAgBSgB8oHHtid0EDaNJGc4J1l1U1FRMWedOYK1Wblce90JpKQkcPHFT1JSUhJ0NEmSKl1FSvJ4oDPQgVAZPgU4rdwxo4CBwPeESvEXQCmwf5ljbic0olzjCjKE5iTnZq2luLAw6ChSYEpLS7n00qfIzMzhpptPJik5gbPOHEFhYVHQ0SRJqlQVKclFwKXAaEIrXTxP6CS8oYTmEY8CRgKvEDrBbw2hIl2rNEpr4vJvUtgtt7xKZmYODzx4HklJ8Qw4cRh5eQVBx5IkqdJU9LLUH4dvZd1a5nE+MOAvXuP2Cr5XREpq2sSpFlIZI0b8h6ysXJ5+ZjBfjhnGsccMZeXKrKBjSZJUKbziXgU1auJIslTe889/xvHH3U23bm347vv72XHH8gvfSJJUM1mSKyipaRMvJCJtxocfjuegA28kMbEB3343nL333inoSJIkbTdLcgU0aJRIXP36XkhE2oIJE2axz97XkpGxjv99fhfHH7930JEkSdouluQKSNp0IRFHkqUtmTt3Ofvucx2TJ8/l7XeGcNVVxwUdSZKkbWZJroCNayQ7J1n6cxkZ2Rx6yM28++53PPDgeYwceTn16lX0/GBJkiKHv70qIKlp6FLZXpJa+mv5+Rs45eThTP11Ibfdfhqdd2zFCX+/h1Wr1gYdTZKkCnMkuQIapTUGINuSLFVIaWkpd9zxOiefdB89e3bkh/Ej2HXX9kHHkiSpwizJFZCUlkZhfgH563KCjiLVKG+/PZa/7X89MTHRjP12OP377xV0JEmSKsSSXAGN0ho7iixto0mT5tBnj6v45ZcFvPefm7jttlOJiooKOpYkSX/KklwBSU3TnI8sbYflyzM56MAbefHFz7nt9tP48KPbaNy4UdCxJEnaIktyBTRq4kiytL0KCgo595yHufCCxzj44F2ZMPEhevXqFHQsSZI2y5JcAY4kS5Xn2WdHs/9+1xMVFcXYb4czaNDhQUeSJOkPLMl/ISY2loSUZEeSpUo0YcIsevf6B2PGTOGZZy9j5MjLadiwftCxJEnaxJL8FxKbpAKQk+HV9qTKlJGRzdFH3cHQO17nnHMP44fxI9hll3ZBx5IkCbAk/6XExqGSvM6SLFW6kpISbr/9NQ4/7BaaNGnEuB8e5KKLjgw6liRJluS/ktg4dCGRnDWZASeRaq///e9Heux2OV999QtPPHkJ77x7A6mpiUHHkiTVYZbkv9CoSbgkZ1iSpaq0cmUWRx91B9dcPZJjj+3D5B//yb77dgs6liSpjrIk/4WN0y0cSZaqXmlpKSNG/Id997mODRsKGfPVPdx995nUqxcbdDRJUh1jSf4LiY1T2bA+n4K8vKCjSHXGhAmz6NXzSl568XNuuPEkxv0wgu7d2wcdS5JUh1iS/0Ji41RHkaUArFu3nkGDHqXfsUNp1iyF8RNGMGTIicTE+M+WJKnq+dvmLzRq0tiSLAXoww/H032XSxk16gfuGTaQr7+5j06dWgQdS5JUy1mS/4IjyVLwMjKyOWnAvZxx+gPstFNrfvzpUa655nhHlSVJVcbfMH8hsUmqayRLEeK1176i+y6D+eyzyQy//1z+b9yD7LZbh6BjSZJqIUvyX3AkWYosS5eu4fjj7uakAffSunUTxk94iHvuOYsGDeoFHU2SVItYkv9Ew6RGxMbFWZKlCPTOO9/SreslvPzS5wy5YQA//vRP/va3XYKOJUmqJSzJf+K3NZKdbiFFoszMHAYNepRDD7mJmJhoxnw1jOdfuJL09JSgo0mSajhL8p/YWJLXebU9KaJ98cXP7Nr9Mu4d9jannfY3ps94ksGDj/bEPknSNvM3yJ/YdElqR5KliLd+fQE33vgyu3a/jB9+mMWjj13ED+NHsPfeOwUdTZJUA1mS/8Sm6RaOJEs1xsyZS+h7xK0MOHEYaWlJfPvd/Yx8/gqnYEiStool+U8kNk6lpKSE3Ky1QUeRtJXeffc7unW9hHuHvc3ppx/AzFlPc/31J1K/flzQ0SRJNYAl+U8kNk4lL2stJcXFQUeRtA1yc/M3TcH48sufGXbvQKbPeIpTTvlb0NEkSRHOkvwnvCS1VDvMnLmE44+7m4MPupGMjHW89vq1fP9/D7DPPl2DjiZJilCW5D/hhUSk2mXMmCns0fsfnD3wIVq3TmPst8N5863r6dixRdDRJEkRxpL8JyzJUu1TWlrKyy9/QZcdL+S2W//FUUf1Ztr0J3nqqcG0atUk6HiSpAhhSf4TiU1SWZfh8m9SbZSXV8Cdd75Bp47n8+QTH3P2OYcwa/YzPPDAuaSlJQUdT5IUMEvyFsTExRGflORIslTLrViRxRVXPMOOnS/k9de+4oor+zFn7rPcccfpJCXFBx1PkhQQS/IWJDZOAbAkS3XEwoWrOO+8f7LLzpfyySeTuOXWU5g77zluvvlkkpMTgo4nSapmluQt2HQhEUuyVKfMmLGYU06+j149r+Tbb6cx9M4zmL9gJEOHnk6TJk7DkKS6wpK8BYmNw5ek9mp7Up00efIc+ve7k567X8Fnn/3Izbecwrz5zzF8+Dk0a5YSdDxJUhWzJG/BxpHkdY4kS3Xajz/O5aQB97LLzoN5//1x/OOq/syd9xwPP3w+rVunBR1PklRFLMlbkJiaAkBupiVZEkydupAzz3iQrjtdzBuvf83FlxzFnLnP8vIrV9Gjxw5Bx5MkVTJL8hYkpKZQXFhEfk5u0FEkRZDZs5dx3nn/pHOnC3js0Q/p339PJk1+hM/+dxd9+/YKOp4kqZJYkrcgISWZ3KysoGNIilALF67i6qtH0rbNuVx/3Qt06dKKjz+5nZ+nPMbZZx9CvXqxQUeUJG0HS/IWJKSmkJOZFXQMSRFu7dpc7r//33Tc4XzOOnMExcUlPP/Clcxf8DxDh57uVfwkqYaqaEnuC8wAZgNDNrO/PvBmeP84oH14+2HARGBK+P7g7charRJTU8i1JEuqoMLCIl599Ut273E5hx92C+PHz+LGm05i3vyRvPX2EA44YJegI0qStkJFSnIM8DhwJNANODV8X9Z5QCbQCXgIuC+8fTVwLNAdGAi8sv2Rq0dCagq5WWuDjiGpBvrf/36kf7876dTxAh4a8R8OOqg7X44Zxs9THuOii44kMbFh0BElSX+hIiW5D6ER4rnABuANoH+5Y/oDL4UfvwMcAkQBk4Gl4e2/Ag0JjTpHvISUZEeSJW2X+fNXcP31L9Km9Tmce87DFBQU8sSTl7B4yYs89thFroohSRGsIiW5FbCozPPF4W1bOqYIWAuUn4h3AjAJKNj6mNUrKjqahslJlmRJlSI/fwMvvvg5e/T+B3vvdQ3vvz+Oc887jEmTH2H8hIe4+OKjvPS1JEWY6jpxb2dCUzAu3ML+C4AJwIS0tOAX549PTiI6OtrVLSRVunHjZjDwrBG0bHEWlw5+kpiYaB5/4mKWLnuJl16+ir/9zbnLkhQJKlKSlwBtyjxvHd62pWNigWQgo8zx7wFnAXO28B7PAL2B3qtXr65ApKqVkJIMQG6mc5IlVY2srFyeeOJjeu5+Bb17XcmLL3xOv359GPPVMGbMfJqbbjqJ9u2bBR1TkuqsipTk8UBnoANQDzgFGFXumFGETswDOBH4AigFUoCPCK2I8e32x60eCRuvtudIsqRqMGnSHAYPfpKWLQYy8KwRLF26hjvvOpO5857jq6/v5YIL+pKS4nQMSapOFSnJRcClwGhgGvAWoZPwhgL9wseMJDQHeTZwFb8tE3cpoRUvbgV+DN/SKyV5FUpISQEgZ01WoDkk1S3r1xfwyitfctCBN9C+3bnceMNLNGnSiKeeHsyy5a/wzrs3cPzxe3uhEkmqBlGlpaVBZ/idCRMmlO6xxx6BZtjzhH6cdPsNDD20H2tXrAo0iyTtvntHzjjjQE497QCaN08lMzOHf7/7HW+9NZYvv/yZoqLioCNKUo1UWlo6kdCU3z9wOGIzElNTAcjNyg44iSTB5MlzmDx5Dtdd9wIHH7wbZ5x5ECedvB/nDTqcjIxs/vPe//H222P54gsLsyRVFkvyZiSkJlOQl0dRQcSvViepDikuLuGzzybz2WeTqV8/jiOO6MmJA/ZlwEm/Feb3//N/vP32t3z++U8WZknaDpbkzUhISXE+sqSIVlBQyKhR4xg1ahz168dx+OG7M+Ck/ThxwH6ce97hrFmzjg8++IEPRv3A6NGTyM3NDzqyJNUoluTNSEhNdmULSTVGQUFhqBB/8MOmwnzCifty7LF9GDjwEAoKCvn885/4YNQ4PvjgB5YuXRN0ZEmKeJbkzUhITSE3yzWSJdU8ZQtzTEw0++7bjX79+tCv/548+dRgnnxqMOPHz+KD8Cj0zz/PDzqyJEUkV7fYjBs/eZd5k3/i9RuHBppDkipT165t6NdvT47t14e99upCdHQ0CxeuYvR/J/Lf/07i889/Ijs7L+iYklRtXN1iK4WmWziSLKl2mTZtEdOmLeK++96hWbMUjjmmD32P7MVJJ+/P+Rf0paiomO++m8bo/07iv/+dxI8/ziXSBlIkqbo4klxObL163DfxKz5+5Ck+f+6lwHJIUnWJjY1hr7260LdvL47o25NevToBsGJFJqNHT+bT0ZP44oufWb48M+CkklS5HEneCvEpyYCXpJZUdxQVFTN27FTGjp3KzTe/Qnp6CocfvjtH9O3JUUf15qyzDgZg6tSFfPnFz3zxxc+MGTOFzMycgJNLUtWxJJeTmJoCQG5mVqA5JCkoK1dm8eqrX/Lqq18SHR1Njx4dOOigXTn4kN0YePYhDL70GEpKSpg8ee6m0vzNN7+6zJykWsWSXE7CxpLsnGRJoqSkhEmT5jBp0hwefPA94uJi2WOPzhx88K4cdPCuXHb5sVxz7d8pLCxiwoTZjP3mV8aOncq3305jzZp1QceXpG1mSS4nITzdImeNc+8kqbzCwiK++24a3303jbvuepMGDeqxzz47ccghu7Hf/jtz+RX9uPa6EwD49deFfDt2Kt+Ei/OCBSsDTi9JFWdJLseRZEmquPz8DXwRnnIBUL9+HL17d2b//bux737dOOnk/bjgwr4ALFq0KjT3+ZupfP/9dKZMmU9xcUmQ8SVpiyzJ5SSmplBSUsL6bP9MKElbq6CgkG+/ncq3304FIDo6mp13brupNO+//86ceuoBAOTlFTBx4mx+GDeTceNmMG7cTBYtWhVkfEnaxJJcTkJqCuuz11FSXBx0FEmq8UpKSpgyZT5TpszniSc+BqBdu3T22qsLe+7ZhT577sjgS4/m6muOB2DZsjWMGzeTH8Klefz4WeTkrA/yW5BUR1mSy0lISXZlC0mqQgsWrGTBgpW8+eY3AMTFxbLrru3Zc88d6bNnF/bcc0eOO24vIFSyp09fzKRJc5gcPoFw8uS5XhlQUpWzJJcTn5JM3trsoGNIUp1RWFjExImzmThx9qbR5tTURPbYozN77tmFnr06csABu3DGGQdt+prZs5cyadLcTcV50qQ5ZGT4b7ekymNJLic+KYm1K50TJ0lByszM4dNPJ/Ppp5M3bWvaNJndd9+Bnj07snvPjvTq1ZGTTtpv0/6FC1cxadIcpvw8n59/nseUKQuYPXsZJSWeHChp61mSy4lPSWLZrDlBx5AklbNq1do/FOeUlAR69PitOPfs2ZFjj92DmJgYANavL2Dq1EVMmbKAX6bM5+efQ/OjV6zICui7kFRTWJLLiU9OIm+ty79JUk2QlZXLmDFTGDNmyqZtDRrUo2vXNnTv3o7u3dvTfdf2HHHE7px99iGbjlm1ai1TpsznlykLmDJlAdOmLWLatEVealvSJpbkMmJiY2mQkOCcZEmqwfLzNzB58hwmT/79XwWbNEmie/d27Lpre7p3b88u3dtx3qDDSUhosOmYFSsymT59CdPDpXnatMVMm7aIxYtXV/e3ISlgluQy4pOTACzJklQLZWRk/2HUOSoqivbt0+natQ1du7Zhp51as1PX1gw4aT8aN2606bicnPVMn76YadMWM33aIqZPX8ysWUuZPXsZ+fkbgvh2JFUxS3IZm0qyV9uTpDqhtLSUefNWMG/eCj7+eMLv9qWnp9C1a2t22ql1qEB3bcOBB+7CmWce9LvjFi1axaxZS5kzexmzZi1l1qxlzJ69lDlzllugpRrMklzGppKc7UiyJNV1K1dmsXJlFl999cvvticmNqRLl1Z06tSCzp1b0jF8f9zxe9O0afLvjl24cBWzZy9l9qxlm0ae585dzrx5K7xIihThLMllxKeE/nFzuoUkaUtyctZvWte5vOTkhE3luVOnFnQK3x//9z8W6IyM7E2j2PM33s9fGb5fQUFBYXV9S5I2w5JcxsaR5FynW0iStsHatbl/WaA7dGi26da+QzN23bU9xx7bhwYN6v3u+KVLMzYV540leuHCVSxatJpFi1aRl1dQXd+WVCdZkstISHYkWZJUNf6sQEdFRdG8eervC3T7dNp3aMY++3Tl5JP3JzY25ndfk5GRHS7Mq1kcLs6/lejVLFmSQWFhUXV9e1KtY0kuo2FyEsWFRRTk5gUdRZJUh5SWlrJs2RqWLVvDd99N+8P+mJho2rRpSps2abRpk0bbtk1p06YprcPP99236+9W4wAoKSlhxYosFi1azcKFq1gcLs5Ll65h6dKN92vIzc2vrm9TqlEsyWXEJyd50p4kKeIUF5cwf35orvKWxMfX30yRTqN1m6bsvHNb+vbtSWJiwz98XXZ23qbivGxZJsvC5blskV62LJP1653eobrFklxGQkqyUy0kSTVSXl4BM2YsZsaMxVs8JikpnhYtGtOy5e9vLVo2oWXLxuy99060bNn4D/OjATIzc8KFeQ0rVmSxauVaVqzIYsWKzPB9FitXrmXlyiw2bHCah2o+S3IZ8UlJrpEsSaq1srPzyM7O+9MiDZCSkkDLcHEue2seLth77dWF9PRkGjWK3+zXZ2bmlCnOWawMPy67beNjT0BUpLIklxGfnETW8i3/KUuSpLogKyuXrKxcpk5d+KfHxcfXJz09hWbNQrf09OTfHjdLJT09mV12aUezQ3b7w5zpjdavL2D16uxNt4yMdWSUeV7+lpGxzou0qFpYksuIT0li6cxZQceQJKlGyMsr+Mu50hvFxcXStGkSzZql/q5Up6Ul0aRJI5qkJZGWlkTbtk3D25K2+Fq5ufl/LM/h+zVr1pGZmUtmZk74cc6mW3FxSWV++6rlLMllxCcnuUayJElVoLCwaNOJgBURExNNamoiaeHyXP7WuMnGx43o2LE5aWlJpKQk/ulrZmfnlSvPuWSWKdJr1uT8rlRvPG7t2jxKS0sr4z+DahBLclhMXBz14+M9cU+SpAhQXFyyaZS4ouLiYklJSaBx40akpiaQmpoYfpz4263xb4932qk1qamh4zd3suJGJSUlrF2bx9q1uX+4X5e9frPb/3DcuvUW7RrGkhy28Wp7lmRJkmqmwsIiVq1ay6pVW/9X4QYN6m0qz43LFOmNRTslJYGk5HiSkxNITo6nVasmdOvWZtPzuLi/rlTZ2Zsv0NllivS6devJycnf9Dj0fP3vnjsnu3pYksMsyZIk1V35+Rs2XdBlWzRsWJ/kMiV6431SUvzvn5fZn56eTOfOLTc9r18/rkLvVVRU/IfiXLZY52ymWJd9npOTT27uxlsBubn5FBUVb9P3XZtZksMSUrwktSRJ2jbr1xewfn0By5dnbvNrxMXFkpjYgEaNGm66JSY23OLzxE3bQl/TtGkSjRrFb3r+Z1NIytuwoZDc3ALy8gr+UKBzc/PJyysgL/eP+/7s+LLH1MSpJpbksE0jyZ64J0mSAlBYWLTppMHKEBsbs4Wi3YCEhI23+sTH1//98/Dj+PjQ6HjLlo037du4PTo6equy5OUV/OG2fv1vj2+68ZW/XHKwulmSw+KTnG4hSZJqj6Ki4kot3WU1aFDvd8V5Y3ku+7xsqd64rUHDesTHh4p5w4ah12jaNJm4uJhKz7i9LMlh8U63kCRJqpD8/A3k528gIyPoJFVn68bKa7H45CSKCgspyMsLOookSZICZkkOi09OchRZkiRJgCV5E0uyJEmSNrIkh8UnJ7myhSRJkoCKl+S+wAxgNjBkM/vrA2+G948D2pfZd0N4+wzgiG0NWtXik5PIy3YkWZIkSRUryTHA48CRQDfg1PB9WecBmUAn4CHgvvD2bsApwM6EivYT4deLOAkpyU63kCRJElCxktyH0EjwXGAD8AbQv9wx/YGXwo/fAQ4BosLb3wAKgHnh1+mz3amrQMOkJPKyLMmSJEmqWEluBSwq83xxeNuWjikC1gJNKvi1ABcAE4AJaWlpFYhUuaKio5kzYRLLZs2p9veWJElS5ImUi4k8E76xevXqar+4d2lJCSMHX1PdbytJkqQIVZGR5CVAmzLPW4e3bemYWCAZyKjg10qSJEkRpSIleTzQGegA1CN0It6ocseMAgaGH58IfAGUhrefQmj1iw7h1/lhu1NLkiRJVagi0y2KgEuB0YRWpnge+BUYSmge8ShgJPAKoRPz1hAqxoSPewuYGn6dwUBx5cWXJEmSKl9UaWm1TwH+UxMmTCjdY489go4hSZKkWq60tHQi0Htz+7ziniRJklSOJVmSJEkqx5IsSZIklWNJliRJksqxJEuSJEnlWJIlSZKkcizJkiRJUjmWZEmSJKkcS7IkSZJUTsRdcQ9YBSwI6L3TgNUBvbeqj59z7ednXDf4OdcNfs51Q1Cfczug6eZ2RGJJDtIEtnBpQtUqfs61n59x3eDnXDf4OdcNEfc5O91CkiRJKseSLEmSJJVjSf69Z4IOoGrh51z7+RnXDX7OdYOfc90QcZ+zc5IlSZKkchxJliRJksqxJIf0BWYAs4EhAWdR1XkeWAn8EnQQVZk2wJfAVOBX4Ipg46iKNAB+AH4i9DnfEWwcVaEYYDLwYdBBVKXmA1OAHwmtchERnG4R+gGcCRwGLAbGA6cS+iWr2uVvQA7wMrBLwFlUNVqEb5OARsBE4Dj8ea5tooAEQj/PccBYQv9D9H9BhlKVuIrQsmBJwDEBZ1HVmU/oc46o9bAdSYY+hEaQ5wIbgDeA/oEmUlX5GlgTdAhVqWWECjLAOmAa0Cq4OKoipYQKMoRKclx4m2qX1sDRwHNBB1HdZEkO/QJdVOb5YvylKtUG7YHdgXEB51DViCH0p9mVwGf4OddGDwPXASUB51DVKwU+JfTXvwsCzrKJJVlSbZQIvAtcCWQHG0VVpBjoQWi0sQ9OoaptjiH0P0ATgw6iarEf0BM4EhhMaHpk4CzJsITQyT4btQ5vk1QzxREqyP8C/h1wFlW9LEIna/YNOIcq175AP0JzVd8ADgZeDTKQqtTG3rUSeI/Q//gGzpIcOlGvM9ABqAecAowKNJGkbRUFjCQ0F3lEwFlUdZoCKeHHDQmdeD09sDSqCjcQGrRqT+j38hfAGUEGUpVJIHSi9cbHhxMhq1BZkqEIuBQYTegX61uElhRS7fM68D3QhdDc8/OCjaMqsC9wJqFRpx/Dt6MCzKOq0YLQ6PHPhAY6PsMlwqSaqhmhFWp+IrS040fAfwNNFOYScJIkSVI5jiRLkiRJ5ViSJUmSpHIsyZIkSVI5lmRJkiSpHEuyJEmSVI4lWZIkSSrHkixJkiSVY0mWJEmSyvl/lKXfml4W4C0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "dark"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "variance_c0_lambda = sympy.utilities.lambdify(s, variance_c0.subs(c, 1.0))\n",
    "variance_c1_lambda = sympy.utilities.lambdify(s, variance_c1.subs(c, 1.0))\n",
    "\n",
    "s_grid_c0 = numpy.linspace(0.0, 1.0, 1000, endpoint=False)\n",
    "s_grid_c1 = numpy.linspace(1.0, 5.0, 1000, endpoint=False)\n",
    "\n",
    "figure, subplot = matplotlib.pyplot.subplots()\n",
    "subplot.plot(\n",
    "    s_grid_c0,\n",
    "    numpy.vectorize(variance_c0_lambda)(s_grid_c0),\n",
    ")\n",
    "subplot.plot(\n",
    "    s_grid_c1,\n",
    "    numpy.vectorize(variance_c1_lambda)(s_grid_c1),\n",
    ")\n",
    "sio.savemat('cnd1.mat',{'cnd1':numpy.asarray(numpy.vectorize(variance_c0_lambda)(s_grid_c0))})\n",
    "sio.savemat('cnd2.mat',{'cnd2':numpy.asarray(numpy.vectorize(variance_c1_lambda)(s_grid_c1))})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}