DL_Kolmogorov.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0.0\n",
      "Num GPUs Available:  0\n"
     ]
    }
   ],
   "source": [
    "from __future__ import absolute_import, division, print_function, unicode_literals\n",
    "\n",
    "try:\n",
    "  # %tensorflow_version only exists in Colab.\n",
    "  %tensorflow_version 2.x\n",
    "except Exception:\n",
    "  pass\n",
    "\n",
    "# TensorFlow and tf.keras\n",
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "tf.random.set_seed(1234)\n",
    "\n",
    "# Helper libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact\n",
    "import datetime, os, sys\n",
    "\n",
    "print(tf.__version__)\n",
    "tf.print(\"Num GPUs Available: \", \n",
    "         len(tf.config.experimental.list_physical_devices('GPU')))\n",
    "\n",
    "# using files from github repository\n",
    "# !rm -r oberwolfach_workshop\n",
    "!git clone https://github.com/juliusberner/oberwolfach_workshop.git\n",
    "    \n",
    "sys.path.append('./oberwolfach_workshop/helper_functions')\n",
    "from NN_constructions import Affine_Comb, Square, R, P\n",
    "from helper_functions import plot_1d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction to Kolmogorov PDEs\n",
    "\n",
    "In this notebook we present a deep learning based method for solving high-dimensional **Kolmogorov equations (PDEs)**, which are defined by\n",
    "\n",
    "$$ \n",
    "\\begin{cases}\n",
    "\\frac{\\partial f}{\\partial t}(t,x) \n",
    "= \n",
    "\\frac12 \\operatorname{Trace}\\!\\big(\n",
    "\\sigma(x)[\\sigma(x)]^{*}\\operatorname{Hess}_x f(t,x)\\big)\n",
    "+ \n",
    "\\big\\langle\n",
    "\\mu(x),\\nabla_x f(t,x) \n",
    "\\big\\rangle  \n",
    "\\\\ \n",
    "f(0,x) = \\varphi(x). \n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Special cases include, e.g. the heat equation and the Black-Scholes equation from financial engineering (pricing of options).\n",
    "\n",
    "Our goal will be to compute the evolution of the solution for time $[0,\\tau]$ on some hypercube $[u,v]^d$, i.e. the function\n",
    "\n",
    "$$\\color{blue}{f^{sol}:=f(t,x)|_{x\\in[u,v]^d, \\ t\\in[0,T]}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple example\n",
    "\n",
    "Let\n",
    "* $\\varphi(x):=\\|x\\|^2_2$\n",
    "* $\\mu(x):=0$\n",
    "* $\\sigma(x):=\\sigma \\in\\mathbb{R}^{d\\times d}$\n",
    "\n",
    "Then it holds that (check!)\n",
    "\n",
    "$$ f(t,x) = \\|x\\|_2^2 + t * \\operatorname{Trace}(\\sigma \\sigma^*).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 1 #space dimension\n",
    "sigma = 0.05*tf.ones([d,d])+0.5*tf.eye(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi(_x):\n",
    "    return tf.reduce_sum(_x**2, axis=1, keepdims=True)\n",
    "\n",
    "def solution(_t, _x):\n",
    "    return phi(_x)+_t*tf.linalg.trace(sigma@tf.transpose(sigma))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAE/CAYAAABin0ZUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hUZeL28e8zkx5SCCWUJPTeMYLYFVSs2H72jn3t3a32teyqq2tX7Ip9wcWCBSs1SIcAIbRQ0oBU0mae948ZeLNICWSSk5m5P9eVK5k5JzP3ycmQm3OeeY6x1iIiIiIiB8bldAARERGRYKYyJSIiItIIKlMiIiIijaAyJSIiItIIKlMiIiIijaAyJSIiItIIKlMiEnaMMV2NMdYYE3GA33+hMWZqoHOJSHBSmRKRFssYc5kx5heHM/yueFlr37XWHu9kLhFpOVSmRERERBpBZUpEGswY08kY84kxptAYs9oYc1O9ZV8YY/5Z7/ZEY8wE/9c9jDHfG2OKjTFFxph3jTHJ9dZNN8Z86n/cYmPMv40x/YAXgVHGmHJjzLY9ZLrMGJNrjCnzZ7rQf7/LGPNnY8xaY0yBMeYtY0zSHh5jjTFmTL3b9xlj3vHf/Mn/eZs/x6hdj5gZYw41xswxxpT4Px9ab9kPxpgHjTG/+jNONca0bfhPXURaOpUpEWkQY4wL+BxYAHQGRgO3GGNO8K9yBXCxMeZYf6EZAdy849uBvwOdgH5AOnCf/3HdwH+BtUBX/2NPtNYuA64FZlhrW1lrd5avepnigWeAE621CcChwHz/4sv8H8cA3YFWwL8PYNOP9H9O9ueYsUuGFGCKP0cb4ElgijGmTb3VLgAuB9oDUcAdB5BDRFoolSkRaaiDgXbW2gestTXW2lzgFeA8AGvtZuA64E3gX8Al1toy/7Ica+031tpqa20hvsJxlP9xR+ArWXdaayustVXW2v0ZJ+UFBhpjYq21m6y1S/z3Xwg8aa3NtdaWA/cC5x3ooPO9OBlYaa1921pbZ619H8gGTq23zuvW2hXW2u3Ah8DQAGcQEQepTIlIQ3UBOhljtu34AP4IpNZb53PADSyvX4iMMan+034bjDGlwDvAjlNd6cBaa23d/gay1lYA5+I7grXJGDPFGNPXv7gTvqNdO6wFInbJGwi7Ps+O5+pc7/bmel9X4jtKJiIhQmVKRBpqPbDaWptc7yPBWntSvXUeBpYBHY0x59e7/xHAAoOstYnARfhO/e143Iw9HDGy+wplrf3aWnsc0BHfEaFX/Is24iuAO2QAdUD+bh6mAoird7vDfmTY9Xl2PNeGfXyfiIQIlSkRaajZQJkx5m5jTKwxxm2MGWiMORjAGHMkvnFBlwCXAs8aY3YcnUkAyoES/3137vK4m4BHjTHxxpgYY8xh/mX5QJoxJmp3gfxHvMb5x05V+5/D61/8PnCrMaabMaYVvkL3wR6OgM3Hdwow0hiTCZxdb1mh/zG77+Hn8gXQ2xhzgTEmwhhzLtAf3zgwEQkDKlMi0iDWWg9wCr7xPquBIuBVIMkYkwi8Bdxgrd1grf0ZeA143RhjgPuB4UAJvsHan+7yuKcCPYF1QB6+U3cA3wNLgM3GmKLdxHIBt+E7OrQF3zis6/zLJgBv43s33mqgCrhxD5v3F6AHsNWf9b16+SrxHXH71X9685Bdfi7F/p/L7UAxcBdwirV2d3lFJAQZa/d5FF1ERERE9kBHpkREREQaQWVKREREpBFUpkREREQaQWVKREREpBFUpkREREQaIdCXVWiwtm3b2q5duzr19CIiIiINNnfu3CJrbbvdLXOsTHXt2pWsrCynnl5ERESkwYwxu142aied5hMRERFpBJUpERERkUZQmRIRERFpBJUpERERkUZQmRIRERFpBJUpERERkUZQmRIRERFpBJUpERERkUZQmRIRERFphJAtU7+t28r5L89kZm6x01FERESkCXi8lvNensGEX1ZTU+d1LEfIlqmfVxQxI7eYf3y9HGut03FEREQkwCbN38DM3C28MX0NxjiXI2TL1PgjupESH0XW2q38sLzQ6TgiIiISQLUeL09/uxKAm0b3ItLtXKUJ2TLVKjqC647qAcA/pi7H69XRKRERkVDx8dw81m2ppHu7eE4f2snRLCFbpgAuHtWF1MRolmws5cvFm52OIyIiIgFQXefh2e98R6VuHdObCAePSkGIl6mYSDc3HtsLgCe/WU6dx7nBaSIiIhIY785cx8aSKvp2SODkQR2djhPaZQrgnMx0MlLiWFVYwWfzNjgdR0RERBqhorqO56blAHDH8X1wuRwcee4X8mUqKsLFrcf5jk49/e1Kqus8DicSERGRAzXhl9UUV9QwPCOZ0f3aOx0HCIMyBXDakM70at+KDdu288Gc9U7HERERkQOwtaKGl3/KBeDOE/pinJwPoZ6wKFNul+H24/sA8Oz3OWyv0dEpERGRYPPiT6soq67jiF5tGdWjjdNxdgqLMgVwwoBUBqclUVhWzZsz1jgdR0RERPZDfmkVb/y6BoA7T+jjbJhdhE2ZMsZwh//o1Is/rqK0qtbhRCIiItJQz3y3kuo6LycO7MDgtGSn4/yPsClTAEf0asuIbilsq6zl1Z9XOx1HREREGmBtcQUfzFmPy8Dtx/d2Os7vhFWZMsbsPDT42s+5bKmocTiRiIiI7MtT36ygzms5c3gaPdsnOB3nd8KqTAEc3DWFo/u0o6LGwws/5DgdR0RERPYie3MpkxZsJMrt4pYxvZyOs1thV6aAnWOn3pqxls0lVQ6nERERkT35x9crsBYuGJlBWus4p+PsVliWqYGdkzhpUAeq67w88/1Kp+OIiIjIbsxdu5Vvl+UTF+XmD8f0dDrOHoVlmQK47bg+uAx8MGc9uYXlTscRERGReqy1PPF1NgBXHNaNdgnRDifas7AtUz3bt+KczHQ8Xss/pi53Oo6IiIjU80tOETNzt5AUG8lVR3Z3Os5ehW2ZArhlTG+iI1x8sWgz89ZtdTqOiIiIsOOolO9Ax3VH9yApNtLhRHsX1mWqQ1IMVxzeDYBHv8zGWutwIhEREfly8WYW5pXQPiGaS0d1dTrOPoV1mQK49qgeJMdFMmv1Fn5YUeh0HBERkbBW6/HuPCp10+hexEa5HU60b2FfppJiI7nB/w6Bx77MxuPV0SkRERGnTJy9jtVFFXRvG8+5B6c7HadB9lmmjDETjDEFxpjFe1h+oTFmoTFmkTFmujFmSOBjNq2LDulC5+RYsjeXMWn+BqfjiIiIhKXy6jr+9Z1vyqK7xvYh0h0cx3wakvINYOxelq8GjrLWDgIeBF4OQK5mFRPp5rbjfNf6+efUFVTVehxOJCIiEn5e+SmXovIahmckc8KADk7HabB9lilr7U/Alr0sn26t3fFWuJlAWoCyNavTh3Wmb4cENmzbzjsz1zodR0REJKwUlFXxys+5ANx7Uj+MMQ4narhAHz8bD3wZ4MdsFm6X4e6xfQH497QcSrbXOpxIREQkfPzr25VU1ng4rn8qB3dNcTrOfglYmTLGHIOvTN29l3WuNsZkGWOyCgtb3jvnju7TjpHdUthWWctLP65yOo6IiEhYWFVYzsQ563EZuHtsH6fj7LeAlCljzGDgVWCctbZ4T+tZa1+21mZaazPbtWsXiKcOKGMM95zoOzo14dfVugiyiIhIM3jiq+V4vJZzD06nZ/sEp+Pst0aXKWNMBvApcLG1dkXjIzlrWEZrThzYgapaL//6Lug3R0REpEWbu3YrXy3ZTEyki1vG9HY6zgFpyNQI7wMzgD7GmDxjzHhjzLXGmGv9q/wVaAM8b4yZb4zJasK8zeKOE/rgdhk+mLOenAJdBFlERKQpWGt59MtlAFx5eHdSE2McTnRgIva1grX2/H0svxK4MmCJWoAe7Vpx7sHpvDdrHY9/lc3Ll2Q6HUlERCTkfLM0nzlrtpISH8U1R7XsixnvTXDMhuWAW0b3IjbSzdSl+czK3eMwMBERETkAdR4vj32VDcCNx/YkIaZlX8x4b1Sm9qB9YgxXH+lryQ9/sQyvLjMjIiISMB/NzWNVYQUZKXFcOLKL03EaRWVqL645qjvtE6JZmFfC5AUbnY4jIiISEipr6njqG9+bvO48oQ9REcFdR4I7fROLi4rgjuN98108/lW2LjMjIiISAK/9vJqCsmoGpyVx8qCOTsdpNJWpfTjroDT6dkhgY0kVE35d7XQcERGRoFZQWsUL/omx7z2xHy5X8Fw2Zk9UpvbB7TL8+eT+ADw/bRVF5dUOJxIREQle/5y6gsoaD8f3T2VUjzZOxwkIlakGOLxXW47p047y6jqe/lYTeYqIiByIJRtL+HDueiJchntP6ud0nIBRmWqgP57UD5eB92evJ6egzOk4IiIiQcVay8NTlmEtXDKqK93axjsdKWBUphqoV2oC543IwOO1PPJFttNxREREgsp3ywqYvqqYpNhIbhrd0+k4AaUytR9uHdOb+Cg332cX8GtOkdNxREREgkKtx8sjX/guG3Pz6F4kx0U5nCiwVKb2Q7uEaK4/xtemH5qyDI8m8hQREdmnd2euJbeogm5t47nokOCeoHN3VKb20/jDu9EpKYZlm0r59Lc8p+OIiIi0aCWVtTz93UoA7j2xb9BP0Lk7obdFTSwm0s2dY30Tef5j6nIqa+ocTiQiItJyPfv9SrZV1jKqexuO65/qdJwmoTJ1AMYN6czgtCTyS6t55SdN5CkiIrI7a4oqeHPGGoyBP53cD2OCf4LO3VGZOgAul+FP/vkxXvxxFZtLqhxOJCIi0vI8+mU2tR7L2cPTGNg5yek4TUZl6gCN7N6GsQM6sL3Ww+NfaaoEERGR+mbmFvPVks3ERrq544Q+TsdpUipTjfDHk/oRFeHi03kb+G3dVqfjiIiItAher+WhKUsBuPaoHqQmxjicqGmpTDVCRps4rjqiGwD3f74Ur6ZKEBER4bN5G1i8oZQOiTFcdWQ3p+M0OZWpRrr+6J60T4hmwfptfDZvg9NxREREHFVeXcdj/uEvd57Qh7ioCIcTNT2VqUaKj47gnhP7AvDYV9mUV2uqBBERCV/PTcuhoKyaIenJnDGss9NxmoXKVACcPrQzQ9OTKSir5vlpOU7HERERccSaogpe+9k3ZdB9p/bH5QrNqRB2pTIVAC6X4W+n9gfg1Z9Xs6640uFEIiIize+hKcuo8Xg5a3gawzJaOx2n2ahMBciwjNacOawzNR4vD3+x1Ok4IiIizeqnFYV8uyyf+Cg3d48N7akQdqUyFUB3je1LXJSbr5fk82tOkdNxREREmkWtx8sD//UdSLhxdC/ah/hUCLtSmQqgDkkx/OGYngA88PlS6jxehxOJiIg0vbdmrCWnoJxubeO5/LCuTsdpdipTATb+8G6ktY5leX4Z789e53QcERGRJlVcXs3T364A4C+n9CM6wu1wouanMhVgMZFu/nyy77p9//xmBdsqaxxOJCIi0nT+MXU5ZVV1HN2nHcf2TXU6jiNUpprACQM6MKp7G7ZV1vL0tyudjiMiItIkFm8oYeKc9US4DH85pb/TcRyjMtUEjDH89dT+uAy8PXMt2ZtLnY4kIiISUNZa7pu8BGvh8sO60qNdK6cjOUZlqon065jIhSO74PFa/vqfJVir6/aJiEjomLxgI1lrt9K2VRQ3ju7ldBxHqUw1oTuO70NKfBSz12xh0vyNTscREREJiPLqOh75YhkAd53Ql8SYSIcTOUtlqgklxUVyz1jfdfse/mIZZVW1DicSERFpvGe+W0l+aTVD05M5+6A0p+M4TmWqiZ19UBpD05MpLKvmXxqMLiIiQW5FfhkTflmNMfDguIFhc/29vVGZamIul+HBcQMxBl6fvoblm8ucjiQiInJArLX8ddJi6ryWC0dmMCgtyelILYLKVDMYlJbEBSMyfIPRJy3WYHQREQlKkxdsZGbuFlLio7jj+PC6/t7eqEw1kztP6EPruEhmrd7C5AUajC4iIsGlrKqWh6f4Bp3fM7YvyXFRDidqOVSmmklyXBR3+wejP/LFMsqr6xxOJCIi0nDPfLeSgjINOt8dlalmdE5mOkPSk8kvreaZ7zQYXUREgsPyzWVM+HUNxsBDp2vQ+a5UppqRbzD6AIyBCb+sZmW+BqOLiEjLtmPQucdruWhkFwZ21qDzXalMNbPBacmcd3AGdV7LXydpZnQREWnZJi/YyKzVGnS+N/ssU8aYCcaYAmPM4j0sN8aYZ4wxOcaYhcaY4YGPGVruOqEPyXGRzMgt1mB0ERFpscqqanlox6DzE/uSFBfeM53vSUOOTL0BjN3L8hOBXv6Pq4EXGh8rtLWOj9o5M/qD/11GyXbNjC4iIi3PU9+spLCsmuEZyZw9XIPO92SfZcpa+xOwZS+rjAPesj4zgWRjTMdABQxV52Smk9mlNUXl1TzxdbbTcURERP7H4g0lvDF9NS4DD2im870KxJipzsD6erfz/PfJXrhchofOGEiEy/DurHXMW7fV6UgiIiIAeLyWP322CK+FSw/tqkHn+9CsA9CNMVcbY7KMMVmFhYXN+dQtUt8OiYw/ohvWwp8+W0ydx+t0JBEREd6dtZYFeSV0SIzhdg0636dAlKkNQHq922n++37HWvuytTbTWpvZrl27ADx18Lt5dC86J8eydFMpb85Y63QcEREJc/mlVTzx1XIA7jutP62iIxxO1PIFokxNBi7xv6vvEKDEWrspAI8bFuKiInhg3AAAnpy6nE0l2x1OJCIi4eyB/y6lrLqO0X3bc8KADk7HCQoNmRrhfWAG0McYk2eMGW+MudYYc61/lS+AXCAHeAW4vsnShqjR/VI5YUAqFTUe7p+81Ok4IiISpqYtL2DKwk3ERrq5f9wAjNGg84bY57E7a+35+1hugT8ELFGYuu+0Afyysoivlmzmu2X5jO6X6nQkEREJI9trPPzlP74pJW89rhdpreMcThQ8NAN6C9ExKZZbj+sNwF8nLaGyRhdCFhGR5vPM9yvJ27qdvh0SuPywbk7HCSoqUy3IZYd2pX/HRDZs284z3+U4HUdERMLE8s1lvPJTLsbAI2cOItKterA/9NNqQSLcLh4+YyDGwKs/55K9udTpSCIiEuK8/jml6ryWC0ZkMDyjtdORgo7KVAszLKM1F43sQp3Xcs8ni/B4dSFkERFpOh9mrSdr7VbatormLv+lzmT/qEy1QHeO7UOHxBjmr9/GWzPWOB1HRERCVEFZFX//0ndJs7+c0o+kWF3I+ECoTLVAiTGRPHj6QACe+Ho5eVsrHU4kIiKh6P7JSynZXsuRvdtx2pBOTscJWipTLdRx/VM5eVBHKms8/Pk/i/HNQCEiIhIYU5dsZsqiTcRFuXnkjIGaU6oRVKZasL+d1p/EmAh+WF7I5AUbnY4jIiIhorSqlr9M8s0pdecJfTSnVCOpTLVg7RNi+PPJ/QG4//OlbKmocTiRiIiEgke/zCa/tJqh6clcMqqr03GCnspUC/d/mWkc2qMNWypqeGiKLjUjIiKNMzO3mPdmrSPSbXj87MG4XTq911gqUy2cMYZHzhhEdISLT3/bwE8rCp2OJCIiQaqq1sO9ny4C4Pqje9I7NcHhRKFBZSoIdG0bzy1jfJea+eNni3SpGREROSDPfLeS1UUV9GrfiuuP6eF0nJChMhUkrjyiG/07JpK3dTtPTl3hdBwREQkySzaW8JL/kjGPnjWY6Ai305FChspUkIh0u3jsrMG4DEz4dTUL1m9zOpKIiASJOo9351U1Lh3VlYO66JIxgaQyFUQGpSVx5RHd8Vq48+MFVNd5nI4kIiJBYMKvq1m0oYTOybHceUIfp+OEHJWpIHPbcb3p3jaeFfnlPPtdjtNxRESkhVtdVMGT3/iGhzx8xkDioyMcThR6VKaCTEykm8fPHowx8MKPq1iUV+J0JBERaaE8XsudHy2gqtbLmcM6c3Sf9k5HCkkqU0Eos2sKVxzWDY/XcsdHC6ip8zodSUREWqA3pq8ha+1W2idE87dTBzgdJ2SpTAWpO47vQ9c2cSzPL+Pf3690Oo6IiLQwq4sqeOLrbAAeOWMQSXGRDicKXSpTQSo2ys3jZw/BGHj+h1Us3qDTfSIi4rPr6b0x/VOdjhTSVKaC2IhuKVw6qit1XsudHy/U6T4REQHg9V9X6/ReM1KZCnJ3je1DRkocyzaV8vwPenefiEi4yy0s54mvlwM6vddcVKaCXFxUBI+fPRiAf3+fw9KNpQ4nEhERp3j8Zyqq63R6rzmpTIWAQ7q34dJRXajzv7uv1qPTfSIi4ej1X1czV6f3mp3KVIi4a2xf0lNiWbqplOenrXI6joiINDOd3nOOylSIiI+O4PGzhgDw7PcrNZmniEgY0ek9Z6lMhZBRPdpwxWHdqPNabv1wPlW1unafiEg4eOXnXJ3ec5DKVIi5a2wferSLJ6fg/x/uFRGR0LV0Yyn/nOr79/6xswfr9J4DVKZCTEykm6fOHYrbZXjtl9VMX1XkdCQREWki1XUebvtwPrUeywUjMzhG195zhMpUCBqclsyNx/YE4M6PFlJaVetwIhERaQpPTl1B9uYyurSJ408n9XM6TthSmQpRfzimJ0PSktiwbTsPfL7U6TgiIhJgs3KLefnnXFwGnjxnKPHREU5HClsqUyEq0u3in+cMJTrCxcdz8/h6yWanI4mISICUVdVy+0cLsBauP7onB3Vp7XSksKYyFcJ6tm/FPSf2BeCPny6iqLza4UQiIhIID/53KXlbtzOgUyI3je7ldJywpzIV4i4d1ZVDe7ShuKKGez9dhLXW6UgiItIIU5ds5sOsPKIiXDx97lCiIvSn3GnaAyHO5TI88X9DSIiO4Jul+Xw0N8/pSCIicoCKyqu599NFANx1Qh96pSY4nEhAZSosdE6O5f5xvknc7p+8hDVFFQ4nEhGR/WWt5Z5PFlFcUcOo7r5JmqVlUJkKE2cM68zJgztSUePh5g/m62LIIiJB5oM56/l2WT4J0RH845whuFzG6UjipzIVJowxPHL6IDolxbBg/Tae/naF05FERKSBcgrKud8/zc0Dpw+gc3Ksw4mkPpWpMJIUF8nT5w3DZeD5H1YxM7fY6UgiIrIP1XUebnp/HttrPZwxrDNnDEtzOpLsQmUqzIzolsIfjumJtXDrB/MpqdTs6CIiLdkTXy1n6aZS0lNieWCcLmLcEqlMhaGbRvdiaHoym0qq+ONnmi5BRKSl+nFFIa/+shq3y/Cv84aREKOLGLdEDSpTxpixxpjlxpgcY8w9u1meYYyZZoyZZ4xZaIw5KfBRJVAi3S7+dd5Q4qPcTFm0SdMliIi0QEXl1dz+4QIAbh3Ti+EZmuW8pdpnmTLGuIHngBOB/sD5xpj+u6z2Z+BDa+0w4Dzg+UAHlcDq0iaeB8YNBOC+yUtYrekSRERaDGstd360gKLyakZ2S+G6o3s6HUn2oiFHpkYAOdbaXGttDTARGLfLOhZI9H+dBGwMXERpKmcO78ypQzpRWePhlonzNF2CiEgL8eb0NUxbXkhSbCRPnTsUt6ZBaNEaUqY6A+vr3c7z31fffcBFxpg84AvgxoCkkyZljOGh0wfSOTmWBXklPPmNpksQEXHask2lPPJlNgCPnjmITpoGocUL1AD084E3rLVpwEnA28aY3z22MeZqY0yWMSarsLAwQE8tjZEUG8nT5w3FZeDFH1fx0wrtFxERp1TV+qZBqKnzcv6IdE4c1NHpSNIADSlTG4D0erfT/PfVNx74EMBaOwOIAdru+kDW2pettZnW2sx27dodWGIJuIO7pnDrmN47p0vIL61yOpKISFi6//MlrCwop0e7eP5yyq7Dk6WlakiZmgP0MsZ0M8ZE4RtgPnmXddYBowGMMf3wlSkd4ggi1x/Tk8N7tqW4ooabJ87D49V0CSIizWnS/A28P3s9UREunjl/GHFREU5HkgbaZ5my1tYBNwBfA8vwvWtviTHmAWPMaf7VbgeuMsYsAN4HLrOavCiouF2Gp84dSttW0czM3cIz3610OpKISNhYVVjOHz9dBMDfTu3PgE5JDieS/WGc6jyZmZk2KyvLkeeWPfs1p4iLXpsFwLvjR3Joz9+drRURkQCqqvVw+nO/kr25jFOHdOKZ84ZijN6919IYY+ZaazN3t0wzoMv/OKxnW248thfWws0fzKewrNrpSCIiIe2B/y4le3MZ3drG88gZA1WkgpDKlPzOzaN7MbJbCoVl1dz24Xy8Gj8lItIkJs3fwHuz1hEV4eLfF+hyMcFKZUp+x+0yPHP+MFLio/h5ZREv/LjK6UgiIiEnt944qb+eonFSwUxlSnYrNTGGJ88ZAsA/py5n9uotDicSEQkdVbUe/vDePCpqPJwyuCMXjsxwOpI0gsqU7NHRfdpz3dE98Fq48f3fKCrX+CkRkUB48L9LWbaplK5t4vj7mYM0TirIqUzJXt12XG8O7tqa/NJqbnxvHnW6fp+ISKNMXrCRd3eOkxqucVIhQGVK9irS7eK5C4bTtlU0M3KL+cdUXb9PRORArcgv4+6PFwLwl1P6M7CzxkmFApUp2af2iTE8d8Ew3C7Diz+u4uslm52OJCISdEqrarn27blsr/Vw5rDOXKRxUiFDZUoaZGT3Ntwzti8Ad3y4gNVFFQ4nEhEJHtZa7vhwAblFFfTtkMDDZ2icVChRmZIGu/KIbpw4sANl1XVc985cKmvqnI4kIhIUXvwxl6lL80mIieCliw8iNsrtdCQJIJUpaTBjDI+fPZju7eLJ3lzGnz5bjC7BKCKyd7/mFPHE19kAPHXOULq0iXc4kQSaypTsl4SYSF686CBiI918Nm8D78xa53QkEZEWa+O27dz4/jzfFDPH9mRM/1SnI0kTUJmS/dY7NYFHzxoEwAOfL2Heuq0OJxIRaXmq6zxc9+5vbKmo4YhebbllTG+nI0kTUZmSAzJuaGcuO7QrtR7L9e9qQk8RkV098PlSFqzfRufkWJ45z/eOaAlNKlNywP54Uj8O6tKaTSVVXP/Ob9TUaUJPERGAj7LW75yY84WLhtM6PsrpSNKEVKbkgEVFuHjhwuF0SIxh9pot3P/5EqcjiYg4bt66rfzpP4sBeOC0AQxOS3Y4kTQ1lSlplPaJMcO51doAACAASURBVLx08UFERbh4d9Y63p211ulIIiKOyS+t4pq351JT5+WiQzI4b4Qm5gwHKlPSaEPSk/n7Gb4B6X+btITZq7c4nEhEpPlV1Xq4+u25FJRVM7JbCn87dYDTkaSZqExJQJx1UBrjD+9Gnddy3Ttz2bBtu9ORRESajbWWP362aOeA8+cvHE6kW39iw4X2tATMvSf25fCebSmuqOGat7PYXuNxOpKISLN47ZfVfPrbBmIj3bxySSZtWkU7HUmakcqUBEyE28Wz5w8jIyWOxRtKuefThZohXURC3k8rCnnki2UA/POcIfTvlOhwImluKlMSUK3jo3jlkkziotxMmr+Rl3/KdTqSiEiTWV1UwQ3v/YbXwk3H9uSkQR2djiQOUJmSgOvTIYEnzxkCwKNfZfPdsnyHE4mIBF5ZVS1XvZVFaVUdx/VP1QznYUxlSprE2IEduXVMb6yFG9+fx9KNpU5HEhEJGI/XcvPE+eQUlNM7tRVPnTsUl2Y4D1sqU9Jkbhrdk3FDO1FZ42H8m3MoKK1yOpKISEA8+N+lfJ9dQHJcJK9ckkmr6AinI4mDVKakyRhjeOyswTsvOXPlW3qHn4gEvzenr+GN6WuIcrt4+eJMurSJdzqSOExlSppUTKSbly4+iLTWsSzMK+G2D+fj9eodfiISnKZlF+y8dNZjZw9iRLcUhxNJS6AyJU2ubatoXr/sYBKiI/hy8Wb+MXW505FERPbb0o2l//+de6N7ccawNKcjSQuhMiXNoldqAs9dOBy3y/D8D6v4eG6e05FERBosv7SK8W/OoaLGw2lDOnHrmF5OR5IWRGVKms2Rvdtx32m+a1Xd++lCZuUWO5xIRGTfKmvquPLNLDaVVJHZpTWPnz0YY/TOPfn/VKakWV18SBeuOKwbtR7LNe/MJbew3OlIIiJ75PFabpk4n0UbSshIieOliw8iJtLtdCxpYVSmpNn96eR+jO7bnm2VtVz6+mwKy6qdjiQisluPfrmMqUvzSYyJYMJlB+uae7JbKlPS7Nwuw7MXDGNIWhLrt2zn8jdmU15d53QsEZH/MeGX1bzy82oiXIYXLzqInu1bOR1JWiiVKXFEXFQEr112MF3a+C6KfP27v1Hr8TodS0QEgP8u3MiDU5YC8PjZgzm0Z1uHE0lLpjIljmnbKpo3Lx9Bm/goflpRyD2fLMJazUElIs6asaqY2z5YgLVw99i+nDlcUyDI3qlMiaO6to3ntcsOJjbSzSe/5fHkNyucjiQiYSx7cylXv51FjcfLpaO6cO1R3Z2OJEFAZUocNzQ9mecuHOYbS/V9Du/MXOt0JBEJQxu3beeyCXMoq6rjxIEd+OupAzQFgjSIypS0CMf2TeWRMwYC8NdJi5m6ZLPDiUQknJRU1nLphNlsLq1iRNcUnjp3KG6XipQ0jMqUtBjnHpzBLWN64bVw4/vzyFqzxelIIhIGqmo9XPVWFisLyunVvhWvXJKpuaRkv6hMSYty8+henD8ineo6L5e/MYelG0udjiQiIazW4+X6d39j9potdEiM4c0rRpAUF+l0LAkyDSpTxpixxpjlxpgcY8w9e1jnHGPMUmPMEmPMe4GNKeHCGMOD4wYydkAHyqrquGTCbFYXVTgdS0RCkNdruf3DBXyfXUDruEjeHj+CTsmxTseSILTPMmWMcQPPAScC/YHzjTH9d1mnF3AvcJi1dgBwSxNklTAR4Xbxr/OHcnjPthSVV3PRq7PYVLLd6VgiEkKstfxt8hImL9hIfJSbNy4fQa/UBKdjSZBqyJGpEUCOtTbXWlsDTATG7bLOVcBz1tqtANbagsDGlHATHeHmpYsPYlhGMhu2befi12azpaLG6VgiEiL+OXUFb89cS1SEi1cuzWRIerLTkSSINaRMdQbW17ud57+vvt5Ab2PMr8aYmcaYsYEKKOErPjqC1y87mD6pCeQUlHPZ67Mpq6p1OpaIBLlXfsrl39NycLsMz10wnEN7aHZzaZxADUCPAHoBRwPnA68YY35X840xVxtjsowxWYWFhQF6agllyXFRvD1+BBkpcSzMK+Gqt7KoqvU4HUtEgtSHc9bz8BfLAPjH/w3muP6pDieSUNCQMrUBSK93O81/X315wGRrba21djWwAl+5+h/W2pettZnW2sx27dodaGYJM+0TY3hn/EjaJ0QzM3cL170zl+o6FSoR2T+TF2zknk8XAnDfqf05Y5guEyOB0ZAyNQfoZYzpZoyJAs4DJu+yzn/wHZXCGNMW32m/3ADmlDCX0SaOd64cSUp8FNOWF3LDe/N0YWQRabAvF23i1g/m47Vw+3G9ueywbk5HkhCyzzJlra0DbgC+BpYBH1prlxhjHjDGnOZf7Wug2BizFJgG3GmtLW6q0BKeeqcm8Pb4ESTFRvLN0nxumTifOhUqEdmHb5bmc+P78/B4LTce25MbR//uxIlIoxhrrSNPnJmZabOyshx5bgluC/O2ceErsyirrmPc0E48eY4u+yAiu/fD8gKufmsuNR4v1xzZnXtO7Kvr7ckBMcbMtdZm7m6ZZkCXoDM4LZk3x48gPsrNpPkbufuThXi9zvynQERarl9WFnH1274idflhXVWkpMmoTElQGp7RmtcvH0FspJuP5+bx50mLceooq4i0PDNzi7nyrTnU1Hm56JAM/npKfxUpaTIqUxK0RnRL4bVLM4mOcPHerHX8ddISHaESEWblFnPFG3OoqvVybmY6D5w2UEVKmpTKlAS1Q3u25ZVLMomKcPH2zLX86T+LVKhEwtj0nCIue30OlTUezhzWmUfOHIRLYyqlialMSdA7snc7Xr3Ed4Tq/dnrufuThXhUqETCzs8rC7n8jTlsr/Vw9kFpPPF/Q/TmFGkWKlMSEo7s3Y7XLzuY2Eg3H83N486PFqhQiYSRacsLGP9mFtV1Xs4fkc7jZw1WkZJmozIlIePQnm154/KDiYty8+m8DdzygeahEgkH3y7N55q35u4cbP7w6Tq1J81LZUpCysjubXjrihG0io7g8wUbuWmiZkoXCWVfLd7Mde/6pj+47NCuPDhuoIqUNDuVKQk5mV1TeGv8CBKiI/hi0Waue2euLo4sEoImzd/ADe/9Rq3HctUR3fjbqZr+QJyhMiUhaXhGa969aiTJcZF8u6yAy16fTVlVrdOxRCRA3p651ncq32u5/uge/PGkfipS4hiVKQlZg9OS+eDqUbRPiGZm7hYufHUWWypqnI4lIo1greW5aTn85T+LsRbuHtuXu8ZqZnNxlsqUhLQ+HRL4+NpDyUiJY2FeCee8NINNJdudjiUiB8Bay9+/zOaJr5djDDxyxiCuO7qH07FEVKYk9GW0iePja0fRJzWBnIJyzn5hBmuKKpyOJSL7weO13PPJIl7+KZcIl+GZ84ZxwcgMp2OJACpTEibaJ8bwwTWHMDQ9mQ3btnP2izNYurHU6Vgi0gDVdR5ufP83PshaT0yki1cuzeTUIZ2cjiWyk8qUhI3kuCjevXIkh/VsQ1F5Nee+NIPpOUVOxxKRvSjZXstlE+bwxaLNJERH8NYVIzmmT3unY4n8D5UpCSvx0RFMuOxgThrUgbLqOi59fTaT5m9wOpaI7Mamku2c8+IMZuQW0y4hmvevPoQR3VKcjiXyOypTEnaiI9z8+/zhXH5YV2o9lpsnzuelH1dhrS4/I9JSLN9cxpnPT2d5fhk92sXz6XWHMrBzktOxRHZLZUrCkstl+NupA/jzyf0A+PuX2dz/+VJdz0+kBZixqpizX5zOppIqMru05pPrDiU9Jc7pWCJ7pDIlYe3KI7rzzPnDiHK7eGP6Gv7w7m+aLV3EQZ8v2MilE2ZTVlXH2AEdeOfKkSTHRTkdS2SvVKYk7J02pBNvXjGChJgIvlqymQtfnUVRebXTsUTCirWWF35YxY3vz9t5nb3nLhxOTKTb6Wgi+6QyJQKM6tGGj689lI5JMcxdu5XTn/uV5ZvLnI4lEhaq6zzc8dFCHvsqG4B7T+zL307tj1sXLJYgoTIl4tenQwKT/nAYQ9KTydu6nbNemM607AKnY4mEtOLyai56dRaf/JZHbKSbFy8azjVH9dDlYSSoqEyJ1NM+MYYPrj6EUwZ3pLy6jvFvzuHVn3P1Tj+RJrAiv4zTn/+VOWu20iExho+uHcXYgR2djiWy31SmRHYRE+nm2fOHccuYXngtPDRlGX/8bBE1dV6no4mEjB+WF3DW89NZv2U7g9OSmHTDYZr6QIKWypTIbhhjuGVMb549fxjRES7en72eSyZoYLpIY1lrefXnXK54Yw5l1XWcPKgjH1w9itTEGKejiRwwlSmRvTh1SCc+uGYU7RKimZm7hdOe/YUF67c5HUskKFXW1HHzxPk8NGUZXgs3HtuTZ88fRmyU3rEnwU1lSmQfhqYn8/kNhzMsI5mNJVX830sz+DBrvdOxRILK2uIKznx+OpMXbCQuys0LFw7n9uP74NI79iQEqEyJNECHpBgmXn0IF47MoKbOy10fL+TP/9E4KpGG+GF5Aac++wvZm8vo1jaeSX84jBMHaaC5hA6VKZEGio5w8/AZg3j8rMFEuV28M3Md578yk/zSKqejibRIXq/l39+v5PI35lBaVceYfu2ZdMNh9EpNcDqaSECpTInsp3MOTueja0ftnODz5Gd+YXpOkdOxRFqUbZU1XP12Fv+YugKA247rzcsXZ5IYE+lwMpHAU5kSOQBD0pP5/MbDOaR7CkXl1Vz42iye/GaFLpQsAsxdu4WT/vUz3y4rIDEmgtcuzeSm0b00PkpClsqUyAFq2yqad8aP5KbRvQB45ruVXPjqTAp02k/ClNdrefHHVZzz0kw2llQxND2ZKTcdwbF9U52OJtKkVKZEGiHC7eK243rzzviRtG3lmz7hxH/9zE8rCp2OJtKstlTUcMWbc3j0y2w8XstVR3Tjw2tGkZ4S53Q0kSanMiUSAIf1bMsXNx/OYT3bUFxRw6Wvz+aJr7Op9ejdfhL6ZuUWc9K/fuaH5YUkx0Xy2qWZ/Onk/kRF6E+MhAf9posESPuEGN66YiS3HdcbAzw3bRVnvzCdVYXlTkcTaRI1dV4e/TKb816ZyebSKg7q0povbjqC0f10Wk/Ci8qUSAC5XYabRvfi/asOoXNyLAvySjj5mZ95e+ZaXSxZQsqK/DJOf+5XXvxxFQa44ZieTLz6EDolxzodTaTZGaf+gc/MzLRZWVmOPLdIcyitquW+SUv4dN4GAI7p047Hzh5M+wRdg0yCl9dreWP6Gh79KpuaOi8ZKXE8ec4QMrumOB1NpEkZY+ZaazN3u0xlSqRpTVm4iT9+toiS7bWkxEfx9zMHccKADk7HEtlvm0uquPPjBfy80jev2jmZafz11AG0io5wOJlI01OZEnHYrn+Exg3txN9OHUBKfJTDyUT2zVrLh1nreWjKMsqq6vSfAglLKlMiLYDXa3lzxhoe+yqbqlovbeKjuO+0AZwyuCPGaDJDaZnWb6nk3k8X8Yt/lv/Rfdvz97MG6XS1hB2VKZEWZG1xBfd8sogZucUAHNc/lYdOH0hqov44Scuxo/w/8fVyKms8tI6L5L7TBnDakE4q/xKW9lamGvRuPmPMWGPMcmNMjjHmnr2sd5YxxhpjdvtkIgJd2sTz3lUjeeSMQbSKjuCbpfmMefJHPpyzXu/4kxYhp6Cc/3tpBvd/vpTKGg+nDO7IN7cdxbihnVWkRHZjn0emjDFuYAVwHJAHzAHOt9Yu3WW9BGAKEAXcYK3d62EnHZkSgU0l2/nTZ4v5PrsAgBFdU3jw9IH06ZDgcDIJR9trPPx72kpe/imXWo+lfUI0D50+kOM1Nkqk0UemRgA51tpca20NMBEYt5v1HgQeA3RhMpEG6pgUy2uXZvL0uUNpEx/F7DVbOOmZn3nov0spr65zOp6ECWstU5dsZsyTP/LctFXUeiznHZzON7cdpSIl0gANKVOdgfX1buf579vJGDMcSLfWTtnbAxljrjbGZBljsgoLde0yEQBjDKcP68z3tx/NxYd0wWstr/6ymtH//IHPF2zUqT9pUuuKKxn/ZhZXvz2XDdu2079jIp9cdyiPnjWYpNhIp+OJBIVGz4BujHEBTwK372tda+3L1tpMa21mu3btGvvUIiElKS6SB08fyOQ/HM6Q9GTyS6u58f15XPTaLFbklzkdT0JMZU0dT3+7guOe+pHvswtIiI7gvlP7M/mGwzioS2un44kElYbMtLYBSK93O81/3w4JwEDgB//AxA7AZGPMafsaNyUivzcoLYnPrjuUiXPW8/jX2fyaU8zYp3/i3IMzuO243rRLiHY6ogQxr9fyyW95/GPqcvJLqwE4Y1hn7j2pr6Y7EDlADRmAHoFvAPpofCVqDnCBtXbJHtb/AbhDA9BFGm9LRQ1PfbOC92avw+O1xEe5ue7oHow/vDuxUW6n40mQmZ5TxENTlrF0UykAAzsn8peT+zOyexuHk4m0fI2eZ8oYcxLwNOAGJlhrHzbGPABkWWsn77LuD6hMiQRUTkEZf/8im+/87/rrmBTDnSf04fShnXG59FZ12bvd/f7cNbYP44bo90ekoTRpp0iI2PXIQt8OCdwypjcnDEjV/D/yO+uKK/nXdyv5bF4eXgvxUW6uP6Yn4w/vRkykjmyK7A+VKZEQ4vVaPp23gX9OXc6mEt9MJAM7J3Lbcb05pk97lSphw7bt/Pv7lXyUlUed1+J2Gc49OJ1bx2jMnciBUpkSCUHVdR4mzl7Pc9NyKCjzDSQelpHMbcf15vCebVWqwlB+aRXPTcth4uz11Hi8uAycMSyNm0b3pEubeKfjiQQ1lSmREFZV6+GdmWt54YdVFFfUAL5Sde1RPTiuX6rGxISBtcUVvPxTLh/NzaOmzosxcOrgTtw8phc92rVyOp5ISFCZEgkDFdV1vDljDS//lMu2yloAerSL55ojezBuWCeiIzRGJtQsyivhxR9X8eXiTXj9/5SfOLADt4zprUsSiQSYypRIGKmsqWPi7PW8+nMuG/1jqlIToxl/eDfOH5FBQoxmtQ5m1lp+ySnixR9X8WtOMQCRbsPpQztzzVHd6dleJUqkKahMiYShWo+Xzxds5KUfc1nun0E9PsrNGcM7c8morvRO1R/dYFJWVcunv23g7ZlrySkoB3z784KRGVxxeDc6JsU6nFAktKlMiYQxay3Tlhfw0o+5zFq9Zef9I7ulcMmorhw/IJVId6OvLCVNZEV+GW/NWMNnv22gosYDQPuEaC49tCsXjexCUpyONIo0B5UpEQFg+eYy3p75v3+YUxOjOTcznTOHp9G1rd7x1RJsr/Hw1ZJNTJy9XgVYpIVQmRKR/7G7U0YAB3VpzVnD0zh5cEeSYnXEozlZa8lau5WPs/KYsmgT5dV1AMRFuTlzeGcuPqSrBpWLOEhlSkR2y1rLzNwtfDw3jy8Xb6LSf7QqKsLF8f1TOWNYZw7v1VbvBGwi1lqW55fxxcJNTF6wkTXFlTuXDctI5qzhaYwb2klvGhBpAVSmRGSfKqrr+HLxZj6Zm8eM3OKd9ydERzCmfyonDuzAkb3b6TIkAbAiv4wpCzcxZdGm/zkymJoYzZnD0zhreBo922t+KJGWRGVKRPZL3tZK/jNvA/9duInszWU774+PcnN0n/Yc07c9R/dpR9tWujRJQ9TUeclas4Xvsgv4PruA1UUVO5e1jotk7MAOnDSoI4f2aItbk6yKtEgqUyJywHILy/ly8Wa+XLyJxRtK/2fZkLQkju7TnqP6tGNQ5yQNiq5n/ZZKpq8q4qcVRfy0opAy/xgogOS4SE7o34GTB3dkVI82+rmJBAGVKREJiHXFlXyfnc+05YXMyC2mps67c1l8lJuDuqZwSPcUDuneJqzKlbWWjSVVzFm9hemripi+qpi8rdv/Z53eqa04tm8qo/u1Z1h6MhFh8rMRCRUqUyIScNtrPMzILWJadiG/rioit7Dif5bHRLoY1DmJIWnJDElPZmh6MmmtY0PiAszl1XUs3lDCvHXbmL9+K/PWbdt5sekdkmIjGdW9DYf2bMMxfdqTnhLnUFoRCQSVKRFpcgWlVcxcvYWZucXMXFVMblHF79ZpHRdJnw4J9ElNoLf/c6/UhBY7DUNVrYe8rZWszC9n2eYysjeVkr25jHVbKn+3blJsJMMzkjm0R1tG9WhD/46Jusi0SAhRmRKRZrelooYFedtYsN7/kVfCloqa3a6bHBdJRkoc6Slxvs+t40hNjKZdgu+jTXw0URGBPS1mraW0qo6C0iryS6vJL60iv6yK9Vu2s7a4grXFlWws2c7u/omMdBt6pyYwLCOZYemtGZaRTLe28SFx1E1Edk9lSkQct2Nc0Yr8MlZsLmN5fhkr8svIKSinqta7z+9PjoskMSaSVtERtIqJICE6gvjoCKIiXES4DG6XIcJlcLkMXq+lxmOp9Xip9XipqfNSXl1HaVUdZdtrKa2qpXR7HTWevT+v22XonBxL93bx9O2QSL+OCfTtkEj3dvFhMx5MRHz2VqYimjuMiIQnY3zFpHNyLMf0ab/zfmstheXVrN9Sybotlawr3k7e1koKy6spLPN9FJVXs62ylm2VtQHN1Co6gvaJ0aQmxJCaGE37xBg6JcXQtW08XdvE07l1rEqTiOyTypSIOMoYQ/uEGNonxHBQl5TdruPxWrZW1lBeVUd5dR1l/s/l1bXUeiwer6XOa/F4vNR5LREuQ2SEi0i3iyi373N8tJvEWN/RrcTYCBJjIjUBqYgEhMqUiLR4bpehbatoTRIqIi2Sjl+LiIiINILKlIiIiEgjqEyJiIiINILKlIiIiEgjqEyJiIiINILKlIiIiEgjqEyJiIiINILKlIiIiEgjqEyJiIiINILKlIiIiEgjGGutM09sTCGwthmeqi1Q1AzP0xJp28NXOG9/OG87hPf2a9vDV3NsfxdrbbvdLXCsTDUXY0yWtTbT6RxO0LaH57ZDeG9/OG87hPf2a9vDc9vB+e3XaT4RERGRRlCZEhEREWmEcChTLzsdwEHa9vAVztsfztsO4b392vbw5ej2h/yYKREREZGmFA5HpkRERESaTNCXKWPM/xljlhhjvMaYPY7kN8aMNcYsN8bkGGPuqXd/N2PMLP/9HxhjoponeWAYY1KMMd8YY1b6P7fezTrHGGPm1/uoMsac7l/2hjFmdb1lQ5t/Kw5MQ7bdv56n3vZNrnd/OOz7ocaYGf7XyEJjzLn1lgXdvt/T67je8mj/vszx79uu9Zbd679/uTHmhObMHQgN2PbbjDFL/fv5O2NMl3rLdvsaCCYN2P7LjDGF9bbzynrLLvW/TlYaYy5t3uSN14Btf6redq8wxmyrtyyo970xZoIxpsAYs3gPy40x5hn/z2ahMWZ4vWXNt9+ttUH9AfQD+gA/AJl7WMcNrAK6A1HAAqC/f9mHwHn+r18ErnN6m/Zz+x8H7vF/fQ/w2D7WTwG2AHH+228AZzu9HU257UD5Hu4P+X0P9AZ6+b/uBGwCkoNx3+/tdVxvneuBF/1fnwd84P+6v3/9aKCb/3HcTm9TgLf9mHqv6+t2bLv/9m5fA8Hy0cDtvwz4926+NwXI9X9u7f+6tdPbFMht32X9G4EJIbTvjwSGA4v3sPwk4EvAAIcAs5zY70F/ZMpau8xau3wfq40Acqy1udbaGmAiMM4YY4BjgY/9670JnN50aZvEOHy5oWH5zwa+tNZWNmmq5rG/275TuOx7a+0Ka+1K/9cbgQJgt5POBYHdvo53Waf+z+RjYLR/X48DJlprq621q4Ec/+MFi31uu7V2Wr3X9UwgrZkzNqWG7Ps9OQH4xlq7xVq7FfgGGNtEOZvC/m77+cD7zZKsGVhrf8J3AGBPxgFvWZ+ZQLIxpiPNvN+Dvkw1UGdgfb3bef772gDbrLV1u9wfTFKttZv8X28GUvex/nn8/oX2sP/w6FPGmOiAJ2w6Dd32GGNMljFm5o7Tm4ThvjfGjMD3P9tV9e4Opn2/p9fxbtfx79sSfPu6Id/bku1v/vH4/re+w+5eA8Gkodt/lv/3+WNjTPp+fm9L1eD8/lO73YDv690d7Pt+X/b082nW/R7RVA8cSMaYb4EOu1n0J2vtpObO09z2tv31b1hrrTFmj2/P9Lf1QcDX9e6+F98f4ih8by29G3igsZkDJUDb3sVau8EY0x343hizCN8f2RYvwPv+beBSa63Xf3eL3vdyYIwxFwGZwFH17v7da8Bau2r3jxC0Pgfet9ZWG2OuwXeE8liHMzW384CPrbWeeveFw753XFCUKWvtmEY+xAYgvd7tNP99xfgOCUb4/xe74/4WZW/bb4zJN8Z0tNZu8v/BLNjLQ50DfGatra332DuObFQbY14H7ghI6AAJxLZbazf4P+caY34AhgGfECb73hiTCEzB95+PmfUeu0Xv+93Y0+t4d+vkGWMigCR8r/OGfG9L1qD8xpgx+Ir2Udba6h337+E1EEx/UPe5/dba4no3X8U3pnDH9x69y/f+EPCETWd/fnfPA/5Q/44Q2Pf7sqefT7Pu93A5zTcH6GV8796KwvcLN9n6RqlNwzeOCOBSINiOdE3Glxv2nf9359L9f4R3jCE6HdjtOyZaqH1uuzGm9Y7TV8aYtsBhwNJw2ff+3/fP8I0p+HiXZcG273f7Ot5lnfo/k//Xzt2zxBFFARh+J01SpElSxSoKAcsIKUJ+gGCKhYCFhQjRxv+wjZ35AykCKdOl20IQzFoKqQwSRN1SkKSzFIu1uCdwlQ27MrMzbnwfGHbm7lyYM2c+zs7HLgLdyHUHWCrS237TwEvgR03LXYWhsRdFMQd8Blr9fv9P1j5wH6htyasxSvzPs8kWcBjj28B8rIcnwDzXr87fdaNs9xRFMUt60Hova/sfcj9MB1iJt/reAOfxQ7HevI/ryfa6BuA96V7oBfAb2I72KWArm+8dcEyqyNtZ+wzpoNoDvgEPm47plvE/A74DJ8AO8DTakN/N6wAAANNJREFUXwNfsvlekCr1Bzf6d4ED0on0K/C46ZiqjB14G/H9jM+1+5R7YBm4BPaz4dWk5n7Qfky6NdmK8UeRy17kdibr245+R8BC07GMIfadOAb+zXMn2v+5D0zSMEL8m8CviHMXmM36rsY20QM+NB1L1bHH9Abw8Ua/ic896QLAWRzHTknPA64D6/F9AXyKdXNA9lZ/nXn3H9AlSZJKuC+3+SRJksbCYkqSJKkEiylJkqQSLKYkSZJKsJiSJEkqwWJKkiSpBIspSZKkEiymJEmSSrgCFS6+0VPEMfMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = 1.\n",
    "plot_1d(lambda _x: solution(t, _x), titles=['exact solution']) #only works for d=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stochastic Interpretation\n",
    "\n",
    "Our main theoretical tool will be the **Feynman-Kac** formula, stating that (under suitable regularity conditions)\n",
    "\n",
    "$$\\color{blue}{f(t,x)=\\mathbb{E}\\big[\\varphi(S_t^x)\\big]},$$\n",
    "\n",
    "where $S_t^x$ is a solution process to the corresponding stochastic differential equations (**SDE**)\n",
    "\n",
    "$$S^x_t = x + \\int_0^t \\mu(S^x_s)\\,ds + \\int_0^t \\sigma(S^x_s)\\,dB_s. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ... back to our example\n",
    "\n",
    "In our example it holds that \n",
    "\n",
    "$$S_t^x = x+ \\sqrt{t}*  \\sigma N,$$\n",
    "\n",
    "where $N\\sim \\mathcal{N}(0,I_d)$ is a $d$-dimensional, uncorrelated standard normal random variable.\n",
    "Accordingly,\n",
    "\n",
    "$$F(t,x) = \\mathbb{E}[\\varphi(x+ \\sqrt{t}*  \\sigma N)].$$\n",
    "\n",
    "This can, of course, also be seen as convolution between $\\varphi$ and the corresponding heat kernel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monte-Carlo Sampling\n",
    "\n",
    "For specific $x\\in\\mathbb{R}^d$, $t\\in \\mathbb{R}$, we can now compute $m$ i.i.d. samples $\\mathcal{S}^{x,i}_{t} \\sim \\mathcal{S}^{x}_{t}$  and approximate the solution to the PDE via **Monte-Carlo sampling**, i.e.\n",
    "\n",
    "$$ \\color{blue}{ f(t,x) = \\mathbb{E}\\big[\\varphi(S_t^x)\\big] \\approx \\frac{1}{m}\\sum_{i=1}^m\\varphi(\\mathcal{S}^{x,i}_{t}) }. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = 0.5*tf.ones((1,d)) #0.5*tf.ones([1, d]) \n",
    "t = 1. #tf.constant([[1.]])\n",
    "m = 50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we compute the realizations efficiently in matrix form\n",
    "def batched_MC(sigma, phi, _t, _x):   \n",
    "    \n",
    "    N = tf.random.normal(shape=tf.shape(_x))\n",
    "    sigma_N = tf.einsum('ij, kj -> ik', sigma, N)\n",
    "    #shape = (d, batchsize)\n",
    "    t_sigma_N = tf.einsum('ij, ki -> ik', tf.sqrt(_t), sigma_N)\n",
    "    #shape = (batchsize, d)\n",
    "    return phi(_x + t_sigma_N)\n",
    "    \n",
    "\n",
    "def repeat(tensor, num):\n",
    "    tensor = tf.reshape(tensor, [1,-1])\n",
    "    return tf.tile(tensor, [num,1]) #row vectors due to tensorflow convention (we want to reuse this later)\n",
    "\n",
    "x_rep = repeat(x,m)\n",
    "t_rep = repeat(t,m)\n",
    "\n",
    "phi_S = batched_MC(sigma, phi, t_rep, x_rep)\n",
    "f_approx = tf.reduce_mean(phi_S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact solution: 0.552500\n",
      "Monte-Carlo Approximation: 0.500079\n"
     ]
    }
   ],
   "source": [
    "print('exact solution: {:.6f}\\nMonte-Carlo Approximation: {:.6f}'.format(solution(t,x)[0][0], f_approx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAE/CAYAAABin0ZUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3gU173/8fdXK60aIBASvQnRDS4g3MDYuBtwwwaXuJDmFNuJnTg3zclNvfHNzzchjh2IW2KbuIArxjG4YBtjUwyY3nsHFZAQ6qvz+2MXsoAAwUoarfbzep59tDszmv0eDYs+OnPmjDnnEBEREZHTE+d1ASIiIiLRTGFKREREJAIKUyIiIiIRUJgSERERiYDClIiIiEgEFKZEREREIqAwJSJRz8zGmdlsr+uIhJn9yswmhZ53MbNiM/N5XZeInJzClEgTYmabzazCzDKOWv6lmTkz61YH7/GxmX0jwn20MLPxZrY1FBo2hF5nnPy765aZDTWzz82s0MwKzOwzMxvc0HWEc85tdc41c84FvKxDRGpHYUqk6dkE3HbohZkNAFK8K+dIZuYHPgTOAK4GWgAXAPnAuaexv/gIamkBTAP+CqQDHYFfA+Wnu08RiT0KUyJNzwvAXWGv7waeD9/AzNLM7HkzyzWzLWb2sJnFhdaNM7PZZvaome0zs01mdk1o3e+Bi4DHQz1Kj4eW9zGz90M9O2vMbOwJ6rsL6ALc6Jxb6Zyrds7tdc791jn379D+fhLqrTpgZivN7Maw2seFeo/+bGb5wK+OfgMzu9DMvgj1Nn1hZhcep5ZeAM65l5xzAedcqXPuPefc0tB+ss1sppnlm1memf3LzFqGvc9mM/uRmS01s4Nm9oyZtTWzd0O1f2BmrULbdgv1Dt5jZjvNbJeZPVRTUWHbxodef2xmvw21+4CZvRfei2dmd4WOY76Z/SJU1+UnOAYiUocUpkSanrlACzPrGxpzcysw6aht/gqkAd2BiwkGnK+GrT8PWANkAH8EnjEzc879HPgUuC90Guo+M0sF3gdeBNqE3u9vZtbvOPVdDkx3zhWfoA0bCIa2NII9RZPMrP1R9W0E2gK/D/9GM0sH3gEeA1oDfwLeMbPWNbzPWiBgZs+Z2TWHgk/47oA/AB2AvkBnjg1vNwFXEAxm1wLvAj8DMgn+H/u9o7YfDvQErgR+fAqh53aCx6gN4AceCrW3H/A34CtAe4I/s4613KeI1AGFKZGm6VDv1BXAKmDHoRVhAeunzrkDzrnNwP8Bd4Z9/xbn3FOhMTvPEfwl3fY47zUK2Oyc+4dzrso59yXwGjDmONu3BnadqHjn3BTn3M5Qr9UrwDqOPAW40zn319D7lR717SOBdc65F0LrXwJWEww6R79PETAUcMBTQK6ZTTWztqH1651z7zvnyp1zuQSD2cVH7eavzrk9zrkdBIPmPOfcl865MuAN4Jyjtv+1c+6gc24Z8A/CTsmexD+cc2tD7Z0MnB1afjPwtnNutnOuAvhlqD0i0kBOe6yBiDRqLwCzgCyOOsVHsLcpAdgStmwLR/Zm7D70xDlXYmYAzY7zXl2B88xsf9iyeOAFM+sCrAzbVzOCY6PacwJmdhfwA6BbaFGzUN2HbDvBt3fgyLbBse07zDm3ChgXet8+BHvxxgO3hULVXwj2kjUn+AfovqN2sSfseWkNr4/+uYXXvgUYcIK2hNsd9rwkbL8dwvcZOl75tdyniNQB9UyJNEHOuS0EB6KPAF4/anUeUEkwBB3ShbDeq5Pt/qjX24BPnHMtwx7NnHPfCbsqrVkoSAF8AFwVOj14DDPrSrCX6D6gtXOuJbCc4Cm349UQbudRbYNats85txr4J9A/tOh/Qu81wDnXArjjqDpOR+ej6toZ4f52AZ0OvTCzZIK9fyLSQBSmRJqurwOXOucOhi8MnbqbDPzezJqHwssPOHZc1fHsITjW6pBpQC8zu9PMEkKPwWbW9zjf/wLBAPZaaOB6nJm1NrOfmdkIIJVggMkFMLOv8p9wUxv/DtVzu5nFm9ktQL9QnUcIvf8PzaxT6HVngqfd5oY2aQ4UA4Vm1hH40SnUcTy/MLMUMzuD4BioVyLc36vAtaFB936CY7oiDXwicgoUpkSaKOfcBufcguOsvh84SHAQ92yCg8efreWu/wLcHLrS7zHn3AGCg6lvJdjLshv4XyDxOHWVExyEvprgwPUiYD7B03jznHMrCY7hmkMwuA0APqtlbTjn8gmO4/ohwVOK/wWMcs7l1bD5AYKD2eeZ2UGCIWp56HshOPh9IFBIcFD70b18p+MTYD3B6SEedc69F8nOnHMrCB7Plwn2UhUDe9H0DiINxpzTOEURkfpmwQlTNwEJzrmqenyfZsB+oKdzblN9vY+I/Id6pkREopyZXRs6dZgKPAosAzZ7W5VI7FCYEhGJftcTPMW6k+AcVrc6nXYQaTA6zSciIiISAfVMiYiIiERAYUpEREQkAp7NgJ6RkeG6devm1duLiIiI1NrChQvznHOZNa3zLEx169aNBQuONwWOiIiISONhZkffpuowneYTERERiYDClIiIiEgEFKZEREREIqAwJSIiIhIBhSkRERGRCChMiYiIiERAYUpEREQkAgpTIiIS9XJLchk3fRx5pXlelyIxSGFKRESi3sSlE1m0ZxETlkzwuhSJQZ7NgC4iIhKpQZMGURGoOPx68prJTF4zGb/Pz8I7FnpYmcQS9UyJiEjUmj56OiOyRpDkSwIgyZfEyKyRzLhphseVSSxRmBIRkaiVmZJJakIq5YFy/D4/5YFyUv2pZCRneF2axBCd5hMRkahWUFbA2N5jGdNrDFPWTtEgdGlw5pzz5I1zcnLcggULPHlvERERkVNhZgudczk1rdNpPhEREZEIKEyJiIiIREBhSkRERCQCClMiIiIiEVCYEhEREYmAwpSIiIhIBBSmRERERCKgMCUiIiISAYUpERERkQgoTImIiIhEQGFKREREJAIKUyIiIiIRUJgSERERiYDClIiIiEgEFKZEREREIqAwJSIiIhIBhSkRERGRCChMiYiIiERAYUpEREQkAgpTIiIiIhFQmBIRERGJgMKUiIiISAQUpkREREQioDAlIiIiEgGFKREREZEIKEyJiIiIREBhSkRERCQCClMiIiIiEVCYEhEREYmAwpSIiIhIBBSmRERERCKgMCUiIiISgZOGKTPrbGYfmdlKM1thZt+vYRszs8fMbL2ZLTWzgfVTroiIiEjjEl+LbaqAHzrnFplZc2Chmb3vnFsZts01QM/Q4zxgQuiriIiISJN20p4p59wu59yi0PMDwCqg41GbXQ8874LmAi3NrH2dVysiIiLSyJzSmCkz6wacA8w7alVHYFvY6+0cG7hEREREmpxahykzawa8BjzgnCs6nTczs3vMbIGZLcjNzT2dXYiIiIg0KrUKU2aWQDBI/cs593oNm+wAOoe97hRadgTn3JPOuRznXE5mZubp1CsiIiLSqNTmaj4DngFWOef+dJzNpgJ3ha7qOx8odM7tqsM6RURERBql2lzNNwS4E1hmZotDy34GdAFwzk0E/g2MANYDJcBX675UERERkcbnpGHKOTcbsJNs44B766ooERERkWihGdBFREREIqAwJSIiIhIBhSkRERGRCChMiYiIiERAYUpEREQkAgpTIiIiIhFQmBIRkaiXW5LLuOnjyCvN87oUiUEKUyIiEvUmLp3Ioj2LmLBkgtelSAyqzQzoIiIijdKgSYOoCFQcfj15zWQmr5mM3+dn4R0LPaxMYol6pkREJGpNHz2dEVkjSPIlAZDkS2Jk1khm3DTD48oklihMiYhI1MpMySQ1IZXyQDl+n5/yQDmp/lQykjO8Lk1iiE7ziYhIVCsoK2Bs77GM6TWGKWunaBC6NDgL3qO44eXk5LgFCxZ48t4iIiIip8LMFjrncmpap9N8IiIiIhFQmBIRERGJgMKUiIiISAQUpkREREQioDAlIiIiEgGFKREREZEIKEyJiIiIREBhSkRERCQCClMiIiIiEVCYEhEREYmAwpSIiIhIBBSmRERERCKgMCUiIiISAYUpERERkQgoTImIiIhEQGFKREREJAIKUyIiEvVyS3IZN30ceaV5XpciMUhhSkREot7EpRNZtGcRE5ZM8LoUiUHxXhcgIiJyugZNGkRFoOLw68lrJjN5zWT8Pj8L71joYWUSS9QzJSIiUWv66OmMyBpBki8JgCRfEiOzRjLjphkeVyaxRGFKRESiVmZKJqkJqZQHyvH7/JQHykn1p5KRnOF1aRJDdJpPRESiWkFZAWN7j2VMrzFMWTtFg9ClwZlzzpM3zsnJcQsWLPDkvUVEREROhZktdM7l1LROp/lEREREIqAwJSIiIhIBhSkRERGRCChMiYiIiERAYUpEREQkAgpTIiIiIhFQmBIRERGJwEnDlJk9a2Z7zWz5cdZfYmaFZrY49Phl3ZcpIiIi0jjVZgb0fwKPA8+fYJtPnXOj6qQiERERkShy0p4p59wsoKABahERERGJOnU1ZuoCM1tiZu+a2Rl1tE8RERGRRq8ubnS8COjqnCs2sxHAm0DPmjY0s3uAewC6dOlSB28tIiIi4q2Ie6acc0XOueLQ838DCWaWcZxtn3TO5TjncjIzMyN9axERERHPRRymzKydmVno+bmhfeZHul8REZHaWp2/mgtevIC1BWu9LkViUG2mRngJmAP0NrPtZvZ1M/u2mX07tMnNwHIzWwI8BtzqnHP1V7KIiMiRfvLpTyiuLOa/Zv2X16VIDDKvck9OTo5bsGCBJ+8tIiJNw4DnBhx33bK7lzVgJdLUmdlC51xOTes0A7qIiEStKaOm0CG1wxHLOqR24LVrX/OoIolFClMiIhK1+rTuQ1J80hHLkuOT6ZXey6OKJBYpTImISFQ7UHGA7LRsHh32KNlp2RRVFHldksSYuphnSkRExDMzx848/PyqrKs8rERilXqmRERERCKgMCUiIiISAYUpERERkQgoTImIiIhEQGFKREREJAIKUyIiEvVyS3IZN30ceaV5XpciMUhhSkREot7EpRNZtGcRE5ZM8LoUiUGaZ0pERKLWoEmDqAhUHH49ec1kJq+ZjN/nZ+EdCz2sTGKJeqZERCRqTR89nRFZI0jyBW8pk+RLYmTWSGbcNMPjyiSWNOkwVVRW6XUJIiJSjzJTMklNSKU8UI7f56c8UE6qP5WM5AyvS5MGUloRoKKq2tMammyY+nLrPoY8MpM3vtzudSkiIlKPCsoKGNt7LC+OeJGxvceSX5rvdUnSQLbvK+GmCZ/zu3dWelpHkx0z9cGqPRwoq+LBV5awfEcRP72mD/G+JpsdRURi1vjh4w8/f/j8hz2sRBrS5xvyuO/FLyk4WEFJRRWFpZWkJSd4UkuTTRc/uqoPv7+xP/FxxjOzN3HXs/MpOFhx8m8UERGRRss5x7OzN3HnM8Hf6xf3yuSte4d6FqSgCYcpgK+c15WX7jmfjGaJfL4hn+sen83KnUVelyUiIiKnoawywA+nLOE301YSqHZ855Jsnh03mLQU74IUNPEwBTC4Wzpv3z+EszqlsX1fKaMnfMbbS3Z6XZaIiIicgp37Sxn79zm8vmgHyQk+Hr/9HH58dR98ceZ1aU0/TAG0T0vmlW9dwE0DO1FWWc39L33JI++uJlDtvC5NRERETmL+pgKue3w2S7cX0qlVMq9950JGndnB67IOi4kwBZCU4OPRMWfyq2v74YszJn6yga/+8wsKSzR9goiISGPknOOFuVu4/am55BVXMKRHa96+byj9OrTwurQjxEyYAjAzxg3JYtLXzyM91c+stblc98Rs1uw+4HVpIiIiEqa8KsBPX1/GL95cTlW14xtDs3juq+fSKtXvdWnHiKkwdcgF2a2Zet8QzujQgi35JdzwxGe8tXiH12WJiIgIwfmjxkycw8tfbCMxPo4/33IWD4/q12inOGqcVTWATq1SePXbFzL6nI6UVgb4/suL+e+3lns+i6qIiJy63JJcxk0fR15pntelSIQ+WZvLqL8Gx0d1bJnMq9++kBvP6eR1WScUs2EKINnv4//GnsXvbuiP3xfHc3O2MPbvc9i5v9Tr0kRE5BRMXDqRRXsWMWHJBK9LkdNUXe0Y/8Faxv1jPvtLKrmkdybvfG8oAzqleV3aSZlz3lzRlpOT4xYsWODJe9dk8bb93PuvRezYX0p6qp/Hbj2HoT11bycRkcZs0KRBVASOnZDZ7/Oz8I6FHlQkp2PfwQoeeGUxn6zNxQwevLwX9w3vQVwjmPbgEDNb6JzLqWldTPdMhTu7c0vevn8oF/XMoOBgBXc+O4/HZ66jWtMniIg0WtNHT2dE1giSfEkAJPmSGJk1khk3zfC4Mqmtpdv3M+qvs/lkbS6tUhJ47qvn8r3LejaqIHUyClNh0lP9/DN0EJ2DR99byzefX6DpE0REGqnMlExSE1IpD5Tj9/kpD5ST6k8lI1lnFho75xwvztvKzRPmsGN/KWd1SuPt+4cyrFem16WdMoWpo/jijB9c0Yt/jBtMWnICH67ey6jHP2X5jkKvSxMRkRoUlBVwXfZ1ZKdlc132deSX5ntdkpxEaUWAH726lJ+9sYyKQDV3nN+Fyd++gE6tUrwu7bQoTB3H8D5tmHb/UAZ0TGNbQSmjJ3zOC3M249UYMxERqdn44eNJjE9kdcFqEuMTGT98vNclyQms33uAG574jFcXbicpIY4/jT2L390wgMR4n9elnTYNQD+JssoAv5m2khfnbQVgxIB2PHLTmbRI8vamiiIiogHo0ebVhdv5xZvLKa0M0D0jlSe+MpC+7RvXbObHowHoEUhK8PE/Nw7gsdvOoVliPP9etptRj81m6fb9XpcmIhLzDg1A9xHs1fDh0wD0RqikooofTl7CQ1OWUFoZ4MZzOvL2/UOjJkidjHqmTsGmvIPc9+IiVuwsIsFn/GxEX8Zd2A2z6LniQESkKVHPVOO3ZvcB7n1xEev3FpOUEMdvruvPmJxOUfe7Uz1TdSQrI5XXvnMhd13QlcqA49dvr+RbLyzU1X4iIh6ZPno6bVPa4rNQz5T5aJvSVj1TjYBzjle+2Mp1j89m/d5ierRpxlv3DmXs4M5RF6RORmHqFCUl+PjN9f3521cG0jwxnvdW7mHEY5/y5dZ9XpcmIhJzMlMyGdZpGNWuGr/PT7Wr5uLOF2tqBI8Vl1fx4CuL+fFryyivqubmQZ2Yet8Qerdr7nVp9UJh6jSNGNCed753EWd2SmPH/lLGTJzDk7M2aJJPEZEGVlBWwNjeY3lxxIuM7T1WUyN4bMXOQq7762zeXLyT5AQf/zfmLB4dcxYp/nivS6s3GjMVoYqqah55dzXPfrYJgIt6ZvB/Y86iTYskjysTERFpONXVjmc/28Qfp6+hIlBN77bNeeIrA+nRppnXpdUJjZmqR/74OH55bT+euiuHVikJfLouj6v/8ikfrNzjdWkiIiINYm9RGXf/Yz6/e2cVFYFqvnJeF968d0iTCVInozBVR67o15bpDwxjaI/gvf2+8fwCHn5zGaUVAa9LExERqTcfrNzD1X/5lE/X5dEqJYEn7xzE728cQLI/eifhPFUKU3WobYsknv/aufx8RF8SfMakucGrGFbtKvK6NBERkTpVWhHg4TeX8Y3nF1BwsIKhPTKY/sAwrjyjndelNTiFqToWF2d8c1h33vjuELpnprJubzHXP/4Zz8zepMHpIiLSJKzcWcS1j89m0tytJPiMn4/oy/NfO5e2MTpeWGGqnvTvmMa0+4dy27ldqAhU89tpKxn3zy/Ye6DM69JEREROS3W14+lPN3LDE5+xfm8x2ZmpvPHdIXxzWHfi4prW3FGnQmGqHqX44/nD6AFMvGMQLVMSmLU2l2vGf8qMFbu9Lk1EROSU7CosPWKQ+e3ndWHa/RfRv2Oa16V57qRhysyeNbO9Zrb8OOvNzB4zs/VmttTMBtZ9mdHt6v7tmP79YVyY3Zr8gxV864WF/HDyEorKNHO6iIg0bs453vxyB1f9eRafrsujZUoCf79zEP8TY4PMT6Q2PVP/BK4+wfprgJ6hxz3AhMjLanrapSUx6evn8ctR/UiMj+O1Rdu5+s+z+Hx9nteliYhEvdySXMZNH0deqf5PrUsFByu498VFPPDKYorKqri0Txvee2AYV8XgIPMTOWmYcs7NAgpOsMn1wPMuaC7Q0sza11WBTUlcnPG1oVm8872LOKtTGjsLy7j96Xn8auoKTaEgIhKBiUsnsmjPIiYs0d/zdeXDVXu48s+z+Pey3aT6fTwyegDP3J2jSalrUKsZ0M2sGzDNOde/hnXTgEecc7NDrz8EfuycO+H05k1lBvTTVRWo5m8fb+CxD9dRVe3onpnKn8aezdmdW3pdmohI1Bg0aRAVgYpjlvt9fhbesdCDiqLfgbJKfjdtFa8s2AbAud3S+b+xZ9E5PcXjyrzVaGZAN7N7zGyBmS3Izc1tyLdudOJ9cXzvsp688d0h9GzTjI25B7lpwuf86b01VAaqvS5PRCQqTB89nRFZI0jyBXtLknxJjMwayYybZnhcWXSauzGfa/7yKa8s2IbfF8fPR/TlpXvOj/kgdTJ1EaZ2AJ3DXncKLTuGc+5J51yOcy4nMzOzDt46+g3olMbb9w/lmxdlUe0cj81cz41/+4zVuzXRp4jIyWSmZJKakEp5oBy/z095oJxUfyoZyRlelxZVSisC/G7aSm57ai7b95VyRocWTPveUL45rDu+GJ7yoLbqIkxNBe4KXdV3PlDonNtVB/uNGUkJPn4+sh8vffN8OrVKZvmOIq7962zGf7CWiir1UomInEhBWQHXZl9Lj7QeXJt9Lfml+V6XFFXmbcznmr/M4unZm4gz43uX9uCN7w6hV9vmXpcWNWozNcJLwBygt5ltN7Ovm9m3zezboU3+DWwE1gNPAd+tt2qbuPO7t2b6A8O48/yuVAYc4z9Yx3WPz2bZ9kKvSxMRabTGDx9PUnwSqwpWkRSfxPjh470uKSoUl1fxy7eWc8uTc9mcX0Lvts15/TsX8oMre+OP1zSUp6JWA9DrQ6wPQD+ZuRvz+fFrS9mSX4IvzvjWsO5877KeJCVoTg8RkUM0AP30fLoul5+8towd+0uJjzPuHd6De4f3UIg6gUYzAF1q7/zurZn+/WF8Y2hwLNXfPt7AyMc+ZeGWfV6XJiLSaBwagJ7oSwQg0ZeoAegnUFhayY9fXcqdz8xnx/5S+ndswdT7hvLgFb0UpCIQ73UBcnzJfh8Pj+rHNQPa81+vLmFD7kFunvg5XxuSxUNX9tbMsyIS88IHoAMagH4CH67aw8/fWM7uojL88XE8cHlP7rmoO/E+hahI6TRflCirDPDXmeuY+MlGAtWOrq1T+J8bBzCkh/7DEJHYpdN8J5dfXM5vp63kzcU7ARjYpSV/vPlMerTRAPNTodN8TUBSgo8fXdWHt+4dQp92zdmSX8JXnp7HDyYvJr+43OvyREQ8oXmmjs85x+QvtnHZnz7hzcU7SUqI4xej+jHl2xcqSNUxneaLMv07pjH1vqE89elGHvtwHa8v2sHM1Xv52Yi+jBnUCTPNByIisePQab6yQBmGURYo02k+YP3eYn72xjLmbwreDe6inhn87ob+dG2d6nFlTZN6pqKQPz6Oe4f3YMYDwxjaI4P9JZX816tLufXJuWzILfa6PBGRBlVQVkB2WjYA2WnZMT3PVFllgD+9v5Zr/jKL+ZsKyGjm5y+3ns3zXztXQaoeacxUlHPO8dbinfx22kryD1bg98XxnUuy+e7wbBLjNUBdRJo2jZn6j8835PHwG8vZmHcQgFsHd+Yn1/ShZYrf48qaBo2ZasLMjBvO6cgHP7iYW3I6UxGo5i8fruOav3zK3I2x+9eZiMSGQ2OmfAT/ePThi7kxUwUHK/jh5CXc/tQ8NuYdpEebZkz+1gU8ctOZClINRD1TTcy8jfn87I1lbMgN/mUyemBHfnpNXzKbJ3pcmYhI3YvlnqnqasfkBdv43+mr2VdSiT8+jvuH9+BbF2drzqh6cKKeKYWpJqi8KsDEjzfyxEfrqQhU0zwxngev6MVdF3TVfCIi0qTkluRy2zu3kVeaR8AF8JmPjOQMXh71cpMehL5k235++dZyloRuN3Zhdmt+d0N/umc287iypkun+WJMYryP71/ek/ceHMbw3pkcKK/iN9NWMuqvs5mnU38i0oRkpmQyrNMwAi6AYQRcgIs7X9xkg9S+gxX89PVl3PC3z1iyvZC2LRJ57LZz+Nc3zlOQ8pCmRmjCumWk8uy4wXy4ai+/nraC1bsPcMuTc7nh7A78bERf2rRI8rpEEZGIHbqab2PhxiZ7NV+g2vHyF1v5fzPWsL+kkvg44+sXZXH/ZT1plqhf5V7Tab4YUVYZYOInG5jw8QbKq6pplhjPA5f35O4Lu5GgU38iEqViYczUl1v38cu3VrBsR/CU3pAerfn1dWdo4s0GpjFTctjW/BJ+M20lH6zaA0Cvts341bVncKFuSyMiUSi3JJdHFzzKzK0zKQuUkeRL4rIul/HQ4Iei/lRf7oFyHp2xhlcWbAOgfVoSD4/sx4gB7TRBswdOFKbUNxhjurRO4em7c/ho9V5+9fYK1u4p5van53FFv7b8bERfsjI0qZuIRI/wGx37ff4mcaPjssoA//hsM098tJ7i8ioSfMY3L+rOfZf2IMWvX9uNkY5KjBrepw0XZLfmmdmb+NtH63l/5R4+XrOXuy7oxvcu7UlaSoLXJYqI1EpBWQFje49lTK8xTFk7hbzSPK9LOi3OOd5dvps/vLuKbQWlAFzWpw0/H9lXg8sbOZ3mE/YWlfHoe2uYsnA7zkGrlAQevKIXt5/bRVMpiIg0gGXbC/nttJXM3xy8l17vts15eFRfLuqZ6XFlcojGTEmtLN8R/DDPC90Ys0ebZjw8si+X9G7jcWUiIk3TnqIy/jh9Da8t2g5AeqqfH17Zi1tyOuuP2UZGYUpqzTnHjBV7+MO7q9iSXwLAxb0yeXhkX3q21ZUjIiJ1obQiwFOfbmTCxxsorQyQ4DO+NiSLey/tQYskDbNojBSm5JSVVwV47vPN/PXD9RworyLOYMygzjxwRU/apyV7XVlDWGMAABy2SURBVJ6ISFSqClQzecF2xn+wlr0HygG4+ox2/HREH7q21gVAjZnClJy2/OJyxn+wjpfmb6Wq2pEYH8dXh2TxnYuzNUhdRKSWgr3+u/njjDVsDN07dUDHNH4+si/nd2/tcXVSGwpTErFNeQd59L01vLN0FwBpyQncOzybuy7oRlKCz+PqREQar3kb8/nDu6tZvG0/AF1bp/Cjq3ozon974uI0X1S0UJiSOrNk234eeXc1c0L3+GuflsSDV/TipoGd8Ok/BRGRw1bvLuKP09cwc/VeADKa+fn+ZT259dwuuvNEFFKYkjrlnGPWujweeXc1q3YVAcGZ1H90VR8u79tGM/OKSEzbVlDC+A/W8fqXwelmUv0+7hmWzTcuyiJV99GLWgpTUi+qqx1Tl+zk0ffWsH1fcIK5szql8eAVvbi4V6ZClYjElJ37S3nio/VMXrCNyoAjwWd85byu3HdpDzKaJXpdnkRIYUrqVXlVgBfnbeWJjzaQVxy8OmVgl5b84IreDOnRWqFKRJq0PUVl/O2j9bw0fxsVgWriDK4/uyMPXt6LLq1TvC5P6ojClDSI0ooAk+ZuYcInGyg4GLyL+7lZ6fzgil66WkVEmpzcA+VM/GQDk+ZuobyqGjMYdWYHvn9ZD3q00bx8TY3ClDSog+VVPDdnM0/O2sj+kkoALsxuzQ+u6EVOt3RvixMRiVDBwQr+/skGnpuzmbLKagCu6d+OBy7vRe92ClFNlcKUeOJAWSX/+GwzT326kQNlVQBc1DOD+y/tyblZClUiEl3yist5ZvYmnv98MwcrAgBc0a8tD1zekzM6pHlcndQ3hSnxVGFpJc/M3sSzszdRXB4MVYO7teLe4T00UF1EGr2d+0t5ctZGXpq/lfKqYE/U8N6ZPHhFL87s1NLj6qShKExJo7DvYAX/+Hwz//xsE0WhnqoBHdO4d3g2V/Zrp8nrROS05Zbk8qNZP+LRix8lIzmjTva5Oe8gEz7ewOtfbqcyEPxdeXnfttw7PJtzurSqk/eQ6KEwJY3KgbJK/jVvK09/uunw1X892jTju5dkc91ZHXSndBE5Zb+d+1umrJnCmN5j+MX5v4hoX6t3F/HERxt4Z+lOqh3EGYw8swPfvSSbvu1b1FHFEm0UpqRRKqsMMHnBNv7+yUZ27A/OU9U5PZlvDcvm5kGddJsaETmpQZMGURGoOGa53+dn4R0LT2lfi7bu428fbeCDVXsASPAZo8/pxLcvySYrQzchjnUKU9KoVVRV8+biHUz8eAMb84I3AG2d6ufOC7py5/ldaa3J7kTkOHJLcnl0waPM3DqTskAZSb4kLutyGQ8NfqhWp/sC1Y73V+7hqU83snDLPgCSEuK4dXAX7hnWnQ4tk+u7CRIlThSmNK+9eM4fH8fYnM7cNLAT7y7fxd8/2ciyHYWM/2AdEz7ewOiBnfjGRVlkZzbzulQRaWQyUzJJTUilPFCO3+enPFBOqj/1pEGqpKKKVxdu55nZm9iSXwJAi6R47ji/K18bmqUZy+WUKExJo+GLM0ad2YGRA9ozb1MBT83ayIer9/LS/K28NH8rl/VpwzeHdee8rHRdASgihxWUFTC291jG9BrDlLVTyCvNO+62e4vKeG7OZv41b+vhefA6pyfz9SFZjMnprHvnyWnRaT5p1NbvLeaZ2Zt4fdH2w5ckD+iYxjcuymLEgPa687qI1Mrq3UU8/ekmpi7eSUUg+H/JOV1a8s2LunPVGe3w6WpiOQmNmZKol19czgtzt/DCnC3kh25V07ZFIref25Xbzu1MmxZJHlcoIo1NVaCa91fu4bk5m5m7sQAAM7iqXzu+OSyLQV01ebDUnsKUNBlllQHe+HIHz87exLq9xQDExxnXDGjPXRd0JadrK50CFIlB4fNMEWjOy/O38q95W9lVWAZAqt/HTYM68bUhWXTTlXlyGhSmpMlxzjFnYz4vzNnCeyv3EKgO/jvu274Fd1/QlevP7kiyX1MriMSK3875LVPWTqGDbzgbV119+FRe98xU7r6gG6MHdqR5UoLHVUo0U5iSJm3n/lJenLeVl7/YSl5x8BRgi6R4xuZ05vbzutBdVwGKNFmDXhhERfWx80wZCfztwvcY0qO1equlTihMSUworwrw7rLdPDdnM19u3X94+blZ6dx2bmeu6d9eE4GKNAHOOZZsL+Tl+VuZunw11a3eJr75CiyuEh9+Luo4nP8e8pM6u62MCGieKYkRifE+bjinIzec05Fl2wuZNHcLby/dyfxNBczfVMAv31rBjed05NbBXejXQbeEEIk2haWVvLV4By/O28rq3QdCS1Pp3KkFhXFVJMT5qayupE2zNAUpaVC16pkys6uBvwA+4Gnn3CNHrR8H/D9gR2jR4865p0+0T/VMSUM4UFbJtKW7eHn+VpZsLzy8/MxOadwyuDPXntWBFhpHIdJoVVc75m7K59WF2/n3sl2UVQbHQrVKSeDmQZ24ZXAXHl/xMKkJqazft54erXpwsPIg44eP97hyaWoiOs1nZj5gLXAFsB34ArjNObcybJtxQI5z7r7aFqUwJQ1t5c4iXvliK298uYOisioAEuPjuPKMdowe2JGLemToJssijcTG3GJeX7SDN77ccfjenQAXZrfmtnO7cOUZbUmM/89p+9/O+hlTNr7NmO7X8oth/+NFydLERRqmLgB+5Zy7KvT6pwDOuT+EbTMOhSmJEmWVAaYv383LX2w9PPcMQEazRK4/uwOjB3akX/sWGrQq0sAKSyp5e+lOXl+0nUVh4x47tkzmxnM6cvOgTsdMa1CXNzoWOZFIx0x1BLaFvd4OnFfDdjeZ2TCCvVgPOue21bCNiOeSEv4ztmr7vhLe/HIHry/awca8gzwzexPPzN5E77bNGT2wI9ef3ZF2aZoQVKS+lFUG+GRtLm8t3sEHK/centIgxe9jxID2jB7YkfOzWhN3nBnKp2/dwaNpKcxMSaYsLo6k6mouKynhocLSGrcXqQ+16Zm6GbjaOfeN0Os7gfPCe6HMrDVQ7JwrN7NvAbc45y6tYV/3APcAdOnSZdCWLVvqriUiETh0ddDri7bz9pKd7Avds8sMBndNZ9RZ7bmmf3sym+vmpyKRqgxUM3t9Hm8v2cn7K/ZwoDx42t0MhmRnMHpgR67u344Ufy3+3j+wm99M/QqvVu4hwTkqzRiT0JZfXPciNG9bzy2RWFLvp/mO2t4HFDjn0k60X53mk8aqoqqaj9fs5fVFO5i5Zi8VoXsCxhmc3701o87swNX925Ge6ve4UpHoEah2zNuYz9tLd/Hu8l2HbzIMcEaHFow6swPXn92BDi2TT3nfD7x4KRn5mxhTUsGUFD95rbMYf/vMuixfJOIwFU/w1N1lBK/W+wK43Tm3Imyb9s65XaHnNwI/ds6df6L9KkxJNDhQVskHq/YwbckuZq3LpTIQ/Lz44owhPTIYNaA9l/VtQ+tm6rESOVploJq5G/OZsWI3M1bsIfdA+eF1Pdo047qzOjDqzPaRT6z78legWVvI+Sos+AcU74Fb/xVh9SJHinjSTjMbAYwnODXCs86535vZb4AFzrmpZvYH4DqgCigAvuOcW32ifSpMSbQpLKlkxsrdTFu6i8/W5x2+hU2cQU63dK46ox1X9mtL5/QUjysV8U5JRRWz1uYyY8UePly15/CVswBdW6cw6sz2XHtWB3q3ba6LPCSqaAZ0kTpWcLCCGSt2M335bj7fkHe4xwqC9we8sl9brjyjra4KlJiw72AFM1fvZcaK3cxal3t4LigI9kBddUZbrjqjHQM6punzIFFLYUqkHh0oq+SjNbm8t2I3H6/Jpbj8P3+Jd2yZzPA+mVzSqw0X9mhduwG1Io2cc44VO4v4eM1eZq7ey+Jt+6kO+1VyVueWhwNUtu6NKU2EwpRIAymvCjBnQz4zVuzh/ZV7yCv+zxgRvy+O87qnc0nvNgzvnUlWRqr+SpeoUVxexex1uXy0OpeP1uxlb9j4pwSfcV5Wa648oy1X9mun6USkSVKYEvFAdbVj2Y5CPlqzl4/X5LJk+37CP25d0lO4pHcmQ3pkcH731qQl67Y20nhUBapZuqOQz9bl8dmGPBZu2XfE6ey2LRIZ3rsNw/u0YUiPDJolqtdVmjaFKZFGIL+4nFnrcvl4TS6frM094tLwOIP+HdO4MDuDIT1ak9M1nWS/7wR7E6lbzjk25BYze10en23IZ+6G/MPzP0Hw3+jALq0Y3qcNw3u3oW97DSCX2KIwJdLIBKodS7bvZ9baXD5fn8+X2478q9/vi2Ng15YMyc5gcFY6Z3duSVKCwpXUnUPhaf6mfXyxuYDPN+Sxp6j8iG26tU5hSI8MhvbI4ILs1rRM0dxqErsUpkQauZKKKuZvKmDOhnw+25DHip1FR5wSTPAZ/TumMbhbOoO7pTOoaytNGiqnpDJQzYqdRXyxqYD5mwtYsLng8Ez/h2Q083NhdjA8XdijNZ1aaZoPkUMUpkSizL6DFczdmM+cjfks2LyPVbuPDFcQvOQ8p2srzunSkjM7taRnm2bE++K8KVganT1FZSzetp8l2/azOPQoqQgcsU2b5okMzkrn3G7pnNc9XXM/iZyAwpRIlCsqq2TRln0s2Bw8JbN4237Kq6qP2CY5wUf/ji04s1NLzurckrM6pdElPUW/HGNAcXkVS7fvZ8m2QpZs28+S7fvZVVh2zHZZGakM7taKwd3SOTcrXf8+RE6BwpRIE1NRVc3ynYUs3LyPJduDvzy3FZQes13LlATO6NCCvu1a0Kd9C/q2b06PNs1IjNf4q2jknGNXYRmrdhWxcmcRq3YHv24pKDmm57J5Yjxndk7jrFC4PqdLS9o0b7pTFqzOX81XZ3yV565+jl7pvbwuR5oghSmRGFBwsOJw78TS7ftZsr3wiHmuDomPM7Izm9G3fXP6tG9B77bNyc5sRsdWyfji1EvRWBQcrGBDbjHr9xazbk8xq3YFw9P+o8Y5QXBMXd/2LQ4Hp7M7t6R7RipxMXQ8R70+ii0HttC1eVemjZ7mdTnSBClMicSg8F6M4OMAq3YVsSn/4DG9GAD++Di6Z6SSndmM7MxUsts0IzuzGd0yUjWHUD2pqKpmx/5StuQfZP3eYjbkFrNh70HW5xZTcLCixu9plZJAv1BvY9/2LejXoQXZmc3wx8fmeLkBzw047rpldy9rwEqkqTtRmNL/kCJNlJnRoWUyHVomc1nftoeXl1YEWLPnwOGQdegX+O6iMlbvPsDq3QeO2VerlAQ6p6fQuVUKndKT6RJ63jk9hQ4tk3Ta8DgC1Y784nJ2FpaxtaCEbQUlbM0vYWtB8LGrsPSI27CES/X76NGm2eFQ27d9c/q1T6Nti0SNcwqTEJdAZXUNvXVxmgRXGo7ClEiMSfb7ODt0KihccXkVG3P/0zuyMS94imlLfgn7SirZV1LI0u2FNe4zPdVPm+aJtEtLom3zJNq2SKRt6HmbFom0SvGTnuonxe9rEkGgoqqafSUVFBwMPvKKy9lVWMbuwjL2FJWxuyj4fO+BcgLHS0sEJ8Ls2DIYTrPbBHsFe7QJPtq1SGoSP6v6VlOQOtFykfqgMCUiADRLjOfMTsFpFsI558gtLmdbQSnbQr0r2/aVsK2glK0FJewuKjscKmrq1Qrnj48jPcVPy5QE0lP9tEr10zI5gWZJ8TTzx5OaGE+zxODX1ETf4ef++Dj8vjj88XEkHP5qJMTF1WpcUHW1oyJQTXllNeWBABVV1cFHIPi1pCJAcVkVxeVVHCivCj2vpLgs+LqwpJKCQ+GpuOKImcFPpnWqn3ZpSXRJTwn26IW+dklPoUPL5Jg9PVdXpoyawphpY45Z/tq1r3lQjcQqhSkROSEzo03zJNo0T2JQ11bHrA9UO/IPlrO3qDzYM3OgjD2FZewpKmfPgTL2FpWzv6SC/IMVlFdVB3ttio69bP90xccZhzpwwseCHXrqnDvuqbTT5YszWh0KhCl+Mpon0q5FUvCRFnq0CPbK6RRo/erTug8p8SmUVJUcXpYSn6Ir+qRBKUyJSER8cf8JW/07pp1w29KKwOHTY4e+FpZWcqCsioPlwUdxeSD4vCLYU3SwvIrKgDvck1QZqKayqjq4LFBNVS2Tkj8+jsTQ41AvV2K8D398HMl+H80T44M9ZIlhj9DrFskJtE4NnqpMT/XTIikhpq6Ua8wGTRpEReDIwfolVSUMmjSIhXcs9KgqiTUKUyLSYJL9PpL9wUHxdcE5d8Q9DQHChxkdeuqLM40/aqKmj57OowseZfqm6VRTTRxxXJN1DQ8Nfsjr0iSGKEyJSNQyM/zxCkmx7OrXrz6iZ6qaat7Z9A7vb31fPVPSYDTyUUREotfxzvJ6M4WixCiFKREREZEIKEyJiEjUmn7TdJLjjxyDlxyfzIybZ3hUkcQihSkREYlamSmZpCelA/+Z9Tw9KZ2M5Awvy5IYowHoIiIS1fqk9yGnbQ5r962lV6teFFcWe12SxBj1TImISFQbP3w8zjlWFazC4Rg/fLzXJUmMUc+UiIhEraMn7Zy6YSpTN0zF7/NragRpMOqZEhGRqOVczXMgHG+5SH1QmBIRkahl1Dxp6/GWi9QHhSkREYlamhpBGgOFKRERiVrhUyMcoqkRpKEpTImISNQaNGkQO4p3HLFsR/EOBk0a5FFFEosUpkREJHrp3nzSCChMiYhI1Jp+03Q6N+98xLIuzbtozJQ0KIUpERGJWpkpmVS7auA/V/AFXEBjpqRBadJOERGJan3S+5BgCWw+sJluLbqR3TLb65IkxihMiYhI1Dp6BvTNRZvZXLSZQZMGaQZ0aTA6zSciIlFLM6BLY6AwJSIiUauyuvKUlovUB4UpERGJWglxCae0XKQ+KEyJiEjUmnHTDHzmO2KZz3y8d/N7HlUksUhhSkREolZmSubhXqhDUyMkxCVoagRpUApTIiIS1YZ0HEK35t1wOLq16MaQjkO8LklijKZGEBGRqKWpEaQxUM+UiIhELU2NII1BrcKUmV1tZmvMbL2Z/aSG9Ylm9kpo/Twz61bXhYqIiBxNUyNIY3DSMGVmPuAJ4BqgH3CbmfU7arOvA/uccz2APwP/W9eFioiIiDRGtemZOhdY75zb6JyrAF4Grj9qm+uB50LPXwUuMzOruzJFREREGqfahKmOwLaw19tDy2rcxjlXBRQCrY/ekZndY2YLzGxBbm7u6VUsIiIS8vB5D9e4/L/P/+8GrkRiWYMOQHfOPemcy3HO5WRmZjbkW4uISBN0S59balx+c++bG7gSiWW1CVM7gM5hrzuFltW4jZnFA2lAfl0UKCIicjKGMbLbyMMTd4o0pNrMM/UF0NPMsgiGpluB24/aZipwNzAHuBmY6XRdqoiINIBldy87/PyRix/xsBKJVScNU865KjO7D5gB+IBnnXMrzOw3wALn3FTgGeAFM1sPFBAMXCIiIiJNXq1mQHfO/Rv491HLfhn2vAwYU7eliYiIiDR+mgFdREREJAIKUyIiIiIRUJgSERERiYDClIiIiEgEFKZEREREIqAwJSIiIhIBhSkRERGRCJhXE5WbWS6wpQHeKgPIa4D3aYzU9tgVy+2P5bZDbLdfbY9dDdH+rs65Gm8s7FmYaihmtsA5l+N1HV5Q22Oz7RDb7Y/ltkNst19tj822g/ft12k+ERERkQgoTImIiIhEIBbC1JNeF+AhtT12xXL7Y7ntENvtV9tjl6ftb/JjpkRERETqUyz0TImIiIjUm6gPU2Y2xsxWmFm1mR13JL+ZXW1ma8xsvZn9JGx5lpnNCy1/xcz8DVN53TCzdDN738zWhb62qmGb4Wa2OOxRZmY3hNb908w2ha07u+FbcXpq0/bQdoGw9k0NWx4Lx/5sM5sT+owsNbNbwtZF3bE/3uc4bH1i6FiuDx3bbmHrfhpavsbMrmrIuutCLdr+AzNbGTrOH5pZ17B1NX4Gokkt2j/OzHLD2vmNsHV3hz4n68zs7oatPHK1aPufw9q91sz2h62L6mNvZs+a2V4zW36c9WZmj4V+NkvNbGDYuoY77s65qH4AfYHewMdAznG28QEbgO6AH1gC9AutmwzcGno+EfiO1206xfb/EfhJ6PlPgP89yfbpQAGQEnr9T+Bmr9tRn20Hio+zvMkfe6AX0DP0vAOwC2gZjcf+RJ/jsG2+C0wMPb8VeCX0vF9o+0QgK7Qfn9dtquO2Dw/7XH/nUNtDr2v8DETLo5btHwc8XsP3pgMbQ19bhZ638rpNddn2o7a/H3i2CR37YcBAYPlx1o8A3gUMOB+Y58Vxj/qeKefcKufcmpNsdi6w3jm30TlXAbwMXG9mBlwKvBra7jnghvqrtl5cT7BuqF39NwPvOudK6rWqhnGqbT8sVo69c26tc25d6PlOYC9Q46RzUaDGz/FR24T/TF4FLgsd6+uBl51z5c65TcD60P6ixUnb7pz7KOxzPRfo1MA11qfaHPvjuQp43zlX4JzbB7wPXF1PddaHU237bcBLDVJZA3DOzSLYAXA81wPPu6C5QEsza08DH/eoD1O11BHYFvZ6e2hZa2C/c67qqOXRpK1zblfo+W6g7Um2v5VjP2i/D3WP/tnMEuu8wvpT27YnmdkCM5t76PQmMXjszexcgn/ZbghbHE3H/nif4xq3CR3bQoLHujbf25idav1fJ/jX+iE1fQaiSW3bf1Po3/OrZtb5FL+3sap1/aFTu1nAzLDF0X7sT+Z4P58GPe7x9bXjumRmHwDtalj1c+fcWw1dT0M7UfvDXzjnnJkd9/LMUFofAMwIW/xTgr+I/QQvLf0x8JtIa64rddT2rs65HWbWHZhpZssI/pJt9Or42L8A3O2cqw4tbtTHXk6Pmd0B5AAXhy0+5jPgnNtQ8x6i1tvAS865cjP7FsEeyks9rqmh3Qq86pwLhC2LhWPvuagIU865yyPcxQ6gc9jrTqFl+QS7BONDf8UeWt6onKj9ZrbHzNo753aFfmHuPcGuxgJvOOcqw/Z9qGej3Mz+ATxUJ0XXkbpou3NuR+jrRjP7GDgHeI0YOfZm1gJ4h+AfH3PD9t2oj30Njvc5rmmb7WYWD6QR/JzX5nsbs1rVb2aXEwzaFzvnyg8tP85nIJp+oZ60/c65/LCXTxMcU3joey856ns/rvMK68+p/Nu9Fbg3fEETOPYnc7yfT4Me91g5zfcF0NOCV2/5Cf6Dm+qCo9Q+IjiOCOBuINp6uqYSrBtOXv8x59JDv4QPjSG6AajxiolG6qRtN7NWh05fmVkGMARYGSvHPvTv/Q2CYwpePWpdtB37Gj/HR20T/jO5GZgZOtZTgVsteLVfFtATmN9AddeFk7bdzM4B/g5c55zbG7a8xs9Ag1VeN2rT/vZhL68DVoWezwCuDP0cWgFXcmTvfGNXm3/3mFkfggOt54QtawrH/mSmAneFruo7HygM/aHYsMe9vka2N9QDuJHgudByYA8wI7S8A/DvsO1GAGsJJvKfhy3vTvA/1fXAFCDR6zadYvtbAx8C64APgPTQ8hzg6bDtuhFM6nFHff9MYBnBX6STgGZet6ku2w5cGGrfktDXr8fSsQfuACqBxWGPs6P12Nf0OSZ4avK60POk0LFcHzq23cO+9+eh71sDXON1W+qh7R+E/g88dJynhpYf9zMQTY9atP8PwIpQOz8C+oR979dC/ybWA1/1ui113fbQ618Bjxz1fVF/7Al2AOwK/T+2neB4wG8D3w6tN+CJ0M9mGWFX9TfkcdcM6CIiIiIRiJXTfCIiIiL1QmFKREREJAIKUyIiIiIRUJgSERERiYDClIiIiEgEFKZEREREIqAwJSIiIhIBhSkRERGRCPx/9vEgzVKOAC4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = [(x,f_approx),(x_rep,phi_S)]\n",
    "plot_1d(lambda _x: solution(t, _x), samples=samples, titles=['Monte-Carlo Sampling'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Corresponding Learning Problem\n",
    "\n",
    "\n",
    "This is very nice (especially because it overcomes the curse dimensionality).\n",
    "\n",
    "But our goal was to numerically compute the solution for **all** $t\\in [0,\\tau]$ and **all** $x\\in[u,v]^d$\n",
    "(and we want to do machine learning)...\n",
    "\n",
    "Note that the unique solution to the statistical learning problem \n",
    "\n",
    "$$\\min_f \\mathbb{E} [ (f(X)-Y)^2 ] $$\n",
    "\n",
    "is given by the **regression** function $E[Y|X=x]$ (assuming $Y$ has finite variance).\n",
    "\n",
    "\n",
    "1. Let$\\color{blue}{(T,X)}\\sim \\mathcal{U}([0,T]\\times [u,v]^d)$ be uniformly distributed on the hypercube $[0,T]\\times [u,v]^d$. \n",
    "2. Define $\\color{blue}{Y:=\\varphi(S^X_T)}$, where $S^X_T$ is the solution to the SDE with initial value $X$ at time $T$.\n",
    "\n",
    "Then, under sufficient regularity assumptions, it the Feynman-Kac formula establishes that\n",
    "\n",
    "$$f(t,x)=\\mathbb{E}\\big[\\varphi(S_t^x)\\big]=\\mathbb{E}\\big[\\varphi(S_T^X)|(T,X)=(t,x)\\big]=E[Y|(T,X)=(t,x)] \\qquad \\forall t\\in[0,\\tau], \\ x\\in[u,v]^d.$$\n",
    "\n",
    "Given $m$ i.i.d. samples $(X^i,Y^i)\\sim (X,Y)$ we can thus approximate\n",
    "\n",
    "$$\\color{blue}{f^{sol}=\\operatorname{argmin}_f \\mathbb{E} [ (f(T,X)-Y)^2 ]\\approx \\operatorname{argmin}_{f\\in \\mathcal{N}} \\tfrac{1}{m} \\sum_{i=1}^m (f(T_i,X_i)-Y_i)^2 },$$\n",
    "\n",
    "where $\\mathcal{N}$ denotes a suitable hypothesis space of neural networks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "u, v = -1., 1. #approximate solution on the hypercube [u,v]^d\n",
    "tau = 2. #approximate evolution of the solution for [0,tau]\n",
    "m = int(1e+6) #number of training samples\n",
    "m_test = 10000 #number of test samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#training data\n",
    "X = tf.random.uniform([m, d], minval=u, maxval=v)\n",
    "T = tf.random.uniform([m, 1], minval=0, maxval=tau)\n",
    "y_train = batched_MC(sigma, phi, T, X)\n",
    "x_train = tf.concat([T,X], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"Kolmogorov_Model\"\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "t_x (InputLayer)             [(None, 2)]               0         \n",
      "_________________________________________________________________\n",
      "dense (Dense)                (None, 100)               300       \n",
      "_________________________________________________________________\n",
      "leaky_re_lu (LeakyReLU)      (None, 100)               0         \n",
      "_________________________________________________________________\n",
      "dense_1 (Dense)              (None, 100)               10100     \n",
      "_________________________________________________________________\n",
      "leaky_re_lu_1 (LeakyReLU)    (None, 100)               0         \n",
      "_________________________________________________________________\n",
      "dense_2 (Dense)              (None, 1)                 101       \n",
      "=================================================================\n",
      "Total params: 10,501\n",
      "Trainable params: 10,501\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "#model (Keras functional API: https://www.tensorflow.org/guide/keras/functional)\n",
    "inputs = keras.Input(shape=(d+1,), name='t_x')\n",
    "hidden = keras.layers.Dense(100)(inputs)\n",
    "#hidden = keras.layers.BatchNormalization()(hidden)\n",
    "hidden = keras.layers.LeakyReLU()(hidden)\n",
    "hidden = keras.layers.Dense(100)(hidden)\n",
    "#hidden = keras.layers.BatchNormalization()(hidden)\n",
    "hidden = keras.layers.LeakyReLU()(hidden)\n",
    "outputs = keras.layers.Dense(1)(hidden)\n",
    "learn_model = keras.Model(inputs = inputs, outputs = outputs, name = 'Kolmogorov_Model')\n",
    "learn_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#test data\n",
    "X = tf.random.uniform([m_test, d], minval=u, maxval=v)\n",
    "T = tf.random.uniform([m_test, 1], minval=0, maxval=tau)\n",
    "y_test = solution(T,X)\n",
    "x_test = tf.concat([T,X], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 1000000 samples, validate on 10000 samples\n",
      "1000000/1000000 [==============================] - 160s 160us/sample - loss: 0.6515 - mae: 0.5189 - val_loss: 0.0040 - val_mae: 0.0479\n"
     ]
    }
   ],
   "source": [
    "learn_model.compile(optimizer='adam', loss='mse', metrics=['mae'])\n",
    "history = learn_model.fit(x_train, y_train, validation_data=(x_test, y_test), epochs=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ba4a307bc8d04d9fb4d7ed8f76737aba",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.0, description='t', max=2.0), Output()), _dom_classes=('widget-inter…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#visualize the evolution of the (learned) solution\n",
    "def evaluate(model, _t,_x):\n",
    "    return model(tf.concat([repeat(_t,_x.shape[0]),_x], axis=1))\n",
    "\n",
    "interact(lambda t: plot_1d([lambda x: evaluate(learn_model, t, x), lambda x: solution(t, x)], \n",
    "                           titles=['model vs. exact solution']), t=(0,tau));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this was just an illustrating example! \n",
    "\n",
    "real power of the deep learning approach: it can **overcome the curse of dimensionality**! \n",
    "\n",
    "* empirical overservations: \n",
    "\n",
    "  Solving a $\\color{blue}{100\\text{-dimensional}}$ Black-Scholes equation from financial engineering (option pricing) on $[90,110]^{100}$ using dense ReLU networks with architecture $(100,200,200,1)$ [[1]](https://arxiv.org/abs/1806.00421):\n",
    "  \n",
    "|Number of descent steps| Relative $\\mathcal{L}^1$ error| Relative $\\mathcal{L}^\\infty$ error | Runtime in seconds|\n",
    "|-----------------------|-------------------------------|-------------------------------------|-------------------|\n",
    "| \t0                   |  1.004285                     | 1.009524                            |  1                |\n",
    "|  100000               |          0.371515             |0.387978                             | 437.9             |\n",
    "| 250000                |   0.001220                    | 0.010039                            |  1092.6           |\n",
    "|500000                 |  0.000949                     | 0.005105                            |   2183.8          |\n",
    "\n",
    "* theoretical explanations: \n",
    "\n",
    "  approximation and generalization results for Kolmogorov equations, that do not suffer from the curse of dimensionality [[2]](https://arxiv.org/abs/1809.07669), [[3]](https://arxiv.org/abs/1903.06652), [[4]](https://arxiv.org/abs/1809.03062), ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can now try higher dimensions, different architectures, \n",
    "# and different Kolmogorov PDEs (see also the SDE solver below)!!\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have any questions, please feel free to contact me:  julius.berner@univie.ac.at"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further Material"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using our Approximation Framework\n",
    "\n",
    "We can also explicitely construct a neural network approximating the PDE solution.\n",
    "\n",
    "Note that this is only the case in our specific simple example, where we know the solution $f(t,x) = \\|x\\|_2^2 + t * \\operatorname{Trace}(\\sigma \\sigma^*)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = max([u,v])\n",
    "trace = tf.linalg.trace(sigma@tf.transpose(sigma)).numpy()\n",
    "\n",
    "approx_param = Affine_Comb([trace,1.],[P('identity'),Square(1e-10,B)])\n",
    "approx_model = R(approx_param)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "64d8567b3df14b0f9bd0b686c4df2af1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.0, description='t', max=2.0), Output()), _dom_classes=('widget-inter…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualizing the learned model, the approximating using the explicit construction, and the exact solution\n",
    "interact(lambda t: plot_1d(lambda x: evaluate(learn_model, t, x), \n",
    "                           lambda x: evaluate(approx_model, t, x),\n",
    "                           lambda x: solution(t, x),\n",
    "                           titles=['learned', 'approx.', 'exact']), t=(0,tau));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Approx. Model #Parameters: 10501\n",
      "Learned Model #Parameters: 2821\n"
     ]
    }
   ],
   "source": [
    "# Note that the learned model is significantly smaller.\n",
    "print('Approx. Model #Parameters: {}\\nLearned Model #Parameters: {}'.format(learn_model.count_params(),\n",
    "      approx_model.count_params()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 1000000 samples, validate on 10000 samples\n",
      "1000000/1000000 [==============================] - 136s 136us/sample - loss: 0.6494 - mae: 0.5180 - val_loss: 0.0023 - val_mae: 0.0364\n"
     ]
    }
   ],
   "source": [
    "# We could also ftry to urther learn the approx. model \n",
    "approx_model.compile(optimizer='adam', loss='mse', metrics=['mae'])\n",
    "history = learn_model.fit(x_train, y_train, validation_data=(x_test, y_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can try to generalize this to the multidimensional setting by \n",
    "# by implementing a network approximating the squared euclidean norm!\n",
    "# \n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Approximation of SDE solution\n",
    "\n",
    "We now show how to approximately solve the SDE for more complicated coefficient functions.\n",
    "\n",
    "Let us partition the time interval\n",
    "$[0,T]$\n",
    "into\n",
    "$N$ \n",
    "subintervals of equal length, i.e. we set\n",
    "\n",
    "$$t_n:=\\frac{nT}{N}, \\quad n\\in\\{0,1,\\dots,N\\}.$$ \n",
    "\n",
    "Note that for sufficiently small step size \n",
    "$\\frac{T}{N}=t_{n+1}-t_n$ \n",
    "one can approximate \n",
    "\n",
    "\n",
    "$$\n",
    "\\color{blue}{ S^x_{t_{n+1}}} = S^x_{t_n} + \n",
    "\\int_{t_n}^{t_{n+1}} \\mu\\big(S^x_s\\big)\\,ds + \n",
    "\\int_{t_n}^{t_{n+1}} \\sigma\\big(S^x_s\\big)\\,dB_s \\approx \n",
    "\\color{blue}{S^x_{t_n} + \\mu\\big(S^x_{t_n}\\big)\\,\\tfrac{T}{N} + \\sigma\\big(S^x_{t_n}\\big)\\underbrace{\\big(B_{t_{n+1}}- B_{t_n}\\big)}_{\\sim \\mathcal{N}\\big(0,\\frac{T}{N}I\\big)}}.\n",
    "$$\n",
    "\n",
    "\n",
    "This defines the so-called **Euler-Maruyama scheme** $\\mathcal{S}^x_n$, which we will implement in the following for affine linear coefficient functions:\n",
    "\n",
    "\n",
    "**input:**\n",
    "\n",
    "* initial value $x\\in\\mathbb{R}^{batchsize \\times d}$\n",
    "* time horizon $T\\in(0,\\infty)$\n",
    "* number of Steps $N\\in\\mathbb{N}$\n",
    "* value function $\\varphi\\colon\\mathbb{R}^d\\to\\mathbb{R}$\n",
    "* affine linear diffusion $\\sigma\\colon\\mathbb{R}^d\\to\\mathbb{R}^{d \\times d}$    \n",
    "  $\\sigma=[\\sigma_{ten},\\sigma_{mat}]$, $\\sigma(x)=[x \\sigma_{ten}[0,:,:],\\dots,x \\sigma_{ten}[d,:,:]]^T + \\sigma_{ten}$ (acting on row vectors)\n",
    "  \n",
    "* affine linear drift $\\mu\\colon\\mathbb{R}^d\\to\\mathbb{R}^{d}$      \n",
    "  $\\mu=[\\mu_{mat},\\mu_{vec}]$, $\\mu(x)=x \\mu_{mat} + \\mu_{vec}$ (acting on row vectors)\n",
    "  \n",
    "**output:**\n",
    "\n",
    "* realizations of $\\varphi$ evaluated at end value of the Euler-Maruyama scheme \n",
    "\n",
    "$$(\\varphi(\\mathcal{S}^{x_i}_{N}))_i\\in\\mathbb{R}^{batchsize}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def EM(x, T, N, phi, mu, sigma): \n",
    "    \n",
    "    batch_size, d = tf.shape(x)\n",
    "    \n",
    "    with tf.name_scope('Euler-Scheme'):\n",
    "        stepsize=tf.divide(T,tf.cast(N,tf.float32), name='Step-Size')\n",
    "        dw_all = tf.random.normal(shape=[N,batch_size,d], mean=0.0, stddev=tf.sqrt(stepsize), name='DW')\n",
    "        realization = x\n",
    "        for dw in dw_all[...,:,:]:\n",
    "            sigma_x = tf.einsum('ij, kjl -> ikl', realization, sigma[0]) + sigma[1] \n",
    "            # einsum[i,k,l]=\\sum_j realization[i,j]*sigma[0][k,j,l]\n",
    "            # equivalent: einsum[i,k,:]=realization[i,:] @ sigma[0][k,:,:]\n",
    "            # resulting dim: [batchsize,d,d] + [d,d] (broadcasting applies to addition)\n",
    "            mu_x = tf.matmul(realization, mu[0]) + mu[1]\n",
    "            realization += (mu_x*stepsize+tf.einsum('ij, ijk -> ik',dw,sigma_x))\n",
    "            # einsum[i,k]=\\sum_j dw[i,j]*sigma_x[i,j,k]\n",
    "            # equivalent: einsum[i,:]=dw[i,:] @ sigma_x[i,:,:]\n",
    "            # resulting dim: [batchsize,d] + [batchsize,d] \n",
    "        return phi(realization)"
   ]
  }
 ],
 "metadata": {
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}