examples/how_to_compute_energies.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f9499194-335f-46c4-9845-6e9f4b29acc0",
   "metadata": {},
   "source": [
    "# How To: compute the energy terms, along with their derivatives.\n",
    "\n",
    "This notebook shows how to compute the required energy and gradient functions. That will be used for the \n",
    "\n",
    "MeanFieldVarGP algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d0405a8d-ff1a-47cb-973f-b18e0fe2cc08",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import sys\n",
    "\n",
    "import dill\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import sympy as sym\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "from sympy.stats import E, Normal\n",
    "\n",
    "from sympy.utilities.iterables import flatten\n",
    "\n",
    "from copy import deepcopy\n",
    "\n",
    "module_path = os.path.abspath(os.path.join('../src/'))\n",
    "\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)\n",
    "# _end_if_\n",
    "\n",
    "from numerical.symbolics import LagrangePolynomial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8887ae9f-482d-47b3-a284-285f37df512e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Help make the output easier to read.\n",
    "sym.init_printing(use_unicode=True)\n",
    "\n",
    "# Objects referred to in the global dictionary are recursively traced and serialized.\n",
    "dill.settings['recurse'] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fffe1703-a2e6-46fa-920f-9237b612bde2",
   "metadata": {},
   "source": [
    "## Dynamical System equations\n",
    "\n",
    "To generate the equations we use the following settings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e756644-c93d-4f05-922d-0ad9c9805cba",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# |------------------|\n",
    "# | Double Well (1D) |\n",
    "# |------------------|\n",
    "# \n",
    "# D = 1\n",
    "# \n",
    "# p = sym.Matrix([sym.symbols(\"theta\", real=True)])\n",
    "# \n",
    "# fx = sym.Matrix((4 * xt) * (p - xt**2), real=True)\n",
    "# \n",
    "\n",
    "# |-------------------------|\n",
    "# | Ornstein-Uhlenbeck (1D) |\n",
    "# |-------------------------|\n",
    "# \n",
    "# D = 1\n",
    "# \n",
    "# p = sym.Matrix([sym.symbols(\"theta mu\", real=True)])\n",
    "# \n",
    "# fx = sym.Matrix([p[0] * (p[1] - xt[0])], real=True)\n",
    "# \n",
    "\n",
    "# |---------------|\n",
    "# | Lorenz63 (3D) |\n",
    "# |---------------|\n",
    "# \n",
    "# D = 3\n",
    "# \n",
    "# p = sym.Matrix([sym.symbols(\"sigma rho beta\", real=True)])\n",
    "#\n",
    "# fx = sym.Matrix([p[0]*(xt[1] - xt[0]),\n",
    "#                 xt[0]*(p[1] - xt[2]) - xt[1],\n",
    "#                 xt[0]*xt[1] - p[2]*xt[2]], real=True);\n",
    "# \n",
    "\n",
    "# |---------------|\n",
    "# | Lorenz96 (nD) |\n",
    "# |---------------|\n",
    "# \n",
    "# D = 4\n",
    "# \n",
    "# p = sym.Matrix([sym.symbols(\"theta\", real=True)])\n",
    "# \n",
    "# fx = sym.Matrix([(xt[1] - xt[3]) * xt[2] - xt[0] + p[0]], real=True)\n",
    "# "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78283361-9825-4f20-967b-4f293326783e",
   "metadata": {},
   "source": [
    "## First define the main variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "37230d60-f326-4a28-a17d-fff9f96b9336",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "⎡   σ⋅(-x0(t) + x1(t))    ⎤\n",
      "⎢                         ⎥\n",
      "⎢(ρ - x2(t))⋅x0(t) - x1(t)⎥\n",
      "⎢                         ⎥\n",
      "⎣ -β⋅x2(t) + x0(t)⋅x1(t)  ⎦\n"
     ]
    }
   ],
   "source": [
    "# State vector dimensions.\n",
    "D = 3\n",
    "\n",
    "# Mean and variance functions: m(t), s(t).\n",
    "mt = sym.symbols(f\"m:{D+1}(t)\", real=True)\n",
    "st = sym.symbols(f\"s:{D+1}(t)\", real=True, positive=True)\n",
    "\n",
    "# Random (state) vector: x(t).\n",
    "xt = sym.Matrix([Normal(f\"x{i}(t)\", mt[i], sym.sqrt(st[i])) for i in range(D)])\n",
    "\n",
    "# System noise (diffusion) coefficient vector: sig := sigma^2.\n",
    "sig = sym.symbols(f\"Sig:{D}\", real=True, positive=True)\n",
    "\n",
    "# Model (drift) parameters vector.\n",
    "p = p = sym.Matrix([sym.symbols(\"sigma rho beta\", real=True)])\n",
    "\n",
    "# Drift function.\n",
    "fx = sym.Matrix([p[0]*(xt[1] - xt[0]),\n",
    "                xt[0]*(p[1] - xt[2]) - xt[1],\n",
    "                xt[0]*xt[1] - p[2]*xt[2]], real=True)\n",
    "\n",
    "# Number of systems dimensions.\n",
    "# NOTE: this can be different from\n",
    "#       the state vector dimensions!\n",
    "dim_D = len(fx)\n",
    "\n",
    "# Pretty printing.\n",
    "sym.pprint(fx)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a8299b9-c2c7-44d8-ad0d-44c671bef574",
   "metadata": {},
   "source": [
    "## Compute the mean/variance Lagrange functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "80c0661e-df68-4350-a4a8-2450b5bfd1ec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done!\n"
     ]
    }
   ],
   "source": [
    "# Store the mean and variance functions.\n",
    "Mfun = []\n",
    "Sfun = []\n",
    "\n",
    "# Store the time derivative of the mean and variance functions.\n",
    "dMt = []\n",
    "dSt = []\n",
    "\n",
    "# Store the mean and variance points.\n",
    "Mpar = []\n",
    "Spar = []\n",
    "\n",
    "# Create one polynomial for each state vector dimension.\n",
    "for i in range(D):\n",
    "    \n",
    "    # Create the mean/variance polynomials.\n",
    "    mt_i, t, mk_i, t_h = LagrangePolynomial(letter=f\"d{i}m\", order=3, fp=\"h\")\n",
    "    st_i, t, sk_i, t_c = LagrangePolynomial(letter=f\"d{i}s\", order=2, fp=\"c\")\n",
    "    \n",
    "    # Store the functions.\n",
    "    Mfun.append(mt_i)\n",
    "    Sfun.append(st_i)\n",
    "    \n",
    "    # Store the Lagrange coefficients.\n",
    "    Mpar.append(mk_i)\n",
    "    Spar.append(sk_i)\n",
    "    \n",
    "    # Get the time derivative functions.\n",
    "    dMt.append(mt_i.diff(t))\n",
    "    dSt.append(st_i.diff(t))\n",
    "    \n",
    "# _end_for_\n",
    "\n",
    "print(\"Done!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04ed8de8-ad92-4bf6-b7c2-552f07f41484",
   "metadata": {},
   "source": [
    "## Compute the expectations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3f31b645-2516-4286-b737-b05bbf52f1b0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAABLCAYAAADODmoeAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAUPUlEQVR4Ae2dXa7dthHHjy8u+li4MdD3OjtwnBX0ZgdJswInO3Dgp+TNSHaQZAVusoMkK3DiHTQLKFD3oo99KPr/8XJkipf6PKS+DgfQoUQNhzN/kaMhpSM++PLLL5+cTqfftKXox6+++uqT1ImaVxGoCFQESiEgv/MPyX6ckq9zD66DE99pH+aQfg8P6n5FoCJQEVgIga9Vz8Oorg91/DF5oeP6Vp7sTcRYDycgIPy4Q9wo5Sawe5pij3g/k8E/K734m11O3C4VV9l9rw95LJzjutp979qIAb6x4vzvAb4FFaXXQ23o97W2z3zKNEGSdB4nPNoe8WM3sil3sZQbt4pruimFEVeao+aOReAHMW55PvAX6fdSHeFHDFJKGP6b0o+0paKkTnusrMrHZZ8pj3o+0HYoks3c6b/X9hft3/YYVwK3w+Lag2PvqRpx9cIz7qQaMsOkX5WmHMA4IQW5vH6PlTqnRVXap/NxTEdrkefvs+dGBYis3oYFvUyGi8/D/Dn7kvFxDjlz6rYyqp8o9Qdt3yrvU23xnIuxulR8Q+1gFm6Sy7XKgmtL4R0fVMeV5+J9ITFMJm6ViAR/Tij3WnlP1DHiDjlkz0cq98Z3qFjsS2W8iDP3eIx92j7R9rn0B6shKonbYXAdAnHM+eq4xqDUw6NGTeRxUrrJaMurzp0+RdzJoeb8SHvgTzlCcEDm70obmVRwdCqN26Xi2tVurrtOXFK+GgURBxHT37TF0QdQfCce7ropYu4j2YlTzFPyVCedn+gG58jdHHJPVZTyaPi1eL4J+P6lPPiZJHc6KTV7WsM68YTknK/PSNrj6wAD5MFPpMYw0+mgNCTq7oryQr4i+17Xs3CbodgSuK2K6wxMihW5eMelRs6TNTqgNQpzYkRQzlmJpy+aooNQvgQx+f2Ftp8knHmWn7T/DRUpRc9/K+WQYZtzbErRFX3+xAnRe3fJ6danYWLO7FGQmbRHcsGHeRY66I1S+LqI9wG7HH1XmZz5OXCbqs8SuK2N61RMivFftONS56Pz8xTsmfbDiWs6O68C9DksuyhEH+YALO/sVHXjUF97QdTBfItzWj7PHNIj5eNUjMjHrphSeTEPx0P20EHfpAoGeeAxtr6g2Pm7BXAbq9QSuK2G61gQluK7aMclkBkeMh/TOC0PvHNYyuep0u3AxcBRJHkor3M4xrGdGDk4URwDehFJUZZOEUc4ODbo5V3S/MIXOjIaexeZ82OIadRpj2e4URrjZWUtBT90HiTZxzVAZkzY/Z7Ofxqf0DHYdL16kgu3RLW9WUvgNhrXXk0PcPLSHRePr1NDGpwCHQBHMpt8+VnvNAV1u06t49AZoRMOKqUj/I0zQ442ZTVDRvaNcA7QmMjyJDnmRBm6ZiHJtLm7ljzlMyTlFY4wymzxpA7Eb9fsLNxSsufmSafsuM3V5Sjlro5iyFQ71JgsIvg1UZa7fOwoEmwuq3T43jU0o2O2Ih/ZRB6dxOX7DoOS2EJ+TBZxhbb22fMUAZLb8Ad1hLKRO8oZhoUy7+fAbYpKS+C2BVynYFKM92IdlxC1IVSrg6kjEm3h1JKRQOJKUN4cQOL02Vk86WwcBdK8s0DHOPJh+MTwkkiM8y7qUMpkvXM6SkMiGozfx+qzp+UMVIeLikKBfp+6LfJJnF4kKwduUxRdArct4DoFk2K8V8Ukb1ywOh0dC4dgnfvkOzudnJcOx3a8N+KfNRxUuV7y+hApxQ6KTom+LYemLJyT5X2u8xZ58T/CtzrG0TjSPnKR8+wup/nts4cyztH78gzl4I8JPEyP+FzxY+lEBz8bt4mKLoHbqrhOxKMo+3VR6dsXToQS/jH4kY4/UMO/naC6vaowochoVjog0VPsBOiUrWGil8jcFq8CMHf3yudZQqN/4Ts1eR9q+6uOY8fTZw8T6TxtRf5Jadf8Ew40dogUWYqy4SYbeQ2FiNpucL8oD+fNqynhH+qXwG1tXJe6foP1PAg+JEiHjRvxoIDK4Dow79fgMFw0sndMZMdse1QWp0Gnfv8cHFSe6HDy5Pw5dZ5btiRuuXA918Y1ywsDbpjcOB9cranIgeomEhk7J7YHs8+xh6e0RCnnElFmKqo8V27J8iVxy4VrSfsXk10dVwaodQdgyEB0QLSxe5prj8oxhL1R2jWEHI2NZPAax64iWOk7qx2oXC9uQ+dHg3ogxuq48l3Mo90R59jjHmzkg3SXkkrgVnGNmkJ1XBEgcw99dMCTPDdxPVfOVspNtcfbjf27ipJy450bt4pr+grVyfk0LjW3IlAR2BgC3onXyfmNXZeqTkWgIjASgTpUHAlUZasIVAS2g0B1XNu5FlWTikBFYCQC1XGNBKqyVQQqAttB4CIclyb1eMdq1NM++ODfziW6r8nR7LlvYZmcilsZXE3qkvju7qmiwOFlPd5QvtXGX1P4a8kr5Sf/rgSYOs+TCL5sMIrEy3szfDK56KN9r5u9cW/OknqTtqC8L7NJe0aBm4lJOBymHWSCZFUxS7RL1dH85WePf7Lmi6K7X9jUOp7S5kue2sch2yKt8R+rrWF2vozoZf4mxvh/k/zhGdz4o/VR6BDtYOsXQ22K/4x+r22NhXA74dnVUFEg4nEZ9jX/YdM+kRfHdOgWef5VFzZtKdQ+oDG0vqAgfYm+sOeeLRTdgj3SYQsLtR6pHXBpN0W6xnyyfPWFcPtA2ZXjkiFEJ6lI5LXy97awKZ9JITqKCftoOHzQMCYcG1FZFzEcjj8MaLwvtXOIhVplx5HagV2fzaRqe/xPdCsL4SZx2Zvjsm8ixcYQpUDNeQHv5oyU9s1TwZ9yhCeVQybfwmpkUkFGQi+ix4djZIpv6/aMMSMXT9c12WM7yIXJKnLWapfXS1jrOz/RAJ3PJqMZO0N80G43C5vmsAWjJadrvskirdjhgtc9J+v14Y+9OEDwPfRCrbLXHP1b2dpFzsn7k4vg5q8DkSAPjGjTzMN2PmTxumVNvA5n9bMZCi2Cb6zXIo5LlTJZfJSFTXPYEl8Hd+wbHp2O1aktejBeGuS9uS/x4cwuaaFW+75/jA84mTPjS7ZGxXHTNWD4zvCeG8hJKREhDw9sUV6yl6BibbNH+eL4puq+SmXmzNNFJIJgDgqiU7JWXvi9JmuIm1/YtIAtDpTghw/w/ah6LCoNTjnsrGOG+bZPAxq6w1PeIhYrt0haALuxdtDmiuEmu3hQwAML57Q8mETLOLKuIa1ny5cUwHesckXx7VJiVMQlUGgk3EHGNpZb8T5TOToS80RMGFMWI+lgIdnQiMnjkOALh0Z9jc+cX+mFTXPZEtrp9oUPTovIKewAIR82gmsX0Umap60dTG5ereNcK1t6EEWkOh7Xcc2FWrfQDkKswCm+0dDOIbBKku8Pc/tUSmaxtpmqLMhbtF1avTiuP2r7p7b/WWacCmQ6TNecTMzeOvZlyXOdQMehMyIfBwXocaeEv3FmnNemrORSYNZA4nkh+O+R5MBP42KBg9EU6HiWLXGFkvucPKVdTisu0jpWuVn2tIREB5IZd0bHoXzmNHioEEbNUen7h+K363sWdsjRRgV2sworW6QdWIXSA1uo8++W59MnPuXGnSSPx6w+lRLo5XHqLHxTsufmSafc7fIP0gVfdbrS9h9tf/b7SopR11AGoFuRQtAgXL4HAMVwetY4OTayRhw6Re7MKV7KPOVHchv+oA5OGSE35Qxz2OLqUL04gveVNk5L+zgGu2ubLiXtsTqWSHNgt5V2wGR86qb7qc9PtZ3SGOfAd4qOS7bL/0oxfJVzXFOUPIeXNfwaR4Egdc4uj0yDYHhJo6ADu7uIUianndNRGhJ3rvj9JRqNObSQl/3WxVUdLoqImXRM3beJ/By2nFTvE8lmQrVxWr4u9ImppD1xXSWPc2C3lXZAu4zbNNeUjTa8BuXAd4req7RLIq7ipI6JA8BJxUMzQKYDty6+snBOlsfngC3yYjGCtzpuOrb2kYucZ9pCIkzvCsUp4+6GvjwRDvwxUd70cOfEl8sW5NABGfqwtiMb/0EkL/UJ5CL2OKMW+pFtubDbQjugDWEPmyPZR55dv1R78pxlklz4TtRulXZ5PVHJuexcXKKnlhNQHhe6NUz0FTC3RSTCE5tXPs8SnMkLf5HI452ZJRc2zWULDRxZz7XF5JxqlInT71r2iwlinB54nZR2zT9xQ4gdPEWWolzYoe/a7cBGAdxk7Bryh3/eOF/cafkLmA1f2UBbY8Ridq65EK43712yu69DvFN9eE/g8zIgDjDlCHoFqAyN4OyFTXsrmXhyC/ZIB6LdyZPzE03Nyl4CN8mkYz9V2hXVZ7Vhy8JK4JuyV/VwY+YG/eAqxXCgvNSj6rHmMe/UFeGMlZGbbwv2EDWnouTctuaUVwI3pide5VRyx7JK4NsLx6EdlzxzkQU6exEteHIL9kgH5uQmR7AFYRkUnRs3ySMaZ5ojnvoY1OWIDLnxHYPRoR2XB2BO5MT801pPhYau29HsGbI31/mcuOG4mLNday4rFyY55eTEd1CvwzsuNS6iAyZQ3cT1ECKeL/VUb6joIuePZs8ioKmSnLhJFv9wYCK+kkcgJ75jQD305PwYACpPRaAisA8E5BwvZnJ+H1ekalkRqAhMQuDwQ8VJaFTmikBFYBcIVMe1i8tUlawIVARCBKrjCtGo+xWBisAuEFjqLz+7AMOU1CTgY+3bXx14grSZ95ZMN6W8o9ZL4mEyc1P69ypc6OTRMJtiD5AesR3Up4pRZ9FF5n9nvKPDl0jtT7O9i7RGIkYfSj4O0r57xT7UWZfn5y8PfN1iFImXd9KQuRnnO0rxDibZwTXhTe1bbfyla7cLAkv3yST7aSeT2gCVqNzu24FsaJ4q1ograDoeGCIU92KhUvflBrHQUUY7i0Bk565kuw6otHnRVfvU07cgbOeLsV7eJSwGy1dDX8pe+2IIOBpmKee8C8xkD/8BZa3NEguvSqz7cz3YHeK/lVdYVMndkVzEowYUvw1Nx7BoKCdUNNLWlxpUN9EXkQSdrUU6x91mtcVt6Vja7CsILd2WOvAY8Afv5r+S2gcvjjeH2RAu0v2htuILr6KHx4mb8qrXcAiTseer43qHFB0zNW/E1yxjZ/au1Pw95tCmLAiLUyMi6yIiwvhjisb7Ujsv7GDHKdFp6v+Br5W/twWBnTNRm1tq4VUu+1HawaJfQAW4XZEaFc6FzeahcurPsIbogYiul8Rj0WBqKGRl0TPVqU8qT1TCvB08e6Yu/bEPas4fDbNz7QEcyThKOzhdY9Clky4oDZ4wGgeBk2LSl0gLp9Kac/C8RDfGq92TfZGVjxoOLm5LAcnpmmvgs79Q6KSQf88peV34cyt6os8mF4PNgZlkYCP09i5J/joH788UwSyHLUnNhzNz2EMttKOuyHVYi41wVMd1dyHo8CzCSue41fZI26/a+OYSnSEcKvJhwnMXt5XI++Q7BfXFC8LiKFNzODRCHC6N+kYpfF2EM8bJrUE5MONrnBDXJyZzZlw3o1KY5bDFdJyS5rCH+tZsB1Ps7eWtjusOHtfg1fHpFM2wUMdcZCbRXXSkY6Ih5lMgHAxPHcPPJFvnGru4rRMU/PDhwtSCsNRlnTNgb3Zp1KFzbU4EO5THMS9KBTAba0N2zArYMuVa5LCH+lZpB1MMHcO7a8elhkQj/kXb2MZ8K95nKhd38nAh2RA3+MIJcfcNJl8vDSmOcGyYxyRoSPDdG+qFDJKJ0yJ6SkVFOER07yKGus2Ttg4mN6fWca7JVv3Yi7yYwPjSF4Itdv1jsBPHOdoAYke1g0T9m8rau+OiM7toaC6q6qh00i6ngnNqHIZ4bd91bB3H5XBQNG7jM7Xgj52ZnTuJ3z2iVppyWg1fakdlHiofPVlM42ySvCbiDIUpn+Ho5G/Nq5xhcRZmyNGGShbVsm8EBlA4L3iXk/iVnFmY5bIlodJZWXPtOavSlQtfrVz/FqpnziKOwEyvrkiJ/FQZOmcr8pFs8ugoLt83Mh3ekY5xCEMLwr4Vj3VOX7JJnrInOY0TjevwnHT4UR3b8+dOcmCGjSkczJk1GIivJGY5bJmKbw57qHPtdjDV7iR/dVzpYRGOgCgGp5OKQJi0DzsJ/F13cZ7gpBa3pcwTncNxxpEWziwkHI51zjCf/VYnkiwXGcVMOsae20T+Ulk5MOMBhXPUkdJE3fE7bCUxy2FLZMLgYQ57qGTtdjBo6BiG6rh0Z1Zn5630hrxDYejFy4GtKEXHXHicVDw0ozGfdL7l0JRFR7M8PgltkRdy6Ijub0XK71sQluiua0iMLk5HyXDOUyn8MVHe9IjPFT2WPrkw28JCsFlsmQF4jjZAtau1gxk2dxa57jxzASfUoSyiYn08JqWZpH/kTf9AeakIhYbLPFbsBHAarWGil8PcFlEVzvGVzyPBaSHrOQcRtZylzuEku5ZKQ++tLwabCzNgouO9EJ7IhHh3bo8LAjvlZQfXlWiatgiVXHgV+dxIn7GzZ7ror0Oo0RDlpIaCm7um0pNXM3CAsVMb1FVl6ORnLW4rGW4IqjR8/WOw7jUZpOuqmOW2/Rx70EXlz24HuW2aIk/6N1+HuJpSsPKuigCR1VwnyxxaV8Q21igizFREObb8GnxrY5bb5nPsQZcc7SC3TbPkXazjkvdmaNf1/tYsMEsWkr6rLm6r+pmLmxztlcRkSPbamA3pN/X8XHuoR2Vp7zdKdxMx9+FzsY5LoIST5n0YbencnDsmc2k82bxUOhpmc+zh2h+qHVys49Kdh7fUU0/fNtvBpS8Rz2EWt10C6KNhNtUeMFYZ5oY2u8jxnHZw0ZPzcwCrZSoCFYF1EPAOmCfoD8LXIfj8bawRf/i95GFGjEc9rghUBBZAQH6HJ8KPu6rCcTH86HpatauhVJeRNb8iUBHYHQI8QeWBQpL+DwBylPnhbDmHAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\sigma \\left(- m0(t) + m1(t)\\right)\\\\- m0(t) m2(t) + m0(t) \\rho - m1(t)\\\\- \\beta m2(t) + m0(t) m1(t)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡      σ⋅(-m0(t) + m1(t))      ⎤\n",
       "⎢                              ⎥\n",
       "⎢-m0(t)⋅m2(t) + m0(t)⋅ρ - m1(t)⎥\n",
       "⎢                              ⎥\n",
       "⎣    -β⋅m2(t) + m0(t)⋅m1(t)    ⎦"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# E[f(x)]\n",
    "Efx = E(fx).simplify(); Efx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1e12af1a-4cd0-41c9-9ae3-80c32c9856c7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABecAAABMCAYAAADneAxPAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2dXdLfNLbu36RS+3IXH1X7vjszCDCCDjMgzQiAGUBxlb5LwQyAEdBkBsAIIJlBegC7qiF1Ls/F6fM8juXI/su2JC/Zkv2oyq9tWR9r/day7HdZf/vef/7znzslERABERABERABEbAi8I9//OMbtPUOlr9i+ReWr5D3GmulnoAYyRVEoA4COhfrsIOkEAEREAEREAEREIGrEniAG9JHUP7FDIDnOP5k5piyRUAEREAEREAERGBEoA90fYc1g/J3WP+E1a9YPuB+bEI93p+8h/UvsXVaKQed+PBiM6NW9JWcIlArAZ2LtVqmTbngT19i+bZN6SW1CIiACIiACIhAKQK4P3iFtjlx7Sbh2L37Xu732P5qsvzoHdemCIiACIiACIiACKwR+GRS4Bn2H+GmI3gzMinb7aIsA/PfYH26wHyv72ZGIW7KEwERSCagczEZmSosEHiJ6xYfSCuJgAiIgAiIgAiIgE+Ak7OmMffnrsADt4E1Z3C99Pa1KQIiIAIiIAIiIAI5BNzrbFjXvc7Gz1tr8wcU+NtaocaP+zxyGDWuvsQXgWoI6FysxhRtC8IHylg+xvI5Fk58UxIBERABERABERCBu9B9Ae8XgKabKHJfjERABERABERABETAigBuMh5i8We8cxY80+9vVst/Ufc7lPgRaxewXq6w01HI81csvIFaTSzH8nMFcWwTo7l2z55PpmQboyfLsXxM2aPKpOhDGVvQaQ+WKdzWmOH4pc5FS3Z72HqPPqyZoL1uVhzW7+whv/oQAREQAREQARFon8C9p0+funfOf4CbCM2cb9+m0kAEREAEREAEqiGAewu+X4+/zlt9Dy/KMJjK8u9iu1hwHm0zaMKfFrIP9vcQCx8IBO+Dermow8coF5VQlq824Idwu3fvL1VCmWhGS+1sOQYZyJ5BJSYX1Kb8QSZviuX/RbtV2SBfk7yaPe8kn2JPqBftV3mS1V0rh1sKM5Q9/FyMtQBkvfQ5FMtpqVwpf0K7X6JfPvj5Yql/HRMBERABERABEbguAdwncNIR/x+457/W5rpEpLkIiIAIiIAIiIA5AdxoMADOn/mvBub7zhnI+B7lGTQvmX5F48/QT/eeP6wZ5HqBNV9HEAqmMyD6JCSQq4tj07qfIY/9LH4IF/VTGYXE2JTX68B3/A869nI5Jv4vITb15VWuxgaeTJs3wY0/TeVrmf6C7SU/zvEpyhflVyzYUirMLYpZ7/Mp49XRiE95DllArcCfeG15hUXBeQuDqg0REAEREAERODmB+yfXT+qJgAiIgAiIgAgcQADBEc4EeAfrlOAE6/xcUtxeLr4eZfgAD7YZROU+A6aj1Jf/HetQ0J5lH2PhTPM/uONS3yYDfZxBGUx926mMgm1tzGQwmQHMIUE2zqInlxsmQ6HMjV7vKmyQqcKoGvShDX/CwlcyfYqFD3tmU69/sk+xQdSlTRb9arZj7wDa+QTLrG96RYttov9duMUwQ5mc8SrIBm0VZ9vL28w5tBOTmvyJ14t/9XYK+okyRUAEREAEREAERMARUHDekdBaBERABERABETAhAACEgxYDz/p5z4W9+75YB/9cQY1S8zS9vvk7PBQH78h/xHkmAZWGaTm7Pa59DEOvEQ9Bk2n6Rkyvp5mch/lkxmF2jHKoywvAm2REwNei7YL1FvLqsIGa0LGHqftsTzBwgdR9KO1tMWn2PasX611XNPxnbnNMoMcNZ2LsSY61TkUq/RSuVr8yZOR4+fwayQvX5siIAIiIAIiIAIiMCKg4PwIh3ZEQAREQAREQAS2EECAhIFcBrN/wzZnkPI1HwxGzs08x6EufYi/DHKGgtx9EZMVA3Gh5PodjkMWzoi/w3pJdpYPBftZj21y9uTQZt9eLiNWL5GoH2fhvlOi8UCbIx7e8d1s4PW56yYYb/IpCos2gn61qyI7d7aV2xwz5Nd2LsaSvew5FAtoqVwpf5r0ydfa8LqmJAIiIAIiIAIiIAKLBB4sHtVBERABERABERCByxNAIINBWwbc/44lFMDle+I5a5iJ70FmmdHrUHCcAcWlxHezzwbBUZ/BKM5SZ3CTwX4mBv6ZPsLChwHfeuX+jTyW57vUu+A51k720StoUMZPXfC0z2D7N4H3vg/qy/ZYnjPuqW8nA9Z+crMn/XZyGfntmm1D9rn34rsZ851der1PYQMLXTINYOFT7DrkV5kipVU7iJ0FtxCzqs5F3xLgzDGX4xgTx0++Com/0qllHOsE2/PPHJMMGUr5ky8KP6bNXx5xof2UREAEREAEREAERCBIQMH5IBZlioAIiIAIiIAIkACCCgzQMvDsAlsMDDFoxIBtF5BHmSGoju13kZ+TGOReCprzg6tfYeE76flu75+x3X1oFmvK9CfWWHWvmOmC99inXJTdyfQeCyCFAiWu7/ffFOn+MhA9esjAXLRLFnzvNwM8j7FmubnE2ZPuwUVXBuWdPHN1Ds+HjHwYQpvwgYfjdRobQC8LXXLsZOFT7PfGr3KEyaxzBDsLbjfMaj0XIRdfM8XxjmPNHdYc9/gQjWNJFeMY5Ng1rTBJlaWIP02EcNcUjqMM1CuJgAiIgAiIgAiIQJCAgvNBLMoUAREQAREQARFAMIRBb84s/Qzb/gdUGRz6DnlDUN6AFtsMtod+HuHYb30fDHTw9TddYL7Pc8Gq95HfBbO8fOowTaG8aRnusy8XYAkdZ4BnLejC+rH9hfo4Ko8PQJ6Dp3vQcRobFPCnFBtZ+BT7O8SvDmRnwe0QZinOwbJgTF35axx/LAu9Gix2XGmeXQKTWNx7MHndC+OuT7GyqZwIiIAIiIAIiMDFCCg4fzGDS10REAEREAERSCDQzZBHYGQIzPd1uyA68i1/rs9AE4NnocT3trvXOTCowqC4nxg4Znr2ZjX8ZTk/wDXXPiu4AIp7jYTLcwEW7k8TZ5dP2UzLkBVlNklkjob4wITrmET5+XDlZUxhlkFZBub5ywB/xv+ZbGClSyxSvxz9bKtPsb0ov4INeQ7TT6eJ/vMejn86PYB98pn7kOVR7Cy4RTEL8AhmFWDr+unGKbT/JzL+ieUnbDcxjh3IxLGLXe/pT7FjdazsKicCIiACIiACInAyAgrOn8ygUkcEREAEREAEDAl8jrb8AK1r+hE2GKRbCjK6srFrtsWAyU3y+umCjJNAFcszCB+Sh+WHgD3bwYKsYD8ugNI9eGChpYR2WJ5Bd75uYrdEHdDZ3DviN8uB9r9kI1iP7N73y0PN28BKF8KwTJDL3KfQZvfLh6mcyOcrmfgRYP8XKNNiN/soT/9j2uQHb5qw+QuZzLnFSIZ+Tdm6PskYC8c0ts8x+HPs82FZ91CyP87iofGSLJgOGccg2yFM3qhs8xc6WPuTO2dsBFQrIiACIiACIiACpyNw/3QaSSEREAEREAEREIHNBBCgcLO9fw809iny/JmcgSJFshicCs0AZ6BwNIMd8jOPQZYuvw+4YLeT2wWwuO+SC3T5enEGa6gs63zIP2h3KO/1wUMusd2oQJmrcNQa8jNg+xDrITCPbQZwnS9QtDPZwEKXVHNZ+BT7PNqv9mZnwe1oZqu+0p9vn2DtgvH8PgUD3o8n5yHHndDYdLpxLIHJKl+vwB7+5OzTxPjvsdGmCIiACIiACIjAzgQUnN8ZuLoTAREQAREQgUYIMHjBNAosIFDyCHkM1gZnSLJCZmI/Lpgx18TfcWAIhrMQ5GEdyjOdwc7XcvBVOJxRz+PdTF+s+YHXLrCOtZ84G53lX3uZlMkFu7zsbnMUnES9bibytBD22bffZqDI8VmQn3blhz6HwHwvFfXy05lsYKGLzyZm28Kn2M/RfrU3OwtuRzOL8Q+OTz+4ghyPsLgPX5OBS1cax2KZODYx6z38yV073LU0Ri6VEQEREAEREAERuCCB+xfUWSqLgAiIgAiIgAisEGBQCEUYCHdB7TvkMbjFQMmT/jg2zdJvaCkUNO866PtmIH4ahGeQkLKNgvbIYlsu7wscdzPov0f+H9gfgs7YZrts5zMsfnqJnblXyLBOFyzr63OGOctPE+s7OabHqtiH3M6uDAR+0y/84C9tTXZOT5Y7hQ2gk4kuGQa08Cl2e5hfHcTOgtthzBL85I/etkMV8P4SO6OHoci70jgWxWQAFrexhz9xjOGY+jpOJJUSAREQAREQARG4KoEHV1VceouACIiACIiACKwS4OxzBmsZZGB6H8sHhYINDJa8g7a5hIIZlIGz4KeBbgaLR6+0wT4T3zXPmeB8Z/OPzPASg3Rf45jT6yPs/w37lMFPfBDAj6OGEj+0yQA227/Deu7d3XxIMA36s0pNiUF4smAQcJr82bpnsoGlLrQ//YQzZd3DrF+RR3Y/Y81AqksWPsW2jvQrM3Y7czuSmbP/2pq/XOFrbViO4yDHN457oTHuKuNYNJPK/In2+R2LkgiIgAiIgAiIgAgsErj39OlT/oz5BRb+sz39p3Sxsg6KgAiIgAiIgAiIgBUB3If8ibY4K38agLfqIrkdyPIKlRjk94PUUe2gDoOYDM4+jKqgQkECZ7PBFn0ICPU3+xXa4C9Hkj8IGzTQTplbuFkwi1WzRrZHs7siE+jM/69/xHruwW2sS6mcCIiACIiACIjACQngHoGTvDjZ6979E+onlURABERABERABNok8E+IzXe515Q4Q370SokE4Tjjc27mfUIzly96Nhts0YfOYOFXfAAWmo1ds7Nt4WbBLJZNjWyPZncpJvgnm7944AQ4/1czsf6jciIgAiIgAiIgAhcjoOD8xQwudUVABERABESgYgIMZHevialFRgRZGFzhDGPOVo5OfXDmMdaaNRlNLVzwbDbI1Yd0UJdBv81+hXb4LuzkX4OELbRPLuRt4lyske3R7C7IhN8weU699zk71IsIiIAIiIAIiEDLBBScb9l6kl0EREAEREAETkQAgQy+Xu8XrPnKjZpSzqzb7sO5NSnRuCxns0GOPjTh1f0qh9vVmblTX+wcibfrUkz4a6tnb7vRlgiIgAiIgAiIgAjME1Bwfp6NjoiACIiACIiACOxPgB9P/Xr/bud7xMMCzjD+AuuoWf19OZZvambyPIHjj5zNBqn60AKoQ/+7tF+lchOzt+eu2L1l4bZKMOl9jrPm9S03B1prERABERABERCBRQL6IOwiHh0UAREQAREQARHYmwCCGpw5z4+wclajkgiIgAiIgAhUTwDXLL7+jB92q+3bKdWzk4AiIAIiIAIicDUCuF8YPgj74GrKU18A4EeR+M5O3kBxVttXyDv9OwGvqjfsW02SDaoxxakEKeFXJdo8FfRMZcQ1Dhw4cdbhHRa+W5sfElQSAREQAREQgVUCuGYc+X8eX2fzZFXIRgsczLZRahJbBPYjoHN0P9bqSQSsCVwuON8PWJzR0P3UHGu+h/JXLB9Yw62pvavqLRvURECylCBQ4twu0WYJ3VtrU1zTLAZez9NqqLQIiIAIiMCVCRx9nUX/p/2119Fsr+zX0l0EYgjoHI2hpDIiUC+BK75zfvqROX6s5xEGM86iP3O6qt412VQ2qMka55GlhF+VaPM8xPM1Edd8dqopAiIgAiIgAmsEdJ1dI5R/XGzz2ammCOxBQOfoHpTVhwgUInDF4DxR+oF49zobP68Q7sOb9XW8kt6Hg/cEkA08GNo0I1DCr0q0aaZwww2Ja8PGk+giIAIiIALVE9B1tpyJxLYcW7UsAhYEdI5aUFQbInAAgSZea9PPaud7Z79fY4QyfKH+L1h3r62Zlkf+w0neo37/90n+obuQkwNrlM4UFOVPofeh0CedX9UGKXqv+d0EqXZ7ApaM0Zb5mFaizVaNX7utcrha6pTT/5nrpLAlB5RfvHYfzeps+vTMo++varfP0f4R6j/FZ1rgezZ9QjY7Os+SMdoyvyc6ms+W/s/I1lKnLWzPWDeFLfVH+arvYXwbpejWkl4THQ8f/67A2Wces90Sk5ZkjWFvVWYvLveePn3K4PQLLB+g05dWCsy1gz74IVZ+qIczt19h4SDy41zfPYikr96jDt8jz4+8BgP0ODYklKEMbP/bIdN4A23zH0F+IIiJ20yUL8i7L5+kMxtEvar0pkxzqdcxmslcO7H56K+431GW2myQyrkvn+R7KTrH2qulctC/uG+lMEZZ8zGtRJtH2PiMtoJOxa8v6CP62lLSrqm6bpVlD3+hjLXwnfLqeSddD/bUp5cv6T4iR6ea7JOj89SuKfvor6rrW4rsFmVb9xcLBqXbKM0Y7ZvfE5VmMtd+zyp6zGuBLWS8zBiTar85P4jNT+0vx18oC+rtfo+IPi/jN7H2DpUDp03jnzjfUm2JCWS9zP+It5ayy+k5Jv0/hDrR4yLK8iEn2793xMz5X9H5M3TefegNaw6uL7D+GEsomE7Fgl+9d3VxfFr3M+Sxn8WPvKI+HxJwln3JwHx38UAfgw59v07nXyDDNOXozDaq0XuqkL8P/XOY+E3kbO/hd5SrGhtkcs7xvSidc4zWSJ09fCuKcT+2mI5pJdo80K6nstWVzvFMXbe62h7+Qhmjzu+tyrj6YMl3kv6A5S/Yfu3yA+uc6wGbKa7PBn/I0am4PgH2N1kbdL5pKyFjj3NgV74J/k9MzfpLgo2LFE3gXIwxZCj+f14ReIFGM8//FtiebowJmO8u036hpqLyMvvL8RfKs+sY3gO4hN9EGXumEHzAYvwT51u+TTDZcQw44vy/tUpGDhiV/H8oi8v9DD2yqwAAnwr8FesuMM+GsM1/DLnPC8Io9eV/xzoUtGfZx1j4ROgP7rjUt8kA1Zcub7ru234H6y+mx4z3+Q8wjTMk9MlZD9TbTGc2jnbZZi16U6S5lMRkrpHYfHDZxe8oT2U2SOLccypyvsXaKqUc5P0Ey+w5ntJWbtme2WnHtF6/zeOkbDX2MCuuaLWpc3yjHyTpOiaevtfbqPi5TcnQ1+q1O12DcQ30wfP4Jyzf4cinWN4ZlxjvoRyvm8nXA7aCusX1QTfJ/pCr0076EN1aStZ5rcGl4z2v4ufAHnzRR5L/k0uvf/I5sIc+Ibuh3xruiZI4l2Tct735/qX3hcPZQo6k878Ftr2MpxhjQufkJC/JfpO6ObtJ/eX6CwVD3T2u+QODFv0GMu86hvSMNo1/rXHeg3FjTHYZA8DE7PzfyYbV3KcMg5q3sWtwHv1y9nhopvhvyH8Eg0z/WWQQm0/95tLHOPCyd4ppmWfI+HqayX2UZ1D/IdZdYJ77WNy750NVtuSxL742aJrIgc4x7XeLzuyjFr2n+vr7qUz8ujnbe/od5avFBqmct/jerM45Bmuozp6+NcuYYxiYmY5pJdo82K5ntNWVzvFUXbe6257+Qllnz++tirA+zufXWJ5g+QK7vOdaS1uuB2y7qD5oP8cftuhUWp81e/B4js4x7c6V2fMcKMo3w//JpHV/mbNrsfwMzkUYQw6eK6b3RMWgxTecev63wPY0Y0yEGVPtF9HkYpHU/rb4CwUpOoZPNL2S30xUX981HP/E+RZ3S0z2HAP2PP9vrZKQU8t9ypzIewfn6SShxCcuTMNxgOOM+Dus52bNu/KhYD/rsc1/YT202bfHYDgD/r/hGJ9i8ucMvCAt9YPD2YntclbA9MHDTYMos0lnNog2atH7Rj8vI5qJV2fL5sgHvIbM/Y5tV2SDaM6QeZPvzenssT7r5m6+NccY+eZjWok2K3CAM9rqSud4tK5Gvrabv1DeufPbSJekZiDLpuvBTvok+cNWnSqxT5LOSUYPF97tHKiE70DhJP4y6FPjRinGaNf8nqgSftHnf0NsrzTGRNvPyN+i+9vqL5QXbQTjD0a6TJu5kt9MdV/chx0sxz9xvqXdEpPdxoCdz/9bqxTKgV6b/h/K4fIgRhc0TEfkLHUKyEA2E4PaTB9hYaD7W6/cv5HH8t8grwueY+2C06NX0KCMnzoAfQbbvwm89318gWNsj+U5456vh+lkwNpPrD99wsX3RLHu6JUyaIMXlSH1/WzSmY2hnbn33nPwZOKJ45KFzmwrW28nSMl1LBMLG6CNI/xukw0s9KYAaGdv3wv53ciVet14Tr7CwrGD358Ifhh5VNF4p5dj0/mNNo7wrRDjqDFtinDFFlltTvuw2L+6raA/HybzmsrE6xRfs9CdM1hXd453Uhb4E6sryrV6v0JqN+e3hT4Z5ih2L2KlT6w/eLpb6HRjH6/94puxOlswRhtVXN8sdMk0zG7+0ut46H3RQZxLMa7m/oW+Z8UW7dRwvU9mC7mD9zDIr2KMyRwfkqvF2g/lNt/D9H63t7+w29VrZK/f6niHcvKbGS9bYZh0jorzLeSWmMzJSq1wbO8xYPX8v6VdJmeJS2KPpe5TZsWICs6jNj+4+hWWn7HN95X+jO3uI6pY8+L6J9ZYda+Y6YL32GfQmQHwd3kA6b03qy640G8OKxewf3/IeRPcHwXQeQzt0vB8rzphPcaaAba59AoHGMgfEso7eYa8mQ0LnYNNQwZeePlggQ80GGxxibps1ZltbdHbybLreoaJhQ2O8LutNrDQO2i/Gc4sa+F7N37nC4G+eaM1vP+ul4U3EbHnpN/c1m0Lxkf41g1jcEzmt2aLnDa3GmSh/mVtBTvwlWi89vK6d4c1r8G82Rpd13jMJZQpeX258T/X7xHrGV1b9RciDPG10CfVPBbXg931mfEHp7uFTiH7uPYPWc/obOEzVVzfANVClxzb7OIvsF8t90VHcC7CGEyT74lyHCShTjG2M+c/RauCLeRbuoepZYxJMKVt0Rn7teovhLN4jYS+UeOd/Gbez9YY4nj0+CfOt5xbYrIi661yyEGd0/+PmMMlCOtNZpFr6UJ/d/eXDvIYFHyE1W99OQaU30NeF5jv89zF9X3kdwEEL989Fe+zulUozz/uttmXC9q7PH9NWGszblk/tr+h7QI6D233G3zA8Rz9uF8huOMWOrOtLL2dEAetR0wK2CDWDw61QQG9p+YccfYOWug963fQix8X5Guk/KAiH+AxWM8LxW6pAOM9fWuWcSzAmmyxJvOVbQXdeU7yl2H+dZXXjOl1Y4rxkHN8KsRO+yNdG/cXIhud3wX0iTWLxfXgCH1G/jBR1kKnkX0m7R+1O9K5gM8cdn0roEuKjYr7C/Sr4r7oQM7FGacYvETZHdiOzn9Ph8PZQnfKEHMPc9gY4/E6anNkv8b9hQxnr5HQLWq8k9/Mu2Isw/kW3h4R57cs3FZLTBJkdeq59WjMcZlYH37N8GTJ3tzAZa7P3bk8gCT/jeV/sfy/Gan43nZ+dJUXTwrIoLifGLxnevZmNfxlOT+owAF7LrkAv/vpPssx7/VcBeQzoPd84TgPMfhHmVOTlc43/YIjTwrO/PeDlK6chc5sK1dvJ8do3dv+V2TG3kDRbp+h3stRQzM7M0ysbHCE31HTXBtY6X1De4azK2fhe0s6c7bENKjozs1Fv4LcrBsK4LMeHxZ+6pTw1uTIn02GkhXjI3xriXFI11Beti1Cjfl5spVPI3sMcI10/gWmfyLjn1h+wrZ/TXXlhjWOl76+RPlfAT8YdHQbM7q2fG5TtSlfK30ctti1xfVgV31m/MHX10KnqX389kfbkIfXp2L3TexsRmcrn6nh+maly8g2kTt7+EvWtRh2r/WeKBLtUGwPxkNnMRstsZ05/52aNbBdu4epYYxxvEZrsG19/B7pw50d/IXdLF0jY8e7Zv2m52w9PrNZl2IZuvJL62Y5w5dLMW6JyZqsN7bfYQxYOv9H8hxow5EcETt7XUv/C7IwHn/H4Pz/wfI/WO5juUmA97rP7IJj2J8GBxiE5w20K+faYPkhYM/jWHiMSk4TL4JMNOpqQjssz8Aef+JvnjxdNuk8FQztfsk8rEOB+Wnx0T7qFNV51NlkB33TtnPvrZqUTttF20EmfZ9sbJMN2A4WtlO931FIK73Zlp/QbpCzX2ZuG3U3+R7q04ZsgwFGPz3qd176mdNt1J8G9bsiyOerrfix5W+ndZb2UZ7+zHQp36LC0H2TLdjGUkL7stUSoIRj9FMsvL6SKWccfY59PtidPiDvWkX+Yed4J4D3B7KY+oHX9KKu6PdU57aVPlN+W/YhU/b1oJQ+aDfb98lii05zLHtdi9w39TIHde77ZZHmr29WuszZKDffwl/QRva1GHVNx9caOVswzrFvK2whZ/D8j9F5L7boZ/Eepj9Okav7/4yyQa5mx++pH0CfQ/0F/UePd71fzN779ser9BsKBflMx2dnS7QbzdDVWVr3HJvkXJBxM2PWmv2mtkf5Q8eAgDylzpNFG07l2LIPptn/DwX6/b/I+x/mBwPygQrM4gkcCqJxsBjNYIewbgDp8nvh2QYD+1RkmtyF2Q/884lQqCzrfsg/aHco7/XBQy6x3aiAv6swWVvo3DUJ+RhIfIj1EJjHNgOLbvYwy1nozHa26s02iqdIJhY22NvvyG6rDSz07mwYydnC9+Z05gz20AO8T/v8Ledop2PmHwvGe/vWHONYBLXaYk3+S9mqvzbwNVAuGM93SPJG5jGPTWFVcI5PRSq2H6lri/5CZnPnt4U+KTaxuB7sok+kP1AWC53m7JPCdnPZSJ0tfKaW65uFLqncS/tLjdfivTmXZpxq85LlzdhGnv+HsoWM/P825h6mljGmpO1HbUfar0V/oZ5z18io8U5+M3KV6U4Uw2ml0L4431JpiUmCrJ2iKL9XDHLu/L8FXiAnlUukCLtfS1OC83+HEryIDgkQ5p4YcADhq3AYkGMggcF6Jn7stAusd3tv//DpNMu/fpvVBdVp5FAaXbRQj053E7Do8/w2Q20t5VnofAf5HqETfuBlCMz3nVJuPzFIuVVntkcWW/RmG0VTAhMLG+ztd2S31QYWeu/te3M68/yfjh08J7hwrDgqWTDe27fmGMcyrNUWa/JfzVb0qx8cFIyXnAnQ/UIE69HDLOzveX3Z6jduVqEAABjDSURBVH9Opax1gq4t+guZzPG10CeFeel7ERN9EvyBulvoNGefFLabyibobMG4luubhS6p3Ev7S43X4r05l2acavOS5U3YJpz/R7ONvYepZYwpafuh7QT7tegv1HPuGhk73slvBm+52YhleFMxkCHOt1BaYhIr613CmHP0NePWIuk50VwSmt6dy4MY4WBYDrYMxE9fI8OLBw0/CrwhiwF4l/cFjnc/XcD6e25j4dN0f1Y92/kbFj+9xM7cT8ooC2HdoR1u8wn9aPY+jyGxvpOjy4j9g/asdGY7dBZ+AJbvyGKizAzC80M5/is5LHRm+9l6s3Lp1LNdZWJog739jgizbWCo996+d6MzdOnOT/CgLF3q82h/jg30+d2TIeO9feuGcSy8Wm2xJv8VbQUmfFI/enAFDl8ib/QzwJ7N6ljaM7a4vmT7Xy9D9ipW14b9hWxu+FrpkwjewleK6hPrD57eFjrd2Mdrv/hmrM5WPoN2Dr++WemSYZxi/gKdqrsvOohzMcYZ9i5WxYpt304r1/s/ADTmHubwMaaY4ScNx9qvYX+hxjfXSOiTMt7JbyZ+w91EhoEWbrLE+QZJ9P9dNYxZsfY7PA50i7loThSXRAleojzHtVDi2GYej44KzqNjGpez4KeBbgoVCorzXfOcKc735P6IxU9U8Gscc8G6j7D/N+xTeT/xQQA/bhdKDHJ/17d/h7Uf4PbL8yHBZ35GwraVzryRYlsMrExTZ1Av00JnNrdFb0+cYpuxTKxsQEX29Dv2t8UGVnrHcqa8Fr4X0plP+pkYiHfnwEPsP8H+9JzvCu70x4oxxd3Tt0KMY5HVaos1+a9oK/7Kig+xyeY1lu4GAPvT620N5zhl3CPF6tqqv5Bh6Pw20wf+w3sqTgxwY8GvyON9yM9Yf4+1SxbXA7ZVUp9Yf7DUKaSPa3+PdazOZj4DpY6+vpnpkuD/tKXFOTDnL+78q+m+6AjOJRnvcT7G9mHFNvb8L+2/MXrH3sOwraPHmBh9LMrE2q9VfyGj0JiXMt7Jb8KelsIw3MI4V5zHPLjXEpNYWWPHHOp/hutxLJe7hPvB3bnce/r0KX8S/wLLBxD0yGAZHWOUIM8rZDDIPw1ij8qFdlCHFzf+w8lAYDNpi85UslW9azLQVW2wRe85v0M+g0EfYs2bb9OENj9Bg8kfhDUVIrExyHzYmFbSFmsYZKs1QvscP9L/qGGLfhBrmS1sezZV3bOcTZ+e8WHjb6wftVxui8+gblX+X9JfoGuR+yK0q3uiQieQ2BYCm9gs7KAxPJFZbPEtbNkH6gfHcOQXGe9i9eplO9xvwCF7fK6BYQxvyHko5y2MY/TLKXM0kxSZa5D1qjaE3pzQzonn9+6nGO2AspwhP/oZf4IMfHrCC0JraYvO1LVVvWuy01VtsEXvOb/jK6t+LGRc/pJnOpO4UFdmzZZgHCtcSVusySBbrRHa5/iR/kcNW/SDWMtsYcs+5sbQ2P6ty51NH/LZolNt9rG2t0V7Z+NbSp9S1+IWx9dSjC382W9DbH0ax2234i/HEcrveQtb9jp3jSw13qVoukW3Ob1S+mfZLWNIDQxj9D2a8xbGMfrllDmaSYrMNch6eRtWPXOe3oQnCPw5AX/6GT17HmX5CgD+XNt8ti5lKp1ydKZMretdmmtK+1e1QY7ec36HfM6i4FP06n6Vk+IL1mUtGcfKJlvEkhqXO8JWYwns986okz2lvBZz2LIn1KvynuVs+vSsL3dPmefNebVyfKZW/y/hL9BV90UT1zqbz0zUO3T3jGzPqNOhTuJ1nsOW1VEveA+D/GrGuxzd5vTykBXfrIlhjLKtco7RLbdMS0xakjXXHjn1SnNB+83MnCe/nCeWfMfS6EM0OYY4sE6OzhS3db0PRH7T9VVtkKP3nN/xpozfqqjqdVk3lt4/w5JxrPSyRSypcbkjbDWWwH7vjDrZU8prMYcte5obQ/OksKt1Nn1IJkenWu1jZ2m7ls7G11ofXYtvfc2a8W0P1805I9sz6lSLh+awpexz18iaxrsc3eb02tNeNTGM0btVzjG65ZZpiUlLsubaI6feblyqnzlPegjucWB6jLX/wbIgWJThk4dfsI6eaR9s6ODMFJ0p6ln0Phj7qPur2iBFb/ndyGWid8Q4GtXhBc9oqzPqdLij9AKksGUVlK/6nuVs+vTML3dP2bvnLqsUn6nd/+Uvu7jM3dl8Zh9qcb2cke0ZdYqzZvlSKWwpDcpXfQ/jE0vRrSW9fB1r2BbnWyu0xKQlWW9Jl8spyQVtDzPnmwjOl8OslkVABERABERABERABERABERABERABERABERABERABERgHwJ+cP7BPl2qFxEQARGoiwAGQn74hO9J7H4yiPVXyHtdl5S20rSuc+vy21pTrYmACIhAHIGrjp1X1TvOK/JLiWs+u7PVlC+czaLSRwREQARE4CgCCs4fRV79ioAIHEag/2fiO6y7119hzfcK/oqlyY9Ix4BsXefW5Y+xkcqIgAiIgDWBq46dV9Xb2n+m7YnrlMh19+UL17W9NBcBERABEbAncN++SbUoAiIgAtUT+GQi4TPsP8I/GpxFf9bUus6ty39Wv5JeIiACdRO46th5Vb1Le6O4libcTvvyhXZsJUlFQAREQAQqJ6CZ85UbSOKJgAjkEegD7Y/72qGPRLvX2bCIe52Nn5fXsVEtJz/Wlh/C9vWrVuceYXU2K2QTI48p04zTecEmZTqObNXJZ3yeRPZuX0z63IXOe3vQlbdYwA+qHvtpjgI6s9nq9aaQpVMBtlVwLaBXaVMc3n4BZof7QgGdDreThQAtcWlJVgvbqA0REAERCBHQB2FDVJQnAiLQNAHc5H0JBf6F9XMsfK88X1vDd8q/DCmGfM7+YZl3se2C1qGiWXlok/+8fNVX5jbTkjwsw9fufNyVjPiDsk7H7lU9a1VQvqjOa/1Pj0OeXW2G/qqzyZSJ1X6qrq7fVJu4erlr9Mdzld+C4Dn4CstDLD8if+68LX6eoP/dUm+noud9rjK9bNFjGPupWZ9cDjn1wKEqv4Y8VY39e/lKjXrLn3IIvKmj8SWdXWlmFucY2qhqvEynXKbGmbmU9ssyFlGrIiACImBDAGPg52iJ///d08x5G6ZqRQREoBIC/QDH2ZddQA/r11gY8OMyF+zmsSIfhEXf3T8aWD9xiHp5XmD9MZZfXL63ZqB9KO/l36E823uBhXX9QPxnyEt5b34xnX15Y7ahBy9Ku9msZ/gN1gNjbJPH0TaJwZVUJlPXO9TLsUmSbIHC9N9n6Ps5j2Hd+TrWU193Vfc4T1xfWWvIzkDoD1j+gu3XK41UqU9vh9TzhapWqc+KDUocrs2vdxv7E/x/D1/ZTe8STuS1eWp/qsxnPOz1blbEzOIcq82/azF8c1wq8stabCg5REAERGCRwP3FozooAiIgAg0RwI0gZ9LeYT2dacsgX3dsqg7K8p8JBoa/nR4z2mdgjoHzIaEvzkBloI4BiVHCMQZFf8faD7z7ZfiqHuryh5+J8myPenAG+mJCmc06o41PYvpaFAQH0cYRNqvOJmucNhxP0pX95Nhkg3xdVfRJv/8r1l1gnpnYpk9z/5DzhDLkJMj9DpafsHyH+p9i4fizmFC2+Hm/KMDywRwfOlwfMDUZo5bRLB/t7VqNX0OezWP/ssbdeZvk/3v4voXeaEP+NDG+BVc2iXaq85mJqqPdSnyhKmYWvoA2mrsP2MMXWuICWavyy9GJqx0REAERqJzA/crlk3giIAIikEKA/zyH3tHOANk0YH/X3/DyRvKLlE4SyzKYzpnu08QZ8+z70eQAA/cMoMwlzv5/iXoMXE4TP2z79TTT399JZ7/Lte0jbFaVTdYAbTyeqiu7S7LJRvlcdf6KIfQrkt+Qz481TwPcRc8TJ1TOmucmlidYOK5Q/phUrT4QPseHatYnxh5WZarxa/gjA1+lr3d36CfV/4v6yl56WznMSjun9KfafGbFBlUcromZ4TlWjX9XYeS3QjTDpSa/fItPWyIgAiLQBgEF59uwk6QUARHIJIAbRQaWuDAAMKQ+/yHWXWCe+1imgfKh/IYNzoDnzMlpcPGmSZRxs8jnZs2zDnUJBTHvUJ8Be75rn2VuUp+/h843fadk9HKWtFk1Nknhklk2Wtel9udsslQn8VjQZ9GGewg1HIcsRc+TRLk3F29AnyQfakCfzTZLaGDw20mdXf26P3+rG/tL+0qtek98IWX38v5U2mdSjNFK2ZLMjM+xKvy7QruekktJv6zQhhJJBERABFYJPFgtoQIiIAIi0ACB/h8EvtaFgTsG4l9h+RQLg+Kj9z2jDIPwnJ3O91vzndBMDNJzdkqX+vY4S921x3xX9iNs/4Yy33rl/o08lue7mYfgObY/QF4ouQcBfiCe7Q91XaW+D8pHXSgPZxLzVR+dDFj7ifVvZtmg/KrOfiN7bPd6mdnMlxlt0760CRMDYXxVkPsOQRU26SQr/Ac6p+h6h/L8JzDKJhS9L7/pPEEb9Gum0aua3mQNf7uAfL9X5Dyx0GWQNm2jan3AJcmHoHoRfdKQbitt4QsV+XXU2G+hcwb1Yr4CfaL0zpA5uYoF29b8aQoJ8s9ek6dlV/aL+cxKvyaHLXwhQ5AizKBL8jk25wfIr+I+IIOtSZWWuMzJmgGiiF9myKEqIiACIlAFAQXnqzCDhBABETAgwIA1g+W8wWcw9n0sv2P5OxYG9vzX2vDDSiw3epc16rKeS/wIJT8S+zMy+O7on7HdvZcea9b9E2usulfMdLPysc9AO9t8lwfmEsoxAEqZKO+oT+SNZGIbKMOAOwOmvJHlDH8GQ+cSH0p0vwaYFIjReVKl+K61zTqBwYevEaLtugcdWNOGDDCGuLg6R9ik63vvP+AxpytFSbEJy1ucJ++xIST/XHiT8zZgz/PZJfp/ifPEQhcnY8q6OX1WfKiUPilMt5a18IVa/Dp27LfQOZV7SV+J1TtV5pzyFmxb86eBE8aL5GvyUPl2o6TP3PZmn2PhC6lSlWKWdI6t+EEt/p3KdnP5lrisyJrKopRfpsqh8iIgAiJQBYEHVUghIURABERgO4EueIcbRwb4hlfYYP8V9vlBQwZnu4S8teD5IxT8rS/OIDrfoet/MNb9E/E+8v2Z7sxn4H4tMdj/HHUHOfsK7Gtp9jBvZP2HDH210Yr1b2RY03nUwn47ZjZzIkNPMmSA2bcLOfu/UHDF/fXuNvE733l7TleKkWIT6/Pkxm9nuJifJ/AXa11mRA9mt6jPkg+Z6xOkViizgC8c5tdEBH0Wr3d9maP8v5ivxOhdyIVGzV7Rn3wA0J82zrkm+83428V8xu+kxHYBX4gVswizlHMswQ8OHS9jgVqVa4lLgqyxeIr4ZWznKicCIiACtRFQcL42i0geERCBXALuFSbT+gxmL31gdVqe+3xvOz+6yn8SePPIoLifGMhgevZmNfxlOT8oPBxwG2iTQS3Ogg/N4n4Px167soE1Zzw/D+T7Wd37of2MLduQk+zY7zSRzXs4/un0APbJb3hFUOC4y7K0mWuze7iB/v9Exj+x/ITtZmwCWcnVzUZzOi2t6S+fod7aQ5uuDZRb8j+WSbGJ1XnS2awT8PYPzwkmX64S54mVLm+kTfvblD4RPlRCn1mikMd6jLLyhRr8epbb5ICVzpNmV3d39ZVVaVBA/jSiZHE/kXxNHklwu7Obz1TsC7dUlnN2Y7YgxpofVD1eFvAFh6olLmuyOp1i1zX4ZaysKicCIiACxQkoOF8csToQAREoTQA3zQwezwVgGVxfCnjfiIf2XPkuKI39adsMwjOY4cq5Nlh+GrB3x+5Q/kvuYB0KzA/lQhuo8w7yqQtf0bJbQr/T2f1d38jnK3b4oVv/FwXRcqGeqc1cx2iXv3KgfSj351ywz4ch0wcsXRXkV2UTyEOfGn7l0Qlp9GdNVxxPskkvK6XbdJ6wHSxsh/+oTRP9nmntlw9dIbSTdZ5Y6dIJYfinNn0gT/b5Qiy5+iwhRZumYxTa4znI1Lxfv1Fj/a+Vzus9xZco4SsxvcufYijFl6FvYYm+Jse3fFvS2mdq9YVbzfNzrJnNSYJ+Fv2gP87qh90HzMnOfMhnep1xfbXEZU1Wp5PFGn1l3ctZ9K02REAEROAoAveP6lj9ioAIiIAhAb7D8+VMe/yncBpcnyl6k826oXYZtBnNYEf/zOPNZJff31hi903CPoPZD7EeAvPYZnCbAXeXOCuFbYTSh8xE+UEXbIfK8h+bqEBmqJMd88xtRpZYPiEjLLQdX+fAf6ge89hUN5ZF3iVsEqlrrk0szhP69Zw/03SD32O75HlioQvlTUlN6BPpQ9S7pD4pXLeWtfCFWvw6loWFzrF9sdxZfCVGZwu2TfkTr7scN7BEXZNjIJ7EZyx8IRJXV+zQ8yzBD5ry7xQDhMq2xCVB1pCqc3mH+uWcUMoXAREQgaMIKDh/FHn1KwIiYEmgm+E4bZA3k8jjseCMl2n5wD4/JusHBe/QJgOIbHc6g52vceGrcDij3vWLrC6g/ggrBj6HwHx34O6OwWE/MagemjXEMqN/5tAW67KfaWLe62lmhfslbMaPhPL7Al0CI87Uch/xHT2wQP5lbJKga65NLM4T2q57ANWbz634KwKeV75PlzxPLHRxsseuq9cnwYeoc0l9YplalLPwhVr8OpaHhc6xfbHcWXwlRmcLtq35U/Q1OQZgX+YMPmPhCwnIDj/PYv2gNf9OsUGobEtcYmUN6TmXd4ZzeU435YuACIhAMoH7yTVUQQREQATqI/AHgkd8hcmQ+mASA+hPsM0bwKSEOgxyMxA/DcLzn6o7HB8F7ZHFwKLL+wLH3Qx6tsObWgaKv+mX77BmHsv5sr1E3twrTShLVxZ1uM0ZaSw/Tazv5Jgeq2nf3GZQjrNwRvzA6EvkjR7OIO8yNknUNdkmffsW58n3tB/aGx5YYZvt8nz7DIufipwnVrr4gkZuV61PzyV2DKPKRfSJZGlSzMoX0M7hfh0LxErn2P76cs37Soy+Vmxb8qeeyx9Yr16TYxh6ZZr2GStf8HjEbB7NLMoPGvTvGPZLZVriEiXrkrKBY0f7ZUAkZYmACIjAcQQeHNe1ehYBERCB7QRwM+9mxn+IbX4YkB+OfL9v+QPkvc7shQFczoKfBroZMBy90qZvn++a5+x4PiT4sc/jikEttsUg8TT5gXke44MAfrAzlKgbg/rdQwis5971zocE02BmqL3D8iB7KZvxlwn8CT11o91pK9pwaq8r2SRK1w02sTpPaDMGcb6GLGyT6SMsf8M+/4HzU6nzxEwXyMzzmL+Coa8z/Yo8nu8/Y82ArZ9q1yfKhzyFSunjdVF808wXIOnRfh0Ly0znBP8/g6/E8DVj25A/kUvsNfnuQj5j5gsNMYv2g8b8mz6+JbXEJVrWhvxyi+1UVwREQATMCdx7+vQpf9r/AguDWNN/gM07VIMiIAIiYEkA4xZno49mRlu2v3db0OUV+mSQfxq4XxUFdfhPH4N/D1cLbyyAPjjDmbP35x4SzPaAOk3ZDPI2YZNZ4BEHZJMISAWLnM3HatAHMmSPUQVNvdh0DdwWBSxwsBWd5U8FjJ/Z5NE+I1/INJxxtaP9gOrU6As1cIk1dUuyxuqkciIgAiKQQgDjICdecgLmvfspFVVWBERABESgOAHOkM992MCZLXMz760F5y8KpjPSrfuopb1WbFILrz3kOJtNpM9br7Eax1oco87mB2+tOr/Vis7yp3kb7n3kaJ+RL+xt8XB/R/sBparRF2rgErbYbW5Lst5KrxwREAERMCSgmfOGMNWUCIjAvgTwhJGvLfkc6+TZ2/tKmtYb9OFP/afvo19spGfBV2eM3u+6WOmAg72czdlMNjnAWVa6PJtNpE83C5FjevXj2Iprbjp8Nj+IgXFFnWO4WJQ5K9uz6mVh87k2zsjsjDrN2S8lvyUuLcmaYgOVFQEREIEYAhgDNXM+BpTKiIAIVE+A71fnrJWzpZyZo3wv9JMGQLRqM9mkPuc6m02kz5tvdLQwjpU8G87mBzGsrqhzDBeLMmdle1a9LGw+18YZmZ1Rpzn7peS3xKUlWVNsoLIiIAIikERAM+eTcKmwCIiACOxDAE9R+f74x1hPPx55I0D/xPUXrJPfU3/TmDJmCcgms2gOO3A2m0ifO41jOJvO5gcxA8QVdY7hYlHmrGzPqpeFzefaOCOzM+o0Z7+U/Ja4tCRrig1UVgREQATWCGD8G2bOKzi/RkvHRUAEREAEREAEREAEREAEREAEREAEREAEREAEREAERMCAgB+cf+C19wIHvN1u8znyrv7z4ikT7YuACIiACIiACIiACIiACIiACIiACIiACIiACIiACIjAIgHE1l+hAN+OEEwMzvM1CF8Fj97dvZzJV7YIiIAIiIAIiIAIiIAIiIAIiIAIiIAIiIAIiIAIiIAIiMA8gW9w6J25w/8f/46NC1K+9FEAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\sigma^{2} \\left(m0(t)^{2} - 2 m0(t) m1(t) + m1(t)^{2} + s0(t) + s1(t)\\right)\\\\m0(t)^{2} m2(t)^{2} - 2 m0(t)^{2} m2(t) \\rho + m0(t)^{2} \\rho^{2} + m0(t)^{2} s2(t) + 2 m0(t) m1(t) m2(t) - 2 m0(t) m1(t) \\rho + m1(t)^{2} + m2(t)^{2} s0(t) - 2 m2(t) \\rho s0(t) + \\rho^{2} s0(t) + s0(t) s2(t) + s1(t)\\\\\\beta^{2} m2(t)^{2} + \\beta^{2} s2(t) - 2 \\beta m0(t) m1(t) m2(t) + m0(t)^{2} m1(t)^{2} + m0(t)^{2} s1(t) + m1(t)^{2} s0(t) + s0(t) s1(t)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡                                                              2 ⎛     2      \n",
       "⎢                                                             σ ⋅⎝m0(t)  - 2⋅m\n",
       "⎢                                                                             \n",
       "⎢     2      2          2                2  2        2                        \n",
       "⎢m0(t) ⋅m2(t)  - 2⋅m0(t) ⋅m2(t)⋅ρ + m0(t) ⋅ρ  + m0(t) ⋅s2(t) + 2⋅m0(t)⋅m1(t)⋅m\n",
       "⎢                                                                             \n",
       "⎢                                    2      2    2                            \n",
       "⎣                                   β ⋅m2(t)  + β ⋅s2(t) - 2⋅β⋅m0(t)⋅m1(t)⋅m2(\n",
       "\n",
       "                  2                ⎞                                          \n",
       "0(t)⋅m1(t) + m1(t)  + s0(t) + s1(t)⎠                                          \n",
       "                                                                              \n",
       "                              2        2                            2         \n",
       "2(t) - 2⋅m0(t)⋅m1(t)⋅ρ + m1(t)  + m2(t) ⋅s0(t) - 2⋅m2(t)⋅ρ⋅s0(t) + ρ ⋅s0(t) + \n",
       "                                                                              \n",
       "          2      2        2              2                                    \n",
       "t) + m0(t) ⋅m1(t)  + m0(t) ⋅s1(t) + m1(t) ⋅s0(t) + s0(t)⋅s1(t)                \n",
       "\n",
       "                   ⎤\n",
       "                   ⎥\n",
       "                   ⎥\n",
       "                   ⎥\n",
       "s0(t)⋅s2(t) + s1(t)⎥\n",
       "                   ⎥\n",
       "                   ⎥\n",
       "                   ⎦"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# E[f(x)^2]\n",
    "Efx2 = E(fx.multiply_elementwise(fx)).simplify(); Efx2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a6f440b4-c4cf-476c-91e8-dcd0b1044a50",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAABLCAYAAACshBFdAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAOGklEQVR4Ae2dW47dNhKGu41GHgceG8j7ODvwZQXp7CCXFTjZQQI/td8azg6SrCCJd+BkBXa8g8kCBojRmMd5GMz/8bA0FJu6ntM6PFIRkEhRxdvPKhZJkdT51dXV47Ozsz90lczrly9fflF64X6OgCNwXAQkm/9UDh6VcqF35xfJix/lhjg1f6YP7nYEHIGqEHil3NzPcvRMz5/jlwr3D5L29xmhPzoCjkClCEheUcgtI7+v5XFLuFtE/uAIOAJlBCRAaEu05o0ueruf6Pq5NuV4T5ly4wg4AtMQ+F3kbyTM3+lCe17r+lXu4vh3WtSHo3bhPhyWHtMGEJAA0+19JPu1FVduNDjPv5pfDbYLdw214Hk4JQT4evRbIcNv5fdYgp5PcBVIl/Fy4V4GZ09lPQhcdhQF7Y3per97u+A9nS1fMNm7TSq2nkx4fKmr1JL+KJpv7jYXpxW7YzZcX4lW/tBDXc24e3WaWxWQLsqhC8XF93q6UsxqfuKCLRQS45glYPQ7H8TXpqVTahP4h6nnMd2r0tyxZWUm87nc6YQHlcJ3fF+Uk3GbY5YBMu6x1BscF3JBqqqEO2G0seDRgiLItviGrvifem4EO2IZhJr4dZVa3QUhry4px2x8lZh2LoUwrf5X6eUx/GoTbgTvyR5A8JmiNJamq47Qu2DfBtcxu41J0Qf+0cU7E+SU7n58qKZ3uJoxt0C3iYx3KeLR/ZXs0ueLAul2vByzWXUNH5kgpxGYwFfDZ1Vp7hSpGW7rMrVaTjEwWhvB/3RGnGsP4phNr2EWqjCUyQ09zvfit2p6h2vS3IBKq3lpqEfNRGV8URPolr9j2xETx2xCRQgzlpt+kB02ZxBUbjQ5n12f81yLWZPmBlM+e70S2NZF57PEEz1X05qSycqMYza9QtDSLxI+e6bnT/VsE7vTY7yDEKsS7ijEpQm1O4BuHVE6ZtPrMWL23fSQy4ZYTbd8Wdg8NUegfgRcuOuvI8+hIzALARfuWbB5IEegfgRcuOuvI8+hIzALARfuWbB5IEegfgRcuOuvI8+hIzALARfuWbB5IEegfgRcuOuvI8+hIzALARfuWbB5IEegfgRcuHvqSCuROOWSLZGDBjroBwmdYNMILMlTq1p+2sc1UfBsyaAJIedOF9cDR3pOb/msL157JzrOZePsauJs7UwzGrfXhYDqmQ0j7BBj78LgzwlED98txlObEG6rBNlskghGbirlD9mf6SrtwQ27ySJ5yxI9lcrPEwmbCjK7gjjmaZ8DJxS8XqPyshvqJ13/kHvrG3Ko62vhEE7+kR34QnbOF1ahi/LUVrrlMGNrO54qAC0OcwJ4y+gdXfF3slPBTWku9UArbPuhwzvRE99vsr8NHiu5qTwcT0Wv5AcViYMvYOJNG2EBj4z+OUGkX5SntiLcCGPpN8VobBiXAx1Sg+CXNuQbDV31ro35/FrmhRGuwRY+HC/Ennh23L1dQ5kOUAZ6gaUeH/iUfk6wOE9tRbjRwLSygxpHNGjkM9ldWpvXNBaliiUc2pvz2qBxs14EuuqX+sc078ULR+Gpi10+pt0j49JyMYnARnXGHcWJqWkxj6dWegiqzWSTB1pGQEQLc0jDX6L5XvaZ7K4xsGnsVJAZU94SXMVBZaG5SDekIz+69G9lh3TkNkP4rpbdaNyuDAHV4yieinTkvjUsy4oTBDr6HYWnJgu3CkZ3la5sOBRBNkzPxMLfY0Hu3IrgchIGAn0mGyF/o4uZyO910ejQguZCJ6+dEQ35pgKgh9YMXe7SOByBZTxNRV3Khq7LkH7Ap4vA/etCQPWJYI/lqQcx9ynfWIFM4FEwZo7CU/cs9TG2AECIPpedMi5aD2FHWJYyjGmvk8T4kwiCaj8jR+jRnH2GyaHXyndoIBJC4rEKSrwbZxhvN09lB+FhFjeng8Acnhpbx0fhqamaG61dEgaqsFhQCQ/+aPbiewJm5kbP6Y8GstfhkWEAdGZoWNCqwU92/lMCowu23iPY0KeNlNHQKqdxm7/ZpNUbv96HMb4FOLatct5XHg5dB8cuViv9A5RxCk/1Nf6m1dOfExyFpxDuv+n6l67/ttDKHgQeTA2T/JK9ehyf32f+4VHhEJSuMW8pyKBfjDOlIw95o5O+b9wKGz5TyS4JdkNXcigM5acVZghwMibiddA6qK3w+5Yxhk+L1clT0OqCFqHNDTyCSedxdj6Fu+I5NE99pGSQ57N7uv6t6+PoltVp6OaW/trxVfQfVZjO2Ge+EDhUAubWJNjO+/930TJebv0IUH7MoiOwZvq61E8hEn2TltxWmRYem0o/Ch5pJtw9DwHV6Riegge66p6EGx6Re0me+o/SQ56DcGOPMWjuNMMwOSBwDY1vx8Q/mkbpMu43YMkXeWl6DnJ/m7wP8eqZfLJyKNfYCHxqEMoHqUfibo23FRdh04bBSPG7sQe360eAutQ1haeYdA2NfVY6ekj5Goij8NRFlrHiYyw0DNswcvSjgN/I3QhWMYIDeiothJl0aVAY+7ZmrfU+CJbsRriiH2GYQLPFKVQkQsyCg3RWnbJ0dWEJEzSywuBurVDSsxnCtxpCe+F2fQioLufwFHsJ2EdAo5AuP/1SJcz/bnMUnhol3MoshccgyGHMKjcz1KxaWkywyYDMO12AiWCRF4T8a7mZJAur0OS2WXN5BYNgI/SW953v7p53n8MntZQgcdMw8LktfF+XnTYKCVlo0Z+nHqfujvjSGBov/C4/sHsjO8f71Io7h6coI404n89M6T3Tc+nnBEfhqfOrqyu6qwgFf+YoCqr8EZynsrs0ml6vx6icfKfuWvzfW1CFpaJheBo/N45AQGApnlI6KB4U0Pm9kdjT1fh5JO0ayEqf/MaWizE9jaEbRyBFYHGeGhRutQBoIsaXmxlDqsx0M/NZ9LSiim6FAydWr3V114vh3HP9CByDpwaFW7Aj3HwCK3bZV1wtczSwTfStGBYv2h4ILMpTg8ItoWYl1+bGjyozk0VMIIbJs6EKjXTQ5xN0Q0H9/UYQWJqnRk2obQR7L6YjcPIIRCUzaULt5AvtBXAEtobAYLd8a4B4eR2BtSDgwr2WmvRyOAIZAi7cGSD+6AisBYGxy09PqryaVHikDPO9eXBZZJyA4IuAz3L31LJj2gPOwCvDLpItxmur09wRSGYLBwUbsCPdqxhuoJrW8Vpl5eQcMKLcrMvHZhly0egdjaVjWkSn31PYsZ+BzUnw4y+6wLET6/7Ypr09Gc0tQNhe+ZOuocPwOxeSKA5WkLGOPl83ziYPTirZxNr5WNZr4ZHuZrIfNJR6MI6pQJtqhC9rJNDUYQGYbA55YBkqV2s349S4x9BXrbkFBBpm9GH4ogXMRQ9+HwNyTTQRo9ZWVfndKI8IOkLcMo5pC47RD8KN3s6Z7HxlJwomvBsd2UzC2oWblm7KYfgctUSr2GVoLfON9EZ7LccLe1ixzRbZ0j6Bt/Kv4jD9lWDPPu/S0JCTi3KBv5MiVy3cU0osIK2lLHUrLSr2IpcY+0zh0V5b+JmA7cc2TMym/JjmvWO6A+RQd+EJtlyjzvvbN92LfSOYGl4FpFti67XZ3E5BEUwmGR7qan4mIPcUw5j8luBGQFmwb90htBPdz9X8TGAsppEOTD9w6zChkYzvNotpBzajvCPPMdYGS/ib8wHQ2PDg0JyRSA5jFhXuyFwvZIeWSzZCHk6pkHvUzwR6ik2XuzRmROABGkbl81jfRAaVQENwMkblgWHGYvogFuymUEATeBpYM5vE1Aq/hx2O7op1A9ZgymkvnIuAwC/SLV9UuFUoxrSMbc2w24zC2tgEoZ9bcOIxBpXzlgnj7Vu+bQ/CIyynZOZgOraMW8V03/oPDaSEG8EOiowI9Yzy4IvPE57v2iwt3Hx+ocBmLuVAqwY/2eHTjL2caKOV0rjz4KQ1FD/jdRj6lMwUTPsaP9PqRz9M/1jgi//uK20+iWKPMfBb6QcaKYZpPCiuvgnflHZv96LCbUKc5JpxdtOyJf4HdcZKQ2gZAqzKTMEUWl2U3wQ5xcIYum9CsqFfI6YRy720quIICqsBqu2AB2kQFjH3FkmlkIhAQLAxtybBdt6T731d6qfEpjSbtOQ2Zk4TgulHMXcaqBb3SEzBoKvsFKXBSO7NYwogEw0LpNDQJcPQMMW3RHMwv0WFW4WeevD7lIIilCWNRByt8Tb5kB+taG4WbVnzxOc8z8CUScfQ2GXpobHyNQCbxDTDZeojmvuWUT3BW7y7856qJb6YcKtwFAzGssIjcI2JhafbOLfbQmvZ1aVCUwWNrPhxs0Kr1LoSfrGWVWntZVSGyZgqDJOXH2TTwAUjN5gwk/t859PcN4dpU/L5DrC1T70hFj3TS2VIyIKsxXqGF/PLMDkknwKY0JryM4EzgcExwWhkaxS6DsMPn9REVzJMYrBgP4Auu+t0UjRazuCl+Grxm4WpMk8jxucz671UdZh+LeBOzYfwhEfRzJzxD88xsfZQF4b/AsxVXLsYJt5XdYaawONTQ74pZBQkCguj+88EMrQc0wyQnkdhxe66xbrdpawofRQYiux8sW55KSN34EdrORdcFq/QS3DTRsAxbeNxMk+rEm61Vown/WcCB2Q/x3QcmMKJeYuu79vjIjkw1aqEO2IzRwMz0cduKTdlBBzTMi6pL/M1VU3Grk641YIyG8nPAVozlmktpO5IB/1is5hp+qfgdkyHa0kYNYcyDFMvQ7GqCbVlIPNUHIF6EYjKKkyopZ/COGYnzzU/q/fuao6KPzsCFSAg2eTr0KOurCDcdEe7ZphLCz264nJ/R8ARWBYBvmQwkVc0/wP409LuEZq0KAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\sigma & \\sigma & 0\\\\\\rho - x2(t) & -1 & - x0(t)\\\\x1(t) & x0(t) & - \\beta\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡   -σ        σ      0   ⎤\n",
       "⎢                        ⎥\n",
       "⎢ρ - x2(t)   -1    -x0(t)⎥\n",
       "⎢                        ⎥\n",
       "⎣  x1(t)    x0(t)    -β  ⎦"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Jacobian: dfx/dx\n",
    "dfxdx = fx.jacobian(xt); dfxdx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b382278d-741d-4477-b8d7-b22476ffc21b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIcAAAAZCAYAAAAbiz05AAAACXBIWXMAAA7EAAAOxAGVKw4bAAADB0lEQVRoBe2b7VEUQRCGD8sAKI1AyQDIQMkAJAI0A/159/fMAMwAyADIQMkAMhAJwfdZZ5Zlbmdvz9sundvuqr75vJme7ne7e+5jazqdXk0mk7fiSAez2ew+NrwchwZk8/c66WnjtJdbAsedBnYanV4duQaEB5zF/OXI9TC64wfD4yWga7WzUeLFnzn+OgYNCAifdc5dlWcqz8Wnqu/mzu7gyGkm6ZcSD8W/xNvJUBFNyf1RguIpLhFY5aOKeWC6FsjDyoJKnjqkQIDwTfwgfiUuFRjVhUPnudUZmsR5qrFmZ6w7OKImWsrwdB0xpDou+bBlWgldeL2vLYIeqy8FTD3NwVGrYlwVgYWkFN7LndwcHBIC10Vs+yBuc8tnmvMpJ6D3r6eBAAJyDcLHF/GdGI+BLd6o/1FlK5kmpNqYTPhH2Bn3DHN1uhbz2cqO5jgwpAhD4nZC6CBvAgivxd/F5FDZfENjEzPPIYFA5o34RPUqQ2ZD1RGKK1T2fs08p8E0ABjQO8DAc1SkNh6EZHv1sKI3R+NS9iE2BwgxwSGU3KtdAyMsUoGC9cW8x4zYQ4sD0L89g5lsyxYeUPafmb2wEzbKUtZzBMNlUZVd8WmAe3VbyCDUABpTYCDGAGdgmX9CQ8iuNUg4CeFtREjptIFJziGhYiwjtqVEMpQTOJ3r7fU0wJeo0ZOnKx2oo9MOWc+RrrRim+QHepZXSFC8BsB5x6CTuQbwHAsUHl7GOiODlefAXYHKWrgg0IX6jlTvdGea4zSMBh6ka8J7TWrzgPIzDezw7OGtJ4WKledgea6tcwkQQwxZ857aRQFD8vIbB25YEeg3QalXKvkC678kyYa83E72VSfxJDHFBlAvO5iBQwIBgraEtJKulBedo9QzkG8AjlzOsdQEJmFl6a4+oQgNODiKMNNqQspjbOsduc83ei/m4OitqqIm7kvazmtqn9OY5Rx9Nvc5NhqQ51gbGEjmnsPGPhuxqoNjI8xoc4jqrwnJ0lyBOj8cSeZ7cwM0IJvzuUjzfyu3vwHQgNnm/ZpkmQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\sigma & -1 & - \\beta\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[-σ  -1  -β]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# E[dfx/dx]\n",
    "Edfxdx = E(dfxdx.diagonal()); Edfxdx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9b6406d-cba5-47f9-8d08-27da8f5cf7e1",
   "metadata": {},
   "source": [
    "## Put all the $E_{sde}$ terms together.\n",
    "\n",
    "These is given by the expression:\n",
    "\n",
    "$E_{sde}(t) = \\frac{1}{2} \\sum_{i=1}^{D} \\frac{1}{\\sigma_{i}^2} \\bigg( E[f_{i}^2]_{q_t} -2 E[f_{i}]_{q_t} \\dot{m_i}(t) + \\dot{m_i}^2(t) + \\frac{(\\dot{s_i}(t) - \\sigma_{i}^2)^2}{4 s_{i}(t)} + (\\sigma_{i}^2 - \\dot{s_i}(t)) E\\big[\\frac{\\partial f_i}{\\partial x_i}\\big]_{q_t}\\bigg)$\n",
    "\n",
    "where $E[\\cdot]_{q_t}$ is the expected value with respect to the Gaussian process.\n",
    "\n",
    "In essence the following cell computes all the terms inside the parenthesis (one for each dimension)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "dd19f1e2-271c-4d04-9bf1-f3748477113a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimating Esde_0 ...\n",
      "Estimating Esde_1 ...\n",
      "Estimating Esde_2 ...\n",
      "Done!\n"
     ]
    }
   ],
   "source": [
    "# Stores the Esde_i.\n",
    "Esde = []\n",
    "\n",
    "# Get the expression for each dimension.\n",
    "for i in range(dim_D):\n",
    "    \n",
    "    # Display info.\n",
    "    print(f\"Estimating Esde_{i} ...\")\n",
    "    \n",
    "    # Update the list.\n",
    "    Esde.append(Efx2[i] - (2.0 * Efx[i] * dMt[i]) + dMt[i]**2 + (dSt[i] - sig[i])**2/(4.0 * Sfun[i]) + (sig[i] - dSt[i]) * Edfxdx[i])\n",
    "    \n",
    "# _end_for_\n",
    "\n",
    "print(\"Done!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4eafe06-4cde-4e6c-a113-9fe0a3c740ce",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Print the energy functions.\n",
    "for n, En in enumerate(Esde):\n",
    "    print(n, \":\")\n",
    "    sym.pprint(En)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5beccbd7-5641-4c91-8e13-dd031f188207",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing dimension: 0 ...\n",
      "  Original expression ops: 755.\n",
      "  Simplify expression ops: 736.\n",
      "Processing dimension: 1 ...\n",
      "  Original expression ops: 1531.\n",
      "  Simplify expression ops: 1506.\n",
      "Processing dimension: 2 ...\n",
      "  Original expression ops: 1103.\n",
      "  Simplify expression ops: 1084.\n",
      "Done! (Count collect = 3/3)\n"
     ]
    }
   ],
   "source": [
    "# Stores the parameterized Esde.\n",
    "Esde_dk = []\n",
    "\n",
    "# Count the succesful simplifications.\n",
    "count_collect = 0\n",
    "\n",
    "# Process each dimension separately.\n",
    "for i in range(dim_D):\n",
    "    \n",
    "    # Display info.\n",
    "    print(f\"Processing dimension: {i} ...\")\n",
    "    \n",
    "    # Initialize temporary result.\n",
    "    tmp_Ei = deepcopy(Esde[i])\n",
    "    \n",
    "    # Replace the mean functions.\n",
    "    for j, mf_j in enumerate(Mfun):\n",
    "        tmp_Ei = tmp_Ei.replace(mt[j], mf_j)\n",
    "    # _end_for_\n",
    "    \n",
    "    # Replace the variance functions.\n",
    "    for j, sf_j in enumerate(Sfun):\n",
    "        tmp_Ei = tmp_Ei.replace(st[j], sf_j)\n",
    "    # _end_for_\n",
    "    \n",
    "    # Atempt a \"collect\" simplification.\n",
    "    collect_expr = tmp_Ei.collect(Mpar[i]).collect(Spar[i])\n",
    "    \n",
    "    print(f\"  Original expression ops: {sym.count_ops(tmp_Ei)}.\")\n",
    "    print(f\"  Simplify expression ops: {sym.count_ops(collect_expr)}.\")\n",
    "    \n",
    "    # Make sure the expressions are equal.\n",
    "    if tmp_Ei.equals(collect_expr):\n",
    "        \n",
    "        # Append the simplified expression.\n",
    "        Esde_dk.append(collect_expr)\n",
    "        \n",
    "        # Increase the success counter.\n",
    "        count_collect += 1\n",
    "    else:\n",
    "        \n",
    "        # Append the non-simplified expression.\n",
    "        Esde_dk.append(tmp_Ei)\n",
    "    # _end_if_\n",
    "    \n",
    "# _end_for_\n",
    "\n",
    "print(f\"Done! (Count collect = {count_collect}/{dim_D})\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df355019-d3b9-468c-8b1a-33b414ebfcd6",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Print the parameterized energy functions.\n",
    "for l, Ei in enumerate(Esde_dk):\n",
    "    \n",
    "    print(f\"Showing dimension {l}:\\n\")\n",
    "    \n",
    "    sym.pprint(Ei)\n",
    "# _end_for_\n",
    "\n",
    "print(\"Finished!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b962ffa-7d55-4def-8dcb-e8f04d8c191b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Compute the gradients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f31bade0-20f7-410e-a611-65950d208abd",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating gradients for dimension: 0 ...\n",
      "Calculating gradients for dimension: 1 ...\n",
      "Calculating gradients for dimension: 2 ...\n",
      "Done!\n"
     ]
    }
   ],
   "source": [
    "# Define the gradient lists.\n",
    "gradMp = []\n",
    "gradSp = []\n",
    "\n",
    "# Process each dimension separately.\n",
    "for i in range(dim_D):\n",
    "    \n",
    "    # Display info.\n",
    "    print(f\"Calculating gradients for dimension: {i} ...\")\n",
    "    \n",
    "    # Reset the lists for each dimension.\n",
    "    gradMi = []\n",
    "    gradSi = []\n",
    "    \n",
    "    # Derivatives with respect to the mean parameters.\n",
    "    for dim_k in Mpar:\n",
    "        for mk in dim_k:\n",
    "            gradMi.append(Esde_dk[i].diff(mk))\n",
    "        # _end_for_\n",
    "    # _end_for_\n",
    "    \n",
    "    # Derivatives with respect to the variance parameters.\n",
    "    for dim_k in Spar:\n",
    "        for sk in dim_k:\n",
    "            gradSi.append(Esde_dk[i].diff(sk))\n",
    "        # _end_for_\n",
    "    # _end_for_\n",
    "    \n",
    "    # Update the lists.\n",
    "    gradMp.append(gradMi)\n",
    "    gradSp.append(gradSi)\n",
    "    \n",
    "# _end_for_\n",
    "\n",
    "print(\"Done!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af10f8be-b9b8-4ff0-a453-6617ee0eaf69",
   "metadata": {},
   "source": [
    "## Generating lambda functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6db82990-9a47-4dc2-934e-06ec4563a5a2",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating Lambda functions 0 ... \n",
      "Generating Lambda functions 1 ... \n",
      "Generating Lambda functions 2 ... \n",
      "Done!\n"
     ]
    }
   ],
   "source": [
    "# Lists to store the results.\n",
    "func_Esde = []\n",
    "grad_Esde_mt = []\n",
    "grad_Esde_st = []\n",
    "\n",
    "for i in range(dim_D):\n",
    "    \n",
    "    # Display information.\n",
    "    print(f\"Generating Lambda functions {i} ... \")\n",
    "    \n",
    "    func_Esde.append(sym.lambdify([t, *flatten([t_h, t_c, Mpar, Spar, sig, p])],\n",
    "                                  Esde_dk[i], modules=[\"scipy\", \"numpy\"], cse=True))\n",
    "    \n",
    "    grad_Esde_mt.append(sym.lambdify([t, *flatten([t_h, t_c, Mpar, Spar, sig, p])],\n",
    "                                     gradMp[i], modules=[\"scipy\", \"numpy\"], cse=True))\n",
    "\n",
    "    grad_Esde_st.append(sym.lambdify([t, *flatten([t_h, t_c, Mpar, Spar, sig, p])],\n",
    "                                     gradSp[i], modules=[\"scipy\", \"numpy\"], cse=True))\n",
    "# _end_for_\n",
    "\n",
    "# Final message.\n",
    "print(\"Done!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87385143-7df7-482c-b268-39d24c8761be",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Save the functions to files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "081b7962-f42d-4e8a-b101-04577a06abc6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Where to save the functions (parent directory).\n",
    "# NOTE: We save the energy and the gradient (lambdified) functions in separate directories\n",
    "# to improve the readability of the code. In practise they can be anywhere.\n",
    "parent_path = Path(\"path/to/save/the/equations/of/MeanFieldVarGP/src/dynamical_systems\")\n",
    "\n",
    "for i in range(dim_D):\n",
    "    \n",
    "    with open(Path(parent_path/f\"energy_functions/L63_Esde_{i}.sym\"), \"wb\") as sym_file:\n",
    "        \n",
    "        dill.dump(func_Esde[i], sym_file)\n",
    "        \n",
    "    # _end_with_\n",
    "    \n",
    "    with open(Path(parent_path/f\"gradient_functions/dL63_Esde_dM{i}.sym\"), \"wb\") as sym_file:\n",
    "        \n",
    "        dill.dump(grad_Esde_mt[i], sym_file)\n",
    "        \n",
    "    # _end_with_\n",
    "    \n",
    "    with open(Path(parent_path/f\"gradient_functions/dL63_Esde_dS{i}.sym\"), \"wb\") as sym_file:\n",
    "        \n",
    "        dill.dump(grad_Esde_st[i], sym_file)\n",
    "        \n",
    "    # _end_with_\n",
    "    \n",
    "# _end_for_\n",
    "\n",
    "print(\"Done!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22eaea1a-3785-4d3f-96f4-193b1d04c1eb",
   "metadata": {},
   "source": [
    "## End of notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0bdba92-650c-4650-b80a-20751f5fec50",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tensorflow_25",
   "language": "python",
   "name": "tensorflow_25"
  },
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}