moment_of_inertia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tensor Field Networks\n",
    "\n",
    "Implementation of moment of inertia demonstration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/sebastian/miniconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "from math import pi, sqrt\n",
    "import numpy as np\n",
    "import torch\n",
    "import random\n",
    "import tensorfieldnetworks.layers as layers\n",
    "import tensorfieldnetworks.utils as utils\n",
    "from tensorfieldnetworks.utils import EPSILON, FLOAT_TYPE\n",
    "xrange = range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_points = 15\n",
    "\n",
    "# radial basis functions\n",
    "rbf_low = 0.\n",
    "rbf_high = 2.\n",
    "rbf_count = 30\n",
    "rbf_spacing = (rbf_high - rbf_low) / rbf_count\n",
    "centers = torch.Tensor(np.linspace(rbf_low, rbf_high, rbf_count))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrix_from_0_2(input_0, input_2):\n",
    "    \"\"\"\n",
    "    inputs:\n",
    "    -shape [N] tensor, L=0 feature at each point\n",
    "    -shape [N, 5] tensor, L=2 feature vector at each point\n",
    "    \n",
    "    returns:\n",
    "    -shape [N, 3, 3] tensor, matrix at each point\n",
    "    \"\"\"\n",
    "    d_xy = input_2[:,0]\n",
    "    d_yz = input_2[:,1]\n",
    "    d_z2 = input_2[:,2]   # -x^2 - y^2 + 2*z^2\n",
    "    d_zx = input_2[:,3]\n",
    "    d_x2y2 = input_2[:,4] # x^2 - y^2\n",
    "    Mxx = - d_z2 / sqrt(3.) + d_x2y2 + input_0\n",
    "    Myy = - d_z2 / sqrt(3.) - d_x2y2 + input_0\n",
    "    Mzz = 2*d_z2 / sqrt(3.) + input_0\n",
    "    \n",
    "    # rows : [N, 3]\n",
    "    row0 = torch.stack([Mxx, d_xy, d_zx], dim=-1)\n",
    "    row1 = torch.stack([d_xy, Myy, d_yz], dim=-1)\n",
    "    row2 = torch.stack([d_zx, d_yz, Mzz], dim=-1)\n",
    "    \n",
    "    # [N, 3, 3]\n",
    "    return torch.stack([row0, row1, row2], dim=-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_inputs(r):\n",
    "    \n",
    "    # rij : [N, N, 3]\n",
    "    rij = utils.difference_matrix(r)\n",
    "\n",
    "    # dij : [N, N]\n",
    "    dij = utils.distance_matrix(r)\n",
    "\n",
    "    # rbf : [N, N, rbf_count]\n",
    "    gamma = 1. / rbf_spacing\n",
    "    rbf = torch.exp(-gamma * (dij.unsqueeze(-1) - centers)**2)\n",
    "    \n",
    "    return rij, dij, rbf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MOIModel(torch.nn.Module):\n",
    "    \n",
    "    def __init__(self, rbf_dimension):\n",
    "        \n",
    "        super().__init__()\n",
    "        self.layer_0 = layers.Filter(0,0, rbf_dimension, nonlin=utils.ssp)\n",
    "        self.layer_2 = layers.Filter(2,2, rbf_dimension, nonlin=utils.ssp)\n",
    "    \n",
    "    def forward(self, masses, rbf, rij):\n",
    "        out_0 = self.layer_0(masses, rbf, rij).squeeze()\n",
    "        out_2 = self.layer_2(masses, rbf, rij).squeeze()\n",
    "        return matrix_from_0_2(out_0, out_2)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_points_and_masses(num_points=num_points, max_coord=rbf_high/4., min_mass=0.5, max_mass=2.0):\n",
    "    \"\"\"\n",
    "    returns:\n",
    "    -shape [num_points, 3] numpy array of points\n",
    "    -shape [num_points] numpy array of masses\n",
    "    \"\"\"\n",
    "    output_points = []\n",
    "    for point in xrange(num_points):\n",
    "        output_points.append(np.array([random.uniform(-max_coord, max_coord) for _ in range(3)]))\n",
    "    \n",
    "    points = np.array(output_points)\n",
    "    masses = np.random.rand(len(output_points)) * (max_mass - min_mass) + min_mass\n",
    "    return points, masses\n",
    "\n",
    "\n",
    "def moi_tensor(points, masses, index):\n",
    "    \"\"\"\n",
    "    inputs:\n",
    "    -shape [N] tensor, L=0 feature at each point\n",
    "    -shape [N, 5] tensor, L=2 feature vector at each point\n",
    "    \n",
    "    returns:\n",
    "    -shape [N, 3, 3] tensor, matrix at each point\n",
    "    \"\"\"\n",
    "    center = points[index]\n",
    "    moi = np.full((3, 3), 0.)\n",
    "    for point, mass in zip(points, masses):\n",
    "        x, y, z = point - center\n",
    "        Ixx = y**2 + z**2\n",
    "        Iyy = x**2 + z**2\n",
    "        Izz = x**2 + y**2\n",
    "        Ixy = -x*y\n",
    "        Iyz = -y*z\n",
    "        Ixz = -x*z\n",
    "        moi += np.array([[Ixx, Ixy, Ixz], \n",
    "                         [Ixy, Iyy, Iyz],\n",
    "                         [Ixz, Iyz, Izz]]) * mass\n",
    "    return moi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def R0_plot(model, x):\n",
    "    gamma = 1. / rbf_spacing\n",
    "    rbf_plot = np.exp(-gamma * (x - centers)**2)\n",
    "    r_layer = model.layer_0.children().__next__().children().__next__()\n",
    "    return r_layer(torch.Tensor(rbf_plot).unsqueeze(0))\n",
    "    \n",
    "def R2_plot(model, x):\n",
    "    gamma = 1. / rbf_spacing\n",
    "    rbf_plot = np.exp(-gamma * (x - centers)**2)\n",
    "    r_layer = model.layer_2.children().__next__().children().__next__()\n",
    "    return r_layer(torch.Tensor(rbf_plot).unsqueeze(0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = MOIModel(rbf_count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch.optim as optim\n",
    "\n",
    "criterion = torch.nn.MSELoss()\n",
    "optimizer = optim.Adam(model.parameters(), lr=0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Step 0: validation loss = 11.978\n",
      "Step 1000: validation loss = 0.006\n",
      "Step 2000: validation loss = 0.006\n",
      "Step 3000: validation loss = 0.000\n",
      "Step 4000: validation loss = 0.006\n",
      "Step 5000: validation loss = 0.004\n",
      "Step 6000: validation loss = 0.056\n",
      "Step 7000: validation loss = 0.000\n",
      "Step 8000: validation loss = 0.004\n",
      "Step 9000: validation loss = 0.001\n",
      "Step 10000: validation loss = 0.003\n"
     ]
    }
   ],
   "source": [
    "max_steps = 10001\n",
    "print_freq = 1000\n",
    "validation_size = 1000\n",
    "\n",
    "x_vals = [index * rbf_high / rbf_count for index in xrange(rbf_count)]\n",
    "\n",
    "\n",
    "\n",
    "# training\n",
    "for step in xrange(max_steps):\n",
    "    shape, rand_masses = random_points_and_masses()\n",
    "    \n",
    "    rand_masses = np.reshape(rand_masses, [-1, 1, 1])\n",
    "    \n",
    "    center_index = 0\n",
    "    # make center point index have zero mass\n",
    "    rand_masses[center_index] = 0.0\n",
    "\n",
    "    truth = moi_tensor(shape, rand_masses, center_index)\n",
    "    \n",
    "    rand_masses = torch.Tensor(rand_masses)\n",
    "\n",
    "    rij, dij, rbf = get_inputs(torch.Tensor(shape))\n",
    "\n",
    "    optimizer.zero_grad()\n",
    "    outputs = model(rand_masses.unsqueeze(0), rbf.unsqueeze(0), rij.unsqueeze(0))[ center_index]\n",
    "\n",
    "    loss = criterion(outputs, torch.Tensor(truth))\n",
    "\n",
    "    loss.backward()\n",
    "    optimizer.step()\n",
    "    \n",
    "    if step % print_freq == 0:\n",
    "        loss_sum = 0.\n",
    "        for example in xrange(validation_size):\n",
    "            shape, rand_masses = random_points_and_masses()\n",
    "            rand_masses = np.reshape(rand_masses, [-1, 1, 1])\n",
    "            rand_masses[center_index] = 0.0\n",
    "            truth = moi_tensor(shape, rand_masses, center_index)\n",
    "            rij, dij, rbf = get_inputs(torch.Tensor(shape))\n",
    "            rand_masses = torch.Tensor(rand_masses)\n",
    "           \n",
    "            outputs = model(rand_masses.unsqueeze(0), rbf.unsqueeze(0), rij.unsqueeze(0))[center_index]\n",
    "            loss = criterion(outputs, torch.Tensor(truth))\n",
    "            \n",
    "            loss_sum += loss.item()\n",
    "        print(\"Step %d: validation loss = %.3f\" % (step, loss_sum / validation_size))\n",
    "\n",
    "\n",
    "        \n",
    "    if step == max_steps - 1:\n",
    "        y_vals = []\n",
    "        for x_val in x_vals:\n",
    "            y_vals.append(R0_plot(model,  x_val))\n",
    "        F0_y_vals = y_vals\n",
    "        \n",
    "        y_vals = []\n",
    "        for x_val in x_vals:\n",
    "            y_vals.append(R2_plot(model,  x_val))\n",
    "        F2_y_vals = y_vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAD8CAYAAABKKbKtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXt8juUbwL/3Zpsx5yFnmzPbLIaUGMmhSCJ0dKgQkkqlEyVKB6Wk8CutVBKVRFGYVGQmyxySY+bUMJvDjB3u3x/X+27vu+O747vD/f18ns/7PPd9P89zvVq7dt3XSWmtMRgMBoOhoHFxtgAGg8FgKJ0YBWMwGAyGQsEoGIPBYDAUCkbBGAwGg6FQMArGYDAYDIWCUTAGg8FgKBScqmCUUn2UUvuUUgeUUlMymX9bKRVhOf5RSsXazCXbzK0sWskNBoPBkBPKWXkwSilX4B/gZuAYsA24S2u9J4v1jwDXaq1HWa4vaq29ikpeg8FgMOQOZ1owHYEDWutDWuurwJfAgGzW3wUsKRLJDAaDwZBvyjnx3fWAKJvrY0CnzBYqpRoBPsAGm+HySqlwIAmYpbVekcW9o4HRABUrVmzfsmXLAhC99HL27FkAatSo4WRJDIbSx759+wBo0aKFkyXJHdu3bz+jta6Z2/ucqWBUJmNZ7dcNA5ZrrZNtxhpqrU8opXyBDUqpSK31wQwP1HohsBAgKChIh4eH51fuUs0nn3wCwPDhw50sicFQ+ggODgZg48aNTpUjtyil/s3Lfc7cIjsGNLC5rg+cyGLtMNJtj2mtT1g+DwEbgWsLXkSDwWAw5BVnKphtQDOllI9Syh1RIhmiwZRSLYBqwBabsWpKKQ/LuTdwA5BpcIDBYDAYnIPTtsi01klKqQnAWsAVWKS13q2Umg6Ea62tyuYu4EttH+7WCliglEpBlOSsrKLPDLmjbdu2zhbBYCi1jBgxwtkiFClOC1N2BsYHYzAYDLlHKbVdax2U2/tMJr/Bjvj4eOLj450thsFQKjlz5gxnzpxxthi54ujRvN9rFIzBjmXLlrFs2TJni2EwlEoGDx7M4MGDnS2GwzzzDDRqlPf7nRmmbDAYDIZixuXLsHcv7NwJhw7l71lGwRgMBkMZ5sgR+PhjCA+H++6D/fshORlq1IAHH4R16yAmJm/PNgrGYDAYyiApKWKldOgASUky5u8P3bvLZ506oJQooMqV8/YOo2AMBoOhDJCSAomJcPo0REbCrl1w6RI0ayZbYgDu7tC7t/19lSrl/Z1GwRjsCArKdSSiwWBwkIcffrjI33noEMydC59/DsHB0Lo1lCsHzZtDQAC0bQv/+x8MHw633Vaw7zZ5MAaDwVAKOX9erJR334WlS2WsSRNYvhxatQIPD8efldc8GGPBGOyIi4sDoEqVKk6WxGAofURFSQH5Bg0a5LAy9yQlwbZtss0VGQn//gtaQ6dOolSSkyEuDho2zJ1yyQ9GwRjsWLFCuh6YasoGQ8Fz3333AQVbTfnyZRg7Fr77Di5ehCeegAYNZDvM3x+qV4erV6FFC7jlFlFARYVRMAaDwVDCSEqCAwdkC2zfPli/XqwTgGrVYPx4iQCz8vTTzpHTKBiDwWAoAezfD/Png5cXlC8PV65AxYoQGCj5Ki+9BLVrS0ixyqzblhMwCsZgMBiKKSkp4kt58014/30Za90aXnkF/PzAxwdcXKBjR7juOujZUyLEigvFSBSDwWAwJCRAdLRsf+3ZI34VW7/JgQPiX7GNw6lZE/r0KXJRc8QoGIMd1113nbNFMBhKLU888USm44mJsGABLF4sSuXRR2UbrHlzsVSaNYPNm0WRDBtWtI76/OBUBaOU6gO8gzQc+1BrPSvd/AjgDeC4Zeg9rfWHlrnhwPOW8Rla60+KROhSTosWLZwtgsFQaunfv7/dtdVSiYyEadPSHPW1asGYMfbhxL//Xry2vxzBaeIqpVyBecDNwDFgm1JqZSadKZdqrSeku7c6MA0IAjSw3XLvuSIQvVRj7VXh7e3tZEkMhtLHrl37WLfOherVm3HxoigYpcDXFwYOhJAQWRcdnTFXpaQpF3CuBdMROKC1PgSglPoSGAA40vq4N/Cz1jrGcu/PQB9gSSHJWmZYvXo1YPJgDIaC5OxZKdfy8surSEmpQVBQM8aNg1tvlax6Ly8ICpItsaFDReGUBpypYOoBUTbXx4BOmawbpJTqCvwDPKa1jsri3nqFJajBYDDkhuRkiI2F3bvFp3LqlESDpaSkAGEcOjSC++6zt0patZKjNOHMjpaZRWqnL4z2PdBYax0ArAOsfhZH7pWFSo1WSoUrpcJPnz6dZ2ENBoMhO65cgXnzoEsXKXX/7ruwYQO4uUmF4tmzwc1tNe7uexg+XCoZl3acacEcA2wL8tQHTtgu0Fqftbn8H/Cazb3B6e7dmNlLtNYLgYUgxS7zI7DBYDCk58wZsVQiI+HJJyXMGKTV8KBB9uHE7dpJdNhbbzlH1qLGmQpmG9BMKeWDRIkNA+62XaCUqqO1Pmm5vA2wdC1gLfCKUqqa5boX8Ezhi2wwGMo6R47Ap59C48bYOeobNIAePeCHH2TduXP2ygVEuZQlnKZgtNZJSqkJiLJwBRZprXcrpaYD4VrrlcBEpdRtQBIQA4yw3BujlHoZUVIA060Of0P+uPHGG50tgsFQLImOhkmTYIkllKhbN+mh0revZNdXqiT5Kt26SURYs2YZn/H8889nHCzFmH4wBoPBkA6tpfMjpDnqT5+Wz2XLZLx1a5krC5h+MIYC4dSpUwBcc801TpbEYCh6Tp2CZ56Rba7y5WHkSNn+atRI6n2NGQNr14oj/447pFaYSy5CpSIiIgAIDAwspG9QvDAKxmDH2rVrAZMHYyg7pKTA0aPw998QESH+lZQUmWvfHrp3lzwVK//9l3dfyqRJk4CC7QdTnDEKxmAwlCm2bYMVK+Cbb6Q5V0wMxMdLToqvL3ToAFu3QtWqEmJsq1yg7Dnq84NRMAaDoUxw5Qr884/0Ttm5U8ZWrYL77oOWLaFpUynP0qCBNPTq0kUUjCHvGAVjMBhKFVqLIlm1CurXF8Wxdy8cPiwZ9j4+aQomOVlyVWwxgZQFh1EwBoOhVDFvHjzyiJw3awZ33y1thDt1EkslMRFq1JA6YDff7FxZSztGwRjs6NGjh7NFMBgc4vhxyZ7v3Vsc73//LZbK/v1pa44cgfvvF6vFto3wRx8VubgAvPLKK855sZMwCsZgR4MGDXJeZDA4kcuXpfLwnj3iI5kxQ5z01mz6YcOkd0rDhmKlXHNN8elRf/311ztbhCLFKBiDHVFRUqTaKBqDs0lKgvBwSWisXBkuXBDrZP9+ifwC2e46dQruuQdatEiL+Nq2rfgoFVs2b94MlB1FYxSMwY4NGzYAJg/G4Fyeew7eew/On4epU6U6sSUHmMqV4brrYPVquOEG6UXfvr39/cVRuQA8++yzgMmDMRgMhkInJkYKRjZsKNeXLsGBA+JPOX9exlavhvHjoWdPcdrXqiU1wDw9oUIF58luyBmjYAwGQ5GzcaOUtg8Pl5IrEyfK1tfx4xJmXKuWrKtZUyoUjxxpf3+NGkUusiEPGAVjMBgKjZQU6eTo45M2dvmybHdZ686uXQv+/pKzEhwsVkqNGhJq3KpV8d3uMuSMUTAGg6HAuXhR8k9+/VWc9X//LYmO+/dDVJSMeXpKcy5fXxg1Km2bzErr1s6R3VBwGAVjsKN3797OFsFQgrhyRSyRdu1EYYBscSUkwB9/SF96gJdegrp1xVnfpYtYKX37ioVStarz5C9q5syZ42wRihSjYAx2mDL9Bkd56CH47DNRJitWiAVy+LAkN54/D7VrSw+VKlXEGnngAWnKZaUsRsKXlTL9VpyqYJRSfYB3kI6WH2qtZ6Wbfxx4EOloeRoYpbX+1zKXDERalh7VWt9WZIKXYg4dOgSAr6+vkyUxFAfOn4fNm6FePfGT2I7HxaX1n3/rLSlrX7GitBL28ZEM+2rVJD/F+FGEdevWAdCzZ08nS1I0OE3BKKVcgXnAzcAxYJtSaqXWeo/Nsh1AkNY6Xin1MPA6MNQyd1lrXbb+HCgCfv31V8AoGIPkoTz6qDjqx46VUOHDh+WIiUnrmVK7NrRtC+PGSdSXUSZZM2PGDMAomKKgI3BAa30IQCn1JTAASFUwWutQm/V/APcWqYQGQykmJUWc73/8IRnxY8akzcXHg7t7mhL59ltRJB4eYqF06CBl7ufMkegvgyEznKlg6gFRNtfHgE7ZrH8A+NHmurxSKhzZPpultV6R2U1KqdHAaICG6cNUDIYyzF9/iXMeREl06yb+k8OHpXjklSvSDrh+fcmcf/BBcdTnpkWwoWzjTAWTmSGtM12o1L1AENDNZrih1vqEUsoX2KCUitRaH8zwQK0XAgsBgoKCMn2+wVAaSU6GDz4QC2XXLti+HVxdxSqJjhb/ibs7XL0Kx47B//4nPpOGDSW50ccHpkwRJ73BkBecqWCOAbZxJPWBE+kXKaV6As8B3bTWV6zjWusTls9DSqmNwLVABgVjMJR2tBYFUbOmfTtfV1eYNUuy4wEWLRKFcvy4KBUQB7yHh9TyuvdeceSXM7GlhgLCmT9K24BmSikf4DgwDLjbdoFS6lpgAdBHax1tM14NiNdaX1FKeQM3IAEAhnxy6623OlsEQy546in4/HM4cQLWrRPL48wZSWaMigJv7zQF8/33Ur4+MFBChOvXh2nTjFO+KFmwYIGzRShSnKZgtNZJSqkJwFokTHmR1nq3Umo6EK61Xgm8AXgBy5T8X2ANR24FLFBKpQAuiA9mT6YvMuQKb29vZ4tgsOH8eYiIgD//FCsjfTvfc+dEuQDMny9l6q2hw56eUiCya1cJIe7Z0z4PxVD0tGjRwtkiFClONYa11j8AP6Qbm2pznmksn9Z6M+Cf2Zwhf+zbtw8oe/8jFAeuXJHtKlveeUfK1QNMmCAtf48fly2xqChJZARpvHXhArRpI9ZJgwZQvbqxToob33//PQD9+/d3siRFQ44KRinlYev7yGrMUDr4448/AKNgiort20WB7NghvU2WLZPxlBTZ6rIto/L992lVhD08ZItr5EjJqA8OliRHQ/Fm9uzZgFEwtmwB2jkwZjAY0pGcDPv2yRbXuXNSITg9P1hs+C1b5PzkSak2nJgoVsk110CTJhJSPHCg1PPy9jbhwobiT5YKRil1DZKr4mlxtluN7cqAafNjMNigteSOpC/l9t9/sm0F0s734Yfh7Fnxm5w8KaXsXVzSQofDwiRMuH17USR168IbbxhlYiiZZGfB9AZGIOHDb9mMXwCeLUSZDIYSRa9eUlH43DnZ1rJuYyUniw+kWjWZu3gRJk9O2/ZydxclMnWqRHZ16iSNtowyMZQWslQwWutPgE+UUoO01l8XoUwGQ7Fi8WL47TfYvRsWLEizSKycPi0KBODLL6UwZHS0HElJYpHUqgXNm0NQEAQEiGVSvbpRJobSjSM+GD+lVJv0g1rr6YUgj8HJ3H777c4WociJjoadO0WB9OhhXzUY4KuvYNUqOQ8Lkz7wVgUSHZ22zt0dfvklzRLp0EEUyYQJYtWYiC7D4sWLnS1CkeKIgrloc14e6AfsLRxxDM6mSimsC6K1WBkHDsh2VatW9vPTpkkOCYi/w99fcklOnxYFYpsd//HH4jcByTOpVUsUSO3a4Ocnn9bGWwZDehqUsSY4OSoYrfVs22ul1JvAykKTyOBUdu/eDUCb9PtAxRytJfIqIcG+/ztIxd/HH5fz8eOlDL2Vq1ftqwGvXCnbVufPp43VqAF33ilNszp2lPIqtWqJ095YJYbcsHTpUgCGDh2aw8rSQV4SLSsApllIKSU8PBwongomJUWaXFWrZj/+448weLCUmB8wQLor2mJbRHvbNsknOXtWjgsX4OhR8ZvUri2tfBs3FgViPapUMYrEUDB88MEHgFEwqSilIkmrcuwK1ASM/8VQ4Fy9KltS6fuLHDwoiuPgQdne+vNP+/latUS5gPhRtmyRhlhWJXLggPhHrJnte/eKVeLrK5/e3vD668bpbjAUNI5YMP1szpOA/7TWSYUkj6GUYs3ziIoS57ctsbGiOE6dgsqVxUqxpUoVURwA+/fLEROTpkSOHZM5Dw9RUmvWiB+kRg2xXgIDYfRoUSTVqxsficFQVDjig/lXKdUO6IJYMr8hrYwNhlT++08SB0+elMKKto5xreUXu1VxxMWJIrFSpYooGRDfx549opBiY+U4d05KyCclSW7Jhx+KknB3FyXSpAl88434XmrUkMPT02xrGQzOxpEtsqnAncA3lqEQpdQyrfWMQpXMUKxYuVJ+8f/3n0RNNWliP9+jh8yDWButW6fNKSW9SqwKZts28XdYFUhsrFT5vXJFHOcLFqT5WSpUkMTEV1+Fpk2hUSNRVjVqSO0to0QMhuKLI1tkdwHXaq0TAJRSs4A/AaNgShBXr4olEBub5new5bXXxFkeG3snMzL5L/vee/Dzz3J+880ZFUzt2mkK5s8/7S2Q2FjxbXh6iuWydKnkh0CaApk5Uxzt3t5yXbWqWDbpqwsbDCWZ5cuXO1uEIsURBXMEyX+xdJnAgxLaOTIpScp1eHpKt7+SREqKbB/Fxsr2U/qaVz/9JFZGbKw4xO+8035+8mSYO1fO58yBRx+1n9+3T5IEoQInT2Z8t23k1tatohguXpQorIsXJTy4Zk2xQNatE8c6pCmQadPslYdRIIaySFnrt+SIgrkC7FZK/Yz4YG4GflNKvQugtZ6YHwGUUn2Ad5AItQ+11rPSzXsAnwLtgbPAUK31EcvcM8ADQDIwUWu9Nrt3RUfDm2/KuYeH/PLz9LT/zOrc01N6bmS1JaO1/JKNj5dfyDVr2s+fOAGbNsl8nTrQt6/9/ObNogDi46Vs+1NP2c+/805aLsfEiXJty/btMG+enNevn1HB2JZ9t5Y1SUxMUxJWhdukyQFCQ13x9PRJVR7x8aKcO3eWbanoaNi4Uf5dKlUSpfLUU/JZqZIoI6NADIaMhISEADBixAinylFUOKJgvrUcVjYW1MuVUq7APERpHQO2KaVWputO+QBwTmvdVCk1DHgNGKqUao20WW4D1AXWKaWaa62Ts3qfx+lj1P5gGvFulWnz4RPEx8svz8uX4ci+K+x491cStDvKw4Og8Z3s7j178irhn/0NWlOxmjv3vdIqVfmULw97dibx8hMxXKIizQMqsOhjhaurOKddXeH33+HBB+VZwcFSet06V64cHD4sdaxAFJktKSnyy9vKyZMSjZWYKFtfiYn2iYG7dkkklXXu6lU4dEiUQ4UK0iFx1qy0zocgiuC++yAm5gLJyUnEx/tQtaooq0qVpNWuVYF4eclR0qxAg8HZGAWTkapaa7u/l5VSj6YfyyMdgQNa60OW534JDABsFcwA4EXL+XLgPSX9kwcAX1oanx1WSh2wPG9LVi+rkfQfD0dPJ57yVGg+UH5zJydDSgrRO09R6/zNJKM4daEOFTp/S0J8ClcvJ5NwOYXj/51kePxrXMWd01dq47ljNFcTUkhMSObSVXA9Gs9sfiAFV2J3VWf7B8NIdnUn2cWNFFc3jv9ziQfZRjyeXNlxDR+90zl1TisXDuxOoBaxXKE8e3eUZ+bM8ri6iuWQlAR//y1RU+XLS6mSjz6y/25xcdC7d1r5kh07ZL2bm3x27y5Vf61jnp5pysL2c9mynSgFw4d3yPwf0WAwGBzEEQUzHNnCsmVEJmN5oR4QZXN9DOiU1RqtdZJSKg6oYRn/I9299dK/QCk1GhgNsscG4IKGTz+1W1fxhIQ4uaKpxAUqb7Hr5Ezt86epimT4ncWbGnU6i+fa8mf8WX2IqnyOxoVo13rU8K5MSgqpx3kVRUW+JgVFXEoDvOLuSJ1LVq7Ex5xgOmtIxpX4eF/KnR1Kkos7VKiAqlwJ1SSG57uF4XpNTVwC/CjXpyflvMrj7qFSlYibmxyurnmPrjJRWQaDoaDIruHYXcDdgI9Syrb2WCXEF1IQZPbrTDu4xpF70VovBBYCtK7RQEfe8ggu7uVoM3JQmoJwccEt5iLHLpbDRSejKlei8pNP2s1XPHKci+fP4pKciGfdevD003bvqfHbb7Dxa4iLo+51TeD5J2V/yrJHVWn1ali7CIDq1wbAyH5pe1iJiRKiteGIPKxxTQh0g6tX4FKMOEl2h8PP38m8nx9c3C57a9b9qshIKfnr7S2myn33pc15ecleoIuL7JEZDAZDEZCdBbMZOAl4A7YFLy8AOwvo/ccA2/Ki9YETWaw5ppQqB1QBYhy8144KPrXxX/xUpnPu10D9L97I8l63po1w+/7LrB/epQup4VcpKRlrjvTtK41F4uKgQQNpDGLLxYuiOOLi4Lrr4P770+a0hpdfhu8sCqZdO9kPsw3jOnQo7bCaMVaUEifQzz/LHttdd4kCqlZNkkqqVxcFZDzyBoOhAMmu4di/wL9A50J8/zagmVLKBziOOO3vTrdmJbJNtwUYDGzQWmuLVfWFUuotxMnfDAgrRFkdJ7OCVg0awL33Zn3P7bfLkRlKiZe9UqW0VPnO6f6zbN2adt6vnygQq/K5eFEKdIF49mNjJaRN2xh8a9bA1q3cU6WK1FUJD09TPpUryzO8vEyxLoMhH/zwww85LypFOJLJf4G0rSd3wA24pLWunPVdjmHxqUwA1iJhyou01ruVUtOBcK31SuAjYLHFiR+DKCEs675CAgKSgPHZRZCVeNq3lyMrpk2D4cNFAbVpI1mLtixeLJZNYiLcfTcMHCiKxlrUy9JRq1xcnISoWTtsgVhDX30lVSKvuQaeew5uukmUT9WqslWXmJgx/M1gMNhRoYxtUSutM7gtsr9BqduBjlrrZwtHpMIjKChIW8vRl0lSUkSZWJN8bLn1Vli/nm0BATBlCh1uuknWnjsnn6NGieIBeOihtFR8pcTCmTNHYqUbNpQY6M6dTZ17gyEd77//PgDjxo1zsiS5Qym1XWsdlPNKe3LdD0ZrvUIpNSW39xmKAS4uGWvEWFm9GlJS2PPhh3D+PB2qVBEF4eMjW2m2Pp2pU+XTav3ExEj2pXX7bd062bLz8JCM01q14P33RSm1by/KypQ0NpRBvvrqK6DkKZi84sgW2R02ly5AEJlEaxlKAS4umTv6lZJM0AsX5NPasMX6GRsryTogW2YTJ0q/YWvP4b/+grWWIgshIRLIUKeOKJ6aNWWrbfNmUT7Nm5utNoOhlOCIBdPf5jwJqU02oFCkMRRvKlWCgICM41WrSjmEqCgpt9yokRxWtmxJs3p8fOQZ0dESWp2QAGfOpNW58faWHCVrO8maNdNq9JgAA4OhROFIP5iRRSGIoYRTrpwoDx+fjHMtW8LXX0sNmypVJMoNZOvtwgX45JO0tdYyBFevpo1FR8OiRWLd9O0Lzz8vys5gMBRrHNkiqw/MBW4greHYo1rrY4Usm6G0UK0a3HGHHLZYAwRatZIgg1274JZb4JlnZNstOlq22T7+WJrFREbKPdY6N/XqiV8nIUHmu3SxL/tsMBicSo5RZJYqyl8Aiy1D9wL3aK1vLmTZCpwyH0VWEsgsSfXRR+Hdd+X8pZekUsGJE3D8uPRMXr8efvtN5ocNg0ceEeVTp47U0DEYDPmiMKPIamqtP7a5DlFKTcrtiwwGh8jMzzJnjpTmCQuDa6+19+8kJKR1QgNRKD/9JOfWVpqbN8u6Ll1gyJC0EGuDwVCoOKJgziil7gWWWK7vouBqkRmKGZs3bwbg+uuvd7IkNiglSiGzSgfly8ONN4rP5q+/xMLx9hYLx2rlrF8Pp06JH2jnTinFY91eq1cvbcvNYChk3rQ0pJo8ebKTJSkaHNkiawi8h5SM0UiNskctpWRKFGaLLGc+sTjchw8f7mRJ8kBCgoRZ2yZ3Xrwofh7rz/mKFeLfOXFC/DYg3dsqVICOHeHVV6FFC5MgaigUgoODAdi4caNT5cgthbZFprU+CtyWJ6kMhqKkfPmMY+7u8OOPsr128qT0kwZROGfPynhsrByrV0vBUS8vaNxYjoYNRdmk71FtMBhyJNeZ/AZDicLdXSpP9+5tP65UWlUDd3fZYrvxRhg8GI4ckWPvXsnrmT9fttJuvhlef13uMxaOwZAjRsEYyja33CLWyx9/SC5P27ZygIy/8oqcHz8u5W/mzZPe040bS86Pl5ckmhoLx2DIgFEwBjvKlSuDPxKentICIT1Vq4p1Y7VwRo2SKtKHD4uFs3u3tD0IDRULZ8wYGDdOSt8YC8eQCZ5lrAafI05+D2AQ0BgbhaS1nl6okhUCxslvyBPx8WLhtGoluTUgPpxz56BPH9i2TcYGDpQyOJUrp1k4rq7QpIlYPQZDCaUw82C+A+KA7cCV3L7AYCjxVKgAPXrYjymV1hHU2mdn2jT5PHIEDh6UkOhFiyS44NprYfZsycUx1o2hjOCIBbNLa+1XoC9VqjqwFLGKjgBDtNbn0q0JBD4AKgPJwEyt9VLLXAjQDVF8ACO01hE5vddYMDmzadMmALp27epkSUoQ8fGwfbsECVjRGv7+W5q/Wf8fmzxZLKBmzaSumq+vjJvq0WWGl19+GYAXXnjByZLkjrxaMI6Up92slPLPg0zZMQVYr7VuBqy3XKcnHrhfa90G6APMUUpVtZl/UmsdaDlyVC4Gxzh8+DCHDx92thgliwoV7JULiJVy9aooGICgILjnHmlxEBkJS5ZIhelKlcSqmTOn6OU2FDnr169n/fr1zhajyHBki6wLMEIpdRjZIlOA1lpnUrfdYQYAwZbzT4CNwNO2C7TW/9icn1BKRQM1gdh8vNdgKDrathVlcvSoFO0MDJQjKUnG5s+XZM/ff4d//5U8nubN5ahb17QnMJR4HFEwfQvhvbW11icBtNYnlVK1sluslOoIuAMHbYZnKqWmYrGAtNaZ+oeUUqOB0QANGzYsCNkNhtzRsKEcVsqVk+0xa5M2kDYEHh5StHPTJgkKiI6W1gWDBkmZHGuAgcFQQnAkk/9fpVRbwLoH8KvW+q+c7lNKrQMySw54LjcCKqXqIJWch2utUyzDzwCnEKWzELF+Mo1q01ovtKwhKCjIdOI0FB/eegvGj5fozL6gAAAgAElEQVRKA716ieVy+TIcOAD//AOLF4sFtH69hEI//risqV7d2ZIbDA7hSD+YR4GHgG8sQ58ppRZqredmd5/Wumc2z/xPKVXHYr3UAaKzWFcZWA08r7X+w+bZJy2nV5RSHwNlo3JcEVDW4vSdTpMmMGFC2rWnJ/j7S0j0Aw+kjTdrBmvWyOHtLfO7dkkztx49TKBACaFGjRrOFqFIcSSKbCfQWWt9yXJdEdiSHx+MUuoN4KzWepZSagpQXWv9VLo17sCPwPda6znp5qzKSQFvAwla68wCBewwUWSGEoPWYsX88INUEFiyRPJu9u+X6LSDB6VszZUrUk1gwwYJJDAh0IZCoDDzYBQSJmwl2TKWH2YBXymlHgCOAncCKKWCgLFa6weBIUBXoIZSaoTlPms48udKqZoWOSKAsfmUx2AoXiglVZ1btEgbq14dOnWS47vv0qpBu7rCqlXS96Z1a4lcq11botiMRWpwIo5YMI8Dw4FvLUO3AyHprYqSgLFgcsYaQnnTTTc5WRJDtuzbJ1Foy5dLBYERI2DPHjh0SLqCnj0L//uf+HbGjJGW1Aan88wzzwDw6quvOlmS3FGY5frfUkptRMKVFTBSa70j9yIaSgLHjh1ztggGR2jRAt5+W6oDXL4sUWft2sn533/DlCli4Xz/vSgbNzexburVM9toTmTLli3OFqFIyVLBKKUqa63PW7Luj1gO61x1rXVM4YtnMBiyxcXFvs6Zp6eUpUlMTBvr1k38OJs3Q5Uqomj+/luKeQ4YIJ8GQyGQnQXzBdAPqUFmu4+mLNe+hSiXwWDID+vXS5TZihVizSQny7banj2ibObMkaABV1fZZhswwFg2hgInSwWjte5n+fQpOnEMBkOBoJSEO/tbqjy5uaX1ugkLk8KcIAomPFyi0lq3lo6e9etLgICHh/PkN5QKcqxFoZTKUDgnszFD6aBy5cpUrlzZ2WIYCpOGDWHWLAlrHjgQhg4VpbJ9O3z0Ebz4olSKHjpUEjwNBUb9+vWpX7++s8UoMrKMIlNKlQcqAKFI3TCr/VwZ+FFr3aooBCxITBSZwZCOxMS0JM0rV6RN9MsvwzeWvGo/P1E6fn5S1NNQJimMKLIxwCSgLuKHsSqY88C8XEtoMBiKH7YVADw8pBinbYBAhw6S7LlmjVQTCAiQ+bp1oVa2JQQNBofyYB7JqSxMScFYMDmzdu1aAHr37u1kSQxOQ2vZLvviC5gxAy5ckOZpO3fK+aefSvXnbt3g/felXI3BISZNmgTAnBLWnqEwM/lTlFJVtdaxlhdVA+7SWr+f25cVRxITEzl27BgJCQnOFqVYUN1SSHHv3r1OlsRQkJQvX5769evj5kjNMqXEPxNk+X1SoQLcfDPcdJNUe37pJRkPDYWlS6WfTUAA1KxZeF+glBARUbZaVzmiYB7SWqduiWmtzymlHgJKhYI5duwYlSpVonHjxigTpsmZM2cA8Pb2drIkhoJCa83Zs2c5duwYPj75CAp1cZG6Z126iKLp0gWaNpV+Nr/+Ku0EmjSREOlRo2RLzVCmcUTBuCillLbspSmlXJEy+aWChIQEo1wMpRqlFDVq1OD06dP5f1i7dqJMDh2CS5ckDPriRcm52bkT5s6FlSslSu2222DZMnAvNb8uDLnEEQWzFilMOR9JsBwLrClUqYoYo1wMpZ0C/xn3tcmz9vKC666TIyQkbfzyZXjzTWkt0LYt+PiYZM4yhiMK5mkkouxhJJLsJ+DDwhTK4DzKlXPkR8JgyIKnnoJPPpEmajNnwsmTUj3gr7+kGvSOHdCoEYwdC9dk1o+wdNO8eXNni1Ck5BhFVprILIps7969tGpV4lJ6DIZcU6Q/6xcuQKVKcp6UJErm119h0iS5dnGBn36SZmnGqin25DWKzJFM/huUUj8rpf5RSh1SSh1WSh3Km5gGg6FMYFUuAOXKSZSZi4soF5CIs02b4L33YMsWiI93jpyGQsWR/ZCPgMeQZMvkHNY6jKVK81KgMVKpeYjW+lwm65KBSMvlUa31bZZxH+BLoDrwJ3Cf1vpqQclXklmxYgWrV68mOjqa8ePH06tXL4fvjY2NBaCqqbBrKGhGjhTFs2ABDBkCXbtKHbS1ayXyDOT6scdkrhRaNqNHjwZg4cKFTpakiNBaZ3sAW3Nak5cDeB2YYjmfAryWxbqLWYx/BQyznM8HHs7pne3bt9fp2bNnT4axoubo0aM6ODhYt2zZUrdu3VrPmTPHbn706NH6t99+0x06dNABAQG6devWeurUqTk+NyYmRo8aNSpXspw+fVqfPn06V/cYSgbF4Wc9lZSUtPNTp7RevVrrgACtJc1T61GjtI6Pd558hUS3bt10t27dnC1GrgHCdR5+z+e4RQaEKqXeUEp1Vkq1sx4FoNsGAJ9Yzj9BOmU6hJKQmB7A8rzcX9woV64cs2fPZu/evfzxxx/MmzePPXv2pM5v3bqV6667jg0bNvDXX38RERHBmjVr+OOPPzJ9XnKyGJozZsxg/PjxRfIdDIZcYWud1K4tEWg2P/O4uUkztW+/hagoUTuGEocjW2SdLJ+2Dh6N/ILPD7W11icBtNYnlVJZFTYqr5QKB5KAWVrrFUANIFZrbdnQ5RhQL7OblVKjgdEADRs2zKfIhUOdOnWoU6cOAJUqVaJVq1YcP36c1q1bs3fvXpo3b46rqyteXl6AVB9ITEy0Cz298847adCgATt27KBHjx5cvHiRvn370q5dQfwtYDAUMtWrw7Ztsn128KBUdA4Ph8hIiUCrWVNK1wwZIu2hq1RxtsQGB3CkZXL3vD5cKbUOyCwW8blcPKah1vqEUsoX2KCUikQKbqYn0z9xtNYLgYUgUWS5eK9TOHLkCDt27KBTJ9HrP/74I3369AHEMmnfvj0HDhxg/PjxqWsAIiMjadWqFaGhobz77rusWLGCuLg4Dhw4wNixY53yXQyGXBEYCB98INaKUtCvH/TqJUomJEQ6cm7eLOHP27dLi4FS6KcpTTgSRTY1s8ORh2ute2qt/TI5vgP+U0rVsbyjDhCdxTNOWD4PARuBa4EzQFWllFVB1gdOOCJTTrz4ovzMKiXn6XniibT52bMzzo8enTafWz/exYsXGTRoEHPmzEntybJ27dpUBePq6kpERATHjh0jLCyMXbt2AVKNICYmhqlT5T/LxIkT2b59O/Pnz8+1cnFzc3OsXpXBUFjYKg13d2jfXvJprLRoIS0EFiwQK+dqyYntCQwMJDAw0NliFBmObJFdsjkvj7RRLohKiCuB4cAsy+d36RdYCmvGa62vKKW8gRuA17XWWikVCgxGIskyvb8kkZiYyKBBg7jnnnu44447AIiPjyc2Npa6devara1atSrBwcGsWbMGPz8/du/eTadOnQokSbKK2XowFEc++gi6d5e/2ubPlyoB4eGwapVEoB07BpUrS55NMS66WdKqKOeXHC0YrfVsm2Mm0nwsU39HLpkF3KyU2g/cbLlGKRWklLJWCmgFhCul/kIan83SWls9gU8DjyulDiA+mY8KQCanoLXmgQceoFWrVjz++OOp46GhoXTvLjuUp0+fTg0hvnz5MuvWraOlpUx6ZGQkAQEBRS+4wVBUVKkC48eLP6ZNG6n0PGaMFNWsVw8+/hheeUX61Pzwg7OlNVjIy5+8FQDfHFflgNb6LHBTJuPhwIOW882Afxb3HwI65leO9Lz4YuZbY1Zmz858a8zKwoW53xr7/fffWbx4Mf7+/qnm8yuvvMKPP/7I4MGDATh58iTDhw8nOTmZlJQUhgwZQr9+/QBRMB07Fsw/xblzkopUrVq1AnmewVBoKCXtn5WSwpsgrQW2bIGzZ6FzZ6nuXIz8NPfeey8An332mZMlKRpyVDAWp7rVOe4K1ASmF6ZQZY0uXbpYc3vseP7553n77bcBCAgIYMeOHZnePzs7jZdLrCHOBkOJYcgQUSJvvy0VnLt3h61b4bPPpOtmw4Zw/Djcfz94ejpV1GPHjjn1/UVNlgpGKeWjtT6M+FysJAH/2YQHGwqRP//809kiGAzFHzc3uPtuuOsuSE6W0jSdO0sLgc2bpSvn1q3w9NPw1luyrWYoErLzwViTGBdprf+1HMeNcjEYDMUSpUS5ALi6SouAu+6SMGeAuDipf7ZqlWyhGQqd7LbIXJRS04DmSqnH009qrd8qPLEMBoOhAHBzg5dfhnfegfLlYeBAaRmwfTs0bw7XXw9Hj8INN0gxTkOBkt2/6DAgAVFClTI5DKUQd3d33E0HQkNpoWJFePxxqQ6wdi0MGJBWTDMqSpI2u3aVYIAiKEDZuXNnOnfuXOjvKS5kacForfcBrymldmqtfyxCmQxOxJrgaTCUKsqVg8aN5dzLSwIBunSRoACAI0ckcdPPT9pCly9fKGK8+uqrhfLc4oojeTBGuRgMhtKHmxtce21a75p+/aQJ2ttvi7UTFwcJCc6VsYRjNh0NdsTExBATE+NsMQyGouGVVySEeckSeOklqfXUvLlEnb3zjlg0/ftDAUV0Dho0iEGDBhXIs0oCRsEY7EhJSSElJcXZYhgMRUelSjBsmJzXrQuDBsHEiVJ08+BBiTrr3FnO88nZs2c5W4Yi2LLLg7kjuxu11t8UvDgGg8FQDKhaFRIT067bt5fEzTZtJCigVlbdRQy2ZBem3D+bOQ0YBWMwGEovM2ZIlYBZs2DaNIk6CwuTBM5WreDCBSms2b+/CXHOguyiyEYWpSCG3LFixQpWr15NdHQ048ePp1evXs4WyWAofQQESKMzkDYB118v/pnNm+G110TJtGgBK1eK78Zgh0PFLpVStwJtkHL9AGitTT0yJ3L77bdz++23c+7cOSZPnlxgCsbDw6NAnmMwlEoqVJAQ5337RLmAWDabNoGHBzRqlO3tN92Uob5vqcaRYpfzkQrK3YEPkR4sYYUslyEbkpOTcXV1BWDGjBmMHz++wJ5dqZLJoTUYcmTAADhwQDpw3n8/nDsnLQMaNYJu3cDbW8Kg0xXXfOGFF5wksHNwxIK5XmsdYEm4fEkpNRvjfyly7rzzTho0aMCOHTu46aabeO6555gyZQp9+/alXbt2zhbPYChb1KkDb74JzzwjlouHh4Qy//47fPqp+GrCwuCpp+Dhh9NybcoYjnimLls+45VSdYFEwCc/L1VKVVdK/ayU2m/5zNB8RCnVXSkVYXMkKKVut8yFKKUO28wVXA9SZ/ZMzobIyEi8vLwIDQ3l+eefZ+7cuaxbt47ly5czf/78AntPWQujNBjyRY0aUhnAzQ06dZLw5q5dYcMGKaj59NOwaJGEPAN9+/alb9++Tha66HDEglmllKoKvAH8iUSQfZj9LTkyBVivtZ6llJpiuX7adoHWOhQIBFFIwAHgJ5slT2qtl1NK6NmzJ6dOncowPnPmTHr37k1MTAxTp05NHZ84cSITJ04scDky60tjMBgcpFw5UTje3pLAWbs2xMTIVlrXriTEx6OLUQO0wiZHBaO1ftly+rVSahVQXmsdl8/3DkBaLwN8AmwknYJJx2DgR611fD7fW2xZt25dlnPbt2+nU6dOlCuXlwakBoOhSOnYUZIyFy8WRePrC7/+CsuXc/uJE+wsXx4mT4YpU2S+FKOy+otVKdVDa70hq4TL/CRaKqVitdZVba7Paa2z7NGrlNoAvKW1XmW5DgE6A1eA9cAUrfWVLO4dDYwGaNiwYft///3Xbn7v3r20atUqr1+lSAgJCWH//v3MnDmz0N915swZALxL+Q9+WaQk/KyXWrSGvXt56/bbqX7iBCMuXZKItHnzYMQIZ0uXI0qp7VrroNzel50Pppvls38mR7+sbrIRaJ1Salcmx4DcCKiUqgP4A2tthp8BWgIdgOpkY/1orRdqrYO01kE1a9bMzauLDZGRkQQEBDhbDIPBkFeUgtatWVm1KqGXLslYfLwkbe7fn+qjKW1kl2g5zfKZp4RLrXXPrOaUUv8ppeporU9aFEh0No8aAnyrtU6t26C1Pmk5vaKU+hiYnBcZSwqzMwsmKCTKF1KZcoPBAP3uvBNatoTwcClHU78+fP45NGwIN90kYc6JiRI0UArIrhZZhi6WtuSzo+VKYDgwy/L5XTZr70IsFlvZrMpJAbcDu/Ihi8EGLy8vZ4tgMJRaJj/5pJwkJ0vuTLVqEt68aZPk0SQnS9TZ9OkwfHhaC+gSSnbSWwO3WyBbUSst1/2BTfl87yzgK6XUA8BR4E4ApVQQMFZr/aDlujHQAPgl3f2fK6VqAgqIAMbmUx6DwWAoOlxd0xz8HTpAYKDkzdx3H5w4AQ8+KFbOBx84V858kt0W2UsASqmfgHZa6wuW6xeBZfl5qdb6LJChZoLWOhx40Ob6CFAvk3U98vN+Q9YYJ7/BUHgEBwcDsHHjRvsJNzfpuGktP+PiAu7u8M03EBwM1asXoZQFhyP2V0Pgqs31VaBxoUhjMBgMZZV69eDwYXj3XTh9Gm65Ja16c7t20uJ550648UYJGigBOKJgFgNhSqlvkSTLgcCnhSqVwWAwlEW8vODZZ9Our7tOcmi2b5cSNMuWiaJ57TWp7FzMcSTRcqZSag3QxTI0Umu9o3DFMhgMBgOVKokl07GjtA4A+O03yZ9p1w6KedSnQ11ytNbbgSXAt8BZpVTDQpXKYDAYDGm4u0tjM1dXSdCsXx/eeUesm6tXc77fSeSoYJRStyml9gOHkWiuw8CPhS1YWUIpxX333Zd6nZSURM2aNenXL8d81kxZuXIls2bNytO9np6eeKYrMZ4fNm7cmPo9cpIrIiKCH374ocDebTAUN4YMGcKQIUNyf2PlyrBggfSh+fJLeOwxaNAA1q8Xn83GjVK1+fDhApc5Pzjig3kZuA5Yp7W+VinVHclNMRQQFStWZNeuXVy+fBlPT09+/vln6tXLEDznMLfddhu33XZbnmUpLHKSKyIigvDwcG655ZZCk8FgcCbjxo3L3wOaNJED4O67pdnZ+vXS3nn9evjwQ6kC/9xz+Za1IHBkiyzRElbsopRysa1yXBoJDg7OcLz//vsAxMfHZzofEhICSIhv+jlH6du3L6tXrwZgyZIl3HVXmg4PCwvj+uuv59prr+X6669n3759ALz11luMGjUKkHIyfn5+xMfHExISwoQJEwAYMWIEDz/8MN27d8fX15dffvmFUaNG0apVK0bY1ECyJlimpKTw1Vdfpc45er8ta9asoWXLlnTp0oVvvkkrWWcr17Jly/Dz86Nt27Z07dqVq1evMnXqVJYuXUpgYCBLly7N8nuHhIRwxx130KdPH5o1a8ZTTz1l9+527drRtm3b1O6Bly5dYtSoUXTo0IFrr72W777LLq/XYCg84uPjiY8vwJq9DRrA4MGwbZtcJyXB7t1w5EjBvSMfOGLBxCqlvJDkys+VUtFAUuGKVfYYNmwY06dPp1+/fuzcuZNRo0bx66+/AtCyZUs2bdpEuXLlWLduHc8++yxff/01kyZNIjg4mG+//ZaZM2eyYMECKlSokOHZ586dY8OGDaxcuZL+/fvz+++/8+GHH9KhQwciIiIIDEz7eyEmJoaLFy/m+f6EhAQeeughNmzYQNOmTRk6dGim33f69OmsXbuWevXqERsbi7u7O9OnTyc8PJz33nsPgPPnz2f6vUGsnR07duDh4UGLFi145JFHKF++PA899BCbNm3Cx8eHmJgYQFoe9OjRg0WLFhEbG0vHjh3p2bNnoVprBkNmWK3zDHkw+cHLC1avlgrN//0nzv+QEGjRAm6+2akVmx1RMAOQpmOPAfcAVYDphSmUM8nuP3yFChWynff29s7zD05AQABHjhxhyZIlGbaI4uLiGD58OPv370cpRWKilGVzcXEhJCSEgIAAxowZww033JDps/v3749SCn9/f2rXro2/vz8Abdq04ciRI3YKIr/3//333/j4+NCsWTMA7r33XhZm0njthhtuYMSIEQwZMoQ77si0YHeW3xukt3mVKlUAaN26Nf/++y/nzp2ja9eu+PhIP7zqluS0n376iZUrV/Lmm28CogSPHj1qKgsbSg9dusCWLXDqlCiUP/6QaLP335cttcWLYeZMaN++SMXKcYtMa31Ja52itU7SWn8CzAP6FL5oZY/bbruNyZMn222PgfTx7t69O7t27eL7778nISEhdW7//v14eXlx4sSJLJ/r4eEBiEKynluvk5LEGFU2iVtXrlzJ9f22KAeSwObPn8+MGTOIiooiMDAw0y6a2X1vWzlcXV1JSkpCa53pu7XWfP3110RERBAREWGUi6F0opS0cnZzk2TMiRNFocyZA2vXQlAQvPBCkYqUpYJRSlVWSj2jlHpPKdVLCROAQ0iFY0MBM2rUKKZOnZpqIViJi4tLdfpb/T3W8UcffZRNmzZx9uxZli/Pe4PP2rVrs3fvXlJSUlJ9QXmhZcuWHD58mIMHDwLiT8qMgwcP0qlTJ6ZPn463tzdRUVFUqlSJC9ZSGWT9vbOic+fO/PLLLxy2RNJYt8h69+7N3LlzU7t17thh0rgMZYCKFSV/JjIybey//+S6iNoDZGfBLEYKXUYi9cF+QopSDtBa56qni8Ex6tevz6OPPpph/KmnnuKZZ57hhhtuIDk5OXX8scceY9y4cTRv3pyPPvqIKVOmEB2dXeeDrJk1axb9+vVj4MCB1K5dO8/foXz58ixcuJBbb72VLl260KhRo0zXPfnkk/j7++Pn50fXrl1p27Yt3bt3Z8+ePalO/qy+d1bUrFmThQsXcscdd9C2bdtU/88LL7xAYmIiAQEB+Pn58UIR/xVnMDiNmjWlvMyAAZL5f+218PXXEm129KisKURlk11Hy0ittb/l3BU4AzS0Fr0siQQFBenw8HC7MdPlzx5rhEtmwQKGko35WXc+Vks8qwjMQuXyZfDwEIWzfr0U1qxQQQICZs+GPn2yrHGW146W2Tn5bRt8JSulDpdk5WJwDKNYDIbCwymKxYo1gTowEFq3lqCABx+EAwekHM3TT0MeE7SzIrstsrZKqfOW4wIQYD1XSp0vUCkMxYbk5GSHtqMMBkPuOXPmTGpLDKfi7i7tAU5amgMrBQkJonQyCdzJK1kqGK21q9a6suWopLUuZ3NeOb8vVkrdqZTarZRKsTQay2pdH6XUPqXUAaXUFJtxH6XUVqXUfqXUUqWUe35lMkjOy7lz55wthsFQKhk8eDCDBw92thhCo0Zw8CCMGQP33y9RZmvXSiHN3bshJUW21fKBM/tx7gLuABZktcDi+5kH3AwcA7YppVZqrfcArwFva62/VErNBx4ASnb7N4PBYChKateG+fPF0a+UbJf99JO0BTh3DpYuhVdfzfPjnaZgtNZ7IceciY7AAa31IcvaL4EBSqm9QA/gbsu6T4AXMQrGYDAYco/193DTpuDrK+2a+/eH6GiwlKPKCw6V63ci9YAom+tjlrEaQKzWOindeAaUUqOVUuFKqfDTp08XqrAGg8FQ4nFxkerNVqVjk9Sc60cVkEiZopRap5TalcnhaB5NZuaNzmY846DWC7XWQVrroJo1azoqepFS0OX6C5rrHeic9+CDD7Jnz54ikCZnRowYkZp0mpNcISEh2VZBMBjKJC1bynbZtGnwyit5fkyhbpFprXvm8xHHgAY21/WBE0hOTlWlVDmLFWMdL5EUdLn+/MqSns2bN+d434cfflgY4uSbnOQKCQnBz8+PunXrFpFEhrLMww8/7GwRHMfLS0r/AzzxRJ4eUdy3yLYBzSwRY+7AMGClluzQUMAajjEcKJAa7J988kmGY5ulFHZiYmKm8xEREYAkKaafc5TsyvVnVW7+yJEj3HjjjbRr14527dqlKoKNGzcSHBzM4MGDadmyJffccw+ZJdQGBwfz2GOP0bVrV1q1asW2bdu45557CAgI4Pnnn09dZy3ln91zg4ODsSaxenl58fTTT9O+fXt69uxJWFgYwcHB+Pr6snLlSsC+dD9Av379UguFOnK/LVprJkyYQOvWrbn11lvtqhlY5UpOTmbEiBH4+fnh7+/P22+/zfLlywkPD+eee+4hMDCQy5cvM336dDp06ICfnx+jR4+2+35PP/00HTt2pHnz5qmVrpOTk5k8eTL+/v4EBAQwd+5cALZv3063bt1o3749vXv35qQ1HNRQphk6dGiWFcZLI05TMEqpgUqpY0BnYLVSaq1lvK5S6gcAi3UyAVgL7AW+0lrvtjziaeBxpdQBxCfzUVF/h4Jk2LBhfPnllyQkJLBz5046deqUOmctN79t2zZCQ0N58sknuXTpErVq1eLnn3/mzz//ZOnSpUycODH1nh07djBnzhz27NnDoUOH+P333zN9r7u7O5s2bWLs2LEMGDCAd955h4iICEJCQjItQOnIcy9dukRwcDDbt2+nUqVKPP/88/z88898++23TJ06Ncd/i9ze/+2337Jv3z4iIyP53//+l6nFFRERwfHjx9m1axeRkZGMHDmSwYMHExQUxOeff05ERASenp5MmDCBbdu2pVqUq1atSn1GUlISYWFhzJkzh5deegmAhQsXcvjwYXbs2MHOnTu55557SExM5JFHHmH58uVs376dUaNG8VwxaQBlcC5RUVFERUXlvLCU4Mwosm+BbzMZPwHcYnP9A5Chj64lsqxjQcs1fPjwLOfc3Nyyna9QoUK289mRXbn+rMrN161blwkTJhAREYGrqyv//PNP6j0dO3akfv36AAQGBnLkyBG6dOmS4b3WDpP+/v60adMGT09PLl++jK+vL1FRUdSoUcNuvSPPdXd3p0+fPqnP9fDwwM3NDX9/f4440Agpt/dv2rSJu+66C1dXV+rWrUuPHj0yrPH19eXQoUM88sgj3HrrrfTq1SvTd4eGhvL6668THx9PTEwMbdq0oX///gCpbQXat2+fKse6desYO3Ys5crJ/0rVq1dn165d7Nq1i5tvvhkQK6dOnTo5fm9D6cfqay3QfjDFGGfmwWjkDGkAABoOSURBVBjSYS3Xv3HjRjvrwVpuvkWLFnbrX3zxRWrXrs1ff/1FSkoK5cuXT53LrJx9ZuS2FL8jz3Vzc0sNP7d9ru0zy5UrR0pKSuo9tqX4Hbk/PTm1CKhWrRp//fUXa9euZd68eXz11VcsWrTIbk1CQgLjxo0jPDycBg0a8OKLL2baIsD2e2fWIkBrTZs2bdiyZUu2MhkMpZ3i7oMpU2RVrj+rcvNxcXHUqVMHFxcXFi9eXKJKvDRu3JiIiAhSUlKIiooiLCwsz8/q2rUrX375JcnJyZw8eZLQ0NAMa86cOUNKSgqDBg3i5Zdf5s8//wSwaxFgVSbe3t5cvHjRofYHvXr1Yv78+akKJyYmhhYtWnD69OlUBZOYmMju3buze4zBUCoxFkwxIqty/S+88AKTJk0iICAArTWNGzdm1apVjBs3jkGDBrFs2TK6d+9eoloA33DDDfj4+KSW7G/Xrl2enzVw4EA2bNiAv78/zZs3p1u3bhnWHD9+nJEjR6ZaTa9aspNHjBjB2LFj8fT0ZMuWLTz00EP4+/vTuHFjOnTokOO7H3zwQf755x8CAgJwc3PjoYceYsKECSxfvpyJEycSFxdHUlISkyZNok2bNnn+jgZDSSTLcv2lEVOuP2eshfi8ndjH21A4mJ915xMcHAyUPB9MYZTrN5RBrCHJBoOh4Hkij/kkJRWjYAx22AYKGAyGgsUakVhWME5+gx2JiYkkJibmvNBgMOSaffv2sW/fPmeLUWQYC8ZgR1xcHGB8MAZDYTBmzBig5Plg8oqxYAwGg8FQKBgFYzAYDIZCwSiYYoCrqyuBgYGph7UMyauvvkrTpk1p0aIFa9euLbD3HTlyBD8/PwDCw8PtaphlxsaNG7NsHXDXXXcREBDA22+/XWDybdy40a6e2Pz58/n0008L7PlWZs6cSZs2bQgICCAwMJCtW7cCMGfOHOLj4wvkHX///TedO3fGw8MjtdSPlTVr1tCiRQuaNm3KrFmzUscPHz5Mp06daNasGUOHDuXq1asAXLlyhaFDh9K0aVM6depkVzansH5WDIZ8obUuM0f79u11evbs2ZNhrKipWLFihrHdu3frgIAAnZCQoA8dOqR9fX11UlJSts/Jad7K4cOHdZs2bTKdO336tD59+rTdWGhoqL711lszrD158qRu2LChQ+/MDdOmTdNvvPFGgT/Xls2bN+vrrrtOJyQkaK3lex8/flxrrXWjRo0y/Bvklf/++0+HhYXpZ5991u47JSUlaV9fX33w4EF95coVHRAQoHfv3q211vrOO+/US5Ys0VprPWbMGP3+++9rrbWeN2+eHjNmjNZa6yVLlughQ4ZorR3/WSkOP+tlnW7duulu3bo5W4xcA4TrPPzONU5+G9asgVOnCvaZ11wDlrqNueK7775j2LBheHh44OPjQ9OmTQkLC6Nz585267y8vHj88cdZu3Yts2fPZsOGDXz//fdcvnyZ66+/ngULFqCUSq3qW6FCBbvilBs3buTNN99k1apVhIWF8eijj3L58mUqVKjAxx9/nKH+mS29evUiOjqawMBA5s6dywsvvMCbb75JUFAQZ86cISgoiCNHjhASEsLKlSuJj4/n4MGDDBw4kNdffx2Qv+KfffZZkpOT8fb25qOPPmL+/Pm4urry2WefMXfuXNavX4+XlxeTJ08mIiKCsWPHEh8fT5MmTVi0aBHVqlUjODiYTp06ERoaSmxsLB999BE33nhjlrKfPHkSb2/v1Ppi1qCGd999lxMnTtC9e3e8vb0JDQ3lp59+Ytq0aVy5coUmTZrw8ccf4+XlRePGjRk6dGhqaZovvviCpk2b2r2nVq1a1KpVK7UVg5WwsDCaNm2Kr68vINW0v/vuO1q1asWGDRv44osvACm++uKLL/Lwww/z3Xff8aKlP8fgwYOZMGECWmuHf1YMzse2DUZZwGyRFQMuX76cuj02cOBAQEqbNGiQ1mutfv36HD9+PMO9ly5dws/Pj61bt9KlS5csy82PHDmSd999N9sCjC1btuTXX38lIiKC6dOn8+yzz2Yr98qVK2nSpAkRERHZ/jIHKZe/dOlSIiMjWbp0KVFRUZw+fZqHHnqIr7/+mr/++otly5bRuHFjxo4dy2OPPZbpc++//35ee+01du7cib+/f2rZfMi8nH5W9OrVi6ioKJo3b864ceP45ZdfAJg4cSJ169YlNDSU0NBQzpw5w4wZM1i3bh1//vknQUFBvPXWW6nPqVy5MmFhYUyYMIFJkyZl+05bsvrv+//2zj6qqjrd459H0VBR0DCHhlFxMCeM0tHQcpYXM8aXTM0xXyoVceWaSi0bNb293Fa6RqtpyhaNxSzf86Y3bYQxNbR0RnwZgcJCuy69YhNBasiLmInIc/84h80BDnBUOBsOv89aZ529z/7tvZ/9nH32c35v3ycvL4+goCBLndn1e3fdx8/Pj8DAQPLy8jy+Vwz2c//993P//Teah7HpYGowLlxPTaM+aNOmjZW0rBx1I+HjTjG4ZcuW/O53v7PW3cnNDx48mIKCAkuja8qUKezYsaPasQoLC5kyZQonT56kRYsW9TofZujQoQQGBgIQERHBN998Q35+PoMHDyYsLAxwSN3XRmFhYaXrmDZtGg8//LC13Z2cfk0EBASQnp7Ovn372LNnDxMnTmTZsmXExsZWKnfo0CGOHTvGoEGDACgpKalUMyhPDDd58mTmzp1bhxcqqOn7re17v559DI2L8t95nz59bLbEO9gSYETkYeBl4HYgSlXT3JT5BbAO+BlQBiSo6nLntpeBx4FzzuL/qY68MT5DaGhopcRE2dnZbtP6+vv707JlS6BmuXl1IynvjhdffJEBAwawcuVKiouLLd0kT3GV4HeVuQf3Mv+e2uUp7uT0a6Nly5ZER0cTHR1NZGQka9eurRZgVJWYmBg++OADt8dwtf9arqWm7zc4OJiCggJKS0vx8/Or9L2X7xMaGkppaSmFhYV06tTJ43vFYD/ltVwzD6ZhyQTGAf+spUwp8AdVvR0YCDwlIhEu299U1T7Ol08FF3Dkhtm4cSOXL18mKyuLEydOEBVVe361muTmg4KCCAwMJCUlBYANGza43b9c/h8cKY2vle7du5Oeng7gkdT9Pffcwz/+8Q+ysrIAh9Q9VJbQdyUwMJCOHTta6YrXr1/vVjnZle+++46hQ4dW+/z48eOcOHHCWs/IyKBbt27Vzj9w4ED279/PyZMnAUdabNfEbps2bbLer6XP4+677+bEiRNkZWVRUlLCxo0bGT16NCLCkCFDLP+tXbuWMWPGAI57ojwN9+bNm7nvvvsQkeu6VwwGb2BLDUZVv4ba//Gpai6Q61y+ICJfAz8HjnnDRrvp3bs3EyZMICIiAj8/P9555x2rplITQUFBNcrNr1692urkHzZsmNv9FyxYwGOPPcaKFStqLFMb8+bNY8KECaxfv95tVsmqdO7cmYSEBMaNG0dZWZmVAvrBBx9k/PjxJCYmWjnuy1m7dq3Vyd+jRw9Wr15d6zlyc3Ot/gxXiouLmT17NgUFBfj5+REeHk5CQgIAM2fOZMSIEYSEhLBnzx7WrFnD5MmTuXz5MgBLlizhtttuAxxDhwcMGEBZWZnbWs73339P//79KSoqokWLFla66Q4dOhAfH8+wYcO4evUqcXFxlpz/q6++yqRJk3jhhRfo27cvM2bMAGDGjBlMmTKF8PBwOnXqxMaNG4Hru1cMBm9gq1y/iOwF5rlrIqtSrjuO2s4dqlrkbCKLBYqANBw1nfy6zmfk+uvG1+T64+Pj6dq1q5Uauj7p3r07aWlpTcZX5l63HyPXX0+IyG4c/SdVeV5VE6/hOAHAFuAZVS1yfrwCWAyo8/0NIK6G/WcCMwG6du3qsf0G32DWrFl2m2AwNFsaLMCo6g2PxRORVjiCywZV/cjl2GdcyvwV2FaLHQlAAjhqMDdqk6/ToUMHu01oMtQ1Us1gqMof//hHu03wKo12mLI4OmhWAl+r6p+rbAtx9tEAPIRj0IChHmjdurXdJhgMPsu9995rtwlexZZRZCLykIhkA/cAH4vIJ87PbxWR8hFhg4ApwH0ikuF8jXRue01EvhKRL4EhgOcTEAy1UlJSYmlfGQyG+uXAgQOVdPZ8HbtGkf0N+Jubz3OAkc7lFMDtMDNVndKgBjZjiooc3VxNpePaYGhKlKtjNLVO/uvFSMUYDAaDoUEwAaYR4E6uf9euXfTr14/IyEj69evHZ599Vq/nDAgIACAnJ4fx48fXWtZV3r8q8+fPp3fv3syfP7/ebMvIyGD79oq5s0lJSZXk7OuLVatWERkZyZ133skdd9xBYqJjcOOaNWvIycmpl3Pk5eUxZMgQAgICqo1oS09PJzIykvDwcObMmWNJvpw/f56YmBh69uxJTEwM+fmOEfiqypw5cwgPD+fOO+/k888/t461du1aevbsSc+ePa3JmAaD7VyPBHNTfTUluf7PP//cko//6quv9NZbb63zOJ7K9dd0TlX3cv21yfu3b9/ekryvL1avXq1PPfVUvR6zKt9++6326NFDCwoKVFX1woULeurUKVV1SKqnpqbWy3mKi4t13759umLFimrXdPfdd+uBAwe0rKxMhw8frtu3b1dV1fnz5+vSpUtVVXXp0qW6YMECVVX9+OOPdfjw4VpWVqYHDx7UqKgoVVXNy8vTsLAwzcvL0/Pnz2tYWJieP3++mi2N4V5v7hi5/uZMI9Lr79u3r7Xcu3dvfvrpJy5fvlxJ0wsck/3i4uJITk5m1qxZXLhwgYSEBEpKSggPD2f9+vW0bduWrKwsHnnkEUpLSxnuYs/p06cZNWoUmZmZnD59mkmTJvHjjz/i5+dHfHx8raNeRo8ezcWLFxkwYACLFi1ix44djBo1yqoRBQQEUFxczN69e3n55ZcJDg4mMzOTfv368f777yMipKam8vTTT3Px4kVuuukmdu3axUsvvcSlS5dISUlh0aJFXLp0ibS0NOLj4/nmm2+Ii4vj3LlzdO7cmdWrV9O1a1diY2Pp0KEDaWlpfP/997z22mu11szOnj1L+/btrZpcQEAAAQEBbN68mbS0NB599FHatGnDwYMHOXbsGM8++yzFxcUEBwezZs0aQkJCiI6Opk+fPhw+fJiioiJWrVpVTaKlXbt2/OY3v7GkZsrJzc2lqKjIkpeZOnUqW7duZcSIESQmJlpt9NOmTSM6OppXX32VxMREpk6diogwcOBACgoKyM3NZe/evcTExFhioTExMezcudMS4jQY7MI0kTUC3Mn1u7Jlyxb69u1bLbiU4+/vT0pKCpMmTWLcuHGkpqZy5MgRbr/9dlauXAnA008/zRNPPEFqaio/+5m7+a9YUi3p6els2rSpzkyXSUlJlhL0xIkTay37xRdfWDIpp06dYv/+/ZSUlDBx4kSWL1/OkSNH2L17N+3ateOVV15h4sSJbo87a9Yspk6dypdffsmjjz5aycbc3FxSUlLYtm0bCxcurNWeu+66iy5duhAWFsb06dP5+9//DjjyrPTv358NGzaQkZGBn58fs2fPZvPmzVZOneeff946zsWLFzlw4AB/+ctfiItzO9fXLd999x2hoaHWuqvE/pkzZyxNuJCQEM6ePWvt406W38j1Nx3eeust3nrrLbvN8BqmBuOKTXr97uT6yzl69CjPPfccycnJNe7v+hDOzMzkhRdeoKCggOLiYktTbP/+/WzZsgVwyPU/99xz1Y5z5coVZs2aRUZGBi1btqwk6nijREVFWQ/U8n6mwMBAQkJCLM00TyZ5Hjx4kI8++si6jgULFljbxo4dS4sWLYiIiODMmTM1HQJw9Hvt3LmT1NRUPv30U+bOnUt6erqV0Kuc48ePk5mZSUxMDABXr161Hv5QIdc/ePBgioqKKCgoICgoqM7r0OuQ2K9pn+s5lsEemotMfzkmwDRisrOzeeihh1i3bh2//OUvayzXrl07azk2NpatW7dy1113sWbNmkrDIet66Lz55psEBwdz+PBhWrVqhb+//zXZ6yrXr6qV5tM0lFy/6/6u53D30HW3b1RUFFFRUcTExDB9+vRqAUZV6d27d42J2qra7+n1hIaGkp2dba27Sux36dKF3NxcQkJCyM3N5ZZbbrH2cSfLHxoaWul7zs7OvuZUCwbvsHv3boBmk3TMNJE1UgoKCnjggQdYunSplezKEy5cuEBISAhXrlypJMs/aNAgS323Nrn+jh07cvHiRdavX8/Vq1evyWZXuf7ExMQ6E5b96le/Iicnh9TUVMv20tLSGuX6wTET2vU6XNM/13aequTk5FQahVWTXH+vXr04d+6cFWCuXLnC0aNHrf3K5fpTUlIIDAy0kqrVRUhICO3bt+fQoUOoKuvWrXMry19Vrn/dunWoKocOHbJqgMOGDSM5OZn8/Hzy8/NJTk6+LjVsQ8OzZMkSlixZYrcZXsPUYBop8fHxnDx5ksWLF7N48WIAkpOTrX+zNbF48WIGDBhAt27diIyMtB6Uy5cv55FHHmH58uWVMmC68uSTTzJ27FiSkpKIiYmpVDPyhMcff5wxY8YQFRXF0KFD69y/devWbNq0idmzZ3Pp0iXatGnD7t27GTJkCMuWLaNPnz4sWrSo0j5vv/02cXFxvP7661Ynf2388MMPbmszV65cYd68eeTk5ODv70/nzp159913AUct8Pe//73Vyb9582bmzJlDYWEhpaWlPPPMM5a0fseOHbn33nutTn53dO/enaKiIkpKSti6dSvJyclERESwYsUKYmNjuXTpEiNGjGDEiBEALFy4kAkTJrBy5Uq6du3Khx9+CMDIkSPZvn074eHhtG3b1rr2Tp068eKLL1pNjS+99FKd2UENBm9gq1y/tzFy/XXja3L927Zt49SpU3UOWLgeoqOj+dOf/kT//tesYm4L5l63HyPXbzD4EKNGjbLbBIOh2WICjMFwnTS1f6EGg7cxAQbqZTSTr+BpJ7WhadGcmsIbM++9957dJniVZh9g/P39ycvL4+abbzZBBmjVqpXdJhjqGVUlLy/vmoedG+qfXr162W2CV2n2AaZ8PsK5c+fsNqVRUD602AQa38Lf37+ScoDBHsoVIx588EGbLfEOzT7AtGrVirCwMLvNaDSUz7+YNm2azZYYDL7HG2+8ATSfAGNXRsuHReSoiJSJSI1D30TktDNzZYaIpLl83klEdonICed7R+9YbjAYDAZPsWsmfyYwDvinB2WHqGqfKmOwFwKfqmpP4FPnusFgMBgaEbYEGFX9WlWP38AhxgDlWZXWAmNv3CqDwWAw1Ce2zuQXkb3APFVNq2F7FpAPKPCeqiY4Py9Q1SCXcvmq6raZTERmAjOdq3fgqD0ZIBj4wW4jGgnGFxUYX1RgfFFBL1Vtf607NVgnv4jsBtwlHnleVRM9PMwgVc0RkVuAXSLyv6rqSbOahTMolQemtOuRO/BFjC8qML6owPiiAuOLClz7wK+FBgswqnrDetSqmuN8PysifwOicPTbnBGREFXNFZEQ4OyNnstgMBgM9UujlesXkXYi0r58GfgtFc1bSUD5ONppgKc1IoPBYDB4CbuGKT8kItnAPcDHIvKJ8/NbRWS7s1gXIEVEjgCHgY9Vdadz2zIgRkROADHOdU9IqLeLaPoYX1RgfFGB8UUFxhcVXJcvmpVcv8FgMBi8R6NtIjMYDAZD08YEGIPBYDA0CD4ZYERkuIgcF5GTIlJtlr+I3CQim5zb/yUi3b1vpXfwwBfPisgxEflSRD4VkW522OkN6vKFS7nxIqK1yRg1dTzxhYhMcN4bR0Xkv71to7fw4DfSVUT2iMgXzt/JSDvsbGhEZJWInBURt3MFxcHbTj99KSK/rvOgqupTL6Al8H9AD6A1cASIqFLmSeBd5/IkYJPddtvoiyFAW+fyE83ZF85y7XEMhT8E9Lfbbhvvi57AF0BH5/otdtttoy8SgCecyxHAabvtbiBfDAZ+DWTWsH0ksAMQYCDwr7qO6Ys1mCjgpKqeUtUSYCMOaRlXXKVmNgNDxTeTwdTpC1Xdo6o/OlcPAb6q6e7JfQGwGHgN+MmbxnkZT3zxOPCOquaDYy6al230Fp74QoEOzuVAIMeL9nkNdUxiP19LkTHAOnVwCAhyzkOsEV8MMD8HvnVZz3Z+5raMqpYChcDNXrHOu3jiC1dm4PiH4ovU6QsR6Qv8QlW3edMwG/DkvrgNuE1E9ovIIREZ7jXrvIsnvngZeMw5tWI7MNs7pjU6rvV54pP5YNzVRKqOxfakjC/g8XWKyGNAf+A/GtQi+6jVFyLSAngTiPWWQTbiyX3hh6OZLBpHrXafiNyhqgUNbJu38cQXk4E1qvqGiNwDrHf6oqzhzWtUXPNz0xdrMNnAL1zWQ6lepbXKiIgfjmpvbVXDpoonvkBE7geeB0ar6mUv2eZt6vJFexxiqHtF5DSONuYkH+3o9/Q3kqiqV1Q1CziOI+D4Gp74YgbwPwCqehDwxyGE2dzw6Hniii8GmFSgp4iEiUhrHJ34SVXKuErNjAc+U2cvlo9Rpy+czULv4QguvtrODnX4QlULVTVYVburancc/VGjtQal7yaOJ7+RrTgGgCAiwTiazE551Urv4Ikv/g0MBRCR23EEmOaYYz0JmOocTTYQKFTV3Np28LkmMlUtFZFZwCc4RoisUtWjIvIKkKaqScBKHNXckzhqLpPss7jh8NAXrwMBwIfOcQ7/VtXRthndQHjoi2aBh774BPitiBwDrgLzVTXPPqsbBg998QfgryIyF0eTUKwv/iEVkQ9wNIkGO/ub/gtoBaCq7+LofxoJnAR+BKbXeUwf9JPBYDAYGgG+2ERmMBgMhkaACTAGg8FgaBBMgDEYDAZDg2ACjMFgMBgaBBNgDAaDwdAgmABjMBgMhgbBBBiDwWAwNAj/D6PBsaT1eqpPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot radial functions vs analytic solution\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "fig, ax1 = plt.subplots(1, sharex=True)\n",
    "\n",
    "ax1.plot(x_vals, [2./3. * r_**2 for r_ in x_vals], \"b:\", lw=3, label=\"$2/3 r^2$\")\n",
    "ax1.plot(x_vals, [-r_**2 for r_ in x_vals], \"r:\", lw=3, label=\"$-r^2$\")\n",
    "ax1.plot([np.sqrt(3 * (rbf_high/4)**2), np.sqrt(3 * (rbf_high/4)**2)], [-4,4], '--', \n",
    "         color=\"#000000\", label=\"Maximum distance\")\n",
    "ax1.plot([0.22, 0.22], [-4,4], '--', color=\"#818181\", label=\"Mean minimum distance\")\n",
    "\n",
    "F0_y_vals_reshape = np.reshape(np.array(F0_y_vals), np.array(F0_y_vals).shape[0])\n",
    "ax1.plot(x_vals, F0_y_vals_reshape, 'b-', label=\"F0 radial function, Step {}\".format(step), alpha=0.5)\n",
    "\n",
    "F2_y_vals_reshape = np.reshape(np.array(F2_y_vals), np.array(F2_y_vals).shape[0])\n",
    "ax1.plot(x_vals, F2_y_vals_reshape,'r-', label=\"F2 radial function, Step {}\".format(step), alpha=0.5)\n",
    "    \n",
    "handles, labels = ax1.get_legend_handles_labels()\n",
    "ax1.legend(handles, labels)\n",
    "\n",
    "ax1.set_xlim(0,1)\n",
    "ax1.set_ylim(-1.5,0.75)\n",
    "ax1.set_ylabel(\"Radial function output\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}