[6]_ODE_PINN_finite_difference.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/ASEM000/Physics-informed-neural-network-in-JAX/blob/main/%5B6%5D_ODE_PINN_finite_difference.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "v77fdC1ZLyg1"
      },
      "outputs": [],
      "source": [
        "#Credits : Mahmoud Asem @Asem000"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "vAR0swbLX_ZI",
        "outputId": "dab17b4c-908d-4213-ac91-c20d0f7ef122"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Requirement already satisfied: optax in /usr/local/lib/python3.7/dist-packages (0.1.1)\n",
            "Requirement already satisfied: numpy>=1.18.0 in /usr/local/lib/python3.7/dist-packages (from optax) (1.21.5)\n",
            "Requirement already satisfied: jaxlib>=0.1.37 in /usr/local/lib/python3.7/dist-packages (from optax) (0.3.2+cuda11.cudnn805)\n",
            "Requirement already satisfied: absl-py>=0.7.1 in /usr/local/lib/python3.7/dist-packages (from optax) (1.0.0)\n",
            "Requirement already satisfied: typing-extensions>=3.10.0 in /usr/local/lib/python3.7/dist-packages (from optax) (3.10.0.2)\n",
            "Requirement already satisfied: jax>=0.1.55 in /usr/local/lib/python3.7/dist-packages (from optax) (0.3.4)\n",
            "Requirement already satisfied: chex>=0.0.4 in /usr/local/lib/python3.7/dist-packages (from optax) (0.1.1)\n",
            "Requirement already satisfied: six in /usr/local/lib/python3.7/dist-packages (from absl-py>=0.7.1->optax) (1.15.0)\n",
            "Requirement already satisfied: toolz>=0.9.0 in /usr/local/lib/python3.7/dist-packages (from chex>=0.0.4->optax) (0.11.2)\n",
            "Requirement already satisfied: dm-tree>=0.1.5 in /usr/local/lib/python3.7/dist-packages (from chex>=0.0.4->optax) (0.1.6)\n",
            "Requirement already satisfied: scipy>=1.2.1 in /usr/local/lib/python3.7/dist-packages (from jax>=0.1.55->optax) (1.4.1)\n",
            "Requirement already satisfied: opt-einsum in /usr/local/lib/python3.7/dist-packages (from jax>=0.1.55->optax) (3.3.0)\n",
            "Requirement already satisfied: flatbuffers<3.0,>=1.12 in /usr/local/lib/python3.7/dist-packages (from jaxlib>=0.1.37->optax) (2.0)\n"
          ]
        }
      ],
      "source": [
        "#Imports\n",
        "import jax \n",
        "import jax.numpy as jnp\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from matplotlib import cm\n",
        "import matplotlib as mpl\n",
        "!pip install optax\n",
        "import optax"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "yoPHsh5lWvyP"
      },
      "outputs": [],
      "source": [
        "import sympy as sp"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7bg4nSbsXVwD"
      },
      "source": [
        "### Generate a a differential equation and its solution using SymPy"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "P9664e-mVMTN"
      },
      "outputs": [],
      "source": [
        "t= sp.symbols('t')\n",
        "f = sp.Function('y')\n",
        "diffeq = sp.Eq(f(t).diff(t,t) + f(t).diff(t)-t*sp.cos(2*sp.pi*t),0)\n",
        "sol = sp.simplify(sp.dsolve(diffeq,ics={f(0):1,f(t).diff(t).subs(t,0):10}).rhs)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 54
        },
        "id": "klgFeU6bcTrC",
        "outputId": "5c0fc009-ce8b-472f-b98b-6da8afd6a66c"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Eq(-t*cos(2*pi*t) + Derivative(y(t), t) + Derivative(y(t), (t, 2)), 0)"
            ],
            "text/latex": "$\\displaystyle - t \\cos{\\left(2 \\pi t \\right)} + \\frac{d}{d t} y{\\left(t \\right)} + \\frac{d^{2}}{d t^{2}} y{\\left(t \\right)} = 0$"
          },
          "metadata": {},
          "execution_count": 4
        }
      ],
      "source": [
        "diffeq"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 60
        },
        "id": "E4Uu2hbiYJtv",
        "outputId": "6444bdb2-9243-46c9-e801-d19960dcb233"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Eq(Subs(Derivative(y(t), t), t, 0), 10)"
            ],
            "text/latex": "$\\displaystyle \\left. \\frac{d}{d t} y{\\left(t \\right)} \\right|_{\\substack{ t=0 }} = 10$"
          },
          "metadata": {},
          "execution_count": 5
        }
      ],
      "source": [
        "sp.Eq(f(t).diff(t).subs(t,0),10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 38
        },
        "id": "29QUbt_2YwlJ",
        "outputId": "0e3cce9d-29c7-4f26-c279-535093e2ad37"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Eq(y(0), 1)"
            ],
            "text/latex": "$\\displaystyle y{\\left(0 \\right)} = 1$"
          },
          "metadata": {},
          "execution_count": 6
        }
      ],
      "source": [
        "sp.Eq(f(t).subs(t,0),1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 81
        },
        "id": "r9KVq1yjYfld",
        "outputId": "9a8437fc-c2fc-4fe7-ad9f-26ca25b974d4"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Eq(y(t), (2*pi*t*exp(t)*sin(2*pi*t) + 8*pi**3*t*exp(t)*sin(2*pi*t) - 16*pi**4*t*exp(t)*cos(2*pi*t) - 4*pi**2*t*exp(t)*cos(2*pi*t) + 16*pi**3*exp(t)*sin(2*pi*t) + exp(t)*cos(2*pi*t) + 12*pi**2*exp(t)*cos(2*pi*t) - exp(t) + 36*pi**2*exp(t) + 336*pi**4*exp(t) + 704*pi**6*exp(t) - 640*pi**6 - 304*pi**4 - 44*pi**2)*exp(-t)/(4*pi**2*(1 + 8*pi**2 + 16*pi**4)))"
            ],
            "text/latex": "$\\displaystyle y{\\left(t \\right)} = \\frac{\\left(2 \\pi t e^{t} \\sin{\\left(2 \\pi t \\right)} + 8 \\pi^{3} t e^{t} \\sin{\\left(2 \\pi t \\right)} - 16 \\pi^{4} t e^{t} \\cos{\\left(2 \\pi t \\right)} - 4 \\pi^{2} t e^{t} \\cos{\\left(2 \\pi t \\right)} + 16 \\pi^{3} e^{t} \\sin{\\left(2 \\pi t \\right)} + e^{t} \\cos{\\left(2 \\pi t \\right)} + 12 \\pi^{2} e^{t} \\cos{\\left(2 \\pi t \\right)} - e^{t} + 36 \\pi^{2} e^{t} + 336 \\pi^{4} e^{t} + 704 \\pi^{6} e^{t} - 640 \\pi^{6} - 304 \\pi^{4} - 44 \\pi^{2}\\right) e^{- t}}{4 \\pi^{2} \\left(1 + 8 \\pi^{2} + 16 \\pi^{4}\\right)}$"
          },
          "metadata": {},
          "execution_count": 7
        }
      ],
      "source": [
        "sp.Eq(f(t),sol)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 37
        },
        "id": "MNVOpPyCW-GU",
        "outputId": "53a70887-e440-43a8-f811-30c862b21866"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0"
            ],
            "text/latex": "$\\displaystyle 0$"
          },
          "metadata": {},
          "execution_count": 8
        }
      ],
      "source": [
        "#verify solution\n",
        "sp.simplify(-t*sp.cos(sp.pi*2*t)+sol.diff(t)+sol.diff(t,t))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NQ61lEQeXgrc"
      },
      "source": [
        "### Constructing the MLP"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Lml6PGLPZgmr"
      },
      "outputs": [],
      "source": [
        "N_b = 1\n",
        "N_c = 100\n",
        "\n",
        "tmin,tmax=0. ,jnp.pi\n",
        "\n",
        "'''boundary conditions'''\n",
        "\n",
        "\n",
        "# U[0] = 1\n",
        "t_0 = jnp.ones([N_b,1],dtype='float32')*0.\n",
        "ic_0 = jnp.ones_like(t_0) \n",
        "IC_0 = jnp.concatenate([t_0,ic_0],axis=1)\n",
        "\n",
        "# U_t[0] = 10\n",
        "t_b1 = jnp.zeros([N_b,1])\n",
        "bc_1 = jnp.ones_like(t_b1) * 10\n",
        "BC_1 = jnp.concatenate([t_b1,bc_1],axis=1)\n",
        "\n",
        "conds = [IC_0,BC_1]\n",
        "\n",
        "#collocation points\n",
        "\n",
        "key=jax.random.PRNGKey(0)\n",
        "\n",
        "t_c = jnp.linspace(tmin,tmax,N_c).reshape(-1,1)\n",
        "colloc = t_c\n",
        "\n",
        "def ODE_loss(t,u):\n",
        "  dt = 0.03173326\n",
        "  u_t = lambda t: (-u(t+2*dt)+8*u(t+dt)-8*u(t-dt)+u(t-2*dt))/(12*dt)\n",
        "  u_tt = lambda t: (-u(t+2*dt) + 16*u(t+dt) -30*u(t) + 16 * u(t-dt) - u(t-2*dt))/(12*dt**2)\n",
        "  return -t*jnp.cos(2*jnp.pi*t) + u_t(t) + u_tt(t)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KoZZJl2TbI_n"
      },
      "outputs": [],
      "source": [
        "def init_params(layers):\n",
        "  keys = jax.random.split(jax.random.PRNGKey(0),len(layers)-1)\n",
        "  params = list()\n",
        "  for key,n_in,n_out in zip(keys,layers[:-1],layers[1:]):\n",
        "    lb, ub = -(1 / jnp.sqrt(n_in)), (1 / jnp.sqrt(n_in)) # xavier initialization lower and upper bound\n",
        "    W = lb + (ub-lb) * jax.random.uniform(key,shape=(n_in,n_out))\n",
        "    B = jax.random.uniform(key,shape=(n_out,))\n",
        "    params.append({'W':W,'B':B})\n",
        "  return params\n",
        "\n",
        "def fwd(params,t):\n",
        "  X = jnp.concatenate([t],axis=1)\n",
        "  *hidden,last = params\n",
        "  for layer in hidden :\n",
        "    X = jax.nn.tanh(X@layer['W']+layer['B'])\n",
        "  return X@last['W'] + last['B']\n",
        "\n",
        "@jax.jit\n",
        "def MSE(true,pred):\n",
        "  return jnp.mean((true-pred)**2)\n",
        "\n",
        "def loss_fun(params,colloc,conds):\n",
        "  t_c =colloc[:,[0]]\n",
        "  ufunc = lambda t : fwd(params,t)\n",
        "  ufunc_t=lambda t:jax.grad(lambda t:jnp.sum(ufunc(t)))(t)\n",
        "  loss =jnp.mean(ODE_loss(t_c,ufunc) **2)\n",
        "\n",
        "  t_ic,u_ic = conds[0][:,[0]],conds[0][:,[1]]  \n",
        "  loss += MSE(u_ic,ufunc(t_ic))\n",
        "\n",
        "  t_bc,u_bc = conds[1][:,[0]],conds[1][:,[1]]  \n",
        "  loss += MSE(u_bc,ufunc_t(t_bc))\n",
        "\n",
        "  return  loss\n",
        "\n",
        "@jax.jit\n",
        "def update(opt_state,params,colloc,conds):\n",
        "  # Get the gradient w.r.t to MLP params\n",
        "  grads=jax.jit(jax.grad(loss_fun,0))(params,colloc,conds)\n",
        "  \n",
        "  #Update params\n",
        "  updates, opt_state = optimizer.update(grads, opt_state)\n",
        "  params = optax.apply_updates(params, updates)\n",
        "\n",
        "  #Update params\n",
        "  # return jax.tree_multimap(lambda params,grads : params-LR*grads, params,grads)\n",
        "  return opt_state,params\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ae1ZDoy0c29c"
      },
      "outputs": [],
      "source": [
        "# construct the MLP of  6 hidden layers of 8 neurons for each layer\n",
        "params = init_params([1] + [20]*4+[1])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jySmbUwic5yk"
      },
      "outputs": [],
      "source": [
        "lr = optax.piecewise_constant_schedule(1e-3,{10_000:5e-3,30_000:1e-3,50_000:5e-4,70_000:1e-4})\n",
        "optimizer = optax.adam(1e-3)\n",
        "opt_state = optimizer.init(params)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kBzGA8OVc8C6",
        "outputId": "ed7565fa-5172-42c8-ce96-39b3cd357d5b"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Epoch=0\tloss=1.026e+02\n",
            "Epoch=1000\tloss=1.134e+01\n",
            "Epoch=2000\tloss=7.139e+00\n",
            "Epoch=3000\tloss=4.549e+00\n",
            "Epoch=4000\tloss=2.824e+00\n",
            "Epoch=5000\tloss=2.343e+00\n",
            "Epoch=6000\tloss=2.147e+00\n",
            "Epoch=7000\tloss=1.904e+00\n",
            "Epoch=8000\tloss=1.563e+00\n",
            "Epoch=9000\tloss=1.016e+00\n",
            "Epoch=10000\tloss=9.470e-01\n",
            "Epoch=11000\tloss=9.235e-01\n",
            "Epoch=12000\tloss=9.041e-01\n",
            "Epoch=13000\tloss=8.905e-01\n",
            "Epoch=14000\tloss=8.731e-01\n",
            "Epoch=15000\tloss=8.661e-01\n",
            "Epoch=16000\tloss=8.594e-01\n",
            "Epoch=17000\tloss=8.494e-01\n",
            "Epoch=18000\tloss=8.424e-01\n",
            "Epoch=19000\tloss=8.380e-01\n",
            "Epoch=20000\tloss=8.348e-01\n",
            "Epoch=21000\tloss=8.259e-01\n",
            "Epoch=22000\tloss=8.208e-01\n",
            "Epoch=23000\tloss=8.104e-01\n",
            "Epoch=24000\tloss=8.035e-01\n",
            "Epoch=25000\tloss=7.996e-01\n",
            "Epoch=26000\tloss=7.991e-01\n",
            "Epoch=27000\tloss=7.855e-01\n",
            "Epoch=28000\tloss=7.851e-01\n",
            "Epoch=29000\tloss=7.818e-01\n",
            "Epoch=30000\tloss=7.784e-01\n",
            "Epoch=31000\tloss=7.729e-01\n",
            "Epoch=32000\tloss=7.824e-01\n",
            "Epoch=33000\tloss=7.703e-01\n",
            "Epoch=34000\tloss=7.798e-01\n",
            "Epoch=35000\tloss=7.726e-01\n",
            "Epoch=36000\tloss=7.629e-01\n",
            "Epoch=37000\tloss=7.736e-01\n",
            "Epoch=38000\tloss=7.748e-01\n",
            "Epoch=39000\tloss=7.663e-01\n",
            "Epoch=40000\tloss=7.684e-01\n",
            "Epoch=41000\tloss=7.751e-01\n",
            "Epoch=42000\tloss=7.682e-01\n",
            "Epoch=43000\tloss=7.782e-01\n",
            "Epoch=44000\tloss=7.733e-01\n",
            "Epoch=45000\tloss=7.665e-01\n",
            "Epoch=46000\tloss=7.643e-01\n",
            "Epoch=47000\tloss=7.625e-01\n",
            "Epoch=48000\tloss=7.619e-01\n",
            "Epoch=49000\tloss=7.627e-01\n",
            "Epoch=50000\tloss=7.688e-01\n",
            "Epoch=51000\tloss=7.621e-01\n",
            "Epoch=52000\tloss=7.633e-01\n",
            "Epoch=53000\tloss=7.616e-01\n",
            "Epoch=54000\tloss=7.713e-01\n",
            "Epoch=55000\tloss=7.645e-01\n",
            "Epoch=56000\tloss=7.589e-01\n",
            "Epoch=57000\tloss=7.626e-01\n",
            "Epoch=58000\tloss=7.668e-01\n",
            "Epoch=59000\tloss=7.682e-01\n",
            "Epoch=60000\tloss=7.605e-01\n",
            "Epoch=61000\tloss=7.642e-01\n",
            "Epoch=62000\tloss=7.570e-01\n",
            "Epoch=63000\tloss=7.593e-01\n",
            "Epoch=64000\tloss=7.541e-01\n",
            "Epoch=65000\tloss=7.565e-01\n",
            "Epoch=66000\tloss=7.604e-01\n",
            "Epoch=67000\tloss=7.644e-01\n",
            "Epoch=68000\tloss=7.555e-01\n",
            "Epoch=69000\tloss=7.600e-01\n",
            "Epoch=70000\tloss=7.641e-01\n",
            "Epoch=71000\tloss=7.505e-01\n",
            "Epoch=72000\tloss=7.549e-01\n",
            "Epoch=73000\tloss=7.627e-01\n",
            "Epoch=74000\tloss=7.575e-01\n",
            "Epoch=75000\tloss=7.568e-01\n",
            "Epoch=76000\tloss=7.699e-01\n",
            "Epoch=77000\tloss=7.498e-01\n",
            "Epoch=78000\tloss=7.505e-01\n",
            "Epoch=79000\tloss=7.544e-01\n",
            "Epoch=80000\tloss=7.511e-01\n",
            "Epoch=81000\tloss=7.579e-01\n",
            "Epoch=82000\tloss=7.470e-01\n",
            "Epoch=83000\tloss=7.527e-01\n",
            "Epoch=84000\tloss=7.583e-01\n",
            "Epoch=85000\tloss=7.445e-01\n",
            "Epoch=86000\tloss=7.476e-01\n",
            "Epoch=87000\tloss=7.633e-01\n",
            "Epoch=88000\tloss=7.500e-01\n",
            "Epoch=89000\tloss=7.467e-01\n",
            "Epoch=90000\tloss=7.447e-01\n",
            "Epoch=91000\tloss=7.495e-01\n",
            "Epoch=92000\tloss=7.487e-01\n",
            "Epoch=93000\tloss=7.542e-01\n",
            "Epoch=94000\tloss=7.583e-01\n",
            "Epoch=95000\tloss=7.612e-01\n",
            "Epoch=96000\tloss=7.489e-01\n",
            "Epoch=97000\tloss=7.502e-01\n",
            "Epoch=98000\tloss=7.474e-01\n",
            "Epoch=99000\tloss=7.563e-01\n",
            "CPU times: user 1min 57s, sys: 256 ms, total: 1min 57s\n",
            "Wall time: 3min 3s\n"
          ]
        }
      ],
      "source": [
        "%%time\n",
        "epochs = 100_000\n",
        "for _ in range(epochs):\n",
        "  opt_state,params = update(opt_state,params,colloc,conds)\n",
        "\n",
        "  # print loss and epoch info\n",
        "  if _ %(1000) ==0:\n",
        "    print(f'Epoch={_}\\tloss={loss_fun(params,colloc,conds):.3e}')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 282
        },
        "id": "eWeNvDsdDEuI",
        "outputId": "1b79ef5d-7f2d-40ee-f4a1-4d0ba1cb0edd"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7f238924fc10>"
            ]
          },
          "metadata": {},
          "execution_count": 36
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deVxU1f/H8ddhExBXoNxS3HPBfWtRcbfc16wMELWvpZmZS2ZppqX+0iyXNHNJs0wzc8sNNc2d3BdMS8UVE1FxQxE4vz8uoiIIwsCdYT7Px2MezcydufO+TH06c+455yqtNUIIIWyPg9kBhBBCpI8UcCGEsFFSwIUQwkZJARdCCBslBVwIIWyUU1Z+mJeXl/bx8cnKjxRCCJu3e/fuS1pr76TPZ2kB9/HxYdeuXVn5kUIIYfOUUqeSe166UIQQwkZJARdCCBslBVwIIWxUlvaBJ+fu3bucPXuW27dvmx1FPAFXV1eKFCmCs7Oz2VGEsFumF/CzZ8+SK1cufHx8UEqZHUekgdaayMhIzp49S/Hixc2OI4TdMr0L5fbt23h6ekrxtiFKKTw9PeVXkxAmM72AA1K8bZB8Z0KYzyoKuBBCZDsxMXDwIMybB0ePZspHSAG3cX5+fqlOjlqyZAmhoaGJj4cNG8a6desyO5oQ9uHOHaNQb9ly/7nmzcHDAypVgjfegOXLM+WjTT+JKTLfkiVLaNmyJeXLlwfg008/NTmREDbmyhUID4eE/4aYOBHWrjVa1idOQHw8lC0Lf/9tbK9ZE6pUgcqVjSJepkymxLL7FvjNmzdp0aIFlStXpmLFiixYsIANGzbQtm3bxNcEBwfTrl07ADw8PBg4cCAVKlSgcePGhISE4OfnR4kSJVi2bBkA33//PW3atMHPz4/SpUszYsQIwGj5fvXVV4n7HTp0KF9//XWqeQDWr19P1apV8fX1JSgoiDt37jxyLB4eHon3Fy1aRGBgINu2bWPZsmUMHDiQKlWqcPz4cQIDA1m0aNFj9+vj48Pw4cOpVq0avr6+/H3vX0whsoPbt+H8eThyBLZvh1Wr4Pvv4d4VyqZNgwYNoFw5yJsX8ueHWrXubz92DM6eNYr00KHw00+wcOH9/Y8cCWPGwKuvQoUKkEnDba2vBe7n9+hznTvD22/DrVvw8suPbg8MNG6XLkHHjg9v27jxsR+3evVqChUqxO+//w5AVFQUuXPn5u233yYiIgJvb29mz55NUFAQYBTYhg0b8sUXX9CuXTs++ugjgoODCQ0NJSAggNatWwMQEhLCoUOHcHd3p2bNmrRo0YKgoCDat29Pv379iI+P5+effyYkJCTVPLdv3yYwMJD169dTpkwZ/P39mTp1Kv369Uvtr8nzzz9P69atadmyJR2T/G1S26+Xlxd79uzhm2++Ydy4ccyYMSPVzxMi02ht3BwcICLCaP1ev/7wzd/fKLbBwTB3rvHctWtw9SpERRnF+qmnYPRoSO6XaNu2RsG+dQtiY6FiRWjSBHx8jFt8PDg6wuTJWX30ybL7Frivry/BwcEMHjyYzZs3kydPHpRSvPHGG8ybN4+rV6+yfft2XnrpJQBcXFxo3rx54nvr16+Ps7Mzvr6+hIWFJe63SZMmeHp64ubmRvv27dmyZQs+Pj54enqyd+9e1q5dS9WqVfH09Ew1z9GjRylevDhlEn6GBQQE8Oeff2b42FPbb/v27QGoXr36Q8cmhMVobRTjvXthxQqjmwIgJARat4YXXjBawU89ZbRid+40tq9cCXXrGg26V16BHj3gvffg9Glj+7lzRp/0iRPGycRChaBOnfuf27Kl0cqeP9/Y17Ztxmtz5TK29+8PmzfDL78Y3SX9+0P79kbxtiLW1wJ/XIvZ3f3x2728Um1xJ1WmTBn27NnDypUr+eijj2jUqBHDhg2jW7dutGrVCldXVzp16oSTk/GncnZ2ThxC5+DgQI4cORLvx8bGJu436TC7e4979OjB999/z4ULFxJb9anladOmTZqO5cHPtMQY7XvH5ujo+NCxCfHELl2C0FAoVsy47d8Pr78Ox48b3Rn3LFwInToZRff0afD0BF9fo1Xt6QlPP228rnFjow86Vy7jZGGuXMYtTx5j+71f5SmpWdO42TjrK+BZ7Pz58+TPn5+uXbuSN2/exG6CQoUKUahQIUaNGpWuERvBwcFcvnwZNzc3lixZwqxZswBo164dw4YN4+7du/z0009pyjNo0CDCwsL4999/KVWqFD/88AP169d/5L1PP/00R44coWzZsvz222/kSmhN5MqVi+vXrz/y+rJly6Zpv0I8sagoGDvWKNT79hn9zQBffAEDBhjFuGRJY7RG0aJQpIhxe/ZZ43Uvvmi8LyWFCxs3O2f3BfzgwYMMHDgQBwcHnJ2dmTp1auK2119/nYiICMqVK/fE+61VqxYdOnTg7NmzdO3alRo1agBGF0yDBg3Imzcvjsn8HEsuj6urK7Nnz6ZTp07ExsZSs2ZNevXq9ch7x4wZQ8uWLfH29qZGjRrcuHEDgC5dutCzZ08mTpyYePISSPN+hUjRjRvw119GF8S2bUY3xccfg6srfP21UaQbNTJGY1SsaJz0A6NYL11qbvbsQGv92BswC7gIHHrgufxAMPBPwj/zpbYfrTXVq1fXSYWGhj7ynLXo3bu3njFjxhO/b/bs2bp3797JbouLi9OVK1fWx44dy2g801nzdycySWzs/fuNG2vt6Hjv1KLW5ctrPXbs/e0xMVmfL5sCdulkampaTmJ+DzRP8twHwHqtdWlgfcLjbKV69eocOHCArl27WmyfoaGhlCpVikaNGlG6dGmL7VeITHP3rtGyHjkS6teHqlXvb3vhBRgyxDgJePkyHD4Mgwbd3y4rVWY6pe+Na3zci5TyAVZorSsmPD4K+Gmtw5VSBYGNWuuyqe2nRo0aOumswSNHjqSri0KYT767bG7iRBg2zOjPVgqqVTO6Qz7/PMOjMe7cuZN4knzjxo2cPHmSS5cucf36dW7evEmRIkV47733AAgMDOTUqVPExsZy9+5dYmNjef7555k4cSJgDJUNDw+/1zuA1pqWLVsyZcoUAMqXL8+1a9cStwG8+uqrjBs3DjDOd8XFxT20/X//+x8jR47k1q1bFCtW7JH8/fv3Z8iQIVy8eJEKFSo8sn3YsGG88847nDhxgtq1a/PFF18Q+LiTqqlQSu3WWtdI+nx6+8Cf1lonjPfhAvD0Yz74TeBNgKJFi6bz44QQmerUKWO699KlMHUqlCpljBbp0AFeesmY1JJkyOvjxMXFJZ7j+fnnn9m2bRunTp0iLCyMsLAwihUrxoEDBwAYMmQIO3bsSHxvzpw5qVu3bmIBj4uLIz4+HhcXF9zd3XFycnpo+G3NmjW5cuUKSqnEkVgPFtWGDRsSHR2duE0pha+vb+L2e3M3HtxeqVIlAJycnOjcufMjx1exYkXAOI+U3PZnE07G5sqVi86dO1OyZMk0/+2eRHpb4Fe11nkf2H5Fa50vtf1ICzx7ke/Oxl26BJMmwbJl90d8PPssTJ9ujLFOo5MnT7Jt2zZCQ0MJDQ3lyJEjREVFEZ4wpvuVV15h1apV+Pj44OPjQ7FixXj22Wfp3bs3AMeOHcPJyQkvLy88PDxwcLD76SmPsHQL/D+lVMEHulAuZiyeECLTxcQY8yScnY0WtYODMd27Vi1jeF/r1o9dsyM+Pp5jx46xe/dudu/ezeeff46rqytTpkxh/PjxODo6Urp0aSpUqED58uWJjY3FycmJuXPn4uLikuISxGUyaZ0Qe5DeAr4MCADGJPxTxgMJYY2uXjVOMi5daqz3cf06NGtmFPD8+eHixfuTX1KwY8cOPv/8c7Zs2cKVK1cAcHNzo0ePHpQvX57evXsTEBBA2bJlcXFxeeT99/q6heWl+ltFKTUf2A6UVUqdVUp1xyjcTZRS/wCNEx7bLKUU77//fuLjcePG8cknnwDwySef4O7uzsWL939kPLhoVHo9uKBUSjZu3Mi2bdsSH0+bNo25c+dm+LNFNhcRcf9+27bGjMdNm4wp58uWwW+/3d/+QPGOiYlh48aNDB8+nIYNGxIcHJz4/N9//027du2YOXMmBw4c4Nq1a4mrWxYvXhxfX99ki7fIXKm2wLXWr6awqZGFs5gmR44cLF68mCFDhuDl5fXIdi8vL8aPH8/YsWOzNNfGjRvx8PDg+eefB5BJNiJ5WhtriSxdatxCQ42Wdd68xvA/Z2ejmySFvuXIyEgCAgLYuHEjN2/exMHBgcqVKycux1CvXj2OHTuWlUck0kjOFmCcaX7zzTeZMGFCstuDgoJYsGABly9fTnEfcXFxBAYGUrFiRXx9fRP3tW/fPurUqUOlSpVo165d4k/QB/n4+HDp0iUAdu3ahZ+fH2FhYUybNo0JEyZQpUoVNm/ezCeffJI49Cml/fr5+TF48GBq1apFmTJl2Lx5c4b+NsLKrV1rTEWvXh1GjYLcuY1hfvfUrWvMjkwo3hEREcyfP5+goCCGDx8OQL58+YiMjCQwMJAlS5Zw+fJl9uzZQ6tWrcw4IvEErG4qvV8yy8l27tyZt99+m1u3bvFyMsvJBgYGEhgYyKVLlx5ZMnVjGhe36t27N5UqVWLQgxMREnh4eBAUFMTXX3+duLZ3Uvv27ePcuXMcOnQIgKtXrwLg7+/PpEmTqF+/PsOGDWPEiBEPrQmeEh8fH3r16oWHhwcDBgwAjLW773ncfmNjYwkJCWHlypWMGDFCrr6TXfz3H/z+uzHcLyDA6B4pVsxoXY8cCS1agLd3sm/94osvmD9/Pnv37gUgf/78BAQEAMZCbNu3b8+ywxCWIy3wBLlz58bf3z9xckBSffv2Zc6cOckuCgVQokQJTpw4wTvvvMPq1avJnTs3UVFRXL16NXGBKEstA5vafmUZ2Gzk7l2jRf3cc1CwIHTvDrt2GScnwbgKzK+/GivveXujtebAgQOMHz+e1157LXFiyj///EPu3Ln57LPPCAkJ4eLFi3z55ZfmHZewCKtrgT+uxezu7v7Y7V5eXmlucSenX79+VKtWjW7duj2yLW/evLz22muJs7uSypcvH/v372fNmjVMmzaNhQsXptglk5STkxPx8fGALANr92Ji4M8/jXWx33gDnJxg9myjP/uTT4yhfpUrGzMjH7Bjxw6mTJnCunXruHDhAgDlypXj0qVLeHt78+2336Y4jE/YLmmBPyB//vx07tyZmTNnJru9f//+fPvtt8kWxUuXLhEfH0+HDh0YNWoUe/bsIU+ePOTLly+xHzql5Vp9fHzYvXs3AL/++mvi8yktA5vW/QobERkJP/xgXHnK29u4AsyQIcbJSaWMJVn/+suY1l6lCreio1mzZg0DBgzg8OHDgLEM8erVq2nQoAGzZs3izJkzhIaG4p3QpSLFO3uyuha42d5//30mp3C5JC8vL9q1a5dsy/rcuXN069YtsSU9evRoAObMmUOvXr24desWJUqUYPbs2Y+8d/jw4XTv3p2PP/74oXMArVq1omPHjixdupRJkyY99J607FdYKa2NayqWLGm0sEeNgq++ggIFjKF+rVoZa47cK7ru7kRFRfHtt98SHBzM5s2buXPnDi4uLlStWpUKFSrQunVr2rZtK7MY7UyaptJbikylz17ku3sCsbGwdasxDnv5cvjnH2Nsdr168O+/xlXPq1dPHC1y/vx5goOD8fDwoEOHDty8eRNPT09Kly5N06ZNadKkCfXq1cPd3d3kAxNZwdJT6YUQaRUaalxh5soVcHExZkG+++79q8+UKgUYo4x+//131q5dm9g18vLLL9OhQwdy5syZeLUmIe6RAi6EJZ08abSwly0zxl+PGgWlS0PHjsYU9qZNIVcu4uPj2bt3L3uXLaNHjx4AfPnll2zYsIG6desSGBhIkyZNHlo1T4q3SMoqCrjWWk6y2Jis7HqzCePHw5w5cPCg8bhcuftjsp2dYfp0wsPDWblwIcHBwaxbt47IyEiUUrRv3578+fMzbdo0vLy8cHNzM+84hE0x/YyHq6srkZGRUhBsiNaayMhIXF1dzY5iDq3hwAHjmo/37NljDPX78kujfzs0lFs9e7Jq1SoiEtYmWZbQ2v7zzz9p0aIF8+bNIzw8PLFl/cwzz0jxFk/E9JOYd+/e5ezZsxYZ/yyyjqurK0WKFMHZni6bdfgw/PgjLFpkFGkHBwgLg2eegbg4tIMDhw4dYs2aNaxZsyZxtMiMGTPo3r07kZGRhIeHU6FCBfnFKZ5ISicxTS/gQtiExYuNq9M4OkLDhkafdtu2XHdzIzIyEh8fH8LDwylUqBBgXMarWbNmNG/enLp160rLWmSIjEIRIq2io42C/cMP0LIl9OkDjRsbY7W7dOFcbCzLly9naUAAGzZsoFmzZixbtoyCBQuycOFC6tSpwzPPPGP2UQg7IAVciHsOHzYuJzZ3rrHWSNGikLCuDLlzw7vvEhAQkLgme8mSJenTp0/i2jMAnTp1MiO5sFNSwIV9i401ZkMC9O4N27YZXSU9e3KqeHEW/PILy+vWZe3atbi5ueHn50e5cuVo3bo15cqVk75sYSop4MI+nT8PU6bArFnGBX2ffhq++YYrzs78uGYNP330UeISq7Vr1+bcuXOUKlUq2YXOhDCL6cMIhchSe/eCvz/4+MDo0VCnDnHXrxMVFQXly/Pv1au888473Lx5k88//5zjx4+zY8cOSiXMlhTCmkgLXNiP8+ehRg1wd4e33iKsfXtmrV/P9w0b0rRpU2bMmEGNGjU4fPhw4vUehbBmUsBF9qW1MaU9JAQ++wwKFYJFi9jk4sKE775jecOGaK1p2rRp4uXDlFJSvIXNkC4Ukf3Ex8Mvv0CVKsZlxxYu5OZ//xmzfdu1Y9Hq1WzdupUPPviAsLAwVq9eTZs2bcxOLcQTk4k8InvZvx9ee81YAbBsWf7r04evz57lm2nTWLp0KfXr1+fy5cu4u7vb71IAwubIRB6RvV2/DrlyQeHC4OFB2MSJfHHkCLMGDuTOnTu0b98eT09PQFb1E9mHFHBh2/btMy4/FhkJO3eClxexW7fyQrFiRERE4O/vz6BBgyhTpozZSYWwOOkDF7bp3DljOGDVqrBzJycbN2bwoEHExsbi5OTEnDlzOHHiBDNmzJDiLbItaYEL27N9u3Hh39hYzvXuzahbt5jxxRc4OTnRsXNnatasSePGjc1OKUSmkwIubIPWEB5uDAWsVo3o119nuNZMmjmT2NhYevbsydChQylcuLDZSYXIMtKFIqxfaKixGmDt2nDrFuTIgfOUKazesYNOnTpx9OhRvvnmGynewu5IARfWKzoahg41xnPv3cvq1q1p8PLLXLt2DScnJ3bu3MncuXMpUaKE2UmFMIV0oQjrFB5uXMn9xAmOtG5N/xs3WP3NN5QsWZKwsDAqVaokF0kQdk8KuLAud+8aFwEuUIDYBg34uHZtxv3yCzlz5mT8+PH06dMHFxcXs1MKYRWkC0VYB62NpV1LlYIzZ0ApHL/7jv1Xr9K1a1eOHTtG//79pXgL8YAMFXCl1HtKqcNKqUNKqflKKZmbLJ7c2bPw8svQvTthBQrQuVcvTp48iVKKpUuXMnv2bJ566imzUwphddJdwJVShYG+QA2tdUXAEehiqWDCTsyeDRUqELtpE//Xpg3lDx5k5aZN7N+/H8C+rnovxBPKaB+4E+CmlLoLuAPnMx5J2JXNmzlYujRBMTHsWrqUtm3bMnHiRLkosBBpkO4WuNb6HDAOOA2EA1Fa67VJX6eUelMptUsptSsiIiL9SUX2oLXR6t6713g8ZQrTatXi1IULLFiwgMWLF0vxFiKN0r2crFIqH/Ar8ApwFfgFWKS1npfSe2Q5WTsXGQk9e8Jvv7G7Qwcchg6latWqXL9+ndu3b+Pt7W12QiGsUkrLyWbkJGZj4KTWOkJrfRdYDDyfgf2J7Cw4GHx9ubNiBUMaNqT2kiUMHjwYgFy5cknxFiIdMlLATwN1lFLuSikFNAKOWCaWyFZWroSmTTng6kpNHx/GbNiAv78/CxcuNDuZEDYtI33gO4FFwB7gYMK+plsol8gOYmKMfzZuzNZevah57hwXr11jxYoVzJo1i7x585qbTwgbl6Fx4Frr4VrrZ7XWFbXWb2it71gqmLBhWsOkSVChAnGXLoGLC7UmTuS9997j4MGDtGjRwuyEQmQLMhNTWNaVK9C+PbpvX+bkzEnF558nMjISZ2dnxowZI33dQliQFHBhOSEhUK0aV5Yvp1OlSgTu3493gQJER0ebnUyIbEkKuLCcTz9lW3Q0Vby9WRoaypgxY/jjjz8oUqSI2cmEyJZkNUKRMVeuwJ07UKAAzJ7NGH9/HI8eZevWrdSqVcvsdEJka1LARfrt3AmvvMJ/RYsSO38+hQsXZva8eTg5OZEnTx6z0wmR7UkXinhyWsPUqVC3LsG3b1P58GG6desGgKenpxRvIbKIFHDxZKKjISiI2Lff5sOiRWl28SJeBQsyYcIEs5MJYXekgIsnc/s2EZs20cTHh9HHj9OjRw9CQkKoUKGC2cmEsDtSwEXabN1qnKzMlw+XLVu4kicPc+fOZfr06bi7u5udTgi7JAVcPJ7WMHYsum5d5nTpwu3bt8lTqBC7d+/mjTfeMDudEHZNCrhI2c2b0Lkz1z/4gFcKFyZwyRJmzZoFgKOjo8nhhBAyjFAk78wZaNOG0H37aO/tzT/nzzN27Fjeeusts5MJIRJIARfJu3aNVadP0ylHDnIqxfr16/Hz8zM7lRDiAdKFIh62Y4fR712hAiU3bqSunx979uyR4i2EFZICLgzx8fDhh9x47jkmBwWhtaZMxYqsWrWKwoULm51OCJEMKeACbtyADh04MXo0z+XLx7tz57Jnzx6zUwkhUiF94PbuwgVo0YINe/fSyd0dDaxevZrq1aubnUwIkQppgdu73bv57vBhmirF08WKERISQpMmTcxOJYRIAyng9ioiwvhnixY8M3curdu0YceOHZQqVcrcXEKINJMCbo9+/pkbxYqxbNQoAJp37szixYvJnTu3ycGEEE9CCrg90Rr+7/849+qr1HN0pOOnn3L69GmzUwkh0klOYtqL+Hh49132T55MCzc3ooClS5dStGhRs5MJIdJJCri9+PVXVk6ezCvOzuT19GTLihVUrlzZ7FRCiAyQLhR70bEjoT16ULpiRXbu3CnFW4hsQAp4dnb1KrptW05v2ABKMeC779i+fTuFChUyO5kQwgKkgGdXFy4QX68e7y1fjm+rVpw4cQKAHDlymBxMCGEp0geeHYWFEdOoEd1OneKn+Hje7dkTHx8fs1MJISxMCnh2ExbGjRdfpON//7EmLo7Ro0czePBglFJmJxNCWJgU8OzmqacYlzs3weHhzJw5k6CgILMTCSEyiRTw7OLYMShQAHLnZuj+/TQJCeGFF14wO5UQIhPJSczs4OBBLjz/PJ3Kl+fixYs4OztL8RbCDkgBt3X79nG2Xj3qR0Wx6soV/vnnH7MTCSGyiBRwW3bgACf9/Kh3/TrhOXKwZu1aaXkLYUcyVMCVUnmVUouUUn8rpY4opZ6zVDCRCq3599VXqXfjBlc9PFj/xx9SvIWwMxltgX8NrNZaPwtUBo5kPJJIE6XIOWMGJapV449Nm6hZs6bZiYQQWSzdBVwplQeoB8wE0FrHaK2vWiqYSMGpU1wYMIDYmBgKPvccG2VdEyHsVkZa4MWBCGC2UmqvUmqGUipn0hcppd5USu1SSu2KuHcVGJE+585xtl49np8wgbf8/QFkgo4QdiwjBdwJqAZM1VpXBW4CHyR9kdZ6uta6hta6hre3dwY+zs5dvMj5+vVpcOYMkW5uvPn++2YnEkKYLCMF/CxwVmu9M+HxIoyCLizt2jX+a9yYhidOcMHVldXBwdLnLYRIfwHXWl8AziilyiY81QgItUgq8RC9ezdtDx/mjIsLK9es4bnnZLCPECLjU+nfAX5USrkAJ4BuGY8kklINGvDlihXccnGhbt26ZscRQliJDBVwrfU+oIaFsogHaU1M376scXKi1YQJPPfSS2YnEkJYGZmJaaXiR4wgYPJkWn/1Ffv27TM7jhDCCkkBt0J6+nT6jhjBz8DYMWOoUqWK2ZGEEFZICri1WbOGEb16MQUY2L8/gwYPNjuREMJKSQG3MvuWL2eE1nTr2pWx48aZHUcIYcXkgg5Wpsrkyaxp1owGzZvLLEshxGNJC9wa3LzJ/hdfZMu33wLQtFUrnJ2dTQ4lhLB20gI3W1wcZ9q14+WtW3E/cYIj3bvj5CRfixAidVIpTBbVrx8tgoO54erK6jVrpHgLIdJMqoWJYmfPptPkyRxxcGDV8uX4+vqaHUkIYUOkgJtFa76fMoVgYMbUqTRu3NjsREIIGyMF3CxKEbRjB4WWLOHljh3NTiOEsEEyCiWr3bnDX+3bc3rHDhycnKR4CyHSTVrgWUlrTgcE0OK333j277/ZdPiwjPUWQqSbtMCz0M0JE2izYAF3cuRg+uLFUryFEBkiLfAsordvp9uAAewHfv/1V5599lmzIwkhbJwU8CzyXffu/KI1Yz/5hJdatDA7jhAiG5ACnkW6BAdzc/Jk+g0bZnYUIUQ2IX3gmSzihx+Ijowkd+HCvDd6tPR7CyEsRgp4Jrq7ciXt/f1pUqUKWmuz4wghshnpQsksZ87wQYcObAF+/PRTaXkLISxOWuCZITaWRU2a8OXt2/Tp2pXXunUzO5EQIhuSAp4Jjg8YQNDRo9QuVYrxM2eaHUcIkU1JAc8E8S1bUqNECRasW4eLi4vZcYQQ2ZT0gVtSdDS4ulK6cWM2HD9udhohRDYnLXBL0Zo1jRrxmo8P169dMzuNEMIOSAG3kPCxY3lj+3YO3b2Lk1zPUgiRBaSAW0D8/v34f/ghNxwc+HntWtzc3MyOJISwA9IHnlF37vB/zZqxTmtmjB9P+YoVzU4khLAT0gLPoJv79vFVRASd69YlqF8/s+MIIeyItMAzKGft2uz5+29cPT1ltqUQIktJCzy9btxgU9++xMfEUKh0afLnz292IiGEnZECnk7ru3TBb9IkJg0ZYnYUIYSdynABV0o5KqX2KqVWWCKQLbj6yy8E/v47ZfPnp3SgsaMAAAsESURBVOfIkWbHEULYKUv0gb8LHAFyW2Bf1i8ykncCAggHti9bhru7u9mJhBB2KkMtcKVUEaAFMMMycazfL23aMC86mo979aLmCy+YHUcIYccy2gL/ChgE5LJAFpvw9Ouv80pMDB9OnGh2FCGEnUt3C1wp1RK4qLXencrr3lRK7VJK7YqIiEjvx5kvPh6Aem+9xc8hITjLdHkhhMky0oXyAtBaKRUG/Aw0VErNS/oirfV0rXUNrXUNb2/vDHycuX5r1oyB1atzOzra7ChCCAFkoIBrrYdorYtorX2ALsAGrXVXiyWzIpcXL+atdetYf+4cjk4y90kIYR2kGqXmxg3e8/cnEli9ZIl0nQghrIZFCrjWeiOw0RL7sjYrX3+duTdv8nFAAFXq1DE7jhBCJJKZmI8Rd+ECfZcvp0L+/Az99luz4wghxEOkC+UxHAsUYNnKldzNk4ccOXKYHUcIIR4iBTwF1/bsIXfVqpRv3tzsKEIIkSzpQknG3aNHeaFGDfq9+KLZUYQQIkVSwJPSmvGtW3NIaxr16GF2GiGESJEU8CROTp3Kp8eO0c7Xl1bdupkdRwghUiQF/AE6Kore/fvj6ODAxGXLzI4jhBCPJQX8AWe2bGFnTAwj+/aliI+P2XGEEOKxZBTKA4q2aMHRM2fI+/TTZkcRQohUSQscID6ezUOHEn/nDl6FC+Mk650IIWyAFHBg38iR+H3+OeNk1IkQwobYfQHXV6/yzmefkd/JiZ5ff212HCGESDO77yv46dVX2XL3LtOHDiVf/vxmxxFCiDSz6xb49b/+YuDq1VT38iJoxAiz4wghxBOx6wJ+5vRpcru7M3nuXBwdHc2OI4QQT8Suu1DKd+jA4bZtpXgLIWySXbbAdXQ0s1q14sbp01K8hRA2yy4L+PIePei+YgVzx483O4oQQqSb0lpn2YfVqFFD79q1K8s+Lzl3Tp6kQsmS5PDwYF9kpFzjUghh9ZRSu7XWNZI+b3d94FO6dOG41qyZNEmKtxDCptlVF8qlTZv4NCSE5j4+NA0IMDuOEEJkiF0V8BsODtQuWJBxP/5odhQhhMgwu+pC8alblzXnz5sdQwghLMI+WuCxsXxVrx6n//jD7CRCCGExdlHANw4ezHubN/PjzJlmRxFCCIvJ9sMI46OiqOnlRYSjI0cvX8bN3T1LP18IITLKbocRzvP3Z09sLPM++kiKtxAiW8nWLfDoEycoU6oUBfLlY2dEBA4OdtFjJITIZuyyBX5bKV7y9eW1QYOkeAshsp1sXcDzFS/O9P37zY4hhBCZIts2S79v354dkyaZHUMIITJNtizg4evX0/u33/h66lSzowghRKbJlgV8VPfuxAAjf/jB7ChCCJFp0l3AlVLPKKX+UEqFKqUOK6XetWSw9Dq+cCHTT52iZ+3alKpe3ew4QgiRaTJyEjMWeF9rvUcplQvYrZQK1lqHWijbk9OaYX364Ax8LAtWCSGyuXQXcK11OBCecP+6UuoIUBgwrYDr+Hgq1qlDRVdXCpYsaVYMIYTIEhaZyKOU8gH+BCpqra8l2fYm8CZA0aJFq586dSrDnyeEEPYkpYk8GT6JqZTyAH4F+iUt3gBa6+la6xpa6xre3t4Z/bgU7Z4yhYVvvkl8TEymfYYQQliTDBVwpZQzRvH+UWu92DKR0iE+ng8++IA+s2YRHR1tWgwhhMhKGRmFooCZwBGt9ZeWi/TkNn/2Getu3GBw587kzJPHzChCCJFl0t0HrpR6EdgMHATiE57+UGu9MqX3ZMpiVvHxNMydm9Dbtzlx5QruuXJZdv9CCGEyiy9mpbXeAqgMpbKAjSNH8sfNm0zo2lWKtxDCrtj8TMwYNzfqennxP5k2L4SwM9l6PXAhhMgOMm0YoVl0bCw/vP46t86dMzuKEEKYwmbXA183bBj+P/1EtLs7b373ndlxhBAiy9lkF4qOi+OFXLk4e/cu/1y7Rg43NwukE0II65StLqm2YdQotkdHMzUoSIq3EMJu2V4LXGsa5MnDsVu3OBEVRY6cOS0TTgghrFS2OYl5IzycOGdnBnToIMVbCGHXbK4LxaNQIf6MjCQ+Ls7sKEIIYSqbaoGHrV9P+PbtADg4OpqcRgghzGVTLfCBXbuy+eJFTl+7hot0nwgh7JzNtMCPLF7Mrxcu0L1uXSneQgiBDRXwMf374wb0mzHD7ChCCGEVbKKAn1i3jh9PneJ/NWviXaqU2XGEEMIq2EQB/2P+fFyAAdL6FkKIRDZRwLvPnMnpkycpVKmS2VGEEMJq2EQBB/Dy8TE7ghBCWBWbKeBCCCEeJgVcCCFslBRwIYSwUVLAhRDCRkkBF0IIGyUFXAghbJQUcCGEsFFSwIUQwkZl6SXVlFIRwKl0vt0LuGTBOFnN1vOD7R+DrecH2z8GW88P5hxDMa21d9Ins7SAZ4RSaldy14SzFbaeH2z/GGw9P9j+Mdh6frCuY5AuFCGEsFFSwIUQwkbZUgGfbnaADLL1/GD7x2Dr+cH2j8HW84MVHYPN9IELIYR4mC21wIUQQjxACrgQQtgoqyvgSqnmSqmjSql/lVIfJLM9h1JqQcL2nUopn6xPmbI05A9USkUopfYl3HqYkTMlSqlZSqmLSqlDKWxXSqmJCcd3QClVLaszPk4a8vsppaIe+PsPy+qMqVFKPaOU+kMpFaqUOqyUejeZ11jt95DG/Fb9PSilXJVSIUqp/QnHMCKZ15hfi7TWVnMDHIHjQAnABdgPlE/ymreBaQn3uwALzM79hPkDgclmZ33MMdQDqgGHUtj+MrAKUEAdYKfZmZ8wvx+wwuycqRxDQaBawv1cwLFk/j2y2u8hjfmt+ntI+Lt6JNx3BnYCdZK8xvRaZG0t8FrAv1rrE1rrGOBnoE2S17QB5iTcXwQ0UkqpLMz4OGnJb9W01n8Clx/zkjbAXG3YAeRVShXMmnSpS0N+q6e1Dtda70m4fx04AhRO8jKr/R7SmN+qJfxdbyQ8dE64JR3xYXotsrYCXhg488Djszz6xSe+RmsdC0QBnlmSLnVpyQ/QIeFn7yKl1DNZE81i0nqM1uy5hJ/Gq5RSFcwO8zgJP8urYrQAH2QT38Nj8oOVfw9KKUel1D7gIhCstU7xOzCrFllbAbcHywEfrXUlIJj7/wcXWWMPxroSlYFJwBKT86RIKeUB/Ar001pfMzvPk0olv9V/D1rrOK11FaAIUEspVdHsTElZWwE/BzzYIi2S8Fyyr1FKOQF5gMgsSZe6VPNrrSO11ncSHs4AqmdRNktJy3dktbTW1+79NNZarwSclVJeJsd6hFLKGaP4/ai1XpzMS6z6e0gtv618DwBa66vAH0DzJJtMr0XWVsD/AkorpYorpVwwTgwsS/KaZUBAwv2OwAadcBbBCqSaP0k/ZWuM/kFbsgzwTxgFUQeI0lqHmx0qrZRSBe71UyqlamH8N2AtDQDAGGECzASOaK2/TOFlVvs9pCW/tX8PSilvpVTehPtuQBPg7yQvM70WOWXlh6VGax2rlOoDrMEY0TFLa31YKfUpsEtrvQzjX4wflFL/Ypys6mJe4oelMX9fpVRrIBYjf6BpgZOhlJqPMULASyl1FhiOcQIHrfU0YCXGCIh/gVtAN3OSJi8N+TsCbymlYoFooIsVNQDueQF4AziY0AcL8CFQFGzie0hLfmv/HgoCc5RSjhj/c1motV5hbbVIptILIYSNsrYuFCGEEGkkBVwIIWyUFHAhhLBRUsCFEMJGSQEXQggbJQVcCCFslBRwIYSwUf8PpGEt+/Jo6qUAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "source": [
        "lam_sol= sp.lambdify(t,sol)\n",
        "\n",
        "dT = 1e-3\n",
        "Tf = jnp.pi\n",
        "T = np.arange(0,Tf+dT,dT)\n",
        "\n",
        "\n",
        "sym_sol =np.array([lam_sol(i) for i in T])\n",
        "\n",
        "plt.plot(T,sym_sol,'--r',label='sympy solution')\n",
        "plt.plot(T,fwd(params,T.reshape(-1,1))[:,0],'--k',label='NN solution')\n",
        "plt.legend()"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "[6] ODE-PINN finite difference.ipynb",
      "provenance": [],
      "include_colab_link": true
    },
    "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.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}