From 58be0b6f8c47f4942a22ff5be3daf3212a001e94 Mon Sep 17 00:00:00 2001
From: Andreagiovanni Reina <areina@ulb.be>
Date: Fri, 18 Mar 2022 18:41:53 +0100
Subject: [PATCH] Updated the demo notebook Variance Suppression with new
 analysis of the stochastic noise.

---
 DemoNotebooks/Variance_suppression.ipynb | 159 +++++++++++++++++++++++
 1 file changed, 159 insertions(+)

diff --git a/DemoNotebooks/Variance_suppression.ipynb b/DemoNotebooks/Variance_suppression.ipynb
index fb430ab..d27348f 100644
--- a/DemoNotebooks/Variance_suppression.ipynb
+++ b/DemoNotebooks/Variance_suppression.ipynb
@@ -259,6 +259,165 @@
     "                          choose_yrange=[0,1], plotProportions=True, ylab=\"Subpopulations $x_i$\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Fokker-Planck equations and noise variance\n",
+    "\n",
+    "MuMoT allows you to easily compute the Fokker-Planck equations for both models and display the first-order and second-order moments of the system's noise."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "posmodel.showFokkerPlanckEquation()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "negmodel.showFokkerPlanckEquation()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Noise equations\n",
+    "\n",
+    "Having the Fokker-Planck equations allows us to derive the equations of motion for the fluctuations.\n",
+    "\n",
+    "### Model without negative social feedback"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "posmodel.showNoiseEquations()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Because MuMoT allows you to access the `sympy` equations, we apply some simplification to the equations in order to find the solution to $\\langle \\eta_{A}^2 \\rangle$ at convergence ($t\\rightarrow \\infty$). We assume that $a\\approx0$ and $\\Phi_U (t\\rightarrow \\infty) \\approx 0$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "noiseEqPos=posmodel.getNoiseEquations()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#in order to handle sympy equations, some imports are necessary\n",
+    "from sympy import (\n",
+    "    collect,\n",
+    "    default_sort_key,\n",
+    "    Derivative,\n",
+    "    lambdify,\n",
+    "    latex,\n",
+    "    linsolve,\n",
+    "    numbered_symbols,\n",
+    "    preview,\n",
+    "    simplify,\n",
+    "    solve,\n",
+    "    Symbol,\n",
+    "    symbols,\n",
+    "    Function\n",
+    ")\n",
+    "from IPython.display import display, Math"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "simplifiedNoiseEqPos={}\n",
+    "\n",
+    "simplifiedNoiseEqPos['etaU2']=noiseEqPos[2][Function('M_2')(Symbol('eta_U')**2)].subs(Symbol('Phi_U'),0).subs(Symbol('a'),0)\n",
+    "display(Math(latex(simplifiedNoiseEqPos['etaU2'])+'=0'))\n",
+    "simplifiedNoiseEqPos['etaU2']=solve( simplifiedNoiseEqPos['etaU2'], Function('M_2')(Symbol('eta_U')**2) )[0]\n",
+    "display(Math(latex( Function('M_2')(Symbol('eta_U')**2) ) + '=' + latex(simplifiedNoiseEqPos['etaU2']) ))\n",
+    "\n",
+    "simplifiedNoiseEqPos['etaA*etaU']=noiseEqPos[2][Function('M_2')(Symbol('eta_A')*Symbol('eta_U'))].subs(\n",
+    "    Function('M_2')(Symbol('eta_U')**2),simplifiedNoiseEqPos['etaU2']).subs(Symbol('Phi_U'),0).subs(Symbol('a'),0)\n",
+    "display(Math(latex(simplifiedNoiseEqPos['etaA*etaU'])+'=0'))\n",
+    "simplifiedNoiseEqPos['etaA*etaU']=solve( simplifiedNoiseEqPos['etaA*etaU'], Function('M_2')(Symbol('eta_A')*Symbol('eta_U')) )[0]\n",
+    "display(Math(latex( Function('M_2')(Symbol('eta_A')*Symbol('eta_U')) ) + '=' + latex(simplifiedNoiseEqPos['etaA*etaU']) ))\n",
+    "\n",
+    "simplifiedNoiseEqPos['etaA2']=noiseEqPos[2][Function('M_2')(Symbol('eta_A')**2)].subs(\n",
+    "    Function('M_2')(Symbol('eta_A')*Symbol('eta_U')),simplifiedNoiseEqPos['etaA*etaU']).subs(Symbol('Phi_U'),0)\n",
+    "display(Math(latex(simplifiedNoiseEqPos['etaA2'])+'=0'))\n",
+    "simplifiedNoiseEqPos['etaA2']=solve( simplifiedNoiseEqPos['etaA2'], Function('M_2')(Symbol('eta_A')**2) )[0]\n",
+    "display(Math(latex( Function('M_2')(Symbol('eta_A')**2) ) + '=' + latex(simplifiedNoiseEqPos['etaA2']) ))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Model with negative social feedback"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "negmodel.showNoiseEquations()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "noiseEqNeg=negmodel.getNoiseEquations()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the system with negative feedback, we assume that $\\Phi_U (t\\rightarrow \\infty) \\approx 0$ and $\\langle \\eta_{A}\\eta_{U} \\rangle \\approx 0$ to find the solution to $\\langle \\eta_{A}^2 \\rangle$:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "simplifiedNoiseEqPos=noiseEqPos[2][Function('M_2')(Symbol('eta_A')**2)].subs(Symbol('Phi_U'),0).subs(\n",
+    "    Function('M_2')(Symbol('eta_A')*Symbol('eta_U')),0)\n",
+    "display(Math(latex(simplifiedNoiseEqPos)+'=0'))\n",
+    "simplifiedNoiseSizePos=solve( simplifiedNoiseEqPos, Function('M_2')(Symbol('eta_A')**2) )[0]\n",
+    "display(Math(latex( Function('M_2')(Symbol('eta_A')**2) ) + '=' + latex(simplifiedNoiseSizePos) ))"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},