[1]_ODE_PINN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/ASEM000/Physics-informed-neural-network-in-JAX/blob/main/%5B1%5D_ODE_PINN.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  Septemeber 2021"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "vAR0swbLX_ZI",
    "outputId": "23ef5ef3-495b-4582-96a2-3108e8dce0bb"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting optax\n",
      "  Downloading optax-0.0.9-py3-none-any.whl (118 kB)\n",
      "\u001b[K     |████████████████████████████████| 118 kB 8.5 MB/s eta 0:00:01\n",
      "\u001b[?25hCollecting chex>=0.0.4\n",
      "  Downloading chex-0.0.8-py3-none-any.whl (57 kB)\n",
      "\u001b[K     |████████████████████████████████| 57 kB 6.1 MB/s  eta 0:00:01\n",
      "\u001b[?25hRequirement already satisfied: jaxlib>=0.1.37 in /usr/local/lib/python3.7/dist-packages (from optax) (0.1.70+cuda110)\n",
      "Requirement already satisfied: jax>=0.1.55 in /usr/local/lib/python3.7/dist-packages (from optax) (0.2.19)\n",
      "Requirement already satisfied: absl-py>=0.7.1 in /usr/local/lib/python3.7/dist-packages (from optax) (0.12.0)\n",
      "Requirement already satisfied: numpy>=1.18.0 in /usr/local/lib/python3.7/dist-packages (from optax) (1.19.5)\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: 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: toolz>=0.9.0 in /usr/local/lib/python3.7/dist-packages (from chex>=0.0.4->optax) (0.11.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) (1.12)\n",
      "Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (from jaxlib>=0.1.37->optax) (1.4.1)\n",
      "Installing collected packages: chex, optax\n",
      "Successfully installed chex-0.0.8 optax-0.0.9\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": 4,
   "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": 19,
   "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": 20,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 54
    },
    "id": "klgFeU6bcTrC",
    "outputId": "cad6a477-37e0-414d-daf6-6e4f05a288cd"
   },
   "outputs": [
    {
     "data": {
      "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$"
      ],
      "text/plain": [
       "Eq(-t*cos(2*pi*t) + Derivative(y(t), t) + Derivative(y(t), (t, 2)), 0)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diffeq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 60
    },
    "id": "E4Uu2hbiYJtv",
    "outputId": "e58180aa-779a-49fa-9771-b9f4763bb037"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left. \\frac{d}{d t} y{\\left(t \\right)} \\right|_{\\substack{ t=0 }} = 10$"
      ],
      "text/plain": [
       "Eq(Subs(Derivative(y(t), t), t, 0), 10)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.Eq(f(t).diff(t).subs(t,0),10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 38
    },
    "id": "29QUbt_2YwlJ",
    "outputId": "3d9edb85-a44f-42cc-8776-bccda97cbc7e"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle y{\\left(0 \\right)} = 1$"
      ],
      "text/plain": [
       "Eq(y(0), 1)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.Eq(f(t).subs(t,0),1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 101
    },
    "id": "r9KVq1yjYfld",
    "outputId": "b6e56baa-b419-49e0-bba9-3706d96b394c"
   },
   "outputs": [
    {
     "data": {
      "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)}$"
      ],
      "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)))"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.Eq(f(t),sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 37
    },
    "id": "MNVOpPyCW-GU",
    "outputId": "c333d114-01de-4132-8ebf-effd19785112"
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 61,
   "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 = jax.random.uniform(key,minval=tmin,maxval=tmax,shape=(N_c,1))\n",
    "colloc = t_c\n",
    "\n",
    "def ODE_loss(t,u):\n",
    "  u_t=lambda t:jax.grad(lambda t:jnp.sum(u(t)))(t)\n",
    "  u_tt=lambda t:jax.grad(lambda t : jnp.sum(u_t(t)))(t)\n",
    "  return -t*jnp.cos(2*jnp.pi*t) + u_t(t) + u_tt(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "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": 69,
   "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": 70,
   "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(lr)\n",
    "opt_state = optimizer.init(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "kBzGA8OVc8C6",
    "outputId": "effacbde-f915-47cd-8f13-a44e6b697062"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch=0\tloss=1.026e+02\n",
      "Epoch=100\tloss=1.500e+01\n",
      "Epoch=200\tloss=7.508e+00\n",
      "Epoch=300\tloss=5.177e+00\n",
      "Epoch=400\tloss=3.458e+00\n",
      "Epoch=500\tloss=2.615e+00\n",
      "Epoch=600\tloss=2.391e+00\n",
      "Epoch=700\tloss=2.254e+00\n",
      "Epoch=800\tloss=2.133e+00\n",
      "Epoch=900\tloss=2.013e+00\n",
      "Epoch=1000\tloss=1.855e+00\n",
      "Epoch=1100\tloss=1.549e+00\n",
      "Epoch=1200\tloss=1.018e+00\n",
      "Epoch=1300\tloss=9.915e-01\n",
      "Epoch=1400\tloss=9.756e-01\n",
      "Epoch=1500\tloss=9.618e-01\n",
      "Epoch=1600\tloss=9.480e-01\n",
      "Epoch=1700\tloss=9.327e-01\n",
      "Epoch=1800\tloss=9.145e-01\n",
      "Epoch=1900\tloss=8.902e-01\n",
      "Epoch=2000\tloss=7.379e-01\n",
      "Epoch=2100\tloss=3.372e-01\n",
      "Epoch=2200\tloss=1.753e-02\n",
      "Epoch=2300\tloss=4.962e-03\n",
      "Epoch=2400\tloss=3.304e-03\n",
      "Epoch=2500\tloss=2.670e-03\n",
      "Epoch=2600\tloss=2.238e-03\n",
      "Epoch=2700\tloss=1.883e-03\n",
      "Epoch=2800\tloss=1.577e-03\n",
      "Epoch=2900\tloss=1.311e-03\n",
      "Epoch=3000\tloss=1.082e-03\n",
      "Epoch=3100\tloss=8.856e-04\n",
      "Epoch=3200\tloss=7.198e-04\n",
      "Epoch=3300\tloss=5.825e-04\n",
      "Epoch=3400\tloss=4.711e-04\n",
      "Epoch=3500\tloss=6.700e-04\n",
      "Epoch=3600\tloss=3.161e-04\n",
      "Epoch=3700\tloss=3.296e-04\n",
      "Epoch=3800\tloss=2.278e-04\n",
      "Epoch=3900\tloss=2.032e-04\n",
      "Epoch=4000\tloss=1.785e-04\n",
      "Epoch=4100\tloss=1.624e-04\n",
      "Epoch=4200\tloss=1.494e-04\n",
      "Epoch=4300\tloss=1.377e-04\n",
      "Epoch=4400\tloss=2.545e-04\n",
      "Epoch=4500\tloss=1.195e-04\n",
      "Epoch=4600\tloss=2.162e-03\n",
      "Epoch=4700\tloss=1.041e-04\n",
      "Epoch=4800\tloss=1.321e-04\n",
      "Epoch=4900\tloss=9.035e-05\n",
      "Epoch=5000\tloss=8.385e-05\n",
      "Epoch=5100\tloss=8.181e-05\n",
      "Epoch=5200\tloss=7.205e-05\n",
      "Epoch=5300\tloss=2.472e-03\n",
      "Epoch=5400\tloss=6.137e-05\n",
      "Epoch=5500\tloss=6.695e-05\n",
      "Epoch=5600\tloss=5.197e-05\n",
      "Epoch=5700\tloss=4.714e-05\n",
      "Epoch=5800\tloss=7.266e-05\n",
      "Epoch=5900\tloss=3.922e-05\n",
      "Epoch=6000\tloss=1.470e-03\n",
      "Epoch=6100\tloss=3.234e-05\n",
      "Epoch=6200\tloss=1.919e-03\n",
      "Epoch=6300\tloss=2.676e-05\n",
      "Epoch=6400\tloss=2.389e-05\n",
      "Epoch=6500\tloss=2.301e-05\n",
      "Epoch=6600\tloss=1.962e-05\n",
      "Epoch=6700\tloss=3.187e-04\n",
      "Epoch=6800\tloss=1.618e-05\n",
      "Epoch=6900\tloss=2.400e-03\n",
      "Epoch=7000\tloss=1.355e-05\n",
      "Epoch=7100\tloss=2.970e-04\n",
      "Epoch=7200\tloss=1.146e-05\n",
      "Epoch=7300\tloss=1.040e-05\n",
      "Epoch=7400\tloss=3.527e-05\n",
      "Epoch=7500\tloss=9.057e-06\n",
      "Epoch=7600\tloss=4.930e-04\n",
      "Epoch=7700\tloss=3.244e-04\n",
      "Epoch=7800\tloss=1.300e-03\n",
      "Epoch=7900\tloss=7.630e-06\n",
      "Epoch=8000\tloss=1.313e-05\n",
      "Epoch=8100\tloss=1.451e-05\n",
      "Epoch=8200\tloss=7.152e-06\n",
      "Epoch=8300\tloss=6.462e-05\n",
      "Epoch=8400\tloss=5.616e-03\n",
      "Epoch=8500\tloss=6.044e-06\n",
      "Epoch=8600\tloss=7.831e-05\n",
      "Epoch=8700\tloss=5.532e-06\n",
      "Epoch=8800\tloss=8.422e-05\n",
      "Epoch=8900\tloss=1.461e-03\n",
      "Epoch=9000\tloss=5.241e-06\n",
      "Epoch=9100\tloss=2.537e-04\n",
      "Epoch=9200\tloss=4.977e-06\n",
      "Epoch=9300\tloss=4.646e-05\n",
      "Epoch=9400\tloss=4.768e-06\n",
      "Epoch=9500\tloss=1.605e-05\n",
      "Epoch=9600\tloss=4.656e-06\n",
      "Epoch=9700\tloss=3.441e-03\n",
      "Epoch=9800\tloss=4.752e-06\n",
      "Epoch=9900\tloss=6.980e-03\n",
      "CPU times: user 17.7 s, sys: 139 ms, total: 17.8 s\n",
      "Wall time: 17.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "epochs = 10_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 _ %(100) ==0:\n",
    "    print(f'Epoch={_}\\tloss={loss_fun(params,colloc,conds):.3e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 282
    },
    "id": "eWeNvDsdDEuI",
    "outputId": "32551eeb-25df-4d2e-8cae-82cc52b41ac5"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f70abfa5750>"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deVwW5f7/8dfFokTgBmguKZmiueVC5ilPUdqquaVp5W6LR1ssLbX6pp3TYqVZlmbmXqaW5ZJpHo8ds8xENNzAyjyauKRgookKeF+/PyB+hijINvfyfj4ePLjve+aeec898OHimplrjLUWERHxPH5OBxARkcJRARcR8VAq4CIiHkoFXETEQ6mAi4h4qIDSXFl4eLiNjIwszVWKiHi8jRs3JltrI3K/XqoFPDIykri4uNJcpYiIxzPG7MnrdXWhiIh4KBVwEREPpQIuIuKhSrUPPC8ZGRkkJSVx6tQpp6PIRQgKCqJGjRoEBgY6HUXEZzlewJOSkggNDSUyMhJjjNNxpACstaSkpJCUlMQVV1zhdBwRn+V4F8qpU6cICwtT8fYgxhjCwsL0X5OIwxwv4ICKtwfSPhNxnlsUcBERb5N+4gRbFy7kw0GD+PHLL0tkHSrgHi4mJibfi6MWLVpEQkJCzvPnn3+e//znPyUdTcQnnDp2jK0LF7J++vSc126qWJGQkBCadOlCr3ff5fN33imRdTt+EFNK3qJFi2jfvj0NGjQA4J///KfDiUQ8h7WWo7/+yqHEROrdfjsA4zp25L/ff8+O33/nfxkZuIDml1zCxv79AYiuU4dW/v40adGCJm3bEnXbbSWSzedb4CdOnKBdu3ZcffXVNGrUiPnz5/PVV1/RqVOnnHlWrlxJ586dAQgJCeGpp56iYcOGtG3bltjYWGJiYqhduzZLliwBYObMmXTs2JGYmBjq1q3LCy+8AGS1fN98882c5T777LO89dZb+eYBWLVqFc2aNaNx48b079+f06dPn7MtISEhOY8XLFhA3759+e6771iyZAlPPfUUTZs25ZdffqFv374sWLDggsuNjIxk1KhRNG/enMaNG7Njx44if9Yi7uJUaioHNm9mx/LlrJ82jeUvvsiHgwblTJ/QtSs3V6xIw6AgKvj5USkykhvuvDNn+rbERPYeO0bzypV5rnVrPho0iJmzZ+dMf33DBl75/nvunTiRhp07ExgcXCLb4X4t8JiYc1+75x4YNAjS0uCsDzFH375ZX8nJ0LXrX6etXn3B1X355ZdUq1aNL774AoDU1FTKlSvHoEGDOHz4MBEREcyYMYP+2X9ZT5w4wc0338zrr79O586dee6551i5ciUJCQn06dOHDh06ABAbG8u2bdsIDg7mmmuuoV27dvTv358uXbowZMgQXC4X8+bNIzY2Nt88p06dom/fvqxatYqoqCh69+7Nu+++y5AhQ/L7NLnuuuvo0KED7du3p2uuzya/5YaHh7Np0yYmTZrE2LFjmTp1ar7rEykprjNnAPDz9+dQQgI/r1nD8ZQUjh85wvHff+d4air93nqLcjVqsOyf/2TerFkcO3mS1FOnSD19mtSMDOJ37SK0Rg3+75ZbGLthwznr6PbKK5QtX55jqamcPnOG+mFh3BwWxhW1ahFZrx7W5cL4+THjp59Ke/Pz5H4FvJQ1btyYoUOHMnz4cNq3b8/f//53AHr16sWHH35Iv379WLduHbOz/7qWKVOG27P/jWrcuDFly5YlMDCQxo0bs3v37pzl3nLLLYSFhQHQpUsXvv32W4YMGUJYWBg//PADv/32G82aNcuZ50J5Nm/ezBVXXEFUVBQAffr0YeLEiQUq4Bfy448/XnC5Xbp0AaBFixZ89tlnRVqXSF6sy0XKzp0kxceTlJDAtd27E3HVVXw7aRKv/+tfpKSlkXLqFCkZGRyxlo3z5nF19+58+uKLDJo795zl3fHAA5SrUYM9iYl8s2cPoQEBlC9ThmqXXkqD4GDOZGQA0HnAAK688krKh4VRPjycilWrcllUFIGXXgrAcytX8lypfhKF434F/EIt5uDgC08PD8+3xZ1bVFQUmzZtYtmyZTz33HO0adOG559/nn79+nHXXXcRFBREt27dCAjI+qgCAwNzTqHz8/OjbNmyOY8zMzNzlpv7NLs/nz/wwAPMnDmTgwcP5rTq88vTsWPHAm3L2essjnO0/9w2f3//v2ybyMU6/NNP7Fi5ktrXXkv16GjWz5zJgIED+eX0ac7+SV1qLe1eeIGTx46x++hRwoKCaBQRQVhoKOGVKlGpZk0A2j36KCsaNiQ0PDzrq0oVQiIiqFi7NgD/mDuXf+RR4P903cMPc93DD5fkJpcK9yvgpWz//v1UqlSJnj17UqFChZxugmrVqlGtWjVefPHFQp2xsXLlSo4cOcIll1zCokWLmJ59hLpz5848//zzZGRk8NFHHxUoz9NPP83u3bvZuXMnderU4YMPPuDGG288571VqlQhMTGRevXqsXDhQkJDQwEIDQ3l+PHj58xfr169Ai1X5GId2bWLcb17E//zz8SnpLA/u/tjUvfu/GPePMIuv5wrK1TgturVqVmzJjWuvJIa9evT4I47ALhlxAg2jxhx3uXX/NvfqPm3v5XKtrgzny/gW7du5amnnsLPz4/AwEDefffdnGn3338/hw8f5qqrrrro5bZs2ZK7776bpKQkevbsSXR0NJDVBXPTTTdRoUIF/P39C5QnKCiIGTNm0K1bNzIzM7nmmmsYOHDgOe8dM2YM7du3JyIigujoaP744w8AevTowYMPPsiECRNyDl4CBV6uyPmcOHyY2Dlz+G7ZMtbGx3Pztdcy7PPPCbzkEl5fu5b6ZcvSpmZNrm7QgEatWtEs++SAOm3asPjgQYfTewFr7QW/gOnAIWDbWa9VAlYCP2d/r5jfcqy1tGjRwuaWkJBwzmvuYvDgwXbq1KkX/b4ZM2bYwYMH5zntzJkz9uqrr7Y//fRTUeM5zp33nZSMMxkZWQ9cLhtToYL1B0v2V8MyZez4Dh1y5j197JhDKb0PEGfzqKkFOY1wJnB7rtdGAKustXWBVdnPvUqLFi3YsmULPXv2LLZlJiQkUKdOHdq0aUPdunWLbbkiJSXz1CnWvf8+/2rblhsrVOD6ihWzJhhDdFQUI6+/nmWjR3Nk5062nT7NkMWLc95bJrsLT0qOySru+cxkTCSw1FrbKPv5j0CMtfaAMaYqsNpaWy+/5URHR9vcVw0mJiYWqotCnKd9591ea9+el774gmOAIetClVuvvpqX1q7F+Pn8JSSlyhiz0Vobnfv1wu6FKtbaA9mPDwJVLrDih4wxccaYuMOHDxdydSJSknavW8db3brRJiyM/dmNrJq1a3NP3bp8PGQIhxITiUtL4+V161S83UiRD2Jaa60x5rzNeGvtFGAKZLXAi7o+ESkev23fzsTBg1kcG8uWkycBaFCmDEmbN1MtOpoeEybQw+GMcmGF/VP6W3bXCdnfDxVfJBEpCeknTrBizBjWTp4MZPVvv/z111QICGDsnXfy8/LlbD91ipYDBjicVAqqsC3wJUAfYEz298UXnl1EnPD7nj0sHzuWxUuWsPzXXzkOdK5alesHDqR6ixYk//ILFbIvfhHPk28L3BgzF1gH1DPGJBljBpBVuG8xxvwMtM1+7rGMMQwdOjTn+dixYxk9ejQAo0ePJjg4mEOH/v8/GWcPGlVYZw8odT6rV6/mu+++y3k+efLknEv6Rc4nZefOnMe3NWzI/e+8w9d799I9KorPn32WOVu25ExX8fZs+bbArbX3nmdSm2LO4piyZcvy2WefMXLkSMLDw8+ZHh4ezrhx43j11VdLNdfq1asJCQnhuuuuA9BFNpIn63Lxw/z5LH7vPRbHxvLTyZMkHzpEcEQErzz/PJeWL0/L/v3x0w2ovY4OJwMBAQE89NBDjB8/Ps/p/fv3Z/78+Rw5cuS8yzhz5gx9+/alUaNGNG7cOGdZ8fHxtGrViiZNmtC5c2d+//33c94bGRlJcnIyAHFxccTExLB7924mT57M+PHjadq0Kd988w2jR49m7NixF1xuTEwMw4cPp2XLlkRFRfHNN98U6bMR9/bFCy9Qs0wZWtx3Hy9+/TXlAgP55513ciY9HYA2Tz9Nq4cfVvH2Um5XwGNiYs75mjRpEgBpaWl5Tp85cyYAycnJ50wrqMGDBzNnzhxSU1PPmRYSEkL//v3PGbv7bPHx8ezbt49t27axdetW+vXrB0Dv3r159dVX2bJlC40bN84ZGzw/kZGRDBw4kCeeeIL4+PicURL/dKHlZmZmEhsby5tvvlng9Yn7O5SYyIwHHqBL9er8+5VXAKhevz7XVK7MjH79OLhtG2tSUxn2xReEVq/ucFopDW5XwJ1Srlw5evfuzYQJE/Kc/thjjzFr1qw8B4UCqF27Nrt27eLRRx/lyy+/pFy5cqSmpnL06NGcAaL69OnDmjVripw1v+WePQzs2UPciuc5lZrKK3fcwXXlynFZgwb0nzaN2IMHSd67F4Cm3bvz2f799J0+nYiGDR1OK6XN7QazWn2B4WCDg4MvOD08PPyC0/MzZMgQmjdvntN6PluFChW47777mDhxYp7vrVixIps3b2bFihVMnjyZjz/++LxdMrkFBATgcrkADQPr69JPnOCbyZNJSUrinvHjKRMSwsR//5vLypZl1I030uHBB2naowcmj4HQxPeoBX6WSpUqcc899zBt2rQ8pz/55JO89957eRbF5ORkXC4Xd999Ny+++CKbNm2ifPnyVKxYMacf+nzDtUZGRrJx40YAPv3005zXzzcMbEGXK57hyK5dzBk8mO41axIREkLbYcN45p13sC4Xfv7+7EhKIi4tjVGrV9Ps/vtVvCWHCnguQ4cOzTmgmFt4eDidO3fO836U+/btIyYmhqZNm9KzZ09eye6jnDVrFk899RRNmjQhPj6e559//pz3jho1iscff5zo6Oi/DDF71113sXDhwpyDmGcryHLFfe1ctQqbPUb2iHbt6DlpEl8nJXFPVBSLR4xgy759OZesh1St6mRUcWMFGsyquGgwK++ifVdwmadP89377/P57Nl8vnkzP6anEzd7Ni169SJx6VKO//Yb0b1762wRydP5BrNyuz5wEW/zw/z5tL33Xo5YSyBwU6VKPHLTTdRs0QKAq9q3dzageCx1oYgUo/+tXcuEbt24JTycl2+9FYB6bdvSoXZtPnniCZJ//ZUVKSk8smABEQ0aOJxWPJ1btMCttefcBFjcW2l2vXmCsXfdxZxVq4jPHtWvfmAg5cuXByA4LIwZZ13eLlJcHG+BBwUFkZKSooLgQay1pKSkEBQU5HQUR1hr2bJwIZPvvz/ntW/i4ggOCGBc+/b8vGIFienpDP7kEwdTii9w/CBmRkYGSUlJxXL+s5SeoKAgatSoQaAPHXTbsXw5H778Mh/HxvJzejr+wP6tW6ncqBEZaWkEBgc7HVG8lNsexAwMDOSKK65wOobIBc197DHue/tt/ICbK1ZkaJcudBo5ksqNGgGoeIsjHO9CEXE3p1JTmfvYY9xRuTIz+vcH4JbBgxnfoQP74uNZeeQID8+dS5UmTRxOKr7O8Ra4iLvYsXw5U557jlk//MARa6np759zsU14vXp/ueO6iDtQARef5srMxC8g69egV9eubE5Lo3ONGjz0yCPc9OSTurBG3JoKuPikg1u3MumRR5i1di2bd+6kQmQkU6dPp2rDhjn92iLuTgVcfEr8ggW8OXIkc3fuJAO4q0oVUpOSqBAZydXduzsdT+SiqICLz9i1Zg3NunXjUuChxo15bOxY6mZfLSniiXQWingt63KxdNQoXrrlFgBq33ADcx95hL2//MLbW7aoeIvHUwtcvI4rM5PPRo7kxYkT2XzyJHUCAxn6++8EVaxIj7ffdjqeSLFRC1y8yqaPPqJJSAjdxo7lZGYmMwcMIOHoUYIqVnQ6mkixUwtcvEJaSgrBYWFUrluXAGP4aNAg7nnjDfyzby8n4o1UwMWjbV6wgGcGD+ZkejqrUlKocc01/HDiRM7dbES8mX7KxSMd2LyZfnXr0qxbN747dIjbW7XClX2vUhVv8RVqgYvH+WbiRO585BHSgWHR0YycN4+KV17pdCyRUqemingE63JxKCEBgObdunFP3bpsX7mS1zZsUPEWn6UWuLi9HcuX82jPnuw+fpztKSlcWrky0376yelYIo5TC1zc1qmjRxl1441cfeedxB05wmN33ZUz8JSIqAUubmpvbCxtWrfm54wM7o+MZNznn1NFg0yJ/IVa4OJWzqSnA1CtWTOiq1Xj3y+/zIf/+5+Kt0geVMDFLViXi9kPP8xVISEcTkjAPzCQj3bv5paRI52OJuK2ilTAjTFPGGO2G2O2GWPmGmN88zblUiT7N22iQ9Wq9JkyhcpBQZxITnY6kohHKHQBN8ZUBx4Doq21jQB/oEdxBRPv92eru2F0NKsOHeLNTp34OjmZyBtucDqaiEco6kHMAOASY0wGEAzsL3ok8RXGz4+ly5bRKDSU6QsWUDd72FcRKZhCt8CttfuAscCvwAEg1Vr779zzGWMeMsbEGWPiDh8+XPik4hWsy8XsgQPZsWwZANM3bODr5GQVb5FCKEoXSkWgI3AFUA241BjTM/d81top1tpoa210RERE4ZOKxzuyaxf31KpFn/feY9KIEQCEXHaZbhwsUkhFOYjZFviftfawtTYD+Ay4rnhiibf56o03aBIVxaKkJMbcdhvj4+KcjiTi8YpSwH8FWhljgo0xBmgDJBZPLPEmi599lrZDhxLi58f3s2cz/Msv8S9TxulYIh6vKH3g64EFwCZga/ayphRTLvECmadPA3Dr0KGMuvFGNu3bR4tevRxOJeI9inQeuLV2lLW2vrW2kbW2l7X2dHEFE89lXS4m3XsvzStU4I+DB7mkUiVGrV5NsI6BiBQrXYkpxSp1717uqVWLwfPmUf3SSzl19KjTkUS8lgq4FJu4OXNoXrs2C5OSeO2OO/ji4EHC69d3OpaI19JohFIsrLU8/cgjZFjLmnff5bqBA52OJOL1VMClSI7u2cOZjAzC6tThg5UrCSpXjrCoKKdjifgEdaFIocV9+CHN69RhwN//DkD16GgVb5FSpAIuhTK1b1+u79WLM9YyYvRop+OI+CR1ochFOX3sGI9eey3v79jBLZUqMXfdOrW6RRyiFrhclNS9e1n288+MbNWK5QcOqHiLOEgFXApkw+zZZJ46ReWGDdn+yy+8vG6dLocXcZgKuFyQdbkY27Ejrfr04Y0uXQAoX6uWw6lEBNQHLhdw8sgR+jdvzrw9e7i7WjX+8f77TkcSkbOoBS552r9pEzdcfjnz9+zh5Vtv5ZO9ewmtXt3pWCJyFhVwydP+7dvZc/Iki555hpErVmD89KMi4m70Wyl/Ef/xxwBE9+rF/w4epMNLLzmcSETORwVcAHCdOcP/3Xgjzbp3Z8mzzwJwaeXKDqcSkQvRQUwhLTmZXk2b8tm+ffSvU4fbs+9XKSLuTS1wH3coIYGbIiNZuG8f49q3Z+qPP1ImNNTpWCJSACrgPm7N9OlsO3GChSNG8OTnn+tgpYgHMdbaUltZdHS0jdPdyN3C0T17qJB9Qc6B+HiqNm3qcCIROR9jzEZrbXTu19Xc8kGfDBtGZGQk6957D0DFW8RD6SCmD7HW8kbnzgxbvJjrQkKIiolxOpKIFIEKuI9wZWYyrFUrxm/cSNdq1Zi9dSuXVKrkdCwRKQJ1ofiIuY8/zviNG3m0SRPm796t4i3iBdQC9xH3vvUWwSEhdHrlFZ1pIuIl9JvsxVL37uXeyEh2f/stfgEBdH71VRVvES+i32YvdSghgZvq1WPBnj3EL13qdBwRKQEq4F5o7/r1/L1pU3acPMmS0aPpNGaM05FEpASoD9zL/LpuHTE33EBKZiYrJ07k+kGDnI4kIiVELXAvU756daLKl2fVzJkq3iJeTi1wL7Hr66+5rGFDytesyZfJyU7HEZFSoBa4F0hYupTrb76ZB6PPGSpBRLyYCriH2/rZZ8R06ADAsxMmOJxGREqTulA8WMKSJdzctStljeGrL74g6vbbnY4kIqWoSC1wY0wFY8wCY8wOY0yiMeZvxRVMLsyVmUmP7t0JMIb/Ll+u4i3ig4raAn8L+NJa29UYUwYILoZMUgB+AQF8NHcufn5+1L31VqfjiIgDCt0CN8aUB24ApgFYa9OttUeLK5jkbW9sLOM7dcK6XDTq1IkG2f3fIuJ7itKFcgVwGJhhjPnBGDPVGHNp7pmMMQ8ZY+KMMXGHDx8uwurkwObNtGndmtGLF5MUG+t0HBFxWFEKeADQHHjXWtsMOAGccztza+0Ua220tTY6IiKiCKvzbSk//0zba69lf0YGyydN4vJWrZyOJCIOK0oBTwKSrLXrs58vIKugSzE7cegQ7Zo145fTp/l87Fiu+8c/nI4kIm6g0AXcWnsQ2GuMqZf9UhsgoVhSyV98+/77/HDiBPNGjOCmoUOdjiMibqKoZ6E8CszJPgNlF9Cv6JEkt9uefZZfbr2VGtdc43QUEXEjRToP3Fobn92/3cRa28la+3txBfN11uVieKtWLB45EkDFW0TOoUvp3dSYO+/ktfXr+XrFCqejiIibUgF3QzMffJBnVqzgvlq1GKvTBUXkPFTA3cxX48bx4NSptK1UiRnbtuEXoOFqRCRvKuBu5j+ffkpU2bIs2LSJMiEhTscRETemAu5mXvr2W9bt3En5WrWcjiIibk4F3A2cPHKEbpdfztZPP8X4+VGuRg2nI4mIB1ABd5grM5NeTZrwaVISv+iApYhcBB0hc9hzN97Ip/v2Ma59ezq9+qrTcUTEg6gF7qD5jz/OK999x4P16vHE4sVOxxERD6MC7hDrcjFr7lyuDw3lnbg4jJ92hYhcHHWhOMT4+bF4926OHzig0wVFpFDU7Ctl6SdOMDQ6mkMJCQQGB1PpyiudjiQiHkoFvJQ9fu21vLFxI2umTXM6ioh4OBXwUvR+795M3r6d4S1b0nXcOKfjiIiHUwEvJRs//JBHPviA28LCeOmbb5yOIyJeQAW8lAx79FGq+Pvz4dq1+Jcp43QcEfECOgullCyIjWX/li2E16uX/8wiIgWgFngJ+89rr5H+xx+E1a1L47vvdjqOiHgRFfAS9N/x47lt+HBebtfO6Sgi4oVUwEvIvk2b6DFsGFFlyjBs3jyn44iIF1IfeAk4k57O/TffzB8uF//95BNCqlZ1OpKIeCEV8BIwpl07vk5NZeaAATTo0MHpOCLipVTAS8BdDz/MybQ0ek+Z4nQUEfFixlpbaiuLjo62cXFxpba+0paRlkZgcLDTMUTEyxhjNlpro3O/roOYxcS6XNxfvz4PXXUV1uVyOo6I+AAV8GIyfcAAPtm7l9q1amlsbxEpFao0xWDHsmU8NnMmbStW5OmlS52OIyI+QgW8iNL/+IP7unYl2Bhmr1qFX4COC4tI6VC1KaLtn3/OrpMnmTl8OFWbNXM6joj4EBXwImp2773satlSd9YRkVKnLpRCSktOZtaDD2JdLhVvEXGECnghjWjblr5Tp7Jpzhyno4iIjypyATfG+BtjfjDG+MzpF/957TXe3ryZIU2b0qJXL6fjiIiPKo4W+ONAYjEsxyMc3bOHfs88Q/0yZXh51Sqn44iIDytSATfG1ADaAVOLJ477G9KmDQfOnGH2e+9xSaVKTscRER9W1LNQ3gSeBkKLIYtHuG/AAJpt2MA1ffs6HUVEfFyhC7gxpj1wyFq70RgTc4H5HgIeAqhZs2ZhV+c463Jh/Py4deRIbnU6jIgIRetCuR7oYIzZDcwDbjbGfJh7JmvtFGtttLU2OiIiogirc9aQ5s15rnVrDVQlIm6j0AXcWjvSWlvDWhsJ9AC+stb2LLZkbuTrCROYsHkzJ9LSNFCViLgNVaN8pKWk8MCwYdQOCODFL790Oo6ISI5iuZTeWrsaWF0cy3I3o+64g50ZGax6/XUurVzZ6TgiIjnUAr+A37Zt450NG3iwXj1uHjbM6TgiIn+hwawuoEqjRsQuWEDN6HPuZCQi4ji1wM8jacMGABrffTfla9VyOI2IyLlUwPPw6/ffU79lS97o2NHpKCIi56UCnofHu3TBBXR56imno4iInJf6wHNZ8uyzLDpwgDG3305k69ZOxxEROS9jrS21lUVHR9u4uLhSW9/FOnHoEA2qVSM0IIAfjhwhMDjY6UgiIhhjNlprzzmbQi3ws2xdsoQ/XC7mjBun4i0ibk8F/CytHniAPe3bE3LZZU5HERHJlw5ikjXS4KKRI3FlZqp4i4jHUAEH5j76KJ3HjOGz4cOdjiIiUmA+fxDzj4MHqVe9OtWCglifmopfgHqVRMS96CDmebzUuTP7XS4+festFW8R8Sg+3YXy88qVvPH99/S58kpaPfCA03FERC6KTxfwI7/+SuPgYF5ZsMDpKCIiF83n+8D/vNeliIi7Ol8fuE9WrvQ//uD19u05ceiQireIeCyfrF5vdu/O0198wdqpU52OIiJSaD5XwA9t386Ly5ZxV5Uq3PrMM07HEREpNJ8r4C90704a8PqMGU5HEREpEp8q4IlLl/Le9u0MbNyYenfc4XQcEZEi8akCbvz8aHfZZYyaP9/pKCIiReZTlx7Wv/NOFh844HQMEZFi4RMt8DMZGTzXujV71693OoqISLHxiQL+4aBBvLR2Ld9+8IHTUUREio3XX4mZlpxMVJUqVL/kEr5PTcX4+5fq+kVEispnRyMc1707+1wu5r36qoq3iHgVr+5COZyYyGtffUXnqlVpPXiw03FERIqVVxdwrOXe+vV5acoUp5OIiBQ7r+5CiWjQgCmJiU7HEBEpEV7bAh/fqRPrp01zOoaISInxygKeuHQpwxYvZt7EiU5HEREpMV5ZwP9v4ECCgWfmzHE6iohIiSl0ATfGXG6M+a8xJsEYs90Y83hxBiusDbNm8em+fQyLiSHiqqucjiMiUmKKchAzExhqrd1kjAkFNhpjVlprE4opW6E88+SThBvDk2p9i4iXK3QBt9YeAA5kPz5ujEkEqgOOFfAz6enENG9OtypVCK1WzakYIiKlolgupTfGRAJrgHofGJcAAAZeSURBVEbW2mO5pj0EPARQs2bNFnv27Cny+kREfEmJ3dTYGBMCfAoMyV28Aay1U6y10dba6IiIiKKu7rzWvP02Hz/xBK7MzBJbh4iIOylSC9wYEwgsBVZYa9/Ib/6SGszKuly0DA0lJT2dH1NTCQwOLvZ1iIg4pdgHszLGGGAakFiQ4l2Svhg9mri0NKb17aviLSI+o9AtcGNMa+AbYCvgyn75GWvtsvO9pyRa4Nbl4prQUH5PT2eHWt8i4oWKvQVurf0WMEVKVQyWjhrFxrQ0pvfrp+ItIj7F46/E9A8I4NawMHq+847TUURESpXX35FHRMTTldhphE6xLhfT+/UjLTnZ6SgiIo7w2PHAF40cyYCZMykTFETPd991Oo6ISKnzyC4U63LRPCSEE5mZJBw7RkBQUDGkExFxT151U+Pl//oX8SdPMmPAABVvEfFZHtcCty4XrStUICktjZ3HjunUQRHxel5zEPP4/v34G8PTXbqoeIuIT/O4LpRyNWqwJjVVg1aJiM/zqBb4/9as4eCWLQD4BXjc3x4RkWLlUQV8aI8eNGvWjIy0NKejiIg4zmMKeMKSJSw8cIAHr79efd8iInhQAR/zxBNcCjw+bZrTUURE3IJHFPBdq1fz0a5dPNyiBWF16zodR0TELXhEAf9qxgzKAEOnTHE6ioiI2/CIAv7ArFns/eknqjVv7nQUERG34REFHFDXiYhILh5TwEVE5K9UwEVEPJQKuIiIh1IBFxHxUCrgIiIeSgVcRMRDqYCLiHgoFXAREQ9VqrdUM8YcBvYU8u3hQHIxxiltnp4fPH8bPD0/eP42eHp+cGYballrI3K/WKoFvCiMMXF53RPOU3h6fvD8bfD0/OD52+Dp+cG9tkFdKCIiHkoFXETEQ3lSAff0sWQ9PT94/jZ4en7w/G3w9PzgRtvgMX3gIiLyV57UAhcRkbOogIuIeCi3K+DGmNuNMT8aY3YaY0bkMb2sMWZ+9vT1xpjI0k95fgXI39cYc9gYE5/99YATOc/HGDPdGHPIGLPtPNONMWZC9vZtMca41W2SCpA/xhiTetbn/3xpZ8yPMeZyY8x/jTEJxpjtxpjH85jHbfdDAfO79X4wxgQZY2KNMZuzt+GFPOZxvhZZa93mC/AHfgFqA2WAzUCDXPMMAiZnP+4BzHc690Xm7wu843TWC2zDDUBzYNt5pt8JLAcM0ApY73Tmi8wfAyx1Omc+21AVaJ79OBT4KY+fI7fdDwXM79b7IftzDcl+HAisB1rlmsfxWuRuLfCWwE5r7S5rbTowD+iYa56OwKzsxwuANsYYU4oZL6Qg+d2atXYNcOQCs3QEZtss3wMVjDFVSydd/gqQ3+1Zaw9YazdlPz4OJALVc83mtvuhgPndWvbn+kf208Dsr9xnfDhei9ytgFcH9p71PIlzd3zOPNbaTCAVCCuVdPkrSH6Au7P/7V1gjLm8dKIVm4Juozv7W/a/xsuNMQ2dDnMh2f+WNyOrBXg2j9gPF8gPbr4fjDH+xph44BCw0lp73n3gVC1ytwLuCz4HIq21TYCV/P+/4FI6NpE1rsTVwNvAIofznJcxJgT4FBhirT3mdJ6LlU9+t98P1toz1tqmQA2gpTGmkdOZcnO3Ar4POLtFWiP7tTznMcYEAOWBlFJJl79881trU6y1p7OfTgValFK24lKQfeS2rLXH/vzX2Fq7DAg0xoQ7HOscxphAsorfHGvtZ3nM4tb7Ib/8nrIfAKy1R4H/ArfnmuR4LXK3Ar4BqGuMucIYU4asAwNLcs2zBOiT/bgr8JXNPorgBvLNn6ufsgNZ/YOeZAnQO/ssiFZAqrX2gNOhCsoYc9mf/ZTGmJZk/Q64SwMAyDrDBJgGJFpr3zjPbG67HwqS3933gzEmwhhTIfvxJcAtwI5cszleiwJKc2X5sdZmGmMeAVaQdUbHdGvtdmPMP4E4a+0Ssn4wPjDG7CTrYFUP5xL/VQHzP2aM6QBkkpW/r2OB82CMmUvWGQLhxpgkYBRZB3Cw1k4GlpF1BsROIA3o50zSvBUgf1fgH8aYTOAk0MONGgB/uh7oBWzN7oMFeAaoCR6xHwqS3933Q1VgljHGn6w/Lh9ba5e6Wy3SpfQiIh7K3bpQRESkgFTARUQ8lAq4iIiHUgEXEfFQKuAiIh5KBVxExEOpgIuIeKj/B6nhTHAnPE76AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "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": {
   "authorship_tag": "ABX9TyPlDK/9ZvMulH91+B+B32BC",
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "[1] ODE-PINN.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3.8.9 64-bit",
   "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.8.9"
  },
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}