Flattening and Full Rank ADVI in PyMC4

advi.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import arviz as az\n",
    "import collections\n",
    "import itertools\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pymc4 as pm\n",
    "import seaborn as sns\n",
    "import tensorflow as tf\n",
    "import tensorflow_probability as tfp\n",
    "from typing import Dict\n",
    "\n",
    "from tensorflow_probability.python.internal import dtype_util\n",
    "\n",
    "tfd = tfp.distributions\n",
    "tfb = tfp.bijectors\n",
    "\n",
    "plt.style.use('arviz-darkgrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.random.normal(12, 2.2, 200)\n",
    "\n",
    "@pm.model\n",
    "def model():\n",
    "    mu = yield pm.Normal('mu', 0, 10)\n",
    "    sigma = yield pm.Exponential('sigma', 1)\n",
    "    ll = yield pm.Normal('ll', mu, sigma, observed=data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['model/mu', 'model/__log_sigma']\n"
     ]
    }
   ],
   "source": [
    "model = model()\n",
    "state, deterministics_names = pm.inference.utils.initialize_sampling_state(model)\n",
    "unobserved_keys = state.all_unobserved_values.keys()\n",
    "print(list(unobserved_keys))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "VarMap = collections.namedtuple('VarMap', 'var, slc, shp, dtyp')\n",
    "\n",
    "\n",
    "class ArrayOrdering:\n",
    "    \"\"\"\n",
    "    An ordering for an array space\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, free_rvs: Dict[str, tf.Tensor]):\n",
    "        self.free_rvs = free_rvs\n",
    "        self.by_name = {}\n",
    "        self.size = 0\n",
    "\n",
    "        for name, tensor in free_rvs.items():\n",
    "            flat_shape = int(np.prod(tensor.shape.as_list()))\n",
    "            slc = slice(self.size, self.size + flat_shape)\n",
    "            self.by_name[name] = VarMap(name, slc, tensor.shape, tensor.dtype)\n",
    "            self.size += flat_shape\n",
    "\n",
    "    def flatten(self):\n",
    "        flattened_tensor = [tf.reshape(var, shape=[-1]) for var in self.free_rvs.values()]\n",
    "        return tf.concat(flattened_tensor, axis=0)\n",
    "    \n",
    "    def split(self, flatten_tensor):\n",
    "        flat_state = dict()\n",
    "        for param in self.free_rvs:\n",
    "            _, slc, shape, dtype = self.by_name[param]\n",
    "            flat_state[param] = tf.cast(tf.reshape(flatten_tensor[slc], shape), dtype)\n",
    "        return flat_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "@tf.function(autograph=False)\n",
    "def logpfn(*values, **kwargs):\n",
    "    if kwargs and values:\n",
    "        raise TypeError(\"Either list state should be passed or a dict one\")\n",
    "        \n",
    "    val = ArrayOrdering(state.all_unobserved_values).split(values[0])\n",
    "    _, st = pm.flow.evaluate_meta_model(model, values=val)\n",
    "    return st.collect_log_prob()\n",
    "\n",
    "def vectorize_logp_function(logpfn):\n",
    "    def vectorized_logpfn(*q_samples):\n",
    "        return tf.vectorized_map(lambda samples: logpfn(*samples), q_samples)\n",
    "\n",
    "    return vectorized_logpfn\n",
    "\n",
    "target_log_prob = vectorize_logp_function(logpfn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _build_posterior():\n",
    "    order = ArrayOrdering(state.all_unobserved_values)\n",
    "    flattened_shape = order.size\n",
    "    dtype = dtype_util.common_dtype(state.all_unobserved_values.values(), dtype_hint=tf.float32)\n",
    "    loc = tf.Variable(tf.random.normal([flattened_shape]), name=\"mu\")\n",
    "    scale_tril = tfb.FillScaleTriL(\n",
    "        diag_bijector=tfb.Chain(\n",
    "            [\n",
    "                tfb.Shift(tf.cast(1e-3, dtype)),  # diagonal offset\n",
    "                tfb.Softplus(),\n",
    "                tfb.Shift(tf.cast(np.log(np.expm1(1.0)), dtype)),  # initial scale\n",
    "            ]\n",
    "        ),\n",
    "        diag_shift=None,\n",
    "    )\n",
    "\n",
    "    cov_matrix = tfp.util.TransformedVariable(\n",
    "        tf.eye(flattened_shape), scale_tril, name=\"sigma\"\n",
    "    )\n",
    "    return tfd.MultivariateNormalTriL(loc=loc, scale_tril=cov_matrix)\n",
    "\n",
    "    \n",
    "approx = _build_posterior()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = pm.adam()\n",
    "\n",
    "@tf.function(autograph=False)\n",
    "def run_approximation():\n",
    "    losses = tfp.vi.fit_surrogate_posterior(\n",
    "        target_log_prob_fn=target_log_prob,\n",
    "        surrogate_posterior=approx,\n",
    "        num_steps=40000,\n",
    "        optimizer=opt\n",
    "    )\n",
    "    return losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /usr/local/lib/python3.7/site-packages/tensorflow_probability/python/math/minimize.py:77: calling <lambda> (from tensorflow_probability.python.vi.optimization) with loss is deprecated and will be removed after 2020-07-01.\n",
      "Instructions for updating:\n",
      "The signature for `trace_fn`s passed to `minimize` has changed. Trace functions now take a single `traceable_quantities` argument, which is a `tfp.math.MinimizeTraceableQuantities` namedtuple containing `traceable_quantities.loss`, `traceable_quantities.gradients`, etc. Please update your `trace_fn` definition.\n"
     ]
    }
   ],
   "source": [
    "mean_field = run_approximation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(2,), dtype=float32, numpy=array([12.167437 ,  0.7683811], dtype=float32)>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "approx.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAAGWCAYAAAAuQ2TEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de1yUVeLH8e8omAFesHQrNKl2Mc2sTa002+6alzazbE3dLNPNLTfT9Zd2sdq0crt4bTc1s3Y3LbU0y6zcSlNT1MB0TUtM0cQbKiKIIDDn98fTACMDMsMMM8/M5/16+QLOcztnHhy+c57nOcdhjDECAACALdQKdgUAAABQdYQ3AAAAGyG8AQAA2AjhDQAAwEYIbwAAADZCeAMAALARwhsAAICNEN4AAABshPAGAABgI1HBroC/ZGVlBfwYDRo0UHZ2dsCPE4oiue1SZLc/ktsuRXb7I7ntUmS3n7YHvu3x8fE+b0vPmxdq1YrclyuS2y5Fdvsjue1SZLc/ktsuRXb7aXtoC/0aAgAAoAThDQAAwEYIbwAAADZCeAMAALARwhsAAICNEN4AAABshPAGAABgI4Q3AAAAGyG8AQAA2AjhDQAAwEYIbwAAADZCeAMAALARwhsAAICNEN6q6KPFRn959JgKCkywqwIAACJYVLArYBcvvWIkFapVS4f63B3s2gAAgEhFz5uXjh+n5w0AAAQP4S3Itmw1OnSIQAgAAKqGy6ZBlJZm9Kc/W8Ft1XJHkGsDAADsgJ63INq0Odg1AAAAdkN4AwAAsBHCWw3IyzPavZv72gAAQPUR3mrAPX806nuv0ebvCXAAAKB6CG814PBh6+uKVYQ3AABQPYQ3AAAAGyG8AQAA2AjhDQAAwEYIbwAAADZCeAMAALARwhsAAICNEN689Pl/g10DAAAQyQhvXtq7V1q23Oixx506doxx2wAAQM0ivPlgzLNGq9dIM2cR3gAAQM0ivFVDdrb/9jV0mFMfLCAMAgCAyhHeQsR3G6WJUwhvAACgcoQ3AAAAGyG81SBHsCsAAABsj/AGAABgI4Q3AAAAGyG8AQAA2AjhDQAAwEYIb37kdBpt/t4oL48hPwAAQGBEBbsC4eSdOdKMmVZwW/6FFBXF86UAAMC/6Hmrhu82SX9/xamcHCuwvflWaY9bVlawagUAAMIZ4a0aDh+WPl4sTZ/JZVIAAFAzCG9+kJER7BoAAIBIQXgDAACwEcIbAACAjRDeAAAAbITwVoMcNTRySG6ukTE8RAEAQDgK2fC2d+9eDRs2TO3bt9fll1+uO++8UwcOHAh2tUJe6gajW3sYjX+Z8AYAQDgKyUF6s7Ky1LdvX1177bV66623VL9+fW3fvl3R0dHBrlrIm/W2Fdo+WSI9/liQKwMAAPwuJMPbG2+8oaZNm2rs2LElZeeff34QawQAABAafLpsumjRIj399NPq1auXWrdurRYtWmjBggWVbrNp0yYNHjxY7dq10+WXX667775bS5Ys8bjusmXL1KpVK/3lL39Rhw4ddOedd2rp0qW+VBUAACCs+BTeJk+erLlz52rv3r1q0qTJaddPTk5W3759lZKSoq5du6pPnz46dOiQhg8frlmzZpVbf8+ePXr33XfVokULvfnmm+rWrZuGDRum9evX+1LdgFv/bbBrAAAAIoVPl03HjRun5s2bKyEhQTNmzNCrr75a4bpFRUUaM2aMHA6HZs+erZYtW0qSHn74Yd11112aMGGCunTpooSEhJJtjDFq06aNhg4dKklq1aqVvv32W82dO1ft27f3pco1rqaeLAUAAJHFp563jh07uoWtyiQnJ2v37t3q0aNHSXCTpHr16mnIkCEqLCzUwoUL3bY5++yzdcEFF7iVXXTRRdq3b58v1a0RO3YYFRcHuxYAACDcBXyokHXr1kmSOnXqVG6Zq+zUy6G//e1vtXv3brey9PR0nXfeeQGqZfXdO5ChOQAAQOAFPLylp6dLkpo3b15uWePGjRUTE6Ndu3a5ld93331KSUnRrFmztGvXLr377rtatmyZ7rnnnkBXFwAAIKQFfKiQ3NxcSdZlUk/i4uKUk5PjVnbZZZdp8uTJmjhxoiZNmqTExERNnjxZV1xxRYXHadCggWrVCmQWPezV2g0aNFR8vKs+1rZ169ZVfHxsyToxMfmSjrttFx8fX51KKjo6W1KRX/Z1Kn/vz24iuf2R3HYpstsfyW2XIrv9tD10heQ4b5J0yy236JZbbqny+tnZ2QGsjfeys48qKkpu98Hl5+crK+tkyc95eeUvtWZlZVW4T2Os++qioip+GqKw0FmlfXkrPj7er/uzm0hufyS3XYrs9kdy26XIbj9tD3zbqxMQA37ZNC4uTpLK9a655ObmVtgrZ3d/fczozrv9dy/cE08Z9bjdKDeX++sAAIhUAQ9viYmJklTuvjZJyszMVF5ensf74cLBuvXS4SP+29/Kb6Tc49LKVf7bJwAAsJeAhzfXuGyrVpVPHK4yu4zdBgAAEGwBD28dOnRQs2bNtHjxYm3durWkPCcnR9OmTVN0dLR69uwZ6GoAAACEBZ8eWJg/f75SUlIkSdu2bSspc43p1rZtW/Xu3ds6QFSUxo0bp0GDBqlfv37q3r27YmNjtXTpUmVkZGjUqFFq2rSpP9oCAAAQ9nwKbykpKeVmRUhNTVVqamrJz67wJklXX3215syZoylTpmjJkiUqKipSUlKSRo4cqW7duvlYdQAAgMjjU3gbP368xo8f79U2bdq00cyZM305nD0xtykAAAiAgN/zBgAAAP8hvAEAANgI4Q0AAMBGQnZ6LNurZBKEDxYaJVZjXGLDBAsAAEQswlsNS91gNHGylb6GD+OpBgAA4B0um9awffuDXQMAAGBnhDcAAAAbIbzVIIefrpJWth/uhwMAILwR3gAAAGyE8AYAAGAjhLdA4UFSAAAQAIQ3AAAAGyG81bDU1NInCnJyglgRAABgS4S3GnQwU/r8v6U/z5zl26OhlT1RytOmAACEN8JbDVq/Ptg1AAAAdkd4C5DjueXLjmTVfD0AAEB4IbwFyMJFXL8EAAD+R3gLkKLiYNcAAACEI8IbAACAjRDeAAAAbITwFiDbtgVu3/6a4B4AANgP4S1AtmwNdg0AAEA4IrwBAADYCOHNhphFAQCAyEV4AwAAsBHCGwAAgI0Q3gAAAGyE8AYAAGAjhDdU6M23nJr6D2ewqwEAAMqICnYF4F/+ehK1qMjorX9Z3/e+0yg+3j/7BQAA1UPPG06rsCjYNQAAAC6ENwAAABshvNlcbq5RURGj9gIAECkIbzaRlVUa0NK2GzmdRkeOGN3aw6j/AMIbAACRggcWbKKgoPT7+R9IRkYXt3BIkvZkVLzdKxOd+m6j9OZ0h844wxHgWgIAgECj580mso66//z+B1Xb7sNFUnq69NVyf9cIAAAEA+HNJoaN8P7S6NIvSrcxDNcGAEBYILzZRF6e99s8N4574QAACDeENwAAABshvIWY3T/TWwYAACpGeAsxL71SPrwdO1b9QLf7Z6P8fKOCAqOtP1hDjVTGX9NsAQAA/yK8hZjvNkovv+r+dME9/U+fpAoKKl/nnTnSvQONHn/KaPAQo/lVfFpVkhhgBACA0EF4C0GLPnb/OfvY6bf5+JPTr7N3r7RuvfX9Bwur3rVGJxwAAKGD8BYmcnOrt33qBqPdu0tjmoPuNgAAQhIzLIQZX+5V25lu9Mhwa8NVy0ltAACEMnreItTevVJ2thXYduwIcmUAAECVEd5srLpPhM5+l7vZAACwG8JbBCssDHYNAACAtwhvIepYjtGhw/SMAQAAdzywEKK63WYFt88WB7kiAAAgpNDzFuIy9ga7BgAAIJQQ3sAUCgAA2AjhLUz4e1Bd5jYFACA0Ed5CXOahqq33xpvep62qBj465gAACB2EtxA3+omKQ9nCD92X7dgRmO4yOuEAAAgdhDcb27LV/ecT+f7bN3ObAgAQmghvIKgBAGAjhLcwkpIa7BoAAIBAI7yFkRkzTUCeEp3+Bne9AQAQKghvYSYQ4W3519L3W4r8v2MAAOA1wlsE8+Zet9xcet8AAAgFzG2Kcua8Z7RyFWENAIBQRHhDuUF4/zmN4AYAQKjisimqPAhvZqZR+q7K187LM9q3j/AHAECg0POGKrujtxXKPnxfOvtszzfM3XaHUUGB9O47UrOmDCAHAIC/0fMWZnb/HPhj7EyveFlBgfWVMecAAAgMwluYyc0Ndg0AAEAgheRl06lTp+q1115zK7vkkku0YMGCINUoTDncvlS+KldAAQAICSEZ3iTp4osv1syZM0t+jooK2aoCAADUmJBNRLVr11bjxo2DXQ38YswzXI8FACAU+BTeFi1apJSUFG3evFnbtm1TYWGhXnzxRfXq1avCbTZt2qSpU6dqw4YNKioqUlJSku677z5169bN4/o7duxQp06ddOaZZ6pdu3YaMWIEYc7PvLkSuifDGbB6AACAqvMpvE2ePFkZGRmKj49XkyZNlJGRUen6ycnJGjRokOrUqaPu3bsrNjZWS5cu1fDhw7V//34NHDjQbf02bdroxRdf1IUXXqgDBw5oypQpGjBggD788EPVqVPHlyqjElt/ZFw2AADswqfwNm7cODVv3lwJCQmaMWOGXn311QrXLSoq0pgxY+RwODR79my1bNlSkvTwww/rrrvu0oQJE9SlSxclJCSUbHPdddeVfN+iRQu1bt1a119/vZYtW6YuXbr4UmVUYs67wTlu+i6jBvWl+HiehgAAoKp8GiqkY8eObmGrMsnJydq9e7d69OhREtwkqV69ehoyZIgKCwu1cOHCSvfRqFEjJSQkaM+ePb5UFyFo3z6j/gOMbruDXj8AALwR8HHe1q1bJ0nq1KlTuWWusvXr11e6j+zsbO3du7fKgRFV42n4j4MH/RSmTrObrT/65zAAAESagD9tmp6eLklq3rx5uWWNGzdWTEyMdu3a5Vb+97//XTfeeKPOPfdc7d+/XxMnTtSvfvUrt8upCIxJU+kJAwAglAU8vOX+MuR/vXr1PC6Pi4tTTk6OW9n+/fs1fPhwHT16VGeddZbat2+vl156SWeeeWaFx2nQoIFq1QpkR+LhAO47OOrWrav4+FiVbduKlaffrl69OMXHV/TgiLWvmJgYxcfXrXAfcbEFkqzfjfj4+CrWOLjsUs9AiOS2S5Hd/khuuxTZ7aftoSskx3mbOHGi19tkZ2cHoCbhLT8/X1lZJ73eLicnV1lZlT9kkJeXp6ysExUuzz1e2sOXlZXldR1qWnx8vC3qGQiR3HYpstsfyW2XIrv9tD3wba9OQAz4PW9xcXGSVK53zSU3N7fCXjkAAAC4C3h4S0xMlKRy97VJUmZmpvLy8jzeDwcAAIDyAh7e2rdvL0latWpVuWWuMtc6qFmGZxMAALCdgIe3Dh06qFmzZlq8eLG2bt1aUp6Tk6Np06YpOjpaPXv2DHQ14MHCRdXfx759RitXGRmSIAAANcKnBxbmz5+vlJQUSdK2bdtKylxjurVt21a9e/e2DhAVpXHjxmnQoEHq16+f2/RYGRkZGjVqlJo2beqPtsBLBQXV30fve6zQ9vxYh667tvr7AwAAlfMpvKWkpJSbFSE1NVWpqaklP7vCmyRdffXVmjNnjqZMmaIlS5aUTEw/cuTICiemh71s2mR03bVMcwUAQKD5FN7Gjx+v8ePHe7VNmzZtNHPmTF8OBxs63UVUrrICAOCbgN/zBgAAAP8hvMFrJ09KBQV0nQEAEAyEtwg35TWn19uMesKox+1GxcUEOAAAahrhLcLNe9+37U7kS+vW+7cuAADg9Ahv8Nn/jS7T8+blg6YOHkwFAMAnhDcAAAAbIbwBAADYCOENfvFtirRrFw8wAAAQaIQ3+MVPP0n9BlQvvG3/yeif053KyfFfCMzMNBr9pFPrvyVYAgDCg08zLACnc7oZFDwtv+8Bq/DoUaMnRvnniYaXXjVakyyt+sZo1XKekgAA2B89bwiIyVONVq7yrbcrLc1/9Th40H/7AgAgFBDeEBDFxdLjT3GpEgAAfyO8AQAA2AjhDQAAwEYIbwi6Q4e4vAoAQFUR3hB0Pe8y2rePAAcAQFUQ3hAS1qeUfu/PeU9PN2QJAAB2Q3gDAACwEcIb4IExRtvSjAoK6LoDAIQWwhuCwp+XRgPh08+kgYON/vIo4Q0AEFoIbwiKUL8X7eNPrApu2RrkigAAcArCGyJa8lqjH7eFeJIEAKAMwhsCKjvb6OPFRsePh15A+nmP0chRRg/8KfTqBgBARaKCXQGEt9FPGv1vs7RmrfTC2Jq/0a2yy7N7MmquHgAA+As9bwio/222vq5YGdx6AAAQLghvCBu7dxstXGRUVFTa3cYFUQBAuOGyKcJG33utqPbRx9JvfmP02F9DfDwSAAB8QHhD6Klm5krbbv27vI3v+wj1oUwAAJGLy6YICXl5/t9n9jH/7xMAgGAjvKHGGGOUtt39njSX1/5ZpiwEer1CfQYIAEDkIryhxsydL90/yOjZ50ylAW1bmn+O53AoJIIgAAD+RHhDjXn3PStJLV8R5Ir84nhusGsAAID3CG8ISc+/6Azo/rOzjZ4dS7ccAMB+CG8ISZ9+Lh067B6uiouNjBePgVZ239rGTb7WDACA4CK8ITiq8ECAKdP5VlBg1LuP0f+Nrl5v2d59Ru/ONco7Ua3dAAAQNIzzBltI3SAdzLT+Vcf9g4yOH5caNqh8PcZ5AwCEKnreUGMOH6nZ43nq3Dt+3Pp6NLtGqwIAgN8Q3hCyiv3wzAI9aACAcEN4Q8h6aGj1kteOnb5vzyC9AIBQRXhDUBzNOv061b2/bfGS6m0PAEAoIrwhKCZO4XomAAC+ILwhpL31L6P8fKNPPvUt7HHPGwAg3DBUCELam28ZvflWsGsBAEDoILzB9l6Z4NSPfprM3oUeOwBAqCK8wfY+/CjYNQAAoOZwzxvCGh1oAIBwQ3gDAACwEcIb4AGD9AIAQhXhDbaVkmo0cLAf5tACAMBGeGABtjVsBHe0AQAiDz1vAAAANkJ4Q8QwXgzexjhvAIBQRXhDWCsbwgYPIZEBAOyP8IaI8cOPwa4BAADVR3gDAACwEcIbAACAjRDeAA8YpBcAEKoIbwhvPKMAAAgzhDcAAAAbIbwhrPk6XhvjvAEAQhXhDQAAwEYIb7Adp5NuMQBA5CK8wXaeH1/18JaxN4AVAQAgCAhvsJ3Plwa7BgAABA/hDbbkzSTzAACEE8IbbOnLrwK7/2AP0puXZ/TFl0Z5eYRUAIC7qGBXAPDFv/4T3qHm2bFGq9dI1/1Oev45pnsAAJQK+Z63Z555Ri1atNA777wT7KoghOzbH9j9B/uq7Oo11tevVwS3HgCA0BPS4W3ZsmX67rvv1KRJk2BXBQAAICSEbHg7dOiQnn32Wb300kuKjo4OdnUQYvzdM3bsWHhfhgUAhA+fwtuiRYv09NNPq1evXmrdurVatGihBQsWVLrNpk2bNHjwYLVr106XX3657r77bi1ZsqTC9R9//HH98Y9/VIsWLXypIuCVsS8Q3gAA9uDTAwuTJ09WRkaG4uPj1aRJE2VkZFS6fnJysgYNGqQ6deqoe/fuio2N1dKlSzV8+HDt379fAwcOdFv/nXfe0YkTJ8qVAy7+7nlbk+zbdu/NM/p8qdGkVx1q0IAHCwAAgedTz9u4ceP01VdfKTk5WX369Kl03aKiIo0ZM0YOh0OzZ8/W2LFjNXr0aC1atEiJiYmaMGGCW/j76aef9M9//lPjx49XrVohe1UXkCS99k+jtO3Sf2bTcwcAqBk+paOOHTsqISGhSusmJydr9+7d6tGjh1q2bFlSXq9ePQ0ZMkSFhYVauHBhSfnGjRt15MgRde7cWa1atVKrVq2UkZGh559/Xrfffrsv1UU4CrGsVFgY7BoAACJFwMd5W7dunSSpU6dO5Za5ytavX19SdvPNN6t169Zu6z3wwAPq1auXevXqFcCawk5OBjgslR2k94cfjDp0COzxAACoqoCHt/T0dElS8+bNyy1r3LixYmJitGvXrpKy+vXrq379+m7rRUdHq3Hjxh73AQRC2XvqBg0xWv6lU1G1g1cfAABcAh7ecnNzJVmXST2Ji4tTTk5OtY/ToEGDAN8jdziA+0YoiI+PL/k+KipbUlHJz3v2FOu3l8d72Mr6vfhgodSsWV0NfuBMP9Wm9PetbL2CJRTqEEyR3P5IbrsU2e2n7aHLFtNjffXV6SeyzM7OroGaIJxlZWWVfF9U5HRb9ulnJ5XYPLfS7SdNydNdvfIDWq9giI+PD3odgimS2x/JbZciu/20PfBtr05ADPjjnHFxcZJUYe9abm5uhb1yQE36aUfFT0GsSeaJBABAaAh4eEtMTJQkt/vaXDIzM5WXl8e9bAgJAwZWHN527CyuwZoAAFCxgIe39u3bS5JWrVpVbpmrzLUOEGyffGp06FDgxyHZvduosDDExjsBANhCwMNbhw4d1KxZMy1evFhbt24tKc/JydG0adMUHR2tnj17BroaQJW8+Hej+wcHNlR9ucyo771GI0cR3gAA3vPpgYX58+crJSVFkrRt27aSMteYbm3btlXv3r2tA0RFady4cRo0aJD69evnNj1WRkaGRo0apaZNm/qjLYBfZGVJCecFbv8LP7RCW0pq4I4BAAhfPoW3lJQUt1kRJCk1NVWpqaV/jVzhTZKuvvpqzZkzR1OmTNGSJUtUVFSkpKQkjRw5Ut26dfOx6gAAAJHHp/A2fvx4jR8/3qtt2rRpo5kzZ/pyOAAAAPyCmd+BIPrhB+57AwB4h/AGBNH+A8GuAQDAbghvgB8VFRk98Cennn/RefqVAQDwAeEN8KOUVOnHbdKnn1e8juFKKQCgGghvgB85CWYAgAAjvAE1zOEIdg0AAHZGeAM82Py9jxtWoeeNy6YAgOogvAEhxpDuAACVILwBIeTECaM+/QhvAICKEd6AKsrPN0rdYFRUVHG4+jb19MGrsnvelq+QMvb6UjsAQKQgvAFV9NQzRo8MN3rz7YoD2tx5p99PpVdF6XQDAJyGT3ObApEoea319T/vSEeznBrwx+A9NpqVZdSwoeTg0VUAiDj0vAE++PgT6fGnqt5NtjPdKOuof7rVFi8xuu0Oo9dn0E0HAJGI8Ab4KG175cuLi42cTqM9e4z+eJ/RbT39E7YmT7X2M+ddv+wOAGAzXDYFAqT/fUZZWVJubrBrAgAIJ4Q3IEB+/jnYNQAAhCMumwIAANgI4Q2oYZUNFfJjGg8hAAAqR3gDatim/1W87P0Pypcd9dNTqgCA8EB4A0Jcj55G/5lNgAMAWAhvgA1Mf4PwBgCwEN4AAABshPAGBNGnn1e/R+1YjtGOnfTMAUCkILwBQfTN6urv47aeRvfeb5S2nQAHAJGA8AbYXHGx9fXblODWAwBQM5hhAahBve9xlis7lmP00F+Mbr7REYQaAQDshp43oAbt21e+bO48o/R0aeas6l32rGzwXwBA+CC8AUFWWBTsGgAA7ITwBgAAYCOEN8BGsrO5NgoAkY7wBtjE+wuMut9udOKE5+UbvjPKyiLcAUC4I7wBNjFpSuXBbE2ydNsdRj/8SIADgHBGeAOCzc9Za9CDRtt/IsABQLgivAFh6H+bg10DAECgEN4AAABshPAGBNnhw75vayoYmZcBewEgfBHegCDLPub7ti9PIKUBQKQhvAFBlrzW920/+thzuTFWr1xBAeEOAMIN4Q0IUy++ZHRTF6MdO08f4E6cMNq7j6AHAHZAeAPC1JJPra+z3z19KOvdx+jue6oW9AAAwUV4A8JRmQz2+dLTr3402/panUu4AICaQXgDwtDadfSgAUC4IrwBYWh1crBrAAAIFMIbEAHGvuDUDz/QGwcA4YDwBkSAz5dKg4aUD2+HDhmtWEmoAwA7iQp2BQAET5/+Rvn5wa4FAMAb9LwBEYzgBgD2Q3gDIojTaZSfz2VSALAzwhsQQR582OjmW42OHiXAAYBdEd6ACLJ1q/U1eV1w6wEA8B3hDQAAwEYIb0AEclRQbriaCgAhj/AGAABgI4Q3IAJ9/l+jHj2dwa4GAMAHDNILRKB16z2Xvz7dKDpauvuuii6sAgCCjZ43AG6mvMaNbwAQyghvAAAANkJ4AwAAsBHCGwAAgI0Q3gD4RW6uUXEx98sBQKAR3gCUU1Bg9MQYpyZPLT+cSGGh0auTnPpmdWlQ23/A6NYeRkOGEt4AINAIbwDKuamL0YqV0vwPyi9b8KG08ENp1BOlQe2rZdZX19ypoaCoiCAJIDwR3gB45WCmd6Fo8RKjhx9xKju75sLUmrVGN9xitPgTAhyA8EN4A+AdL/PQ+JeMNm6S3vpXzQWpx580MkYa/zLhDUD4IbwBqDZHFSZkOH488PVwqUp9AMCuCG8AKnX/oGy9N88oL8/qxTJlOrNGPWE90FCVsFSTfWBkNwDhjLlNAVRq3foirVtvzYc64WWH5r1fuuyb1dLXK4w2fHf6aGYClN4KC42io93jmoOPpQDCWEi+xb399tvq3r27fvvb36pdu3a69957tXHjxmBXC4ho69ZLd9xVfuiQJ582+mZ16c+z3jbatctDUgtAePvnNKduuMUoLc1957XoegMQxkIyvJ133nkaPXq0Fi1apPfee0/NmzfXAw88oKysrGBXDYhomYdOv86st436DSif1ALR8zbnPevrjDdP2TnhDUAYC8nw1rlzZ1177bU6//zz9etf/1qjRo1STk6O0tLSgl01AD5yBvCmt1ODIQ8sAAhnPt3ztmjRIqWkpGjz5s3atm2bCgsL9eKLL6pXr14VbrNp0yZNnTpVGzZsUFFRkZKSknTfffepW7dulR7r5MmTmjt3rho0aKCkpCRfqgsgBATqnjdJSl4r7d1ndN65VmojvAEIZz6Ft8mTJysjI0Px8fFq0qSJMjIyKl0/OTlZgwYNUp06ddS9e3fFxsZq6dKlGj58uPbv36+BAweW2+bbb7/V4MGDlZ+fr7PPPluzZs1Sw4YNfakugBBxLMdo7PNGt3Z26KYb/Zuwxjxj9OYMwhuA8OfTZdNx48bpq6++UnJysvr06VPpukVFRRozZowcDodmz56tsWPHltzPlpiYqEF+DWIAACAASURBVAkTJngMf61bt9aHH36o9957T7/73e/06KOP6siRI75UF0AIMEaa9ZbRmmTpmedKu+EOHPTPhPZ795V+X6vMO1vqBqOhw5zamc6AvQDCg0/hrWPHjkpISKjSusnJydq9e7d69Oihli1blpTXq1dPQ4YMUWFhoRYuXFhuu7p166p58+a67LLL9Pzzz6tWrVpasGCBL9UFEAK+WiYdOeWZo08+NbrzbqPrbjLqdL1TP+2oOGD9tMPo+y0VL3eWeRC2bMfbI8ONvttozboAAOEg4A8srFu3TpLUqVOncstcZevXrz/tfowxOnnypH8rByCgvlzmHpjy8kq/P37c6MW/uy8f/lejnByr7FiOe4/cgIFGDz5klHXUcwhzFpd+7+my6eHDVa83k9oDCGUBD2/p6emSpObNm5db1rhxY8XExGjXrl1u5S+//LJSUlKUkZGhLVu26KmnntL+/fvVpUuXQFcXgB898zf3EJS8tvT7114vH5COZEldbzN6/Emnut1m9cid6lAFw5W4Pc3qIbxVNY59vtTo+puNli0nwAEITQGfYSE3N1eSdZnUk7i4OOXk5LiVHTx4UCNGjNDhw4fVsGFDtW7dWrNnz9ZFF11U4XEaNGigWrUCmUW9+NgOQPHx8ars/83HiyveduU3pd8fy6mv5ufXLtnX4cOxurL9GVqTfFJS6XuHMa5jSrVrH9Gpcc2h0uWVGfuCdZwxzxp9v7FRmbZEpkhuuxTZ7aftoSskp8d6+eWXvd4mOzs7ADUB4Ct/Daq9cWO26tcr7UobOSpXV1+Vp0EPus/2UFxc9pgeBgn2oU5ZWVmKj49XVlaWJk1xKipKGvpQSA6PGRCutkeqSG4/bQ9826sTEAP+LhQXFydJ5XrXXHJzcyvslQNgX52uLz+Vli8ef6r8fW7vzi0fztweWPBw2TQ/3/P+TSUD0LnuuTt82Oj9BdJ786QTJ7icCiC4Ah7eEhMTJancfW2SlJmZqby8PI/3wwGAy2093QPTPzzcL+fKYD/vMRXeF3eqn3YY9bzL6KPFngNZ525GaWlFKirzMITTP5kUAHwW8PDWvn17SdKqVavKLXOVudYBgOrYs8fonv4V94w9+bRTHyywlv+4zWjAQKPDh6WXXjFK3WBUWOi+bUGB9PKEPE+7AoCgCXh469Chg5o1a6bFixdr69atJeU5OTmaNm2aoqOj1bNnz0BXA0AE6FNJcJOkr1dIE6dY6zzwJ/d1HxluNGFy5Zdjq6qgwOh/m/0z+DAAnMqnBxbmz5+vlJQUSdK2bdtKylxjurVt21a9e/e2DhAVpXHjxmnQoEHq16+f2/RYGRkZGjVqlJo2beqPtgBAlVR0P56nJ2DXJBdq46bSm+hOvUUuO9to7AtG993rUOtLrPXGPGO0Oln60yCH7u3vt2oDgCQfw1tKSkq5WRFSU1OVmppa8rMrvEnS1VdfrTlz5mjKlClasmRJycT0I0eOPO3E9Hb25wcden06n7wBu3tuXOn/49QNUvt2Rrd0NaoTbZWdLJSS1xot/0KKinJodbJV/v4HRvf2r3yi1cOHjT5bKnXrKsU3ZFJWAKfnMJU9amUjgX6sd9gIp1J+yab//dShW7pW/LL9+iLp6accuvACh09P3I1+zKHUVKOlX/haWwDB0PN2qX9fh+76g/X+EBsrjX3WoXPPlXbulNK2G/W526FN/5P+b7TRC2MdevvfRtvSpMsvk16bHFrDkETycBFSZLeftof2UCGEtyoa9YRT36y2vl+1vJbSdxkNHWZ0xhnSgQNSYqL08osOnXuu+yfn7zYaffGV0YeLqn6sVcutN/BpM5xa8qn0u2utT/u7fy5d5/bbpNt/79DAwdbpmzvboZwcadAQ6+ek30jb0srv+0+DHOpzt3Rj57A47UBYee5Zh85qJGVmSuefL114gdWTJ0nrvzU65xypWVOHjh0zql+/6r10xcVGtWpJjl/GUNn8vdErE4weGerQFb+teD+R/Adciuz203bCW40I9Au9Z4/R0Eelu++S+vbx/tJGdrbRtDeMund1qFlTqfvt1svet4+1vx6/DIUw/nmHOl1T9v4aU/KGe+CgUcyZ0g8/Wp/So6MdSl5rFBenknttnE6jWrUccjqNTp6U3pljFBXl0MxZ1v6X/deh6GiH3ptn9No/jca/4NBFF0iHDktNmkh33m2td+010l+GOvSXR416/t6h+wc01Oo1R3UkSxr3glGPbtKo/3Po8BHp7LMcyssz6tzN2nb4MIfaXCrdP6j8r9ZjIx166RWr/PrrpOVfly476yzv5p8EIkGrltaHsQ8/sn7+w93S3HnS4Aes++mWfS0t+sio3z0OtbxYeulVoy6dHbqmg/XBMj5e6j/AqHlz6ZW/Wx8Mb+riVEGBtb+3Zjq0YqVR3z4OnXmm+3tbJP8BlyK7/bSd8FYjauKFbtiwoY4ePeqXfeXkGH23Ubr6KiuEBZIxRiP+zyg2Rhr3XOWXZf42zqnvvpNm/9uhmJjSepX9ZS4oMDrjjPJ1PnDQaPP30vW/k2rXdui7jUZNmkir10gpqUaPj3KoTrR0863Wr9yXnzv02VIpNdVozJMOFRZKi5dYwfHbFOnMM6VLW0tHjpT2KL4xzaG4OFU4HETP26UPF0m1a1l/lNaul778yuiHH63lHTtIzc+X3p3rvl3tWlKxF1e4B97n0Ky3w+K/DsLIXXdK73/gXtY0QdqTYX2/anktHcsx6nZb+d/dPneXnz0ikv+AS5HdftpOeKsRNfVCR8Ivs6v3rix/tv1YjlEthxQXV/XQ+sMPVhBs1MjaJivLaPfP0iWtrHCXtl06flzqfIvnfc6Y6dT/NksTXrZ6HqfPdOo/71jLhvzJof59HTLG6NobrP8OgwY61CJJiomRLmvjUHx8vLb+cKTkXibXpe30XUb9B5T+F/r4Q4dWrJQKC6VPlhilbbfKH7jfoe5drfucftpRWq8vPnMoKsoaT+zBh4zSfxnL+sr20rr15dvR+WZxLyR81qSJdPBgxcs/XuhQQYH0l0eNBt7nUN97GikrK0sFBUZ16pRedo0UkfKe7wltJ7zVCMJbYIVj240x2r9fOuec0j9Ka9cZrVxl9JeHHW69i672L/nUKDZOuu7a0mVpaUaPPW40+AGHunV1/+M2731raqcHB5X2aOTlGZ15puc/hHl5RhkZUkys9Ie+1n/NJ0c79N58o2fHOHRBorXN9p+MMjOl3OPS38Z6/i/csKF0akfx+3Md+sfrVk/kgP4O/epXVvnMWUbfb6niC1fBvhF+Bg2sq/POLdBzz1u/Y88/59BlbaSGEfJUbDi+71UVbSe81QjCW2BFctul4LT/lYlO1YuTHhxc+aXu3FyjXbuty8HGSHXqqKSXpLjYKDu7tMeyIpmZRnf09vxW0L5dlJzOIr083qGiIpVcTi8qMvrXf6xtvlkj/TLkY5U8NtKhkyelSVPC4u0n4tx+m3ThBQ7d2cuh1WuMtv5g9dSFW89cJL/v0XbCW40gvAVWJLddioz2r11nFBsrnTghbdxkdFcvh+rWlc49t5FXbd+3z6hxY+mr5db4aI3PljLLzDU6eYJDba+w/sgbY2SMNVfp3PnSTTdIXy4rXffUnxFaPlvs0K09rD8hL45z6NpOhLdwQdsJbzWC8BZYkdx2KbLb74+2G2OUlyft2y/9+iLPf+BdT1bPm2/UoIHUpbO13rLlRmOeLX2benCwQ9PfKP+2NelVhx79q1V+3e+kp5906KYuYfH2ZgvXXycNGezQG7OM/tDboVYtHSXh/NR7aO2C//e0PdDH8RXhzQv8Mkdm26XIbn8otL3sUDT/muXQ8ePS408aFZyU8vOlBvWlTz4qf3nZl0Gy4R8rlzk0YKDRjp1W4O56q1TLIdWta1163/qDNW5dvbjSsexCTSj87gcLbQ/t8ObT9FgAUJNiYhx6cZx1+fWiC60/9J985FBRkdGKVVKbS6u2n45XSz9nSGfUkWZOd+j6m8Pis2tIWv61tGOn9f30N4xmvmkNxxMVJT38Z2nyVOu1b9hAmveu3IYmAlA5et68wCeRyGy7FNntt3Pbn3rGqeVfS48MdSjhPOmK30pnnCE5HNYDHRu+M/rLo57fApcucSgnt3TgagROt67SE6NCa2owyd6/+9VF20O75y30/rcAgJ88O8ahd/7lUO87pWs6WjMI1KpV+lTkby93aNXyWhr9mPVzYnPpjp7StH9Yg1T/qgm9QTVhyafS/v1Gm783evpvTu0/YAXmYznWUDuS9UT05KlO7f6ZMA1w2RRA2IqKciix+enX697VepDigkR5nD3EpaLBk1F9e/dJjwy3gtlXy4zOamR0+Ii17I/9jObOk04WSp8tNfr0Y/dzdPKkUZ06BG1EDsIbgIjncDh0cQvPy2b/u74+/eyY7h9gDdy8eo01KLPLxFcc+nmPNGESPULV4QpuLq7gJkn/mV36fU6O+3abvzca8rDReecZXXapNYZgoKccBIKN8AYAlbj8smg1P7/0DpOOHRz69GPpyaeNbr7JofbtHGrfTurV06Gp/3Bq7vyK93VmXelEfg1UOszN/8DI6ZTq15fe/8AKfXv3Wv8aNDAa+hDhDeGNBxa8wA2ckdl2KbLbH8ltl7xvf0XDk3z4gUPfrJZefjUs3nJDmmvu4eqK5N992s4DCwAQMdq1tb4+/ZRDl7Syvr/lZunssxzq3jV49YokG74z2pZm9OZbTuXmGqVtN5ox06mCgvLB+dBho7Q0qzw/3+jIkdJ1nE6joiLCNkIPPW9e4JNIZLZdiuz2R3LbJe/bX1holLFXSmzuUHa20eo11owPp45jdurMEagZf+wnxcU51PoSadP/VDJbx6svOfTM34xyj0u/vkh66UWHHnu8lrKyivXmDIcKC6Vzz42cy7GR/P/eDj1vhDcv8MscmW2XIrv9kdx2KbDt73G7U0ezre97/l7663CHho0wSt0QkMOhmt5/z6Gdu6TFnxhdeIH09Qop4Tzp1lsdim8oxcVKtWtb627YKPXoJn2/RTr7LOm886xy1zA1ZRUXG9WuXVqek2MUF+d53ZoSyf/v7RDeeGABAELAyBHWXSy33yalbgiLz9Rh564+pefl6xXW1x07pZXfeD5fr0zwVFrRufVU7v3vQZPG0sHM06/X6RopJUV6YZxDZ54pHTsmffGl0dIvrOWvv3ZSq75x6qorHUpJNTpxQtqWJrVr61BmplGLFlZgPXTYGh8xOloqKLAuWUdFOfTby6VataTjx6X9+6WEBCvYnjghFRdLDRtayw8clJxOa9aT+HhrAO3iYmub+g2kkwVSvXpSnTrWttHR0sGD1v6OHpXi4koD8/79UpMmUlGRdOiQdNZZUvJaa37dli0d+lUTz6/F0WwpO9vaV6N4qXnzYtWv5/VLX6MIbwAQQlx/iABfVCW4SdKqb6yvw0d6Doh/HmqNyTL7Xffl32385eePKwuWofjhw5s6HdXrrzl0aevQvUxOeAOAILnwQpW7RHpNR6n5+VYvwAtjHZo4xajdFQ61bi1deIG04EPprEbSs2ONiourdpyoKKs3AvDk3HOkffurt4+E86SMvaU/n9VIOqOuNXyLS+OzrfmJJemcX0lOIxXkSwUnpfxfhtA591xp3z7re9fvbaN4a18yVj3rRFvbNmhgjft31lnW19xca7vYWGt/jeJLjydJMTFWj19ZB/Zbc+6WbUeTJlE695wq/ucKEu558wL3AERm26XIbn8kt10KbPsPHTJ6822jXj0d+s2vvfuUf+iQ0dYfrD9ul7a2ymJirJkKmjW1ZoooKrIuYUmSMUY//2z9YSsuti5Z1aljzfVau7Z16So/3/qXm2v94TyYGauLLjyumBhr0Nzdu6VLWlnrFhZal7jyC6xLXtYxrD+i9epZf1yNsfZ/xi/HKSiwjpl3wiozxqpLUZF1Sezss63LcA3qW8dw/fGuU8c6Xp06Vrk1N631r1Ytaz+S9bOLMe4/l1XVe8ki+Xeftof2PW+ENy/wyxyZbZciu/2R3HYpstsfyW2XIrv9tD20wxvjvAEAANgI4Q0AAMBGCG8AAAA2QngDAACwEcIbAACAjRDeAAAAbITwBgAAYCOENwAAABshvAEAANgI4Q0AAMBGCG8AAAA2QngDAACwEcIbAACAjRDeAAAAbMRhjDHBrgQAAACqhp43AAAAGyG8AQAA2AjhDQAAwEYIbwAAADZCeAMAALCRqGBXINRt2rRJU6dO1YYNG1RUVKSkpCTdd9996tatW7Cr5tGNN96ojIwMj8uuvPJK/ec//3ErO3nypGbMmKGPPvpI+/btU4MGDXTDDTfo0Ucf1VlnneVxPx999JH+/e9/a/v27YqOjtYVV1yhRx55RJdcconH9f39Gi5atEgpKSnavHmztm3bpsLCQr344ovq1auXx/Vzc3M1depULV26VJmZmWrSpIm6dOmioUOHKjY2ttz6TqdTs2fP1rx587Rr1y7FxMSoY8eOGj58uJo1a+bxGCtXrtT06dP1/fffy+Fw6JJLLtFDDz2kDh06eFx/586dmjRpkpKTk3XixAklJiaqT58+uueee+RwOPzW/qlTp+q1116rcF9ffvmlmjZtWuPt8facSNKBAwf06aefasWKFdqxY4cOHTqkBg0a6IorrtCgQYN02WWXVfs4oXzuvW1/OJ37goICTZgwQZs3b9auXbuUnZ2t+vXrq1mzZurdu7d+//vfKzo6ulrHCeVz7237w+ncV2TGjBl69dVXJUlz587V5ZdfXq1jhfL594ShQiqRnJysQYMGqU6dOurevbtiY2O1dOlSZWRkaNSoURo4cGCwq1jOjTfeqGPHjmnAgAHlliUkJLj9gXc6nRo8eLBWrVqlyy+/XO3bt9euXbv03//+V02bNtW8efPUqFEjt328/vrrmjRpkhISEtS5c2cdP35cn3zyiQoLC/X222+rbdu2busH4jV0BdT4+HjFxMQoIyOjwvCSl5envn37auvWrerUqZNatmyprVu3atWqVbr00ks1e/ZsnXHGGW7bPPXUU5o/f75+85vf6LrrrtPBgwf16aefKjY2VnPnzlViYqLb+osWLdJjjz2mRo0alQTSJUuWKCsrS5MmTdKtt97qtv727dvVp08f5efnq2vXrmrSpIm+/vprpaWlqX///hozZozf2u96E7/jjjuUkJBQbvmAAQNUv379Gm2PL+dEkl555RW98cYbOv/883XllVeqUaNG2rVrl7744gsZY/Tqq6+6fSAIt3PvbfvD6dwfOXJE119/vdq0aaPExEQ1atRI2dnZWrlypTIyMtSpUye98cYbqlWrls/HCeVz7237w+nce7Jt2zbdeeedioqKUl5eXrnwFm7n3yMDjwoLC83NN99sWrdubbZs2VJSfuzYMdO5c2dzySWXmD179gSxhp7dcMMN5oYbbqjSuu+//75JSkoyI0aMME6ns6R8zpw5JikpyYwZM8Zt/Z07d5pWrVqZzp07m2PHjpWUb9myxbRu3dp07drVFBcXl5QH6jX85ptvSrabPn26SUpKMh988IHHdSdPnmySkpLMyy+/7Fb+8ssvm6SkJDNt2jS38jVr1pikpCTTr18/U1BQUFK+fPlyk5SUZAYOHOi2/tGjR027du3MVVddZfbt21dSvm/fPnPVVVeZq666yuTk5Lht069fP5OUlGSWL19eUlZQUGD69u1rkpKSTGpqqt/aP2XKFJOUlGSSk5Mr3WdNtsfbc+Ly+eefm7Vr15YrX79+vbnkkktM+/bt3c5ZuJ17b9sfTue+uLjYrW0uhYWFpn///iYpKcksW7bM5+OE+rn3tv3hdO5PdfLkSXPHHXeY3r17m5EjR5qkpCSzYcOGah0r1M+/J9zzVoHk5GTt3r1bPXr0UMuWLUvK69WrpyFDhqiwsFALFy4MYg2rb/78+ZKkESNGuHXZ9unTR82aNdPHH3+s/Pz8kvIFCxaoqKhIf/7zn1WvXr2S8pYtW6pHjx766aeflJKSUlIeqNewY8eOHj9NnsoYo/nz5ysmJkYPPfSQ27KHHnpIMTExJa+Bi+vnYcOGqU6dOiXl1113na688kqtWrVKe/fuLSn/7LPPdOzYMfXv31/nnHNOSfk555yj/v37KysrS1988UVJ+c6dO7V+/XpdddVVuu6660rK69Spo2HDhkmS5s2b55f2+yLQ7fHlnLh07txZV155Zbnydu3a6aqrrlJ2drZ+/PFHn48T6ufem/b7IpTPfa1atdzOiUtUVJRuueUWSdKuXbt8Pk6on3tv2u+LUD73p5o2bZrS0tL0wgsvqHbt2uWWh+P594TwVoF169ZJkjp16lRumats/fr1NVqnqjp58qQWLFigadOm6Z133tHGjRvLrVNQUKCNGzfqggsuKBcEHA6HOnbsqLy8PG3evLmk3PWaXHPNNeX253pNXOuU/T5Yr2F6eroOHjyoK664QjExMW7LYmJidMUVV+jnn3/Wvn37SsrXrl1bsuxU1157rSTv21jV9du2bauYmJiAvCbr16/XjBkzNHPmTH3xxRc6fvy4x/UC3R5fzklVREVFuX2NtHN/avvLCudz73Q6tXLlSklSUlKSz8ex67n31P6ywu3cf//995o2bZqGDh2qX//61x7XiZTzzwMLFUhPT5ckNW/evNyyxo0bKyYmplqfdAIpMzNTjz/+uFvZpZdeqgkTJuj888+XJO3evVtOp7PcdXwXV3l6erratWtX8n1MTIwaN25cbn3X61T2NQn2a+jad2VtXLVqldLT03XuuecqLy9PmZmZSkpK8viJzts2ert+7dq11bRpU23fvl1FRUUe/xD7aurUqW4/169fX08++aR69uzpVh7o9nh7Tqpi7969Wr16tRo3blzyByySzr2n9pcVTuf+5MmTmj59uowxOnr0qNasWaMdO3aoV69eJTeJh/O5r0r7ywq3cz9q1ChdfPHFGjRokMd1ytYzHM9/WYS3CuTm5kqS2+XBsuLi4pSTk1OTVaqSXr16qW3btkpKSlJMTIzS09P11ltvadGiRbrvvvv00UcfudU9Li7O435c5a7XwfX9qQ8wnLp+2dck2K+ht22s6vpVbaMvr0lsbKycTqeOHz+uBg0aeFzHGxdffLFeeOEFXXnllWrSpIkyMzO1fPlyTZkyRaNHj1a9evV000031Vh7fPm9q0xhYaEee+wxnTx5UiNHjix5842Uc19R+6XwPPeFhYVuT1E6HA4NHDhQf/3rX0vKwvncV6X9Unie+8mTJys9PV0LFizwGLJcwvn8l0V4CzNDhw51+7lly5Z66aWXJFlPx8yfP1/3339/MKqGIHDdD+PStGlT9e/fXxdddJHuv/9+TZo0ye1N3E6cTqdGjx6t9evX6+677y7XmxDuTtf+cDz3sbGx+vHHH+V0OnXw4EF99dVXmjhxor777ju98cYbFf4BDhdVbX+4nfsNGzZo1qxZGjp0qMfe5UjEPW8V8JSey8rNza0wRYeiP/zhD5Kk1NRUSaWfACr6pOMqL/tmWFlPmadPFsF+Db1tY1XXr2obfXlNjh8/LofD4fWYR97q0KGDzj//fG3bts2tvYFujy+/d544nU498cQTWrx4sX7/+9/rb3/7m9vycD/3p2t/Zex+7iXrBv5zzjlHffv21XPPPafU1FS9/vrrPh3Hbudeqrz9lbHjuS8qKtLo0aPVokUL/elPf6q8gT4cy47nXyK8Vch1vdzTPVmZmZnKy8vzeP06VMXHx0uyxr+RpGbNmqlWrVol1+JP5Sove99AYmJiyf0Bp3K9TmVfk2C/hq59V7WNrvv59uzZo+Li4nLre9tGb9cvLi7Wnj171LRpU7/e71YR1+/EiRMnqlQ/f7TH23PiidPp1OOPP66FCxeqR48eGj9+fMn4Vi7hfO6r0v7Tseu59+TUG8TD+dx74ukG+crY7dzn5eUpPT1dW7duVevWrdWiRYuSf67RCv7whz+oRYsW+uKLLyLm/BPeKtC+fXtJ0qpVq8otc5W51rGDTZs2SVLJk6V169ZVmzZttHPnznIzMhhjtHr1asXExKh169Yl5a72fvPNN+X273pNyg5lEOzXMDExUU2aNFFqampJaHXJy8tTamqqmjZt6naD7JVXXlmy7FSup7rK1rkqbazqa5KSkqK8vLwa+b3Ky8tTWlqaYmJiSt7MT1c/f7THl3NSliu4fPjhh+rWrZteeuklj/e/hOu5r2r7K2PXc1+RgwcPSip90jZcz31FTm1/Zex47uvUqaO77rrL4z9XKLrxxht11113KSEhIXLOv1ejwkWQwsJCc9NNN1U6wOzPP/8cxBqWt337dpOXl+ex/JprrjFJSUlm3bp1JeXeDtK7Y8cOrwfpDfRrGAqD9LZt29avgzWmpKT4pf05OTlmx44d5cpPnDhhRowYYZKSkszo0aNrvD3VGah11KhRJikpyTzyyCOmsLCwglfFt+OE+rn3pv3hdu7T0tI8vrfl5eWZBx54wCQlJZnXX3/d5+OE+rn3pv3hdu4r4/r/EIxBeoP5vm+MMUyPVQm7TY81depUvfXWW2rfvr3OO+88nXnmmUpPT9eKFStUWFioBx98UCNGjChZ39P0WLt379bSpUuVkJCg+fPnh+T0WPPnzy8ZDHjbtm36/vvvdcUVV5R0U7dt21a9e/eWZH3Suueee/TDDz+oU6dOatWqlbZs2VIyTco777yjunXruu3/1GlSMjMztWTJEsXGxuq9997TBRdc4LZ+ZdOkTJw4UV27dnVbPy0tTffcc4/y8/PVrVs3NW7c2KtpUqra/j179ujmm2/WpZdeqosuukhnn322Dh8+rNWrV2v//v1KSkrSv//9b7dP4DXRHl/OiVQ65U9MTIzuvfdejz0NN998c8mA0OF27r1pfzie+7feektt27ZVQkKC4uLidODAAa1YIYpe7QAAAc5JREFUsUJHjx5Vu3bt9Oabb5ZsG47nvqrtD7dzX5nRo0dr4cKFHqfHCqfz7wnh7TQ2bdqkKVOmuE2qfv/994fkxPTr1q3TnDlztHXrVh06dEj5+fmKj49XmzZt1LdvX48DBLompl+0aJH27dunhg0b6vrrr9ejjz6qs88+2+NxPvroI/3rX/9ym5h+2LBhlU5M78/X0PUftiJ33HGHxo8fX/JzTk5OyQTFhw4dUuPGjXXrrbfq4Ycf9niDrNPp1DvvvONxgmLXOHmnWrFihaZPn64tW7ZIklq3bq0///nP6tixo8f1d+zYoUmTJmnt2rXKy8srmaC4b9++p52guKrtz83N1YQJE7Rp0yZlZGTo2LFjOuOMM3TRRRepS5cu6t+/f4VvloFuj7fnpCrtllRujtdwOvfetD/czv3//vc/zZs3Txs2bNCBAweUl5enuLg4tWjRQt27dy+Z57I6xwnlc+9N+8Pt3FemovDmy7FC+fx7QngDAACwER5YAAAAsBHCGwAAgI0Q3gAAAGyE8AYAAGAjhDcAAAAbIbwBAADYCOENAADARghvAAAANkJ4AwAAsBHCGwAAgI0Q3gAAAGyE8AYAAGAjhDcAAAAb+X9JfUcA+Uge7AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(mean_field)\n",
    "plt.yscale('log')"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}