{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classroom problems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#using Pkg\n",
    "#Pkg.add(\"StatsBase\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Riemann integrals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have some function $f(t)$ we want to integrate over the range of $t\\in[0, T]$. \n",
    "\n",
    "First, we partition the interval into $\\{0=t_0,t_1,t_2,\\dots, t_n, \\dots, t_{N-1},t_N=T\\}$.\n",
    "\n",
    "Consider\n",
    "$$S(T) = \\int_0^T f(t)dt = \\lim_{N\\to\\infty} \\sum_{n=0}^N f(\\tau_n)\\Delta t_n,\\quad \\Delta t_n = t_{n+1}-t_{n},\\quad \\tau_n\\in[t_n, t_{n+1}].$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check how approximations to the integral of $f(t) = t$ converge to the same value no matter where you fix $\\tau_n$ in the interval $[t_n, t_{n+1}]$. Choose $T=1, \\Delta t_n = 0.2$. \n",
    "\n",
    "*Hint*: you can start by computing the integral with the left boundary $\\tau_n = t_n$, and then compare it with the integral using the midpoint $\\tau_n = \\frac{1}{2}(t_n+t_{n+1})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Δt = ####\n",
    "T = ####\n",
    "tt = ####\n",
    "f = function(t) return #### end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# S_left: S(T) approximated with the left point \\tau_n=t_n\n",
    "\n",
    "S_left = zeros(length(tt))\n",
    "S_left[1, :] .= 0\n",
    "for n in 1:length(tt)-1\n",
    "    S_left[n+1] = S_left[####] + f(####)*Δt\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# S_mid: S(T) approximated with the midpoint \\tau_n=0.5*(t_n+t_{n+1})\n",
    "\n",
    "S_mid = zeros(length(tt))\n",
    "S_mid[1, :] .= 0\n",
    "for n in 1:length(tt)-1\n",
    "    S_mid[n+1] = S_mid[####] + f(####)*Δt\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "plot(tt, S_left, label = \"Left point\", color=:black, linestyle=:dot, legend=:topleft, size=(600,300), xlabel=\"t\", ylabel=\"S(T)\", title = \"Time discretisation: \".*\"$(Δt)\")\n",
    "plot!(tt, S_mid, label = \"Mid point\", color=:black)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now increase $N$, that is, choose smaller $\\Delta t_n, n=0,1,\\dots, N$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Δt = ####\n",
    "tt = ####\n",
    "\n",
    "S_left = zeros(length(tt))\n",
    "S_left[1, :] .= 0\n",
    "for n in 1:length(tt)-1\n",
    "    S_left[n+1] = ####\n",
    "end\n",
    "\n",
    "S_mid = zeros(length(tt))\n",
    "S_mid[1, :] .= 0\n",
    "for n in 1:length(tt)-1\n",
    "    S_mid[n+1] = ####\n",
    "end\n",
    "\n",
    "using Plots\n",
    "plot(tt, S_left, label = \"Left point\", color=:black, linestyle=:dot, legend=:topleft, size=(600,300), xlabel=\"t\", ylabel=\"S(T)\", title = \"Time discretisation: \".*\"$(Δt)\")\n",
    "plot!(tt, S_mid, label = \"Mid point\", color=:black)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Stochastic integrals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$dX(t) = f(X(t), t)dt + g(X(t), t)dW, \\quad X:[0,\\infty)\\to \\mathbb{R}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $\\{0=t_0,t_1,t_2,\\dots, t_n, \\dots, t_{N-1},t_N=T\\}$, choose a point $\\tau_n \\in [t_n, t_{n+1})$ and solve the stochastic differential equation according to\n",
    "\n",
    "$$X(t_n+\\Delta t_n) = X(t_n)+f(X(\\tau_n), \\tau_n)\\Delta t_n + g(X(\\tau_n), \\tau_n)\\eta_{t_n}, \\quad X(0) = x_0,\\quad\\eta_{t_n}\\sim\\mathcal{N(0, \\Delta t_n)}\\quad\\text{ and }\\quad n=1,2,\\dots,N,.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now consider the case where $f(X(t), t) = 0$, $g(X(t), t) = 1$ and $X_0 = 0$. Here, the time interval is $t\\in[0, 50]$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Itô integral ($\\tau_n =t_n$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the expectation of following sum\n",
    "\n",
    "$S_N^{ito} = \\sum_{n=1}^N W_{t_n}(W_{t_{n+1}}-W_{t_n}),\\quad t_N=T=50$,\n",
    "\n",
    "which approximates the stochastic integral\n",
    "$S(T) = \\int_0^T W_{t}dW_t$. \n",
    "\n",
    "Hint: it is enough to run 200 realisations of your stochastic integral in order to get an estimate for the expectation. You are free to choose the time discretisation, but you can use for example $\\Delta t = 0.5$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Δt = \n",
    "T = \n",
    "tt = \n",
    "\n",
    "x = [0.0]\n",
    "Nx = length(x) \n",
    "\n",
    "xx = zeros(length(tt), Nx)\n",
    "xx[1, :] .= x\n",
    "\n",
    "f = function(x, t) return [####] end\n",
    "g = function(x, t) return ones(####,####) end "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ns = \n",
    "Exp_S_ito = zeros(size(x))\n",
    "\n",
    "for ns in 1:Ns\n",
    "\n",
    "    ww = zeros(length(tt), Nx)\n",
    "    ww[1, :] .= 0\n",
    "    for n in 1:length(tt)-1\n",
    "        ww[n+1, :] .= ww[n, :] .+ f(####)*Δt .+ g(####)*(sqrt(Δt)*randn(Nx, 1))[]\n",
    "    end\n",
    "    \n",
    "    xx = zeros(length(tt), Nx)\n",
    "    xx[1, :] .= x\n",
    "    for n in 1:length(tt)-1\n",
    "        xx[n+1, :] .= xx[n, :] .+ ####.*(ww[n+1,:]-ww[n,:])\n",
    "    end\n",
    "    \n",
    "    Exp_S_ito += xx[end, :]/Ns \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Exp_S_ito"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratonovich integral ($\\tau_n =\\frac{1}{2}(t_n+t_{n+1})$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the following sum\n",
    "\n",
    "$S_N^{str} = \\sum_{n=1}^N W_{\\frac{1}{2}(t_n+t_{n-1})}(W_{t_n}-W_{t_{n-1}})$,\n",
    "\n",
    "which approximates the stochastic integral\n",
    "$S_t = \\int_0^T W_{t}dW_t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Exp_S_strato = zeros(size(x))\n",
    "\n",
    "for ns in 1:Ns\n",
    "  \n",
    "    ww = zeros(length(tt), Nx)\n",
    "    ww[1, :] .= 0\n",
    "    for n in 1:length(tt)-1\n",
    "        ww[n+1, :] .= ww[n, :] .+ f(####)*Δt .+ g(####)*(sqrt(Δt)*randn(Nx, 1))[]\n",
    "    end\n",
    "    \n",
    "    xx = zeros(length(tt), Nx)\n",
    "    xx[1, :] .= x\n",
    "    for n in 1:length(tt)-1\n",
    "        xx[n+1, :] .= xx[n, :] .+ ####.*(ww[n+1,:]-ww[n,:])\n",
    "    end\n",
    "\n",
    "    Exp_S_strato += xx[end, :]/Ns\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Exp_S_strato"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider now a smaller value of $\\Delta t$. Recompute the expectation of the stochastic integral.\n",
    "\n",
    "*Hint*: Decrease the number of realisations to 25, so that you can plot each realisation by the end of your estimation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ns = ####\n",
    "Δt = ####\n",
    "tt = ####\n",
    "\n",
    "exp_plot = plot(legend=:topleft, size=(600,300), ylims=(-25,100), title = \"Estimation with \".*\"$(Ns)\".*\" realisations\" )\n",
    "\n",
    "Exp_S_ito = zeros(size(x))\n",
    "for ns in 1:Ns\n",
    "    ww = zeros(length(tt), Nx)\n",
    "    ww[1, :] .= 0\n",
    "    for n in 1:length(tt)-1\n",
    "        ww[n+1, :] .= #### \n",
    "    end\n",
    "    \n",
    "    xx = zeros(length(tt), Nx)\n",
    "    xx[1, :] .= x\n",
    "    for n in 1:length(tt)-1\n",
    "        xx[n+1, :] .= ####\n",
    "    end\n",
    "    plot!(tt, xx, colour=:black, linestyle=:dot, label=false)\n",
    "\n",
    "    Exp_S_ito += xx[end, :]/Ns \n",
    "end\n",
    "    \n",
    "Exp_S_strato = zeros(size(x))    \n",
    "for ns in 1:Ns\n",
    "  \n",
    "    ww = zeros(length(tt), Nx)\n",
    "    ww[1, :] .= 0\n",
    "    for n in 1:length(tt)-1\n",
    "         ww[n+1, :] .= ####\n",
    "    end\n",
    "        \n",
    "    xx = zeros(length(tt), Nx)\n",
    "    xx[1, :] .= x\n",
    "    for n in 1:length(tt)-1\n",
    "        xx[n+1, :] .= ####\n",
    "    end\n",
    "    plot!(tt, xx, colour=:darkorange,label=false)\n",
    "\n",
    "    Exp_S_strato += xx[end, :]/Ns\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter!([tt[end]],[Exp_S_ito], markersize=6, markercolour=:black, label=\"E[S(T)] with Ito integral\")\n",
    "scatter!([tt[end]],[Exp_S_strato], markersize=6, markercolour=:darkorange, label=\"E[S(T)] with Stratonovich integral\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. General code for solving Itô integrals with the Euler-Maruyama method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "de_EM_solver = function(f, g, x₀, p, T, Δt)\n",
    "\n",
    "    \"\"\"\n",
    "    f = f(x, p, t) and g = g(x, p, t), where x is the process variable, p stands for a vector of parameters and t is time.    \n",
    "    Outputs:\n",
    "        tt - a partitioned time interval [0, T]\n",
    "        xx - a solution to the stochastic integral dX_t = f(X_t, p, t)dt + g(X_t, p, t)dW\n",
    "    \"\"\"\n",
    "    \n",
    "    x = copy(x₀)\n",
    "    t = 0.0\n",
    "\n",
    "    Nx = length(x) \n",
    "    tt = 0.0:Δt:T\n",
    "    xx = zeros(length(tt), Nx)\n",
    "\n",
    "    xx[1, :] .= x\n",
    "    \n",
    "    Nw = size(g(x, p, 0))\n",
    "    \n",
    "    for n in 1:length(tt)-1\n",
    "        t⁻ = tt[n]\n",
    "        x⁻ = xx[n, :]\n",
    "        xx[n+1, :] = x⁻ + f(x⁻, p, t⁻).*Δt + sqrt(Δt)*g(x⁻, p, t⁻)*randn(Nw)\n",
    "    end\n",
    "    \n",
    "    return tt, xx\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 20\n",
    "Δt = ####\n",
    "Nₛ = ####;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (3.1)   $X(t+dt) = X(t) + dW, \\quad X:[0,\\infty)\\to \\mathbb{R}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "\n",
    "x₀ = [0.0]\n",
    "p = []\n",
    "f1(x, p, t) = ####  # In vector form\n",
    "g1(x, p, t) = ####  # In matrix form\n",
    "\n",
    "p = plot(legend = false, xlabel = \"t\", ylabel = \"X\")\n",
    "for nₛ in 1:Nₛ    \n",
    "    tt, xx = de_EM_solver(f1, g1, x₀, p, T, Δt)\n",
    "    plot!(tt, xx, linetype=:steppre)\n",
    "end\n",
    "p\n",
    "plot!(0:Δt:T, t -> sqrt.(t), colour=:black, linewidth=2)  # +1 standard deviation\n",
    "plot!(0:Δt:T, t -> -sqrt.(t), colour=:black, linewidth=2) # -1 standard deviation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (3.2) $X(t+dt) = X(t) + dt + dW, \\quad X:[0,\\infty)\\to \\mathbb{R}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f2(x, p, t) = #### # In vector form\n",
    "g2(x, p, t) = #### # In matrix form\n",
    "\n",
    "p = plot(legend = false, xlabel = \"t\", ylabel = \"X\")\n",
    "for nₛ in 1:Nₛ    \n",
    "    tt, xx = de_EM_solver(f2, g2, x₀, p, T, Δt)\n",
    "    plot!(tt, xx, linetype=:steppre)\n",
    "end\n",
    "p\n",
    "plot!(0:Δt:T, t -> t+sqrt.(t), colour=:black, linewidth=2) # +1 standard deviation\n",
    "plot!(0:Δt:T, t -> t-sqrt.(t), colour=:black, linewidth=2) # -1 standard deviation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (3.3) $X(t+dt) = X(t) +  \\begin{bmatrix}\\sqrt{2D} & 0 \\\\ 0 & \\sqrt{2D} \\end{bmatrix}dW, \\quad X:[0,\\infty)\\to \\mathbb{R}^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x₀ = zeros(2)\n",
    "D = 1.0\n",
    "p = [D]\n",
    "\n",
    "f3(x, p, t) = #### # In vector form\n",
    "g3 = function(x, p, t)\n",
    "    D = p[1]\n",
    "    return ####  # In matrix form\n",
    "end\n",
    "\n",
    "plots = plot(legend = false, xlims = (-15,15), ylims=(-15,15),xlabel = \"X1\", ylabel = \"X2\")\n",
    "for nₛ in 1:Nₛ    \n",
    "    tt, xx = de_EM_solver(f3, g3, x₀, p, T, Δt)\n",
    "    plot!(xx[:, 1], xx[:, 2], linetype=:steppre)\n",
    "    scatter!(xx[end:end, 1], xx[end:end, 2], markersize = 7, markercolour = :black)\n",
    "end\n",
    "plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. System with multiple favourable states (Homework from week 3 revisited)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\frac{dx}{dt} = -κ₁*x^3+κ₂*x^2-κ₃*x+κ₄,\\quad x_0 = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stoichiometric matrix\n",
    "Products = [2; \n",
    "            3; \n",
    "            0; \n",
    "            1]\n",
    "Reactants = [3; \n",
    "             2; \n",
    "             1; \n",
    "             0]\n",
    "SM = (Products-Reactants)'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Δt = 0.001\n",
    "κ₁ = 0.0015\n",
    "κ₂ = 0.36\n",
    "κ₃ = 37.5\n",
    "κ₄ = 2200\n",
    "V = 0.5\n",
    "\n",
    "p = [κ₁,κ₂,κ₃,κ₄, V]\n",
    "x₀ = zeros(1)\n",
    "T = 40\n",
    "\n",
    "fode = function(x, p, t) \n",
    "    κ₁,κ₂,κ₃,κ₄,V = p\n",
    "    return ####\n",
    "end\n",
    "gode = function(x, p, t)\n",
    "    Nx = length(x)\n",
    "    Nw = length(p)-1\n",
    "    return zeros(####, ####)\n",
    "end\n",
    "\n",
    "tt, xx_ode = de_EM_solver(fode, gode, x₀, p, T, Δt);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(tt, xx_ode, colour=:black, ylabel=\"X(t)\", xlabel=\"t\", label = \"ODE\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fsde = function(x, p, t) \n",
    "    κ₁,κ₂,κ₃,κ₄,V = p\n",
    "    return ####\n",
    "end\n",
    "gsde = function(x, p, t)\n",
    "    κ₁,κ₂,κ₃,κ₄,V = p\n",
    "    Nw = length(p)-1\n",
    "    return SM*sqrt(diagm([####; ####; ####; ####]))\n",
    "end\n",
    "tt, xx_sde = de_EM_solver(fode, gsde, x₀, p, T, Δt);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(tt, xx_sde, label=\"SDE\", xlabel=\"t\", ylabel=\"X(t)\")\n",
    "plot!(tt, xx_ode, colour=:black, linewidth=5, label=\"ODE\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, change the volume $V$ to $20V$. What is the effect of chaning the volume $V$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = ####\n",
    "p = [κ₁,κ₂,κ₃,κ₄,κ₅, V]\n",
    "\n",
    "tt, xx_ode = de_EM_solver(fode, gode, x₀, p, T, Δt)\n",
    "tt, xx_sde = de_EM_solver(fsde, gsde, x₀, p, T, Δt);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(tt, xx_sde, label=\"SDE\", xlabel=\"t\", ylabel=\"X(t)\")\n",
    "plot!(tt, xx_ode, colour=:black, linewidth=5, label=\"ODE\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "plot(tt, xx_sde, label=\"SDE\", xlabel=\"t\", ylabel=\"X(t)\")\n",
    "plot!(tt, xx_ode, colour=:black, linewidth=5, label=\"ODE\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1e3                     #number of realisations\n",
    "n_upperbound_A = 100\n",
    "n_upperbound_B = 150\n",
    "pp_sde = zeros(n_upperbound_A+1, n_upperbound_B+1)   #vector of probabilities\n",
    "\n",
    "for n in 1:N   \n",
    "     tt, xx = de_EM_solver(fsde, gsde, x₀, p, T, Δt)\n",
    "  \n",
    "    nA, nB = Integer.(round.(xx[end, :]))    \n",
    "\n",
    "     if (nA <= n_upperbound_A && nB <= n_upperbound_B) \n",
    "        pp_sde[nA+1, nB+1] += 1\n",
    "     end\n",
    "end\n",
    "\n",
    "pp_sde /= N;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap(0:n_upperbound_B, 0:n_upperbound_A, pp_sde, c=:coolwarm,\n",
    "    xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"Empirical density with the \\nSDEs @ t = 20sec\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "\n",
    "a = heatmap(0:n_upperbound_B, 0:n_upperbound_A, pp,\n",
    "    c=:coolwarm,\n",
    "    xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"Solution to the \\nChemical Master equation @ t = 20sec\")\n",
    "\n",
    "b = heatmap(0:n_upperbound_B, 0:n_upperbound_A, pp_sde, c=:coolwarm,\n",
    "    xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"Empirical density with the \\nGillespie algorithm @ t = 20sec\")\n",
    "\n",
    "layout = @layout [a b]\n",
    "Plots.plot(a, b; layout = layout,  colorbar = false, size = (800, 400),  margin = 30Plots.px)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.2",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}