tests/test_rbp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9435e6ca",
   "metadata": {},
   "source": [
    "# Test radiative backpropagation (RBP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fdacec0-a831-4ca4-8ba3-89e574199421",
   "metadata": {},
   "source": [
    "Testing is done by comparing the gradients computed by RBP to the ones computed by Dr.JiT using conventional automatic differentiation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2e87e056",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "gpu_num = 0 # Use \"\" to use the CPU\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = f\"{gpu_num}\"\n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6da01ce7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import mitsuba as mi\n",
    "import drjit as dr\n",
    "\n",
    "from instant_rm import MapTracer, PathlossMapRBPTracer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53f900fd",
   "metadata": {},
   "source": [
    "# Simple Reflector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "eef5996e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RBP Grad:  [ 0.16182931 -0.04538707  0.13174717 -0.04447149]\n",
      "Ref. Grad:  [ 0.1618286  -0.04538754  0.13174748 -0.04447149]\n",
      "Relative Error:  4.104959e-06\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f9960257210>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhwAAAGdCAYAAABZ+qqcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJrUlEQVR4nO3dfVxUZd4H/s+MCGMKgy7CMEVKrc8PUFgTrpauswL2MmnbVog71Ejv2x+0umT5kAKld5SmpelKekfUL01zf0Xd6k0hhmxJmKibuuqqi4obgynCAC2P5/z+cDk5MjwMcw4zcj7vXudlc851rvnO4QBfrus616URRVEEERERkYK0rg6AiIiIej4mHERERKQ4JhxERESkOCYcREREpDgmHERERKQ4JhxERESkOCYcREREpDgmHERERKQ4D1cHIAdBEPDDDz/A29sbGo3G1eEQEZGDRFFEdXU1jEYjtFrl/hauq6tDQ0OD0/V4enpCp9PJEJF69IiE44cffkBQUJCrwyAiIieVlpbirrvuUqTuuro6BA/qB8uVZqfrMhgMKCkpYdLhgB6RcHh7ewMAJmAaPNDbxdG4N613P1nqEaprZKnHXXkYAlwdQitCba2rQ2hFrvtArvsS6Nn3pjteJ7liahIbUVDzsfTzXAkNDQ2wXGlGSfEg+Hh3vRXFWi0gOOwiGhoamHA4oEckHC3dKB7oDQ8NE472aDWestQj9PDr7KGV5zrJSdA0ujqEVuS6D+S6L4GefW+643WSMyYA3dIt7uOtdSrhoK7pEQkHERFRZzWLApqdWLa0WRTkC0ZFmHAQEZGqCBAhoOsZhzPnqplibUqbNm3C4MGDodPpYDKZcOjQoXbL79q1C8OHD4dOp8OYMWOwd+9epUIjIiIVE2T4jxynSMKxc+dOJCcnIzU1FUeOHEFISAgiIiJw5coVu+UPHjyI2NhYJCQk4OjRo4iOjkZ0dDROnDihRHhERETdTu4/xEVRREpKCgIDA9GnTx+YzWacPXvWpkxFRQXi4uLg4+MDX19fJCQkoKbm5wHDdXV1mD17NsaMGQMPDw9ER0fbjSU/Px/3338/vLy88Mtf/hJZWVkOf35FEo5169Zh7ty5mDNnDkaOHImMjAzccccdyMzMtFt+/fr1iIyMxAsvvIARI0Zg5cqVuP/++7Fx40YlwiMiIhVrFkWnN0cp8Yf46tWrsWHDBmRkZKCoqAh9+/ZFREQE6urqpDJxcXE4efIkcnNzsXv3bhQUFGDevHk/X4vmZvTp0wd/+MMfYDab7cZSUlKCRx99FJMnT8axY8ewcOFCPPvss/jiiy8cugYaUezClWtHQ0MD7rjjDvz5z3+2yZRmzZqFyspKfPbZZ63Oufvuu5GcnIyFCxdK+1JTU5GdnY2//vWvrcrX19ejvr5eem21WhEUFIRJmMGnVDqglemRM6G6WpZ63JVHoMHVIbQi1LjjY7Hy3Ady3ZdAz7433fE6yRVTk9iA/dXbUFVVBR8fH1nqvJXVaoVer8fF00anH4sdNPwHh2I1mUx44IEHpD+kBUFAUFAQnnvuOSxZsqRV+ZkzZ6K2tha7d++W9j300EMIDQ1FRkYGRFGE0WjE888/j0WLFgEAqqqqEBAQgKysLMTExODUqVMYOXIkvvvuO4wbNw4AkJOTg2nTpuHy5cswGo027zl79mxUVlYiOzvbZv/ixYuxZ88em2QnJiYGlZWVyMnJ6dTnBxRo4bh69Sqam5sREGA7j0FAQAAsFovdcywWi0Pl09PTodfrpY2TfhERUXezWq02281/CN+soaEBxcXFNi0IWq0WZrMZhYWFds8pLCxs1eIQEREhlS8pKYHFYrEpo9frYTKZpDKFhYXw9fWVkg0AMJvN0Gq1KCoq6vTn7CiWzrotH0ReunQpqqqqpK20tNTVIRER0W1CgIhmJ7aWp1SCgoJs/vhNT0+3+35K/CHe8m9HZfz9/W2Oe3h4YMCAAW2+ryOxWK1W/Otf/+p0PbI/Fuvn54devXqhvLzcZn95eTkMBvvN1AaDwaHyXl5e8PLykidgIiJSFbkeiy0tLbXpUuHvpfbJ3sLh6emJsLAw5OXlSfsEQUBeXh7Cw8PtnhMeHm5THgByc3PbLE9ERORqPj4+NltbCYcSf4i3/NtRmVsHpTY1NaGioqLN93UkFh8fH/Tp06fT9SjSpZKcnIytW7fi/fffx6lTpzB//nzU1tZizpw5AID4+HgsXbpUKr9gwQLk5ORg7dq1OH36NNLS0nD48GEkJSUpER4REalYdz+losQf4sHBwTAYDDZlrFYrioqKpDLh4eGorKxEcXGxVGb//v0QBAEmk6nT8cvVKKDITKMzZ87Ejz/+iJSUFFgsFoSGhiInJ0fqA7p06ZLN8sPjx4/H9u3bsXz5cixbtgxDhgxBdnY2Ro8erUR4RESkYsK/N2fOd1RycjJmzZqFcePG4cEHH8Rbb73V6g/xO++8UxoHsmDBAjzyyCNYu3YtHn30UezYsQOHDx/Gli1bANxYc2bhwoVYtWoVhgwZguDgYKxYsQJGo1F6QnTEiBGIjIzE3LlzkZGRgcbGRiQlJSEmJsbmCZW//e1vaGhoQEVFBaqrq3Hs2DEAQGhoKADgv/7rv7Bx40a8+OKLeOaZZ7B//358/PHH2LNnj0PXQPbHYl2h5VEnPhbbMT4W2zl8LLZz3O3RSqBn35vueJ1ux8diT58KgLcTj8VWVwsYPqLc4Vg3btyINWvWSH+Ib9iwQWppmDRpEgYPHmwzodauXbuwfPlyXLhwAUOGDMHq1asxbdo06bgoikhNTcWWLVtQWVmJCRMm4E9/+hOGDh0qlamoqEBSUhL+93//F1qtFk888QQ2bNiAfv1+XuV38ODBuHjxYqt4b04P8vPz8cc//hF/+9vfcNddd2HFihWYPXt2pz87wIRDdZhwdA4Tjs5xt19aQM++N93xOt2OCcfJU/5OJxyjRlxRNNaeiIu3qYxcP2Tk/IXcVNb5x7PaI+cPY7nI9dkA90wW5boP5LxO7sgtE1iZ7gPZ6hEbZamnM5pFOLlarHyxqAkTDiIiUhVXjOGg23TiLyIiIrq9sIWDiIhURYAGzdA4dT45jgkHERGpiiDe2Jw5nxzHLhUiIiJSHFs4iIhIVZqd7FJx5lw1Y8JBRESqwoTDNdilQkRERIpjCwcREamKIGogiE48peLEuWrGhIOIiFSFXSquwS4VIiIiUhxbOIiISFWaoUWzE39vN8sYi5ow4SAiIlURnRzDIXIMR5cw4SAiIlXhGA7X4BgOIiIiUhxbOIiISFWaRS2aRSfGcHAtlS5hwkFERKoiQAPBiQZ+Acw4uoJdKkRERKQ4tnBQlzSVWVwdQitCdbVb1iUXd4zJHe+DXkPvlaUeTXWtLPUAgFAjT13ueA/cjjho1DWYcBARkao4P4aDXSpdwS4VIiIiUhxbOIiISFVuDBp1YvE2dql0CRMOIiJSFcHJqc35lErXsEuFiIiIFMcWDiIiUhUOGnUNJhxERKQqArSc+MsFmHAQEZGqNIsaNDux4qsz56oZx3AQERGR4tjCQUREqtLs5FMqzexS6RImHEREpCqCqIXgxKBRgYNGu4RdKkRERKQ4tnAQEZGqsEvFNZhwEBGRqghw7kkTQb5QVEX2LpX09HQ88MAD8Pb2hr+/P6Kjo3HmzJl2z8nKyoJGo7HZdDqd3KERERGRi8iecBw4cACJiYn49ttvkZubi8bGRkydOhW1tbXtnufj44OysjJpu3jxotyhERERSRN/ObOR42TvUsnJybF5nZWVBX9/fxQXF+Phhx9u8zyNRgODwSB3OERERDacn9qcCUdXKD6Go6qqCgAwYMCAdsvV1NRg0KBBEAQB999/P1599VWMGjXKbtn6+nrU19dLr61Wq3wBU6dovb1lq0uorpalHneMibqfprr91lRX0PbrK0s9vC/pdqZomiYIAhYuXIhf/epXGD16dJvlhg0bhszMTHz22Wf48MMPIQgCxo8fj8uXL9stn56eDr1eL21BQUFKfQQiIuphBGic3shxiiYciYmJOHHiBHbs2NFuufDwcMTHxyM0NBSPPPIIPvnkEwwcOBDvvPOO3fJLly5FVVWVtJWWlioRPhER9UAtXSrObOQ4xbpUkpKSsHv3bhQUFOCuu+5y6NzevXvjvvvuw7lz5+we9/LygpeXlxxhEhGRyjg/DwcTjq6Q/aqJooikpCR8+umn2L9/P4KDgx2uo7m5GcePH0dgYKDc4REREZELyN7CkZiYiO3bt+Ozzz6Dt7c3LBYLAECv16NPnz4AgPj4eNx5551IT08HALzyyit46KGH8Mtf/hKVlZVYs2YNLl68iGeffVbu8IiISOUEUQPBmYm/uDx9l8iecGzevBkAMGnSJJv97733HmbPng0AuHTpErTanxtXrl+/jrlz58JisaB///4ICwvDwYMHMXLkSLnDIyIilROc7FLhPBxdI3vCIXZiFb38/Hyb12+++SbefPNNuUMhIiIiN8G1VIiISFWcX56eLRxdwYSDiIhUpRkaNDsxl4Yz56oZ0zQiIiJSHFs4iIhIVdil4hpMOIiISFWa4Vy3SLN8oagK0zQiIiJSHFs4iIhIVdil4hpMOIiISFWcXYCNi7d1DRMOIiJSFdHJJeZFPhbbJUzTiIiISHFs4SAiIlVhl4prMOEgl9N6e8tSj1BdLUs91Hlyfe3k1FRmcXUIrch1ndzxet+O33dcLdY1mKYRERGR4tjCQUREqtLs5PL0zpyrZkw4iIhIVdil4hpM04iIiEhxTDiIiEhVBGid3rpi06ZNGDx4MHQ6HUwmEw4dOtRu+V27dmH48OHQ6XQYM2YM9u7da3NcFEWkpKQgMDAQffr0gdlsxtmzZ23KVFRUIC4uDj4+PvD19UVCQgJqampsynz//feYOHEidDodgoKCsHr16laxvPXWWxg2bBj69OmDoKAg/PGPf0RdXZ1Dn58JBxERqUqzqHF6c9TOnTuRnJyM1NRUHDlyBCEhIYiIiMCVK1fslj948CBiY2ORkJCAo0ePIjo6GtHR0Thx4oRUZvXq1diwYQMyMjJQVFSEvn37IiIiwiYRiIuLw8mTJ5Gbm4vdu3ejoKAA8+bNk45brVZMnToVgwYNQnFxMdasWYO0tDRs2bJFKrN9+3YsWbIEqampOHXqFN59913s3LkTy5Ytc+gaaERRFB06ww1ZrVbo9XpMwgx4aHq7OhxV4ON5BPA+6Cx3vE5yket6N4mNyMdnqKqqgo+Pjyx13qrld8X8v/wWXv26/ruivqYRmyd+4lCsJpMJDzzwADZu3AgAEAQBQUFBeO6557BkyZJW5WfOnIna2lrs3r1b2vfQQw8hNDQUGRkZEEURRqMRzz//PBYtWgQAqKqqQkBAALKyshATE4NTp05h5MiR+O677zBu3DgAQE5ODqZNm4bLly/DaDRi8+bNeOmll2CxWODp6QkAWLJkCbKzs3H69GkAQFJSEk6dOoW8vDwplueffx5FRUX4+uuvO33d2MJBRESq0jJo1JkNuJHA3LzV19fbfb+GhgYUFxfDbDZL+7RaLcxmMwoLC+2eU1hYaFMeACIiIqTyJSUlsFgsNmX0ej1MJpNUprCwEL6+vlKyAQBmsxlarRZFRUVSmYcfflhKNlre58yZM7h+/ToAYPz48SguLpa6gP7xj39g7969mDZtWieu9s+YcBARkaqI/14ttqub+O+ZRoOCgqDX66UtPT3d7vtdvXoVzc3NCAgIsNkfEBAAi8X+RHUWi6Xd8i3/dlTG39/f5riHhwcGDBhgU8ZeHTe/x1NPPYVXXnkFEyZMQO/evXHvvfdi0qRJDnep8LFYIiJSlWZo0OzEAmwt55aWltp0qXh5eTkdmzvKz8/Hq6++ij/96U8wmUw4d+4cFixYgJUrV2LFihWdrocJBxERURf4+Ph0agyHn58fevXqhfLycpv95eXlMBgMds8xGAztlm/5t7y8HIGBgTZlQkNDpTK3DkptampCRUWFTT323ufm91ixYgWefvppPPvsswCAMWPGoLa2FvPmzcNLL70ErbZznSXsUiEiIlURRGfHcTj2fp6enggLC7MZdCkIAvLy8hAeHm73nPDwcJvyAJCbmyuVDw4OhsFgsCljtVpRVFQklQkPD0dlZSWKi4ulMvv374cgCDCZTFKZgoICNDY22rzPsGHD0L9/fwDATz/91Cqp6NWrF4Abj+Z2FhMOIiJSFWfGb7RsjkpOTsbWrVvx/vvv49SpU5g/fz5qa2sxZ84cAEB8fDyWLl0qlV+wYAFycnKwdu1anD59GmlpaTh8+DCSkpIAABqNBgsXLsSqVavw+eef4/jx44iPj4fRaER0dDQAYMSIEYiMjMTcuXNx6NAhfPPNN0hKSkJMTAyMRiOAG+MzPD09kZCQgJMnT2Lnzp1Yv349kpOTpVimT5+OzZs3Y8eOHSgpKUFubi5WrFiB6dOnS4lHZ7BLhYiISGEzZ87Ejz/+iJSUFFgsFoSGhiInJ0caoHnp0iWbVoTx48dj+/btWL58OZYtW4YhQ4YgOzsbo0ePlsq8+OKLUtdGZWUlJkyYgJycHOh0OqnMtm3bkJSUhClTpkCr1eKJJ57Ahg0bpON6vR5ffvklEhMTERYWBj8/P6SkpNjM1bF8+XJoNBosX74c//znPzFw4EBMnz4d//3f/+3QNeA8HNQl7jivgDvOv9DT8T7oHHe8TnK5HefhePqrWHj28+z4hDY01DTg/538kaKx9kRs4SAiIlXp6myhN59PjuMYDiIiIlIcWzhuIWfTpzs27VLnyHUfaPv1laUeOTWV2Z9oqCvk+nxCTa0s9fR0cv5Mkese9wi0/1inw4QGQL5bs/236uLAz5vPJ8cx4SAiIlUR8PP05F09nxzHNI2IiIgUxxYOIiJSFREap1opRLZwdAkTDiIiUpWbV3zt6vnkOCYcRESkKhw06hqyX7W0tDRoNBqbbfjw4e2es2vXLgwfPhw6nQ5jxozB3r175Q6LiIiIXEiRNG3UqFEoKyuTtq+//rrNsgcPHkRsbCwSEhJw9OhRREdHIzo6GidOnFAiNCIiUjnnFm5zrjtGzRRJODw8PGAwGKTNz8+vzbLr169HZGQkXnjhBYwYMQIrV67E/fffj40bNyoRGhERqZzw70GjzmzkOEUSjrNnz8JoNOKee+5BXFwcLl261GbZwsJCmM1mm30REREoLCxs85z6+npYrVabjYiIiNyX7AmHyWRCVlYWcnJysHnzZpSUlGDixImobmOGPIvFIq2W1yIgIAAWS9tTzqWnp0Ov10tbUFCQrJ+BiIh6LnapuIbsCUdUVBSefPJJjB07FhEREdi7dy8qKyvx8ccfy/YeS5cuRVVVlbSVlpbKVjcREfVsTDhcQ/HHYn19fTF06FCcO3fO7nGDwYDy8nKbfeXl5TAY2p6f38vLC15eXrLGSURERMpR/GHimpoanD9/HoGBgXaPh4eHIy8vz2Zfbm4uwsPDlQ6NiIhUiC0criF7wrFo0SIcOHAAFy5cwMGDB/H444+jV69eiI2NBQDEx8dj6dKlUvkFCxYgJycHa9euxenTp5GWlobDhw8jKSlJ7tCIiIiYcLiI7F0qly9fRmxsLK5du4aBAwdiwoQJ+PbbbzFw4EAAwKVLl6DV/pznjB8/Htu3b8fy5cuxbNkyDBkyBNnZ2Rg9erTcoREREZGLyJ5w7Nixo93j+fn5rfY9+eSTePLJJ+UOhYiIqBURzi0xL8oXiqpwLRUiIlIVLt7mGkw4iIhIVZhwuAYTjlsIbUxQ1hVab2956unXV5Z6AECoqZWtLrnIec3l4o4xyUWu+xIAmsranqCPfibX/STn16566khZ6tEfvChLPdTzMeEgIiJVYQuHazDhICIiVWHC4RqKT/xFRERExBYOIiJSFVHUQHSilcKZc9WMCQcREamKAI1T83A4c66asUuFiIiIFMcWDiIiUhUOGnUNJhxERKQqHMPhGuxSISIiIsWxhYOIiFSFXSquwYSDiIhUhV0qrsGEg4iIVEV0soWDCUfXcAwHERERKY4tHEREpCoiAFF07nxyHBMOIiJSFQEaaDjTaLdjlwoREREpji0cRESkKnxKxTV6VMLhYQiAh9bT1WFIhJpaWeppKrPIUg8AaL29ZalHqK6WpR5Avpi0/frKUg8g39dOzuvUk8l1D8jJHe9xOWP61y/kaeDuK9PPpyaxUZZ6OkMQNdBwHo5uxy4VIiIiUlyPauEgIiLqiCg6+ZQKH1PpEiYcRESkKhzD4RrsUiEiIiLFsYWDiIhUhS0crsGEg4iIVIVPqbgGEw4iIlIVDhp1DY7hICIiIsWxhYOIiFTlRguHM2M4ZAxGRZhwEBGRqnDQqGuwS4WIiIgUxxYOIiJSFfHfmzPnk+OYcBARkaqwS8U12KVCREREipM94Rg8eDA0Gk2rLTEx0W75rKysVmV1Op3cYREREd0gyrCRw2TvUvnuu+/Q3NwsvT5x4gR+85vf4Mknn2zzHB8fH5w5c0Z6rdGwuYqIiBTiZJcK2KXSJbInHAMHDrR5/dprr+Hee+/FI4880uY5Go0GBoNB7lCIiIha4UyjrqHoGI6GhgZ8+OGHeOaZZ9pttaipqcGgQYMQFBSEGTNm4OTJk0qGRURERN1M0adUsrOzUVlZidmzZ7dZZtiwYcjMzMTYsWNRVVWFN954A+PHj8fJkydx11132T2nvr4e9fX10mur1QoAaLKUA5reTsXsEeh+LS1ab2/Z6hKqq2WrSy7uGBPdvjSB/rLU49Gvryz1AIDoLU9dvWT6bABQ9UidLPUY/lemn5lCA2CRp6qO8CkV11C0hePdd99FVFQUjEZjm2XCw8MRHx+P0NBQPPLII/jkk08wcOBAvPPOO22ek56eDr1eL21BQUFKhE9ERD2RqHF+I4cplnBcvHgR+/btw7PPPuvQeb1798Z9992Hc+fOtVlm6dKlqKqqkrbS0lJnwyUiIiIFKZZwvPfee/D398ejjz7q0HnNzc04fvw4AgMD2yzj5eUFHx8fm42IiKgzWgaNOrN1xaZNmzB48GDodDqYTCYcOnSo3fK7du3C8OHDodPpMGbMGOzdu/eWzyEiJSUFgYGB6NOnD8xmM86ePWtTpqKiAnFxcfDx8YGvry8SEhJQU1NjU+b777/HxIkTodPpEBQUhNWrV7eKpbKyEomJiQgMDISXlxeGDh3aKp6OKJJwCIKA9957D7NmzYKHh+0wkfj4eCxdulR6/corr+DLL7/EP/7xDxw5cgT/8R//gYsXLzrcMkJERNQpLpiHY+fOnUhOTkZqaiqOHDmCkJAQRERE4MqVK3bLHzx4ELGxsUhISMDRo0cRHR2N6OhonDhxQiqzevVqbNiwARkZGSgqKkLfvn0RERGBurqfx+fExcXh5MmTyM3Nxe7du1FQUIB58+ZJx61WK6ZOnYpBgwahuLgYa9asQVpaGrZs2SKVaWhowG9+8xtcuHABf/7zn3HmzBls3boVd955p0PXQJFBo/v27cOlS5fwzDPPtDp26dIlaLU/5znXr1/H3LlzYbFY0L9/f4SFheHgwYMYOXKkEqERERF1u3Xr1mHu3LmYM2cOACAjIwN79uxBZmYmlixZ0qr8+vXrERkZiRdeeAEAsHLlSuTm5mLjxo3IyMiAKIp46623sHz5csyYMQMA8MEHHyAgIADZ2dmIiYnBqVOnkJOTg++++w7jxo0DALz99tuYNm0a3njjDRiNRmzbtg0NDQ3IzMyEp6cnRo0ahWPHjmHdunVSYpKZmYmKigocPHgQvXvfeDBj8ODBDl8DRVo4pk6dClEUMXTo0FbH8vPzkZWVJb1+8803cfHiRdTX18NisWDPnj247777lAiLiIhIekrFmQ240Tpw83bz05M3a2hoQHFxMcxms7RPq9XCbDajsLDQ7jmFhYU25QEgIiJCKl9SUgKLxWJTRq/Xw2QySWUKCwvh6+srJRsAYDabodVqUVRUJJV5+OGH4enpafM+Z86cwfXr1wEAn3/+OcLDw5GYmIiAgACMHj0ar776qs0kn53BtVSIiEh9ZOhOCQoKsnliMj093e5bXb16Fc3NzQgICLDZHxAQAIvF/rPAFoul3fIt/3ZUxt/f9lFqDw8PDBgwwKaMvTpufo9//OMf+POf/4zm5mbs3bsXK1aswNq1a7Fq1Sq7sbeFq8USERF1QWlpqc1DC15eXi6MRjmCIMDf3x9btmxBr169EBYWhn/+859Ys2YNUlNTO10PEw4iIlIVuSb+6uxTkn5+fujVqxfKy8tt9peXl7e5rIfBYGi3fMu/5eXlNk91lpeXIzQ0VCpz66DUpqYmVFRU2NRj731ufo/AwED07t0bvXr1ksqMGDECFosFDQ0NNt0x7WGXChERqUs3P6Xi6emJsLAw5OXlSfsEQUBeXh7Cw8PtnhMeHm5THgByc3Ol8sHBwTAYDDZlrFYrioqKpDLh4eGorKxEcXGxVGb//v0QBAEmk0kqU1BQgMbGRpv3GTZsGPr37w8A+NWvfoVz585BEASpzN///ncEBgZ2OtkAmHAQEZHqaGTYHJOcnIytW7fi/fffx6lTpzB//nzU1tZKT63cOmXEggULkJOTg7Vr1+L06dNIS0vD4cOHkZSUdOMTaDRYuHAhVq1ahc8//xzHjx9HfHw8jEYjoqOjAdxohYiMjMTcuXNx6NAhfPPNN0hKSkJMTIw0A/hTTz0FT09PJCQk4OTJk9i5cyfWr1+P5ORkKZb58+ejoqICCxYswN///nfs2bMHr776KhITEx26BuxSISIiUtjMmTPx448/IiUlBRaLBaGhocjJyZEGaN46ZcT48eOxfft2LF++HMuWLcOQIUOQnZ2N0aNHS2VefPFF1NbWYt68eaisrMSECROQk5MDnU4nldm2bRuSkpIwZcoUaLVaPPHEE9iwYYN0XK/X48svv0RiYiLCwsLg5+eHlJQUm7k6goKC8MUXX+CPf/wjxo4dizvvvBMLFizA4sWLHboGGlG8/RfatVqt0Ov1mIQZ8HCjxduEmlrZ6pJLT14oTc5F7uTSk6+3nOT82sm1eJumWr7vX7kWb5PT31PlmaF52KIyWeppEhqwz7IFVVVVis0e3fK7ImhzGrR9dB2f0AbhX3UonZ+maKw9EVs4iIhIXbo4W6jN+eQwjuEgIiIixbGFg4iI1MXZJea5PH2XMOEgIiJVcWbF15bzyXFMOG7RVGZ/mllXcsfBkHLGJNfAyp4+QFOuAc1yDmZ2y6+dG94H2hr3+x4+N/n/k6WeiLJQWeppEhs7LkS3NSYcRESkLhw06hJMOIiISF04hsMl+JQKERERKY4tHEREpCoa8cbmzPnkOCYcRESkLhzD4RJMOIiISF04hsMlOIaDiIiIFMcWDiIiUhd2qbgEEw4iIlIXJhwuwS4VIiIiUhxbOIiISF3YwuESTDiIiEhd+JSKS7BLhYiIiBTHFg4iIlIVzjTqGkw4iIhIXTiGwyXYpUJERESKY8JBREREimOXChERqYoGTo7hkC0SdelRCYfWux+0Gk+n6hCqq2WKBtB6e8tWl1zcMSa5uONnk/N+aiqzyFaXXOS65j39+04u2n59XR1CKx6BBnkqEhqA7rrF+VisS7BLhYiIiBTXo1o4iIiIOsSnVFyCCQcREakLEw6XYJcKERERKY4tHEREpCqcadQ1HG7hKCgowPTp02E0GqHRaJCdnW1zXBRFpKSkIDAwEH369IHZbMbZs2c7rHfTpk0YPHgwdDodTCYTDh065GhoREREHRNl2MhhDicctbW1CAkJwaZNm+weX716NTZs2ICMjAwUFRWhb9++iIiIQF1dXZt17ty5E8nJyUhNTcWRI0cQEhKCiIgIXLlyxdHwiIiIyA05nHBERUVh1apVePzxx1sdE0URb731FpYvX44ZM2Zg7Nix+OCDD/DDDz+0agm52bp16zB37lzMmTMHI0eOREZGBu644w5kZmY6Gh4REVH72MLhErIOGi0pKYHFYoHZbJb26fV6mEwmFBYW2j2noaEBxcXFNudotVqYzeY2z6mvr4fVarXZiIiIOqNlDIczGzlO1oTDYrkxTVxAQIDN/oCAAOnYra5evYrm5maHzklPT4der5e2oKAgGaInIiIipdyWj8UuXboUVVVV0lZaWurqkIiI6HbRMrW5Mxs5TNbHYg2GG3Pql5eXIzAwUNpfXl6O0NBQu+f4+fmhV69eKC8vt9lfXl4u1XcrLy8veHl5yRM0ERGpCyf+cglZWziCg4NhMBiQl5cn7bNarSgqKkJ4eLjdczw9PREWFmZzjiAIyMvLa/McIiKiruIYDtdwuIWjpqYG586dk16XlJTg2LFjGDBgAO6++24sXLgQq1atwpAhQxAcHIwVK1bAaDQiOjpaOmfKlCl4/PHHkZSUBABITk7GrFmzMG7cODz44IN46623UFtbizlz5jj/CYmIiMjlHE44Dh8+jMmTJ0uvk5OTAQCzZs1CVlYWXnzxRdTW1mLevHmorKzEhAkTkJOTA51OJ51z/vx5XL16VXo9c+ZM/Pjjj0hJSYHFYkFoaChycnJaDSQlIiJyGrtUXMLhhGPSpEkQxbavtkajwSuvvIJXXnmlzTIXLlxotS8pKUlq8SAiIlKMs90iTDi65LZ8SoWIiIhuLz1q8TahugaCprdTdWi9vWWKBhCqq2WpxyPQ/tM6XdFUZn9uk55Aruvd07njPd7TaQL9ZalHzj+sI08/Kks9Ytk/ZamnSWyUpZ5OYZeKS/SohIOIiKhDTDhcgl0qREREpDi2cBARkao4O5cG5+HoGrZwEBERkeKYcBAREZHi2KVCRETqwkGjLsGEg4iIVIVjOFyDCQcREakPk4ZuxzEcREREpDi2cBARkbpwDIdLMOEgIiJV4RgO12CXChERESmOLRxERKQu7FJxCSYcRESkKuxScQ12qRAREZHi2MJBRETqwi4Vl2DCQURE6sKEwyXYpUJERESKYwvHLYTqaleH0EpTmcXVIShK6+0tSz1yfu3kiskdaQL9ZatLrr9YtP36ylSTm36/yHRvynlfVmaOkqWe/t5WWerRig1AN/345aBR12ALBxERqYsow9YFmzZtwuDBg6HT6WAymXDo0KF2y+/atQvDhw+HTqfDmDFjsHfvXtuPIYpISUlBYGAg+vTpA7PZjLNnz9qUqaioQFxcHHx8fODr64uEhATU1NTYlPn+++8xceJE6HQ6BAUFYfXq1W3GtGPHDmg0GkRHRzv24cGEg4iI1MYFCcfOnTuRnJyM1NRUHDlyBCEhIYiIiMCVK1fslj948CBiY2ORkJCAo0ePIjo6GtHR0Thx4oRUZvXq1diwYQMyMjJQVFSEvn37IiIiAnV1dVKZuLg4nDx5Erm5udi9ezcKCgowb9486bjVasXUqVMxaNAgFBcXY82aNUhLS8OWLVtaxXThwgUsWrQIEydOdPwCgAkHERGR4tatW4e5c+dizpw5GDlyJDIyMnDHHXcgMzPTbvn169cjMjISL7zwAkaMGIGVK1fi/vvvx8aNGwHcaN146623sHz5csyYMQNjx47FBx98gB9++AHZ2dkAgFOnTiEnJwf/8z//A5PJhAkTJuDtt9/Gjh078MMPPwAAtm3bhoaGBmRmZmLUqFGIiYnBH/7wB6xbt84mnubmZsTFxeHll1/GPffc06VrwISDiIhUpWUMhzMbcKN14Oatvr7e7vs1NDSguLgYZrNZ2qfVamE2m1FYWGj3nMLCQpvyABARESGVLykpgcVisSmj1+thMpmkMoWFhfD19cW4ceOkMmazGVqtFkVFRVKZhx9+GJ6enjbvc+bMGVy/fl3a98orr8Df3x8JCQkdXt+2MOEgIiJ1kalLJSgoCHq9XtrS09Ptvt3Vq1fR3NyMgIAAm/0BAQGwWOwPcrZYLO2Wb/m3ozL+/raDxD08PDBgwACbMvbquPk9vv76a7z77rvYunWr3Vg7i0+pEBERdUFpaSl8fHyk115eXi6MRhnV1dV4+umnsXXrVvj5+TlVFxMOIiJSFbkei/Xx8bFJONri5+eHXr16oby83GZ/eXk5DAaD3XMMBkO75Vv+LS8vR2BgoE2Z0NBQqcytg1KbmppQUVFhU4+992k5dv78eVy4cAHTp0+XjguCAOBGa8mZM2dw7733dngNAHapEBGR2nTzUyqenp4ICwtDXl6etE8QBOTl5SE8PNzuOeHh4TblASA3N1cqHxwcDIPBYFPGarWiqKhIKhMeHo7KykoUFxdLZfbv3w9BEGAymaQyBQUFaGxstHmfYcOGoX///hg+fDiOHz+OY8eOSdtjjz2GyZMn49ixYwgKCur0dWALBxERkcKSk5Mxa9YsjBs3Dg8++CDeeust1NbWYs6cOQCA+Ph43HnnndI4kAULFuCRRx7B2rVr8eijj2LHjh04fPiw9LiqRqPBwoULsWrVKgwZMgTBwcFYsWIFjEajNEfGiBEjEBkZiblz5yIjIwONjY1ISkpCTEwMjEYjAOCpp57Cyy+/jISEBCxevBgnTpzA+vXr8eabbwIAdDodRo8ebfNZfH19AaDV/o4w4SAiInVxwVoqM2fOxI8//oiUlBRYLBaEhoYiJydHGqB56dIlaLU/dzqMHz8e27dvx/Lly7Fs2TIMGTIE2dnZNr/kX3zxRdTW1mLevHmorKzEhAkTkJOTA51OJ5XZtm0bkpKSMGXKFGi1WjzxxBPYsGGDdFyv1+PLL79EYmIiwsLC4Ofnh5SUFJu5OuSiEUXxtp+k1Wq1Qq/XYxJmwEPT29XhkIM4tXn3knNqc7HM/qRFjurxU5vLRM778voMmaY2/+ykLPU0iQ3YX70NVVVVnRoX0RUtvytG/j+vopeXruMT2tBcX4e//WmZorH2RBzDQURERIpjlwoREakLl6d3CSYcRESkKlwt1jUc7lIpKCjA9OnTYTQaodFopDnbAaCxsRGLFy/GmDFj0LdvXxiNRsTHx0tztrclLS0NGo3GZhs+fLjDH4aIiKhDLlotVu0cTjhqa2sREhKCTZs2tTr2008/4ciRI1ixYgWOHDmCTz75BGfOnMFjjz3WYb2jRo1CWVmZtH399deOhkZERERuyuEulaioKERFRdk9ptfrkZuba7Nv48aNePDBB3Hp0iXcfffdbQfi4dHmjGtERESyYitFt1N8DEdVVRU0Go00UUhbzp49C6PRCJ1Oh/DwcKSnp7eZoNTX19usyme1WuUMWTbu+LinO5Lr83kEypewuuOjlXJ9PkGmR1kB+b52cj4WK9d1kvMecMefBb4ne/bPlfZwDIdrKPpYbF1dHRYvXozY2Nh2n1U2mUzIyspCTk4ONm/ejJKSEkycOBHVbXxzpaen26zQ58jUqkRERNT9FEs4Ghsb8fvf/x6iKGLz5s3tlo2KisKTTz6JsWPHIiIiAnv37kVlZSU+/vhju+WXLl2KqqoqaSstLVXiIxARUU/EQaMuoUiXSkuycfHiRezfv9/hmdh8fX0xdOhQnDt3zu5xLy+vHrkMMBERKY9dKq4hewtHS7Jx9uxZ7Nu3D7/4xS8crqOmpgbnz5+3WXKXiIiIbl8OJxw1NTXSErUAUFJSgmPHjuHSpUtobGzE7373Oxw+fBjbtm1Dc3MzLBYLLBYLGhoapDqmTJmCjRs3Sq8XLVqEAwcO4MKFCzh48CAef/xx9OrVC7Gxsc5/QiIiopuxS8UlHO5SOXz4MCZPniy9Tk5OBgDMmjULaWlp+PzzzwEAoaGhNud99dVXmDRpEgDg/PnzuHr1qnTs8uXLiI2NxbVr1zBw4EBMmDAB3377LQYOHOhoeERERO1il4prOJxwTJo0Ce0tMNuZxWcvXLhg83rHjh2OhkFERES3Ea6lQkRE6sLF21yCCQcREakLEw6XYMJBRESqwjEcrqHoTKNEREREAFs4iIhIbdil4hJMOIiISFU0oghNJ56obO98chy7VIiIiEhxbOEgIiJ1YZeKSzDhICIiVeFTKq7BLhUiIiJSHFs4iIhIXdil4hJMOG7hEWiQrS6hpla2utyN1ttbtrqE6mrZ6urJmsosrg5BMXJ+NrnuTXf8WSDn9x3OXZKlGrm+fwWxUZZ6OoNdKq7BLhUiIiJSHFs4iIhIXdil4hJMOIiISFXYpeIaTDiIiEhd2MLhEhzDQURERIpjCwcREakOu0W6HxMOIiJSF1G8sTlzPjmMXSpERESkOLZwEBGRqvApFddgwkFEROrCp1Rcgl0qREREpDi2cBARkapohBubM+eT45hwEBGRurBLxSXYpUJERESKYwsHERGpCp9ScQ0mHEREpC6c+MslmHAQEZGqsIXDNZhw3KKpzOLqEFrRenu7OoRWhOpq2eqS6/MJNbWy1COnnv61k4tHoEG2ukTvvrLU0/T387LUA/Ts+0Cuz6YVGwD3uzVJRkw4iIhIXfiUiksw4SAiIlVhl4pr8LFYIiIiUhxbOIiISF34lIpLMOEgIiJVYZeKazjcpVJQUIDp06fDaDRCo9EgOzvb5vjs2bOh0WhstsjIyA7r3bRpEwYPHgydTgeTyYRDhw45GhoRERG5KYcTjtraWoSEhGDTpk1tlomMjERZWZm0ffTRR+3WuXPnTiQnJyM1NRVHjhxBSEgIIiIicOXKFUfDIyIiap8ow0YOc7hLJSoqClFRUe2W8fLygsHQ+efq161bh7lz52LOnDkAgIyMDOzZsweZmZlYsmSJoyESERG1iV0qrqHIUyr5+fnw9/fHsGHDMH/+fFy7dq3Nsg0NDSguLobZbP45KK0WZrMZhYWFds+pr6+H1Wq12YiIiMh9yZ5wREZG4oMPPkBeXh5ef/11HDhwAFFRUWhubrZb/urVq2hubkZAQIDN/oCAAFgs9mf9TE9Ph16vl7agoCC5PwYREfVUguj8Rg6T/SmVmJgY6f/HjBmDsWPH4t5770V+fj6mTJkiy3ssXboUycnJ0mur1cqkg4iIOoczjbqE4hN/3XPPPfDz88O5c+fsHvfz80OvXr1QXl5us7+8vLzNcSBeXl7w8fGx2YiIiDpDg5/HcXRpc/UHuE0pnnBcvnwZ165dQ2BgoN3jnp6eCAsLQ15enrRPEATk5eUhPDxc6fCIiIioGziccNTU1ODYsWM4duwYAKCkpATHjh3DpUuXUFNTgxdeeAHffvstLly4gLy8PMyYMQO//OUvERERIdUxZcoUbNy4UXqdnJyMrVu34v3338epU6cwf/581NbWSk+tEBERyaZlplFnNnKYw2M4Dh8+jMmTJ0uvW8ZSzJo1C5s3b8b333+P999/H5WVlTAajZg6dSpWrlwJLy8v6Zzz58/j6tWr0uuZM2fixx9/REpKCiwWC0JDQ5GTk9NqICkREZGz+FisaziccEyaNAliO9ndF1980WEdFy5caLUvKSkJSUlJjoZDREREtwGuFktEROrioplGHV3CY9euXRg+fDh0Oh3GjBmDvXv32n4MUURKSgoCAwPRp08fmM1mnD171qZMRUUF4uLi4OPjA19fXyQkJKCmpsamzPfff4+JEydCp9MhKCgIq1evtjm+detWTJw4Ef3790f//v1hNpu7tPwIEw4iIlIVjSg6vTnK0SU8Dh48iNjYWCQkJODo0aOIjo5GdHQ0Tpw4IZVZvXo1NmzYgIyMDBQVFaFv376IiIhAXV2dVCYuLg4nT55Ebm4udu/ejYKCAsybN086brVaMXXqVAwaNAjFxcVYs2YN0tLSsGXLFqlMfn4+YmNj8dVXX6GwsBBBQUGYOnUq/vnPfzp0DTRie/0jtwmr1Qq9Xo9JmAEPTW+n6tJ6e8sUFaDt11eWeprK7E+A1lPIdc2F6mpZ6pGTR2Dnp/jvSE++D+S8TnIRamrlq0ume1POn0/upklswP7qbaiqqlJsqoOW3xUTJ6XCw0PX5Xqamurwl/yXHYrVZDLhgQcekB6YEAQBQUFBeO655+wu4TFz5kzU1tZi9+7d0r6HHnoIoaGhyMjIgCiKMBqNeP7557Fo0SIAQFVVFQICApCVlYWYmBicOnUKI0eOxHfffYdx48YBAHJycjBt2jRcvnwZRqMRmzdvxksvvQSLxQJPT08AwJIlS5CdnY3Tp0/b/SzNzc3o378/Nm7ciPj4+E5fN7ZwEBGRuggybA7oyhIehYWFNuUBICIiQipfUlICi8ViU0av18NkMkllCgsL4evrKyUbAGA2m6HValFUVCSVefjhh6Vko+V9zpw5g+vXr9uN7aeffkJjYyMGDBjgyGWQf6ZRIiIid9bVbpGbzwfQah0vLy8vmycyW7S3hEdbrQgWi6XdJT9a/u2ojL+/v81xDw8PDBgwwKZMcHBwqzpajvXv379VbIsXL4bRaGyVEHWELRxERERdEBQUZLOuV3p6uqtDUtxrr72GHTt24NNPP4VO51i3FFs4iIhIXWRaS6W0tNRmDIe91g2ga0t4GAyGdsu3/FteXm4zk3d5eTlCQ0OlMrcOSm1qakJFRYVNPfbe5+b3aPHGG2/gtddew759+zB27Fi7cbeHLRxERKQuMs00euuaXm0lHF1ZwiM8PNymPADk5uZK5YODg2EwGGzKWK1WFBUVSWXCw8NRWVmJ4uJiqcz+/fshCAJMJpNUpqCgAI2NjTbvM2zYMJvulNWrV2PlypXIycmxGRPiCCYcRESkKk4t3NbFWUo7WsIjPj4eS5culcovWLAAOTk5WLt2LU6fPo20tDQcPnxYmiBTo9Fg4cKFWLVqFT7//HMcP34c8fHxMBqNiI6OBgCMGDECkZGRmDt3Lg4dOoRvvvkGSUlJiImJgdFoBAA89dRT8PT0REJCAk6ePImdO3di/fr1Niuyv/7661ixYgUyMzMxePBgWCwWWCyWVvN5dIRdKkRERArraAmPS5cuQav9uQ1g/Pjx2L59O5YvX45ly5ZhyJAhyM7OxujRo6UyL774ImprazFv3jxUVlZiwoQJyMnJsRlbsW3bNiQlJWHKlCnQarV44oknsGHDBum4Xq/Hl19+icTERISFhcHPzw8pKSk2c3Vs3rwZDQ0N+N3vfmfzmVJTU5GWltbpa8B5OG7BeTi6H+fh6JyefB9wHo7O4Twczmn5XfFI+HKn5+E4ULhK0Vh7IrZwEBGRqmiEG5sz55PjOIaDiIiIFMcWDiIiUpebnjTp8vnkMCYcRESkLjLNw0GOYZcKERERKY4tHEREpCpyraVCjmHCQURE6sIxHC7BLhUiIiJSHFs4iIhIXUQAzsylwQaOLulRCYfWux+0Gk+n6pBztsqePLugO14nOck2+6mMs1X25Nk45bxOcpFrpmB35W7fd4LY2HEhmXAMh2v0qISDiIioQyKcHMMhWySqwjEcREREpDi2cBARkbrwKRWXYMJBRETqIgDQOHk+OYxdKkRERKQ4tnAQEZGq8CkV12DCQURE6sIxHC7BLhUiIiJSHFs4iIhIXdjC4RJMOIiISF2YcLgEu1SIiIhIcWzhICIideE8HC7BhIOIiFSFj8W6BhMOIiJSF47hcAmHx3AUFBRg+vTpMBqN0Gg0yM7Otjmu0WjsbmvWrGmzzrS0tFblhw8f7vCHISIiIvfkcAtHbW0tQkJC8Mwzz+C3v/1tq+NlZWU2r//v//4PCQkJeOKJJ9qtd9SoUdi3b9/PgXmw8YWIiBQgiIDGiVYKgS0cXeHwb/WoqChERUW1edxgMNi8/uyzzzB58mTcc8897Qfi4dHqXCIiItmxS8UlFH0stry8HHv27EFCQkKHZc+ePQuj0Yh77rkHcXFxuHTpUptl6+vrYbVabTYiIiJyX4r2W7z//vvw9va22/VyM5PJhKysLAwbNgxlZWV4+eWXMXHiRJw4cQLe3t6tyqenp+Pll19utV+oroGg6S1b/O5CqK52dQitaO18XXoSd7zm7hiTXOS8n+S6TnJeb7k+nzveA3J9Nq3YAHTbx3OyhQNs4egKRVs4MjMzERcXB51O1265qKgoPPnkkxg7diwiIiKwd+9eVFZW4uOPP7ZbfunSpaiqqpK20tJSJcInIqKeqKVLxZmNHKZYC8df/vIXnDlzBjt37nT4XF9fXwwdOhTnzp2ze9zLywteXl7OhkhERETdRLEWjnfffRdhYWEICQlx+NyamhqcP38egYGBCkRGRESqJojOb+QwhxOOmpoaHDt2DMeOHQMAlJSU4NixYzaDPK1WK3bt2oVnn33Wbh1TpkzBxo0bpdeLFi3CgQMHcOHCBRw8eBCPP/44evXqhdjYWEfDIyIiap8oOL+RwxzuUjl8+DAmT54svU5OTgYAzJo1C1lZWQCAHTt2QBTFNhOG8+fP4+rVq9Lry5cvIzY2FteuXcPAgQMxYcIEfPvttxg4cKCj4REREZEbcjjhmDRpEsQOBszMmzcP8+bNa/P4hQsXbF7v2LHD0TCIiIi6hvNwuASn8yQiInURRDj1aCvHcHQJEw4iIlIXtnC4hKLzcBAREREBbOEgIiK1EeFkC4dskagKEw4iIlIXdqm4BLtUiIiISHFs4SAiInURBABOTN4lcOKvrmDCQURE6sIuFZdglwoREREpji0cRESkLmzhcAkmHEREpC6cadQlelTC4WEIgIfW09VhSJrKLK4OQTFCdbWrQ2hF6+3t6hBUxx2vuVwxyXmPu+P3i1zk+myC2ChLPeS+elTCQURE1BFRFCA6scS8M+eqGRMOIiJSF1F0rluEYzi6hAkHERGpi+jkGA4mHF3Cx2KJiIhIcWzhICIidREEQOPEOAyO4egSJhxERKQu7FJxCXapEBERkeLYwkFERKoiCgJEJ7pU+Fhs1zDhICIidWGXikuwS4WIiIgUxxYOIiJSF0EENGzh6G5MOIiISF1EEYAzj8Uy4egKdqkQERGR4tjCQUREqiIKIkQnulREtnB0CVs4iIhIXUTB+a0LNm3ahMGDB0On08FkMuHQoUPtlt+1axeGDx8OnU6HMWPGYO/evbYfQxSRkpKCwMBA9OnTB2azGWfPnrUpU1FRgbi4OPj4+MDX1xcJCQmoqamxKfP9999j4sSJ0Ol0CAoKwurVqx2OpTOYcBARkaqIguj05qidO3ciOTkZqampOHLkCEJCQhAREYErV67YLX/w4EHExsYiISEBR48eRXR0NKKjo3HixAmpzOrVq7FhwwZkZGSgqKgIffv2RUREBOrq6qQycXFxOHnyJHJzc7F7924UFBRg3rx50nGr1YqpU6di0KBBKC4uxpo1a5CWloYtW7Y4FEtnaMQe0DZktVqh1+thNsyDh9bT1eFImsosrg5BVbTe3rLVJVRXy1ZXTybnNXc3vAe6V5PYiHx8hqqqKvj4+CjyHi2/KyZpHoeHpneX62kSG5EvfupQrCaTCQ888AA2btwIABAEAUFBQXjuueewZMmSVuVnzpyJ2tpa7N69W9r30EMPITQ0FBkZGRBFEUajEc8//zwWLVoEAKiqqkJAQACysrIQExODU6dOYeTIkfjuu+8wbtw4AEBOTg6mTZuGy5cvw2g0YvPmzXjppZdgsVjg6Xnj9+eSJUuQnZ2N06dPdyqWzuoRYzhacqYmocHFkdhqEhtdHYKqaEX5vv4Cv3adIuc1dze8B7pXE25c7+74G7hJrHdqAbaWWK1Wq81+Ly8veHl5tSrf0NCA4uJiLF26VNqn1WphNptRWFho9z0KCwuRnJxssy8iIgLZ2dkAgJKSElgsFpjNZum4Xq+HyWRCYWEhYmJiUFhYCF9fXynZAACz2QytVouioiI8/vjjKCwsxMMPPywlGy3v8/rrr+P69evo379/h7F0Vo9IOKr//ZdI/pUs1wZCrsU/SLsfrznJrLq6Gnq9XpG6PT09YTAY8LXF8fEHt+rXrx+CgoJs9qWmpiItLa1V2atXr6K5uRkBAQE2+wMCAqRWhFtZLBa75S0Wi3S8ZV97Zfz9/W2Oe3h4YMCAATZlgoODW9XRcqx///4dxtJZPSLhMBqNKC0thbe3NzQaTZvlrFYrgoKCUFpaqliTnRIYd/e6XeMGbt/YGXf3cse4RVFEdXU1jEajYu+h0+lQUlKChgbnW+ZEUWz1+8Ze6wb9rEckHFqtFnfddVeny/v4+LjNN5kjGHf3ul3jBm7f2Bl393K3uJVq2biZTqeDTqdT/H1u5ufnh169eqG8vNxmf3l5OQwGg91zDAZDu+Vb/i0vL0dgYKBNmdDQUKnMrYNSm5qaUFFRYVOPvfe5+T06iqWz+JQKERGRgjw9PREWFoa8vDxpnyAIyMvLQ3h4uN1zwsPDbcoDQG5urlQ+ODgYBoPBpozVakVRUZFUJjw8HJWVlSguLpbK7N+/H4IgwGQySWUKCgrQ2Nho8z7Dhg1D//79OxVLp4kqUlVVJQIQq6qqXB2KQxh397pd4xbF2zd2xt29bte4b2c7duwQvby8xKysLPFvf/ubOG/ePNHX11e0WCyiKIri008/LS5ZskQq/80334geHh7iG2+8IZ46dUpMTU0Ve/fuLR4/flwq89prr4m+vr7iZ599Jn7//ffijBkzxODgYPFf//qXVCYyMlK87777xKKiIvHrr78WhwwZIsbGxkrHKysrxYCAAPHpp58WT5w4Ie7YsUO84447xHfeecehWDpDVQlHXV2dmJqaKtbV1bk6FIcw7u51u8Ytirdv7Iy7e92ucd/u3n77bfHuu+8WPT09xQcffFD89ttvpWOPPPKIOGvWLJvyH3/8sTh06FDR09NTHDVqlLhnzx6b44IgiCtWrBADAgJELy8vccqUKeKZM2dsyly7dk2MjY0V+/XrJ/r4+Ihz5swRq6urbcr89a9/FSdMmCB6eXmJd955p/jaa6+1ir2jWDqjR8zDQURERO6NYziIiIhIcUw4iIiISHFMOIiIiEhxTDiIiIhIcT0u4ZB7+V+lpaen44EHHoC3tzf8/f0RHR2NM2fOtHtOVlYWNBqNzdbdE9mkpaW1imH48OHtnuPqaw0AgwcPbhW3RqNBYmKi3fKuvNYFBQWYPn06jEYjNBpNq3ULxE4sTW2Po98jcsbd2NiIxYsXY8yYMejbty+MRiPi4+Pxww8/tFtnV+43OeMGgNmzZ7eKITIyssN6XXm9Adi93zUaDdasWdNmnd1xvUl9elTCocTyv0o7cOAAEhMT8e233yI3NxeNjY2YOnUqamtr2z3Px8cHZWVl0nbx4sVuivhno0aNsonh66+/brOsO1xrAPjuu+9sYs7NzQUAPPnkk22e46prXVtbi5CQEGzatMnu8c4sTX0rR79H5I77p59+wpEjR7BixQocOXIEn3zyCc6cOYPHHnusw3odud/kjrtFZGSkTQwfffRRu3W6+noDsIm3rKwMmZmZ0Gg0eOKJJ9qtV+nrTSrk8IO0buzBBx8UExMTpdfNzc2i0WgU09PT7Zb//e9/Lz766KM2+0wmk/if//mfisbZnitXrogAxAMHDrRZ5r333hP1en33BWVHamqqGBIS0uny7nitRVEUFyxYIN57772iIAh2j7vDtRZFUQQgfvrpp9JrQRBEg8EgrlmzRtpXWVkpenl5iR999FGb9Tj6PSJ33PYcOnRIBCBevHixzTKO3m/Oshf3rFmzxBkzZjhUjzte7xkzZoi//vWv2y3T3deb1KHHtHC0LP9781K9nVn+9+bywI0ld9sq3x2qqqoAAAMGDGi3XE1NDQYNGoSgoCDMmDEDJ0+e7I7wbJw9exZGoxH33HMP4uLicOnSpTbLuuO1bmhowIcffohnnnmm3UX/3OFa36qjpant6cr3SHeoqqqCRqOBr69vu+Ucud+Ukp+fD39/fwwbNgzz58/HtWvX2izrjte7vLwce/bsQUJCQodl3eF6U8/SYxKO9pb/bWsJXbmW3JWLIAhYuHAhfvWrX2H06NFtlhs2bBgyMzPx2Wef4cMPP4QgCBg/fjwuX77cbbGaTCZkZWUhJycHmzdvRklJCSZOnIjqavvrlbvbtQaA7OxsVFZWYvbs2W2WcYdrbU9nlqa+VVe+R5RWV1eHxYsXIzY2tt1FxBy935QQGRmJDz74AHl5eXj99ddx4MABREVFobm52W55d7ze77//Pry9vfHb3/623XLucL2p5+kRq8X2FImJiThx4kSHfaXh4eE2i+aMHz8eI0aMwDvvvIOVK1cqHSYAICoqSvr/sWPHwmQyYdCgQfj444879deTO3j33XcRFRXV7nLY7nCte6rGxkb8/ve/hyiK2Lx5c7tl3eF+i4mJkf5/zJgxGDt2LO69917k5+djypQp3RKDszIzMxEXF9fhwGd3uN7U8/SYFg4llv/tTklJSdi9eze++uor3HXXXQ6d27t3b9x33304d+6cQtF1zNfXF0OHDm0zBne61gBw8eJF7Nu3D88++6xD57nDtQZsl6a+WXvXtCvfI0ppSTYuXryI3Nxch5dI7+h+6w733HMP/Pz82ozBna43APzlL3/BmTNnHL7nAfe43nT76zEJhxLL/3YHURSRlJSETz/9FPv370dwcLDDdTQ3N+P48eMIDAxUIMLOqampwfnz59uMwR2u9c3ee+89+Pv749FHH3XoPHe41kDnlqa+VVe+R5TQkmycPXsW+/btwy9+8QuH6+jofusOly9fxrVr19qMwV2ud4t3330XYWFhCAkJcfhcd7je1AO4etSqnJRY/ldp8+fPF/V6vZifny+WlZVJ208//SSVuTXul19+Wfziiy/E8+fPi8XFxWJMTIyo0+nEkydPdlvczz//vJifny+WlJSI33zzjWg2m0U/Pz/xypUrdmN2h2vdorm5Wbz77rvFxYsXtzrmTte6urpaPHr0qHj06FERgLhu3Trx6NGj0tMcnVma+te//rX49ttvS687+h5ROu6GhgbxscceE++66y7x2LFjNvd8fX19m3F3dL8pHXd1dbW4aNEisbCwUCwpKRH37dsn3n///eKQIUNsVlx1t+vdoqqqSrzjjjvEzZs3263DFdeb1KdHJRyiKP/yv0oDYHd777332ox74cKF0mcMCAgQp02bJh45cqRb4545c6YYGBgoenp6infeeac4c+ZM8dy5c23GLIquv9YtvvjiCxFAq2WcRdG9rvVXX31l995oia8zS1MPGjRITE1NtdnX3veI0nGXlJS0ec9/9dVXbcbd0f2mdNw//fSTOHXqVHHgwIFi7969xUGDBolz585tlTi42/Vu8c4774h9+vQRKysr7dbhiutN6sPl6YmIiEhxPWYMBxEREbkvJhxERESkOCYcREREpDgmHERERKQ4JhxERESkOCYcREREpDgmHERERKQ4JhxERESkOCYcREREpDgmHERERKQ4JhxERESkOCYcREREpLj/H/8P5G/jNaH2AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Setup scene\n",
    "\n",
    "scene_file = \"../scenes/simple_backscattering_reflector/simple_backscattering_reflector.xml\"\n",
    "tx_position = np.array([1.0, 0.0, 0.5])\n",
    "tx_orientation = np.array([0.0, 0.0, 0.0])\n",
    "tx_pattern = 'iso'\n",
    "fc = 3.5e9\n",
    "tx_slant_angle = 0.0\n",
    "\n",
    "# Path loss map used for evaluation\n",
    "cm_center = np.array([-1., 0, 0.6])\n",
    "cm_orientation = np.array([0., 0., 0.])\n",
    "cm_size = np.array([2., 2.])\n",
    "cm_cell_size = np.array([0.1, 0.1])\n",
    "num_samples = int(1e7)\n",
    "max_depth = 1\n",
    "\n",
    "scene = mi.load_file(scene_file)\n",
    "params = mi.traverse(scene)\n",
    "\n",
    "params['mat-itu_metal.eta_r'] = 1.\n",
    "params['mat-itu_metal.eta_i'] = 1.0\n",
    "params['mat-itu_metal.s'] = 0.7\n",
    "params['mat-itu_metal.lambda'] = 0.7\n",
    "params.update()\n",
    "\n",
    "# Define Loss and set trainable variables\n",
    "\n",
    "def func_obj(x):\n",
    "    return -1e3*dr.sum(x)\n",
    "\n",
    "variables =[params['mat-itu_metal.eta_r'],\n",
    "            params['mat-itu_metal.eta_i'],\n",
    "            params['mat-itu_metal.s'],\n",
    "            params['mat-itu_metal.lambda']\n",
    "           ]\n",
    "for v in variables:\n",
    "    dr.enable_grad(v)\n",
    "    \n",
    "# Compute gradients with RBP\n",
    "    \n",
    "# Reset grads\n",
    "for v in variables:\n",
    "    dr.set_grad(v, 0.0)\n",
    "\n",
    "rbp_tracer = PathlossMapRBPTracer(scene, fc, tx_pattern, tx_slant_angle,\n",
    "                                  cm_center, cm_orientation, cm_size, cm_cell_size,\n",
    "                                  num_samples, max_depth)\n",
    "\n",
    "cm_test,_ = rbp_tracer(tx_position, tx_orientation, func_obj, loop_record=True)\n",
    "g_test = np.array([dr.grad(v).numpy()[0] for v in variables])\n",
    "print(\"RBP Grad: \", g_test)\n",
    "\n",
    "# Compute gradients with conventional automatic differentiation\n",
    "\n",
    "# Reset grads\n",
    "for v in variables:\n",
    "    dr.set_grad(v, 0.0)\n",
    "\n",
    "ref_tracer = MapTracer(scene, fc, tx_pattern, tx_slant_angle,\n",
    "                       cm_center, cm_orientation, cm_size, cm_cell_size,\n",
    "                       num_samples, max_depth)\n",
    "cm_ref,_,_,_ = ref_tracer(tx_position, tx_orientation, loop_record=False)\n",
    "l = func_obj(cm_ref)\n",
    "dr.backward(l)\n",
    "g_ref = np.array([dr.grad(v).numpy()[0] for v in variables])\n",
    "print(\"Ref. Grad: \", g_ref)\n",
    "\n",
    "print(\"Relative Error: \", np.linalg.norm(g_ref-g_test)/np.linalg.norm(g_ref))\n",
    "\n",
    "cm_ref_db = 10.0*np.log10(cm_ref.numpy()).T\n",
    "cm_test_db = 10.0*np.log10(cm_test.numpy()).T\n",
    "ale = np.abs(cm_ref_db - cm_test_db)\n",
    "plt.figure()\n",
    "plt.imshow(ale)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fea85354",
   "metadata": {},
   "source": [
    "# Double Reflector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e2933551",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RBP Grad:  [ 1.8727098e-04 -7.4322161e-04  8.3550689e-04 -6.5419845e-05]\n",
      "Ref. Grad:  [ 1.8727082e-04 -7.4322108e-04  8.3550601e-04 -6.5419837e-05]\n",
      "Relative Error:  9.075952e-07\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f98f43b8e10>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAGsCAYAAAD9ro91AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0nUlEQVR4nO3de3RU5b3/8c9wm6AmQSWQpEYuIncICBICWqGkhByKYHuQZtFDRMRzXOEcJUcPpr8qeDkG61FsCyuoS4guq6CnXFyKcCA1UAuRa1ahtZTQQIKQICxJSCwJndm/P2xGB3KZyd4z2ZP9fnU9a7H3PM+zv7Md+uV59uVxGYZhCAAA2Fan9g4AAAC0jGQNAIDNkawBALA5kjUAADZHsgYAwOZI1gAA2BzJGgAAmyNZAwBgcyRrAABsjmQNAIDNkawBALa3a9cuzZgxQ4mJiXK5XNq0aVPIj/n555/rJz/5iW688UZ1795dI0aM0P79+0N+3KaQrAEAtldXV6fk5GStWrUqLMf78ssvNXHiRHXt2lUfffSR/vSnP+nFF1/U9ddfH5bjX8nFQh4AgEjicrm0ceNGzZo1y7evvr5e/+///T+98847unDhgoYPH67nn39ekyZNatMxHn/8cf3+97/X7373O2uCNomRNQAg4i1atEh79uzRunXr9Ic//EGzZ8/WtGnTdOzYsTb19/7772vs2LGaPXu2evXqpdGjR+u1116zOOrAMbIGAESUK0fW5eXl6t+/v8rLy5WYmOirl5aWpnHjxum5554L+hhRUVGSpJycHM2ePVv79u3Tww8/rNWrVysrK8uS7xGMLmE/IgAAFjp8+LA8Ho8GDhzot7++vl433nijJOnPf/6zhgwZ0mI/S5Ys0fLlyyVJXq9XY8eO9SX60aNH68iRIyRrAADaora2Vp07d9aBAwfUuXNnv8+uu+46SVL//v312WeftdhPY2KXpISEBA0dOtTv8yFDhug3v/mNRVEHh2QNAIhoo0ePlsfj0dmzZ3XnnXc2Wadbt24aPHhwwH1OnDhRR48e9dv3l7/8RX369DEVa1uRrAEAtldbW6vS0lLfdllZmUpKSnTDDTdo4MCBmjt3rubNm6cXX3xRo0eP1hdffKHCwkKNHDlS06dPD/p4ixcv1oQJE/Tcc8/p3nvv1d69e/Xqq6/q1VdftfJrBYwbzAAAtldUVKTJkydftT8rK0sFBQW6fPmynn32Wb355pv6/PPP1bNnT40fP15PPfWURowY0aZjfvDBB8rNzdWxY8fUr18/5eTkaOHChWa/SpuQrAEAsDmeswYAwOZI1gAA2FyHuMHM6/Xq9OnTio6Olsvlau9wAABBMgxDFy9eVGJiojp1Ct048tKlS2poaDDdT7du3XwvTgmHDpGsT58+raSkpPYOAwBgUkVFhW666aaQ9H3p0iX163OdKs96TPcVHx+vsrKysCXsDpGso6OjJUl36J/URV3bORogNP66/HZL+un/+D5L+gGs9Hdd1ifa4vv/81BoaGhQ5VmPyg70UUx020fvNRe96jfmpBoaGkjWwWic+u6iruriIlmjY+pk0f8p8HcEtvSP55LCcSkzJrqTqWTdHjpEsgYAIFAewyuPiYeWPYbXumACRLIGADiKV4a8anu2NtO2rUI2D7Bq1Sr17dtXUVFRSklJ0d69e1us/95772nw4MGKiorSiBEjtGXLllCFBgBwMK8F/wu3kCTr9evXKycnR0uXLtXBgweVnJys9PR0nT17tsn6u3fvVmZmphYsWKBDhw5p1qxZmjVrlo4cORKK8AAAiCghSdYvvfSSFi5cqPnz52vo0KFavXq1rrnmGq1Zs6bJ+r/4xS80bdo0PfbYYxoyZIieeeYZ3XbbbVq5cmUowgMAOJjHMEyXcLM8WTc0NOjAgQNKS0v75iCdOiktLU179uxpss2ePXv86ktSenp6s/Xr6+tVU1PjVwAACETjNWszJdwsT9bnzp2Tx+NR7969/fb37t1blZWVTbaprKwMqn5eXp5iY2N9hReiAAA6ssh60OwfcnNzVV1d7SsVFRXtHRIAIEJ4ZchjorTHyNryR7d69uypzp07q6qqym9/VVWV4uPjm2wTHx8fVH232y23221NwAAAR+HRLX39cvMxY8aosLDQt8/r9aqwsFCpqalNtklNTfWrL0nbt29vtj4AAE4Skpei5OTkKCsrS2PHjtW4ceP08ssvq66uTvPnz5ckzZs3T9/5zneUl5cnSXr44Yd111136cUXX9T06dO1bt067d+/X6+++moowgMAOJjZO7rb427wkCTrOXPm6IsvvtCTTz6pyspKjRo1Slu3bvXdRFZeXu63BNqECRP09ttv62c/+5l++tOf6tZbb9WmTZs0fPjwUIQHAHAw7z+Kmfbh5jKMdvgngsVqamoUGxurSZrJIgXosEpXjLeknwGLiy3pB7DS343LKtJmVVdXKyYmJiTHaMwVf/6st6JNLORx8aJXg4dUhTTWK/FucACAozTe1W2mfbiRrCOAVSMqqWOPqradLrGsr/TEUZb1ZZWO/N+O33hgOE/W8BgyueqWdbEEimQNAHCUSLxmHZEvRQEAwEkYWQMAHMUrlzxymWofbiRrAICjeI2vi5n24cY0OAAANsfIGgDgKB6T0+Bm2rYVyRoA4CiRmKyZBgcAwOZI1gAAR/EaLtMlGH379pXL5bqqZGdnB9wH0+AAAEcJ9zT4vn375PF4fNtHjhzR97//fc2ePTvgPkjWAACEUFxcnN/28uXLdcstt+iuu+4KuA+SNQDAUTzqJI+Jq8CNY+Samhq//W63W263u8W2DQ0Neuutt5STkyOXK/AROtesAQCOYpi8Xm3845p1UlKSYmNjfSUvL6/VY2/atEkXLlzQfffdF1TMjKwBAI5i1TXriooKv/WsWxtVS9Lrr7+ujIwMJSYmBnVMkjUAAG0QExPjl6xbc/LkSe3YsUMbNmwI+lgkawCAo3iMTvIYJq5Zt/Hd4GvXrlWvXr00ffr0oNuSrAEAjuKVS14Tt2x5FXy29nq9Wrt2rbKystSlS/CplxvMAAAIsR07dqi8vFz3339/m9ozskablK4Yb1lfAxYXW9JPeuIoS/qRrPt+Vn23js7K89SR/9sdn7Pasr7SF4+yrK9I0x7vBp86daoMo+1ra5KsAQCOYv6adfgXtGYaHAAAm2NkDQBwlK9vMGv7NLiZtm1FsgYAOIrX5OtG23I3uFlMgwMAYHOMrAEAjhKJN5iRrAEAjuJVp7C/FMUskjUAwFE8hksew8Rz1ibathXXrAEAsDlG1gAAR/GYvBvcwzQ4AACh5TU6yWviBjMvbzADAABXYmQNAHAUpsEBALA5r8zd0e21LpSAWT4NnpeXp9tvv13R0dHq1auXZs2apaNHj7bYpqCgQC6Xy69ERUVZHRoAABHJ8mS9c+dOZWdnq7i4WNu3b9fly5c1depU1dXVtdguJiZGZ86c8ZWTJ09aHRoAAL6Xopgp4Wb5NPjWrVv9tgsKCtSrVy8dOHBA3/3ud5tt53K5FB8fb3U4AAD4Mf+60Q6QrK9UXV0tSbrhhhtarFdbW6s+ffrI6/Xqtttu03PPPadhw4Y1Wbe+vl719fW+7ZqaGusCtqEBi4vbO4Sr2DEmBGbb6RLL+kpPHGVZX1bpyL9NO55vhEdI/3ng9Xr1yCOPaOLEiRo+fHiz9QYNGqQ1a9Zo8+bNeuutt+T1ejVhwgSdOnWqyfp5eXmKjY31laSkpFB9BQBAB9O4nrWZEm4hTdbZ2dk6cuSI1q1b12K91NRUzZs3T6NGjdJdd92lDRs2KC4uTq+88kqT9XNzc1VdXe0rFRUVoQgfANABNU6DmynhFrJp8EWLFumDDz7Qrl27dNNNNwXVtmvXrho9erRKS0ub/NztdsvtdlsRJgDAYcw/Zx3+ZG35EQ3D0KJFi7Rx40b99re/Vb9+/YLuw+Px6PDhw0pISLA6PAAAIo7lI+vs7Gy9/fbb2rx5s6Kjo1VZWSlJio2NVffu3SVJ8+bN03e+8x3l5eVJkp5++mmNHz9eAwYM0IULF/TCCy/o5MmTeuCBB6wODwDgcF7DJa+Zl6K0wxKZlifr/Px8SdKkSZP89q9du1b33XefJKm8vFydOn0zqP/yyy+1cOFCVVZW6vrrr9eYMWO0e/duDR061OrwAAAO5zU5Dd4hnrM2AliNpKioyG97xYoVWrFihdWhAADQIfBucACAo5hfIrMDjKwBALAzj1zymHhW2kzbtmI9awAAbI6RNQDAUZgGBwDA5jwyN5XtsS6UgDENDgCAzTGyBgA4CtPgAADYHOtZAwBgc4bJZS4NHt0CAABXYmQNAHAUpsFhe6UrxlvSz4DFxZb0Y1cd+fulJ45q7xCuYtXvUurY/+1gjUhcdYtpcAAAbI5kDQBwFM8/lsg0U4L1+eef6yc/+YluvPFGde/eXSNGjND+/fsDbs80OADAUcI9Df7ll19q4sSJmjx5sj766CPFxcXp2LFjuv766wPug2QNAEAIPf/880pKStLatWt9+/r16xdUH0yDAwAcxatOposk1dTU+JX6+vomj/f+++9r7Nixmj17tnr16qXRo0frtddeCypmkjUAwFE8hst0kaSkpCTFxsb6Sl5eXpPH++tf/6r8/Hzdeuut2rZtmx566CH9x3/8h954442AY2YaHACANqioqFBMTIxv2+12N1nP6/Vq7Nixeu655yRJo0eP1pEjR7R69WplZWUFdCySNQDAUay6wSwmJsYvWTcnISFBQ4cO9ds3ZMgQ/eY3vwn4mCRrAICjGCZX3TKCbDtx4kQdPXrUb99f/vIX9enTJ+A+SNYAAEfxyCWPicU4gm27ePFiTZgwQc8995zuvfde7d27V6+++qpeffXVgPvgBjMAAELo9ttv18aNG/XOO+9o+PDheuaZZ/Tyyy9r7ty5AffByBoA4Chew9z7vb1G8G1+8IMf6Ac/+EGbj0myBgA4itfkNWszbduKaXAAAGyOkTUAwFG8cslr4gYzM23bimQNAHCUb7+FrK3tw41pcAAAbI6R9RVKV4y3rK8Bi4st68sqdowJgbHqt2nlb8COMQGticQbzEjWAABH8crk60bb4Zo10+AAANgcI2sAgKMYJu8GN7gbHACA0LJq1a1wIlkDABwlEm8ws/yIy5Ytk8vl8iuDBw9usc17772nwYMHKyoqSiNGjNCWLVusDgsAgIgVkn8eDBs2TGfOnPGVTz75pNm6u3fvVmZmphYsWKBDhw5p1qxZmjVrlo4cORKK0AAADtc4DW6mhFtIknWXLl0UHx/vKz179my27i9+8QtNmzZNjz32mIYMGaJnnnlGt912m1auXBmK0AAADtf4ulEzJdxCkqyPHTumxMRE9e/fX3PnzlV5eXmzdffs2aO0tDS/fenp6dqzZ0+zberr61VTU+NXAADoqCxP1ikpKSooKNDWrVuVn5+vsrIy3Xnnnbp48WKT9SsrK9W7d2+/fb1791ZlZWWzx8jLy1NsbKyvJCUlWfodAAAdF9PgkjIyMjR79myNHDlS6enp2rJliy5cuKB3333XsmPk5uaqurraVyoqKizrGwDQsUVisg75o1s9evTQwIEDVVpa2uTn8fHxqqqq8ttXVVWl+Pj4Zvt0u91yu92WxgkAgF2F/GGx2tpaHT9+XAkJCU1+npqaqsLCQr9927dvV2pqaqhDAwA4UCSOrC1P1o8++qh27typEydOaPfu3brnnnvUuXNnZWZmSpLmzZun3NxcX/2HH35YW7du1Ysvvqg///nPWrZsmfbv369FixZZHRoAABGZrC2fBj916pQyMzN1/vx5xcXF6Y477lBxcbHi4uIkSeXl5erU6Zt/I0yYMEFvv/22fvazn+mnP/2pbr31Vm3atEnDhw+3OjQAACKS5cl63bp1LX5eVFR01b7Zs2dr9uzZVocCAMBVDJlb5tKwLpSA8W5wAICjsJAHAAA2R7LuAAYsLrasr9IV4y3px8qYrGLVd5Ps+f3syI7nyY4xWYXfOOyEZA0AcBRG1gAA2FwkJuvwr6ANAACCwsgaAOAohuGSYWJ0bKZtW5GsAQCOYnZN6g6znjUAALAOI2sAgKNE4g1mJGsAgKNE4jVrpsEBALA5RtYAAEdhGhwAAJuLxGlwkjUAwFEMkyNrrlkDANDBLFu2TC6Xy68MHjw4qD4YWQMAHMWQZBjm2gdr2LBh2rFjh2+7S5fg0i/JGgDgKF655ArzG8y6dOmi+Pj4Nh+TaXAAANqgpqbGr9TX1zdb99ixY0pMTFT//v01d+5clZeXB3UskjUAwFEa7wY3UyQpKSlJsbGxvpKXl9fk8VJSUlRQUKCtW7cqPz9fZWVluvPOO3Xx4sWAY2YaPIQGLC5u7xBCpiN/NwAdm9dwyWXBc9YVFRWKiYnx7Xe73U3Wz8jI8P155MiRSklJUZ8+ffTuu+9qwYIFAR2TZA0AQBvExMT4JetA9ejRQwMHDlRpaWnAbZgGBwA4imGYL2bU1tbq+PHjSkhICLgNyRoA4ChWXbMO1KOPPqqdO3fqxIkT2r17t+655x517txZmZmZAffBNDgAACF06tQpZWZm6vz584qLi9Mdd9yh4uJixcXFBdwHyRoA4Cjhfjf4unXr2nysRiRrAICjWHU3eDiRrAEAjmL2JjGzN5i1BTeYAQBgc4ysAQCO8vXI2sw1awuDCRDJGgDgKOG+wcwKTIMDAGBzjKwBAI5iqG1rUn+7fbiRrAEAjsI0OAAAsBwjawCAs0TgPLjlI+u+ffvK5XJdVbKzs5usX1BQcFXdqKgoq8MCAOBrZhfx6AhvMNu3b588Ho9v+8iRI/r+97+v2bNnN9smJiZGR48e9W27XOE/EQAAZ4jEN5hZnqyvXEVk+fLluuWWW3TXXXc128blcik+Pt7qUAAA6BBCes26oaFBb731lnJyclocLdfW1qpPnz7yer267bbb9Nxzz2nYsGHN1q+vr1d9fb1vu6amxtK4gdIV4y3pZ8DiYkv6sauOfJ6sjKkjn6dIxN3gV9i0aZMuXLig++67r9k6gwYN0po1a7R582a99dZb8nq9mjBhgk6dOtVsm7y8PMXGxvpKUlJSCKIHAHRIjdedzZQwC2myfv3115WRkaHExMRm66SmpmrevHkaNWqU7rrrLm3YsEFxcXF65ZVXmm2Tm5ur6upqX6moqAhF+AAA2ELIpsFPnjypHTt2aMOGDUG169q1q0aPHq3S0tJm67jdbrndbrMhAgAcKBJvMAvZyHrt2rXq1auXpk+fHlQ7j8ejw4cPKyEhIUSRAQAczbCghFlIkrXX69XatWuVlZWlLl38B+/z5s1Tbm6ub/vpp5/W//3f/+mvf/2rDh48qJ/85Cc6efKkHnjggVCEBgBAxAnJNPiOHTtUXl6u+++//6rPysvL1anTN/9G+PLLL7Vw4UJVVlbq+uuv15gxY7R7924NHTo0FKEBABwuEu8GD0mynjp1qoxmJvWLior8tlesWKEVK1aEIgwAAJrWHktnmcBCHgAA2BwLeQAAHIVpcAAA7C4CV90iWQMAHMb1j2KmfXhxzRoAAJtjZA0AcBamwQEAsLkITNZMgwMAYHOMrAEAzmJ2mUse3QIAILQicdUtkjUANKF0xXjL+hqwuNiSfuwYE8KDZA0AcJYIvMGMZA0AcJYIvGbN3eAAANgcI2sAgKO4jK+LmfbhRrIGADhLBF6zZhocAOAsjdeszRQTli9fLpfLpUceeSTgNiRrAADCZN++fXrllVc0cuTIoNqRrAEAzmJYUNqgtrZWc+fO1Wuvvabrr78+qLYkawCAs7RTss7Oztb06dOVlpYWdFtuMAMAoA1qamr8tt1ut9xud5N1161bp4MHD2rfvn1tOhYjawCAs1g0sk5KSlJsbKyv5OXlNXm4iooKPfzww/r1r3+tqKioNoXMyBoA4CwWvcGsoqJCMTExvt3NjaoPHDigs2fP6rbbbvPt83g82rVrl1auXKn6+np17ty5xUOSrAEAaIOYmBi/ZN2cKVOm6PDhw3775s+fr8GDB2vJkiWtJmqJZA0AcJhwv8EsOjpaw4cP99t37bXX6sYbb7xqf3NI1gAAZ4nAN5iRrAEACLOioqKg6nM3OAAANsfIGgDgKC6ZvGZtWSSBI1lfoXTF+PYO4SoDFhe3dwghZdU5t/I8dfRzjtZZ+Ruw42/c0Sx6dCucmAYHAMDmGFkDAJyFu8EBALC5CEzWTIMDAGBzjKwBAI4S7jeYWSHokfWuXbs0Y8YMJSYmyuVyadOmTX6fG4ahJ598UgkJCerevbvS0tJ07NixVvtdtWqV+vbtq6ioKKWkpGjv3r3BhgYAQOvaaT1rM4JO1nV1dUpOTtaqVaua/PznP/+5fvnLX2r16tX69NNPde211yo9PV2XLl1qts/169crJydHS5cu1cGDB5WcnKz09HSdPXs22PAAAOhwgk7WGRkZevbZZ3XPPfdc9ZlhGHr55Zf1s5/9TDNnztTIkSP15ptv6vTp01eNwL/tpZde0sKFCzV//nwNHTpUq1ev1jXXXKM1a9YEGx4AAC1zwsi6JWVlZaqsrFRaWppvX2xsrFJSUrRnz54m2zQ0NOjAgQN+bTp16qS0tLRm29TX16umpsavAAAQiMZr1mZKuFmarCsrKyVJvXv39tvfu3dv32dXOnfunDweT1Bt8vLyFBsb6ytJSUkWRA8AgD1F5KNbubm5qq6u9pWKior2DgkAECkaXzdqpoSZpY9uxcfHS5KqqqqUkJDg219VVaVRo0Y12aZnz57q3Lmzqqqq/PZXVVX5+ruS2+2W2+22JmgAgLM4/aUo/fr1U3x8vAoLC337ampq9Omnnyo1NbXJNt26ddOYMWP82ni9XhUWFjbbBgCAtorEa9ZBj6xra2tVWlrq2y4rK1NJSYluuOEG3XzzzXrkkUf07LPP6tZbb1W/fv30xBNPKDExUbNmzfK1mTJliu655x4tWrRIkpSTk6OsrCyNHTtW48aN08svv6y6ujrNnz/f/DcEACDCBZ2s9+/fr8mTJ/u2c3JyJElZWVkqKCjQf/3Xf6murk4PPvigLly4oDvuuENbt25VVFSUr83x48d17tw53/acOXP0xRdf6Mknn1RlZaVGjRqlrVu3XnXTGQAApkXgNHjQyXrSpEkyjOYjdblcevrpp/X00083W+fEiRNX7Vu0aJFvpA0AQMiYncqO9GvWAADAeizkcYUBi4vbOwS0UemK8Zb1xe8gMJynwHCebMYJ0+AAAES0CEzWTIMDAGBzjKwBAI7iiPWsAQBAeJGsAQCwOabBAQDOEoE3mJGsAQCOEonXrEnWAADnaYeEawbXrAEAsDlG1gAAZ+GaNQAA9haJ16yZBgcAwOYYWQMAnIVpcAAA7I1pcAAAYDlG1gAAZ2EaHAAAm4vAZM00OAAAIZSfn6+RI0cqJiZGMTExSk1N1UcffRRUH4ysHWbb6RJL+klPHGVJP5I0YHGxJf1Y9d0kKX3xKMv6QutKV4y3rC+rfk92xHmyRrhvMLvpppu0fPly3XrrrTIMQ2+88YZmzpypQ4cOadiwYQH1QbIGADhLmKfBZ8yY4bf93//938rPz1dxcTHJGgCAJlmUrGtqavx2u91uud3uFpt6PB699957qqurU2pqasCH5Jo1AABtkJSUpNjYWF/Jy8trtu7hw4d13XXXye1269/+7d+0ceNGDR06NOBjMbIGADiKVdesKyoqFBMT49vf0qh60KBBKikpUXV1tf73f/9XWVlZ2rlzZ8AJm2QNAHAWi6bBG+/uDkS3bt00YMAASdKYMWO0b98+/eIXv9Arr7wSUHumwQEACDOv16v6+vqA6zOyBgA4Srgf3crNzVVGRoZuvvlmXbx4UW+//baKioq0bdu2gPsgWQMAnCXMj26dPXtW8+bN05kzZxQbG6uRI0dq27Zt+v73vx9wHyRrAABC6PXXXzfdB8kaAOAsEfhucJI1AMBRXP8oZtqHG3eDAwBgc4ysAQDOwjQ4AAD2Fu5Ht6wQ9DT4rl27NGPGDCUmJsrlcmnTpk2+zy5fvqwlS5ZoxIgRuvbaa5WYmKh58+bp9OnTLfa5bNkyuVwuvzJ48OCgvwwAAK0yLChhFnSyrqurU3JyslatWnXVZ1999ZUOHjyoJ554QgcPHtSGDRt09OhR3X333a32O2zYMJ05c8ZXPvnkk2BDAwCgQwp6GjwjI0MZGRlNfhYbG6vt27f77Vu5cqXGjRun8vJy3Xzzzc0H0qWL4uPjgw0HAIDgtcPo2IyQX7Ourq6Wy+VSjx49Wqx37NgxJSYmKioqSqmpqcrLy2s2udfX1/u9U/XKNUXhTOmJo9o7BLTRgMXF7R1CROA8WcMR16yDcenSJS1ZskSZmZktrkySkpKigoICbd26Vfn5+SorK9Odd96pixcvNlk/Ly/Pbw3RpKSkUH0FAADaXciS9eXLl3XvvffKMAzl5+e3WDcjI0OzZ8/WyJEjlZ6eri1btujChQt69913m6yfm5ur6upqX6moqAjFVwAAdEQReINZSKbBGxP1yZMn9dvf/jbg9T4b9ejRQwMHDlRpaWmTn7vd7hYX+QYAoDlMg+ubRH3s2DHt2LFDN954Y9B91NbW6vjx40pISLA6PAAAIk7Qybq2tlYlJSUqKSmRJJWVlamkpETl5eW6fPmy/vmf/1n79+/Xr3/9a3k8HlVWVqqyslINDQ2+PqZMmaKVK1f6th999FHt3LlTJ06c0O7du3XPPfeoc+fOyszMNP8NAQD4NidMg+/fv1+TJ0/2befk5EiSsrKytGzZMr3//vuSpFGjRvm1+/jjjzVp0iRJ0vHjx3Xu3DnfZ6dOnVJmZqbOnz+vuLg43XHHHSouLlZcXFyw4QEA0KJInAYPOllPmjRJhtF8pC191ujEiRN+2+vWrQs2DAAAHIN3gwMAnIWFPAAAsDmSNQAA9haJ16xD+gYzAABgHiNrAICzMA0OAIC9uQxDrgCeXGqpfbgxDQ4AgM0xsgYAOAvT4AAA2Bt3gwMAAMsxsgYAOAvT4LC79MRR7R3CVUpXjLeknwGLiy3pB+Fn1W9A4neA1jENDgAALMfIGgDgLEyDAwBgb5E4DU6yBgA4SwSOrLlmDQCAzTGyBgA4TntMZZtBsgYAOIthfF3MtA8zpsEBALA5RtYAAEeJxLvBGVkDAJzFsKAEIS8vT7fffruio6PVq1cvzZo1S0ePHg2qD5I1AAAhtHPnTmVnZ6u4uFjbt2/X5cuXNXXqVNXV1QXcB9PgAABHcXm/LmbaB2Pr1q1+2wUFBerVq5cOHDig7373uwH1QbIGADiLRS9Fqamp8dvtdrvldrtbbV5dXS1JuuGGGwI+JNPgAAC0QVJSkmJjY30lLy+v1TZer1ePPPKIJk6cqOHDhwd8LEbWAABHsepu8IqKCsXExPj2BzKqzs7O1pEjR/TJJ58EdUySNQDAWSx6KUpMTIxfsm7NokWL9MEHH2jXrl266aabgjokyRoA4Cjhfs7aMAz9+7//uzZu3KiioiL169cv6GOSrB2mdMV4S/oZsLjYkn4Aid8TOrbs7Gy9/fbb2rx5s6Kjo1VZWSlJio2NVffu3QPqgxvMAADOEuaXouTn56u6ulqTJk1SQkKCr6xfvz7gPhhZAwAcpT2mwc1iZA0AgM0xsgYAOEsELpFJsgYAOIojVt3atWuXZsyYocTERLlcLm3atMnv8/vuu08ul8uvTJs2rdV+V61apb59+yoqKkopKSnau3dvsKEBANAhBZ2s6+rqlJycrFWrVjVbZ9q0aTpz5oyvvPPOOy32uX79euXk5Gjp0qU6ePCgkpOTlZ6errNnzwYbHgAALQvz3eBWCHoaPCMjQxkZGS3Wcbvdio+PD7jPl156SQsXLtT8+fMlSatXr9aHH36oNWvW6PHHHw82RAAAmuWIafBAFBUVqVevXho0aJAeeughnT9/vtm6DQ0NOnDggNLS0r4JqlMnpaWlac+ePU22qa+vV01NjV8BAKCjsjxZT5s2TW+++aYKCwv1/PPPa+fOncrIyJDH42my/rlz5+TxeNS7d2+//b179/a95eVKeXl5fiudJCUlWf01AAAdldcwX8LM8rvBf/zjH/v+PGLECI0cOVK33HKLioqKNGXKFEuOkZubq5ycHN92TU0NCRsAEBiL1rMOp5C/FKV///7q2bOnSktLm/y8Z8+e6ty5s6qqqvz2V1VVNXvd2+12+1Y7CXbVEwCAs7n0zXXrNpV2iDnkyfrUqVM6f/68EhISmvy8W7duGjNmjAoLC337vF6vCgsLlZqaGurwAACwvaCTdW1trUpKSlRSUiJJKisrU0lJicrLy1VbW6vHHntMxcXFOnHihAoLCzVz5kwNGDBA6enpvj6mTJmilStX+rZzcnL02muv6Y033tBnn32mhx56SHV1db67wwEAsEzjG8zMlDAL+pr1/v37NXnyZN9247XjrKws5efn6w9/+IPeeOMNXbhwQYmJiZo6daqeeeYZud1uX5vjx4/r3Llzvu05c+boiy++0JNPPqnKykqNGjVKW7duveqmMwAAzIrER7eCTtaTJk1qcQWRbdu2tdrHiRMnrtq3aNEiLVq0KNhwAADo8Hg3OADAWSLwbnCSNQDAUVyGIZeJ685m2rYVydphBiwubu8QrmLHmBBe206XWNbXLev/zZJ++F3CTkjWAABn8f6jmGkfZiRrAICjROI0eMhfigIAAMxhZA0AcBbuBgcAwObMvoWMu8EBAAitSHyDGdesAQCwOUbWAABnYRocAAB7c3m/LmbahxvT4AAA2BwjawCAszANDgCAzUXgc9ZMgwMAYHOMrAEAjhKJ7wYnWQMAnCUCr1kzDQ4AgM0xsgYAOIshc2tSs5CHORv/clgx0eYmC9ITR1kTTAdXumK8ZX0NWFxsWV+ITFb+vRsgfk9oGdesAQCwO0Mmr1lbFknAuGYNAIDNMbIGADgLd4MDAGBzXgtKkHbt2qUZM2YoMTFRLpdLmzZtCqo9yRoAgBCrq6tTcnKyVq1a1ab2TIMDABylPe4Gz8jIUEZGRpuPSbIGADiLRdesa2pq/Ha73W653W4zkTWLaXAAANogKSlJsbGxvpKXlxeyYzGyBgA4i0Uj64qKCsXExPh2h2pULZGsAQBOY1GyjomJ8UvWocQ0OAAANsfIGgDgLF5JLpPtg1RbW6vS0lLfdllZmUpKSnTDDTfo5ptvbrU9yRoA4Cjt8ejW/v37NXnyZN92Tk6OJCkrK0sFBQWttidZAwCcpR1eNzpp0iQZJo4Z9DXr1l6Z5nK5miwvvPBCs30uW7bsqvqDBw8O+ssAANARBT2ybnxl2v33368f/vCHV31+5swZv+2PPvpICxYs0I9+9KMW+x02bJh27NjxTWBdGPQDAELAa0guEyNrbwSsZ93aK9Pi4+P9tjdv3qzJkyerf//+LQfSpctVbQEAsByrbvmrqqrShx9+qAULFrRa99ixY0pMTFT//v01d+5clZeXN1u3vr5eNTU1fgUAgI4qpHPNb7zxhqKjo5ucLv+2lJQUFRQUaNCgQTpz5oyeeuop3XnnnTpy5Iiio6Ovqp+Xl6ennnrqqv3JG+5Xp6goUzEPULGp9nZXumJ8e4cQEaw6TwMWd+zfE2DV3xXvpUvS45st6at1JkfW6mAj6zVr1mju3LmKaiWBZmRkaPbs2Ro5cqTS09O1ZcsWXbhwQe+++26T9XNzc1VdXe0rFRUVoQgfANARNU6DmylhFrKR9e9+9zsdPXpU69evD7ptjx49NHDgQL8HyL8tlCubAABgNyEbWb/++usaM2aMkpOTg25bW1ur48ePKyEhIQSRAQAczWuYL2EWdLKura1VSUmJSkpKJH3zyrRv3xBWU1Oj9957Tw888ECTfUyZMkUrV670bT/66KPauXOnTpw4od27d+uee+5R586dlZmZGWx4AAC0zPCaL2EW9DR4IK9MW7dunQzDaDbZHj9+XOfOnfNtnzp1SpmZmTp//rzi4uJ0xx13qLi4WHFxccGGBwBAhxN0sg7klWkPPvigHnzwwWY/P3HihN/2unXrgg0DAIC2icDnrHlNGADAWbyGTD1+FQlvMAMAIKJF4Mg6pM9ZAwAA8xhZAwCcxZDJkbVlkQSMZA0AcBamwQEAgNUYWQMAnMXrlWTixSbeCHgpCgAAEY1pcAAAYDVG1gAAZ4nAkTXJGgDgLLzBDN9WumK8Jf0MWFxsST9W9mXVd7MrK885YJVtp0vaO4SrpCda08/fjcsqb72aY5GsAQCOYhheGSaWuTTTtq1I1gAAZzEMc1PZXLMGACDEDJPXrHl0CwAAXImRNQDAWbxeyWXiujPXrAEACDGmwQEAgNUYWQMAHMXwemWYmAbn0S0AAEKNaXAAAGA1RtYAAGfxGpIrskbWJGsAgLMYhiQzj24xDQ4AAK7AyBoA4CiG15BhYhrcYGQNAECIGV7zpQ1WrVqlvn37KioqSikpKdq7d2/AbUnWAABHMbyG6RKs9evXKycnR0uXLtXBgweVnJys9PR0nT17NqD2JGsAAELspZde0sKFCzV//nwNHTpUq1ev1jXXXKM1a9YE1L5DXLNuvH7gvXTJdF9/Ny6b7qORFfFI1sZkFau+m2TP7wfYUc3F8L85qzVW/f39u77uJxzXg/9u1JtajKMx1pqaGr/9brdbbrf7qvoNDQ06cOCAcnNzffs6deqktLQ07dmzJ7CDGh1ARUVF4+toKBQKhRLBpaKiImS54m9/+5sRHx9vSZzXXXfdVfuWLl3a5HE///xzQ5Kxe/duv/2PPfaYMW7cuIBi7xAj68TERFVUVCg6Oloul6vZejU1NUpKSlJFRYViYmLCGKE5xB1ekRq3FLmxE3d42TFuwzB08eJFJSYmhuwYUVFRKisrU0NDg+m+DMO4Kt80Naq2SodI1p06ddJNN90UcP2YmBjb/ECDQdzhFalxS5EbO3GHl93ijo2NDfkxoqKiFBUVFfLjfFvPnj3VuXNnVVVV+e2vqqpSfHx8QH1wgxkAACHUrVs3jRkzRoWFhb59Xq9XhYWFSk1NDaiPDjGyBgDAznJycpSVlaWxY8dq3Lhxevnll1VXV6f58+cH1N5Rydrtdmvp0qUhva4QCsQdXpEatxS5sRN3eEVq3JFszpw5+uKLL/Tkk0+qsrJSo0aN0tatW9W7d++A2rsMox3emwYAAALGNWsAAGyOZA0AgM2RrAEAsDmSNQAANtfhknWwS5C99957Gjx4sKKiojRixAht2bIlTJF+LS8vT7fffruio6PVq1cvzZo1S0ePHm2xTUFBgVwul18J90P+y5YtuyqGwYMHt9imvc+1JPXt2/equF0ul7Kzs5us357neteuXZoxY4YSExPlcrm0adMmv88Nw9CTTz6phIQEde/eXWlpaTp27Fir/ZpZps9s3JcvX9aSJUs0YsQIXXvttUpMTNS8efN0+vTpFvtsy+/Nyrgl6b777rsqhmnTprXab3ueb0lN/t5dLpdeeOGFZvsMx/lGcDpUsg52CbLdu3crMzNTCxYs0KFDhzRr1izNmjVLR44cCVvMO3fuVHZ2toqLi7V9+3ZdvnxZU6dOVV1dXYvtYmJidObMGV85efJkmCL+xrBhw/xi+OSTT5qta4dzLUn79u3zi3n79u2SpNmzZzfbpr3OdV1dnZKTk7Vq1aomP//5z3+uX/7yl1q9erU+/fRTXXvttUpPT9elFhZZMbtMn9m4v/rqKx08eFBPPPGEDh48qA0bNujo0aO6++67W+03mN+b1XE3mjZtml8M77zzTot9tvf5luQX75kzZ7RmzRq5XC796Ec/arHfUJ9vBMnEO9FtZ9y4cUZ2drZv2+PxGImJiUZeXl6T9e+9915j+vTpfvtSUlKMf/3Xfw1pnC05e/asIcnYuXNns3XWrl1rxMbGhi+oJixdutRITk4OuL4dz7VhGMbDDz9s3HLLLYbX623yczuca8MwDEnGxo0bfdter9eIj483XnjhBd++CxcuGG6323jnnXea7SfYvyNWx92UvXv3GpKMkydPNlsn2N+bWU3FnZWVZcycOTOofux4vmfOnGl873vfa7FOuM83WtdhRtaNS5ClpaX59rW2BNmePXv86ktSenp64EuWhUB1dbUk6YYbbmixXm1trfr06aOkpCTNnDlTf/zjH8MRnp9jx44pMTFR/fv319y5c1VeXt5sXTue64aGBr311lu6//77W1wAxg7n+kplZWWqrKz0O6exsbFKSUlp9py25e9IOFRXV8vlcqlHjx4t1gvm9xYqRUVF6tWrlwYNGqSHHnpI58+fb7auHc93VVWVPvzwQy1YsKDVunY43/hGh0nW586dk8fjueptML1791ZlZWWTbSorK4OqH2per1ePPPKIJk6cqOHDhzdbb9CgQVqzZo02b96st956S16vVxMmTNCpU6fCFmtKSooKCgq0detW5efnq6ysTHfeeacuXrzYZH27nWtJ2rRpky5cuKD77ruv2Tp2ONdNaTxvwZzTtvwdCbVLly5pyZIlyszMbHFBiWB/b6Ewbdo0vfnmmyosLNTzzz+vnTt3KiMjQx6Pp8n6djzfb7zxhqKjo/XDH/6wxXp2ON/w56jXjdpddna2jhw50uq1odTUVL+Xv0+YMEFDhgzRK6+8omeeeSbUYUqSMjIyfH8eOXKkUlJS1KdPH7377rsB/avdDl5//XVlZGS0uCSfHc51R3X58mXde++9MgxD+fn5Lda1w+/txz/+se/PI0aM0MiRI3XLLbeoqKhIU6ZMCUsMZq1Zs0Zz585t9SZJO5xv+OswI+u2LEEWHx9vaskyKy1atEgffPCBPv7446CW+5Skrl27avTo0SotLQ1RdK3r0aOHBg4c2GwMdjrXknTy5Ent2LFDDzzwQFDt7HCuJfnOWzDn1Ipl+qzSmKhPnjyp7du3B71MY2u/t3Do37+/evbs2WwMdjrfkvS73/1OR48eDfo3L9njfDtdh0nWbVmCLDU11a++JG3fvj3gJcusYBiGFi1apI0bN+q3v/2t+vXrF3QfHo9Hhw8fVkJCQggiDExtba2OHz/ebAx2ONfftnbtWvXq1UvTp08Pqp0dzrUk9evXT/Hx8X7ntKamRp9++mmz59SKZfqs0Jiojx07ph07dujGG28Muo/Wfm/hcOrUKZ0/f77ZGOxyvhu9/vrrGjNmjJKTk4Nua4fz7XjtfYebldatW2e43W6joKDA+NOf/mQ8+OCDRo8ePYzKykrDMAzjX/7lX4zHH3/cV//3v/+90aVLF+N//ud/jM8++8xYunSp0bVrV+Pw4cNhi/mhhx4yYmNjjaKiIuPMmTO+8tVXX/nqXBn3U089ZWzbts04fvy4ceDAAePHP/6xERUVZfzxj38MW9z/+Z//aRQVFRllZWXG73//eyMtLc3o2bOncfbs2SZjtsO5buTxeIybb77ZWLJkyVWf2elcX7x40Th06JBx6NAhQ5Lx0ksvGYcOHfLdNb18+XKjR48exubNm40//OEPxsyZM41+/foZf/vb33x9fO973zN+9atf+bZb+zsS6rgbGhqMu+++27jpppuMkpISv998fX19s3G39nsLddwXL140Hn30UWPPnj1GWVmZsWPHDuO2224zbr31VuPSpUvNxt3e57tRdXW1cc011xj5+flN9tEe5xvB6VDJ2jAM41e/+pVx8803G926dTPGjRtnFBcX+z676667jKysLL/67777rjFw4ECjW7duxrBhw4wPP/wwrPFKarKsXbu22bgfeeQR33fs3bu38U//9E/GwYMHwxr3nDlzjISEBKNbt27Gd77zHWPOnDlGaWlpszEbRvuf60bbtm0zJBlHjx696jM7neuPP/64yd9GY3xer9d44oknjN69extut9uYMmXKVd+pT58+xtKlS/32tfR3JNRxl5WVNfub//jjj5uNu7XfW6jj/uqrr4ypU6cacXFxRteuXY0+ffoYCxcuvCrp2u18N3rllVeM7t27GxcuXGiyj/Y43wgOS2QCAGBzHeaaNQAAHRXJGgAAmyNZAwBgcyRrAABsjmQNAIDNkawBALA5kjUAADZHsgYAwOZI1gAA2BzJGgAAmyNZAwBgcyRrAABs7v8DU8AinBtbNfYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Setup scene\n",
    "\n",
    "scene_file = \"../scenes/double_backscattering_reflector/double_backscattering_reflector.xml\"\n",
    "tx_position = np.array([-2.0, 0.0, 5.0])\n",
    "tx_orientation = np.array([0.0, 0.0, 0.0])\n",
    "tx_pattern = 'iso'\n",
    "fc = 3.5e9\n",
    "tx_slant_angle = 0.0\n",
    "\n",
    "# Path loss map used for evaluation\n",
    "cm_center = np.array([2., 0, 0.0])\n",
    "cm_orientation = np.array([0., 0., 0.])\n",
    "cm_size = np.array([2., 2.])\n",
    "cm_cell_size = np.array([0.1, 0.1])\n",
    "num_samples = int(1e7)\n",
    "max_depth = 2\n",
    "\n",
    "scene = mi.load_file(scene_file)\n",
    "params = mi.traverse(scene)\n",
    "\n",
    "params['mat-itu_metal.eta_r'] = 1.1\n",
    "params['mat-itu_metal.eta_i'] = 1.0\n",
    "params['mat-itu_metal.s'] = 0.7\n",
    "params['mat-itu_metal.lambda'] = 0.7\n",
    "params.update()\n",
    "\n",
    "# Define Loss and set trainable variables\n",
    "\n",
    "def func_loss(x):\n",
    "    return -1e3*dr.sum(x)\n",
    "\n",
    "variables =[params['mat-itu_metal.eta_r'],\n",
    "            params['mat-itu_metal.eta_i'],\n",
    "            params['mat-itu_metal.s'],\n",
    "            params['mat-itu_metal.lambda']\n",
    "           ]\n",
    "for v in variables:\n",
    "    dr.enable_grad(v)\n",
    "    \n",
    "# Compute gradients with RBP\n",
    "    \n",
    "# Reset grads\n",
    "for v in variables:\n",
    "    dr.set_grad(v, 0.0)\n",
    "\n",
    "rbp_tracer = PathlossMapRBPTracer(scene, fc, tx_pattern, tx_slant_angle,\n",
    "                                  cm_center, cm_orientation, cm_size, cm_cell_size,\n",
    "                                  num_samples, max_depth)\n",
    "\n",
    "cm_test,_ = rbp_tracer(tx_position, tx_orientation, func_loss, loop_record=True)\n",
    "g_test = np.array([dr.grad(v).numpy()[0] for v in variables])\n",
    "print(\"RBP Grad: \", g_test)\n",
    "\n",
    "# Compute gradients with vanilla AD\n",
    "\n",
    "# Reset grads\n",
    "for v in variables:\n",
    "    dr.set_grad(v, 0.0)\n",
    "\n",
    "ref_tracer = MapTracer(scene, fc, tx_pattern, tx_slant_angle,\n",
    "                       cm_center, cm_orientation, cm_size, cm_cell_size,\n",
    "                       num_samples, max_depth)\n",
    "cm_ref,_,_,_ = ref_tracer(tx_position, tx_orientation, loop_record=False)\n",
    "l = func_loss(cm_ref)\n",
    "dr.backward(l)\n",
    "g_ref = np.array([dr.grad(v).numpy()[0] for v in variables])\n",
    "print(\"Ref. Grad: \", g_ref)\n",
    "\n",
    "print(\"Relative Error: \", np.linalg.norm(g_ref-g_test)/np.linalg.norm(g_ref))\n",
    "\n",
    "cm_ref_db = 10.0*np.log10(cm_ref.numpy()).T\n",
    "cm_test_db = 10.0*np.log10(cm_test.numpy()).T\n",
    "ale = np.abs(cm_ref_db - cm_test_db)\n",
    "plt.figure()\n",
    "plt.imshow(ale)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "026304b1",
   "metadata": {},
   "source": [
    "# Floor-wall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "84a2a33a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RBP Grad:  [-2.3387974e-02 -4.9913162e-04  1.0319018e-01 -9.5513858e-02\n",
      " -2.5012636e-01 -5.4635536e-02 -6.8191779e-01]\n",
      "Ref. Grad:  [-2.3388032e-02 -4.9913180e-04  1.0319018e-01 -9.5513962e-02\n",
      " -2.5013155e-01 -5.4637268e-02 -6.8189418e-01]\n",
      "Relative Error:  3.2644723e-05\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_1863661/1668180513.py:68: RuntimeWarning: divide by zero encountered in log10\n",
      "  cm_ref_db = 10.0*np.log10(cm_ref.numpy()).T\n",
      "/tmp/ipykernel_1863661/1668180513.py:69: RuntimeWarning: divide by zero encountered in log10\n",
      "  cm_test_db = 10.0*np.log10(cm_test.numpy()).T\n",
      "/tmp/ipykernel_1863661/1668180513.py:70: RuntimeWarning: invalid value encountered in subtract\n",
      "  ale = np.abs(cm_ref_db - cm_test_db)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f98f42869d0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAGsCAYAAAD9ro91AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3IklEQVR4nO3df3QV9Z3/8dcNkASUXERKQiRYbBDwB6goEOxWVCoiy5GV46Eez0JZ6p52AweIq924Ciu2ja1HQWsErQrH3aWxbBc4uhbLxgJfV1AIzSla5UiXbbJCgvZsEsiaBHLn+wfl1rvkfiZ3PndyZ7jPxzlz2szcz8xnPncuHz+fmfe8I47jOAIAAIGVk+kKAAAAMzprAAACjs4aAICAo7MGACDg6KwBAAg4OmsAAAKOzhoAgICjswYAIODorAEACDg6awAAAo7OGgDQp3bv3q05c+aouLhYkUhEW7du9fV4//AP/6BIJJKwjBs3ztdjphudNQCgT7W3t2vixImqrq7us2NeeeWVOnbsWHx5++23++zY6dA/0xUAAGSXWbNmadasWUm3d3Z26u///u/105/+VC0tLbrqqqv0wx/+UNOnT/d8zP79+6uoqMhz+UxjZA0ACJQlS5Zoz549qqmp0W9+8xvdfffduv322/Xxxx973ufHH3+s4uJiXXbZZbr33nvV0NCQxhr7L0KKTABApkQiEW3ZskVz586VJDU0NOiyyy5TQ0ODiouL45+bMWOGJk+erB/84AcpH+MXv/iFTp48qbFjx+rYsWN69NFH9cknn+j999/X4MGD03UqvmIaHAAQGAcPHlR3d7cuv/zyhPWdnZ26+OKLJUkfffSRxo8fb9zPd7/7XT3++OOSlDDlPmHCBE2ZMkWXXnqpfvazn2nx4sVpPgN/0FkDAALj5MmT6tevn+rq6tSvX7+EbRdeeKEk6bLLLtOHH35o3M/Zjr0nQ4YM0eWXX67Dhw/bV7iP0FkDAALj2muvVXd3t44fP64/+7M/6/Ezubm5VqFXJ0+e1O9+9zv95V/+ped99DU6awBAnzp58mTCqPbIkSOqr6/X0KFDdfnll+vee+/VggUL9OSTT+raa6/Vp59+qtraWk2YMEGzZ89O+Xh/+7d/qzlz5ujSSy/V0aNHtWrVKvXr10/33HNPOk/LVzxgBgDoUzt37tTNN998zvqFCxdq48aNOnXqlL73ve/plVde0SeffKJhw4Zp6tSpevTRR3X11VenfLxvfOMb2r17t/7whz/oS1/6kr761a/q+9//vr7yla+k43T6BJ01AAABR5w1AAABR2cNAEDA+faAWXV1tZ544gk1NTVp4sSJ+vGPf6zJkye7lovFYjp69KgGDx6sSCTiV/UAAD5xHEcnTpxQcXGxcnL8GxN2dHSoq6vLej+5ubnKz89PQ4185PigpqbGyc3NdV5++WXngw8+cO677z5nyJAhTnNzs2vZxsZGRxILCwsLS8iXxsZGP7oYx3Ec5/PPP3eKhvdLSz2Lioqczz//3Le6poMvD5hNmTJFN9xwg5599llJZ0bLJSUlWrp0qf7u7/7OWLa1tVVDhgzRV3WH+mtAuquWPW640rx93wdJN/1hkXkG5OIN73mpUTBZtJNV2Uwx1dnmXG34dVyb9vfYTr7+drx+dxZcz+f99h7Xn+7u1NsHnlRLS4ui0agfVVNbW5ui0aiO1F2qgsHeR+9tJ2IaPen3am1tVUFBQRprmF5pnwbv6upSXV2dKisr4+tycnI0Y8YM7dmz55zPd3Z2qrOzM/73iRMn/lixAeofobP2rL/LlI6hbfvlmsueV9+LRTtZlc0UU51tztWGX8e1aX+P7eTrb8frd2fB9Xz6dxu398WtzILBOVaddVik/Qw/++wzdXd3q7CwMGF9YWGhmpqazvl8VVWVotFofCkpKUl3lQAA56luJ2a9hEHG/3OksrJSra2t8aWxsTHTVQIAhERMjvUSBmmfBh82bJj69eun5ubmhPXNzc09Jv7Oy8tTXl5euqsBAMgCMcVkMza2K9130t5Z5+bmatKkSaqtrY3nJ43FYqqtrdWSJUvSfTgks/c3xs2f/XVZ0m3DXjj32YIEUyd4qZFrnayOadj34TVTk24rXbHXa438Ox+b/box7dvr9+qyX9O1JknDTF+BX21hcT2Zyrr+doLIdD6/6fkBsrhk7eScsqgQeuJLnHVFRYUWLlyo66+/XpMnT9batWvV3t6uRYsW+XE4AECW6nYcdVsENdmU7Uu+dNbz58/Xp59+qpUrV6qpqUnXXHONtm/ffs5DZwAA2LC975y196zPWrJkCdPeAACkAfmsAQChFZOjbkbWAAAEV7ZMg2c8zhoAAJgxsgYAhFa2PA3uSyIPG2dfzj5dd55f76AGessm3tnF4fmDkm4rffV/kxf0M/bbcL6fTbjAl0P69S4Bt/oa45b9bGMTH2L+TzuntFPbfE2Ocbav+OjDQg22eDf4iRMxjRvfHPhEHkyDAwAQcEyDAwBCq9vyaXCbsn2JzhoAEFrdzpnFpnwY0FkDAEIr9sfFpnwYcM8aAICAY2QNAAitmCLqVsSqfBiErrNufaPUuD16x2HP+7ZKGxk2mUrRmAFB/F7d0kaaGOvsEnLkOSWoj+FkJjahTsY29ikd6DBlpp2shPz3HnPOLDblw4BpcAAAAi50I2sAAM7qtpwGtynbl+isAQChlS2dNdPgAAAEHCNrAEBoxZyIYo7F0+AWZfsSnTUAILSYBgcAAIEQupG1TRy1m9DFUvsUK2rFrzq57NevVIpGFnUyxg/bcPtePcbXu6Z+tIj99qudTHWyinO3iaUOYkxzyN+50K0cdVuMO7vTWBc/ha6zBgDgLMfynrXDPWsAAPzFPWsAABAIjKwBAKHV7eSo27G4Zx2Sd4PTWQMAQiumiGIWk8QxhaO3ZhocAICAC93Iekdss3H713Pu7qOa9A1jekeL9IE24RoZSTnpUqdhpsyPhnP1LVWlLEJ8bMJl3ELnTOFZAUwlasXQFr6dj1v7ZyK0MQThVzay5QGz0HXWAACcZX/PmmlwAACQBoysAQChdeYBM4tEHkyDAwDgr5jl60Z5GhwAAKQFnTUAILTOPmBms6SiqqpKN9xwgwYPHqzhw4dr7ty5OnTokGu5zZs3a9y4ccrPz9fVV1+tN954I6Xjhm4a/HwLzXLjOcTEJkTEhW918jNkKQmrEJ4AZhizOR+bLFUt45NPJQ750HxP0FR22AverwljO5lC/VwcXjM16bbSV//X+46DGGLl8Vrs7uqQNmzzo0bniCmnT1+KsmvXLpWXl+uGG27Q6dOn9dBDD+m2227Tb3/7W11wQc/t8c477+iee+5RVVWV/vzP/1ybNm3S3LlzdeDAAV111VW9Om7oOmsAAM7qdiLqtsiclWrZ7du3J/y9ceNGDR8+XHV1dfra177WY5mnn35at99+ux544AFJ0mOPPaYdO3bo2Wef1fr163t1XKbBAQBZr62tLWHp7OzsVbnW1lZJ0tChQ5N+Zs+ePZoxY0bCupkzZ2rPnt7PhNFZAwBCq/uPT4PbLJJUUlKiaDQaX6qqqlyPHYvFtHz5ct14443G6eympiYVFhYmrCssLFRTU1Ovz5NpcABAaMWcHMUs3mAW++MbzBobG1VQUBBfn5eX51q2vLxc77//vt5++23Px+8tOmsAQNYrKChI6KzdLFmyRK+//rp2796tkSNHGj9bVFSk5ubmhHXNzc0qKirq9fGYBgcAhFa6psF7y3EcLVmyRFu2bNFbb72l0aNHu5YpKytTbW1twrodO3aorKz3iYQYWQMAQium1J/o/r/lU1FeXq5NmzZp27ZtGjx4cPy+czQa1cCBAyVJCxYs0CWXXBK/771s2TLddNNNevLJJzV79mzV1NRo//79euGFF3p9XDrr81WG0iz6Ui7T+/bhmKZYXbe4ZJt4aGMctkVKyRbD+RhTuUoaZmhGm7SdA+Z+mnyj27+RhrZw+36MLOLrjWX92q9L2WTfwWnnlPmYIbZu3TpJ0vTp0xPWb9iwQd/85jclSQ0NDcrJ+dOIfdq0adq0aZMefvhhPfTQQxozZoy2bt3a6xhric4aABBi9i9FSX0a3M3OnTvPWXf33Xfr7ru9v9SLzhoAEFr2+azD8ehWOGoJAEAWY2QNAAgt8lkDABBw2TINTmcNAAgtL7HS/7d8GNBZn6dcQ3hM4TQZSFXpyq86WezXrY1NSld4C6GyPa7n78elnDE1pEUb26SyjN5x2HthU50nWLS/Da/fnc1vJ4hpO7MUnTUAILRiTkQxm5eiWJTtS3TWAIDQillOg9vEaPelcNQSAIAsxsgaABBa9ikywzFmpbMGAIRWtyLqtoiVtinbl8LxnxQAAGQxRtYAgNBiGhyh5paW0LeUen6Us+VTKsvSVy3a2IJNjLwx5aTbNeOVRez9ZxMuSLrNrb42ZU28phmV7FJ+ZuS7C4Fu2U1ld6evKr4Kx39SAACQxRhZAwBCi2lwAAACjkQeAAAEnGOZItMhdAsAAKQDI2sAQGgxDZ7E7t279cQTT6iurk7Hjh3Tli1bNHfu3Ph2x3G0atUq/eQnP1FLS4tuvPFGrVu3TmPGjEnpOH9YNFn9cvPPWe8W3pBVLFI/WoWB+JUG0yc24TKlKww5Gl3awRz2lTylpCnkyM0wWXw3mQjnc9m3ze/d1BaH5w8yljWm/DRw++6Mvy2XNjaVNYapuV0TIU+DmS1Zt1L+T4r29nZNnDhR1dXVPW7/0Y9+pGeeeUbr16/Xu+++qwsuuEAzZ85UR0eHdWUBAMhGKY+sZ82apVmzZvW4zXEcrV27Vg8//LDuvPNOSdIrr7yiwsJCbd26Vd/4xjfsagsAwBd0W6bItCnbl9JayyNHjqipqUkzZsyIr4tGo5oyZYr27Ol5Oquzs1NtbW0JCwAAvXF2GtxmCYO0dtZNTU2SpMLCwoT1hYWF8W3/V1VVlaLRaHwpKSlJZ5UAAAi9jI//Kysr1draGl8aGxszXSUAQEjElGO9hEFaQ7eKiookSc3NzRoxYkR8fXNzs6655poey+Tl5SkvLy+d1QAAZIluJ6Jui6lsm7J9Ka2d9ejRo1VUVKTa2tp459zW1qZ3331X3/nOd1La18Ub3lP/yIB0Vg9fYJM5yCrExy8ew39MYV2SXTYjY9iXyQRznYxcwnCMYTwWoYC+ZdYyNaHFMYd86PIPtMdwppbxjssnMpAJzE/JvoPTHdK+bX1bl/Ncyp31yZMndfjw4fjfR44cUX19vYYOHapRo0Zp+fLl+t73vqcxY8Zo9OjReuSRR1RcXJwQiw0AQDpkS5x1yp31/v37dfPNN8f/rqiokCQtXLhQGzdu1IMPPqj29nb99V//tVpaWvTVr35V27dvV37+uS84AQDAhmOZdcs5X99gNn36dDlO8qmeSCSi1atXa/Xq1VYVAwDATbci6rZIxmFTti+F4z8pAADIYiTyAACEVsyxu+8cc3smMCDorAEAoRWzvGdtU7YvhaOWAABkMUbWYeZXSkOb9Jo+xXua0k1K5rhZU6yun/GpXtOQusXb2qTQNDJ8727tb0wp6Rb77TEc3Y0pDeaQD81lvV7jbrH1bnH9Rh5/71ZpO90kq5Nzyvs+UxRTRDGLh8RsyvYlOmsAQGhlyxvMmAYHACDgGFkDAEIrWx4wo7MGAIRWTJavGw3JPetw/CcFAABZjJE1ACC0HMunwZ2QjKzprM9XHlP89UYm0vF5TjfpIz/Taxr3a9P+HsN/bEKS3EKzvIYzubWvMZzPpQ1NdXILYzMpfdVQZ5s0pIayxrSovTmuSQBSZJJ1CwCAgMuWB8zCUUsAALIYI2sAQGgxDQ4AQMBly+tGmQYHACDgGFkDAEKLaXAAAAKOzhqh5hoDnIFYaTc2qTe9pqO0STfpV51ceYy3dWNqC7dYXZvzaRnvJN1mipW22e8wlxSy5vMxtJPb78otda0fXK4JU9y4MfWpad99mCIzW9BZAwBCi5E1AAABly2dNU+DAwAQcIysAQCh5cguVjr5kw3BQmcNAAitbJkGp7MGAIQWnTXSxiYkySur/dqEl1iEDtnU2WtZt/SNgWRoY5u0naY2dEsLaUyh6XI9+ZX+tNRmt4Y621ynxlDBCebvzhSKVmoKrXP5TbqGZyEQ6KwBAKHFyBoAgIDLls6a0C0AAAKOkTUAILQcJyLHYnRsU7Yv0VkDAEKLfNYAACAQsmpkHcZMVJnglonKczu5hYRZhH35lYnKtF/XdvKYico1E5ghxMcqE5gh1Mk1vMemnQxhR8ZMYC7t27oyeZ2jdxw2lvV8zbhl8zJ9P66ZwAwbba5xE7f9Jqvz6Q5p37b016cH2fKAWVZ11gCA80u23LNmGhwAgIBjZA0ACC2mwQEACLhsmQanswYAhJZjObIOS2fNPWsAAFKwe/duzZkzR8XFxYpEItq6davx8zt37lQkEjlnaWpq6vUx6awBAKHlSHIci8XDMdvb2zVx4kRVV1enVO7QoUM6duxYfBk+fHivyzIN3gd8i9/2KQbYtb4eU2j6GT/sV1zs4fmDkm4rfdUljtpUJ1MKRpf4Ydd2NDClWRzmU2pUmzSkphhsN66x1AZe09q6/u5s0nb6FUttIdn5dnf1k/b1TR1iiijSx28wmzVrlmbNmpVyueHDh2vIkCEpl5MYWQMAoLa2toSls7Mz7ce45pprNGLECH3961/Xf/zHf6RUls4aABBaZ58Gt1kkqaSkRNFoNL5UVVWlrY4jRozQ+vXr9fOf/1w///nPVVJSounTp+vAgQO93gfT4ACA0Io5EUXSEGfd2NiogoKC+Pq8vDzrup01duxYjR07Nv73tGnT9Lvf/U5r1qzRP/7jP/ZqH3TWAICsV1BQkNBZ+23y5Ml6++23e/15OmsAQGidfarbpnwm1NfXa8SIEb3+PJ01ACC0MvEGs5MnT+rw4T9FGxw5ckT19fUaOnSoRo0apcrKSn3yySd65ZVXJElr167V6NGjdeWVV6qjo0Mvvvii3nrrLf3yl7/s9TGzqrM+31Jg2qQPtAoD8VjWNWzFEDp0eM1UY1FjCkdTfV3OpdRUZ7dQJ49pO21SuVq1k4FN2J0b4/kYQufcmL47v9rYNZzPJjzO637dfq8W/xYk+3fm9OkOz/sMg/379+vmm2+O/11RUSFJWrhwoTZu3Khjx46poaEhvr2rq0v333+/PvnkEw0aNEgTJkzQv//7vyfsw01WddYAgPNLJkbW06dPl2OYP9+4cWPC3w8++KAefPDBlI/zRXTWAIDQStfT4EFHZw0ACK2wPmCWKl6KAgBAwDGyBgCE1pmRtc096zRWxkd01gCA0MrEA2aZwDQ4AAABx8i6D3hNt2d1TLdUfabUgxZxl27xqyam2HCv8cGu3OJeLWK0jd+7IQbYLUbeGOe7wiWY3WOcb90/rDNun3nXAk/7tTHkQ5cRkeFcTalCJanFJV49KZv3F7h8N57/HbG5xkPAkbec1F8sHwZ01gCA0GIaHAAABAIjawBAeGXJPDidNQAgvCynwRWSaXA6awBAaPEGMwAAEAgpjayrqqr0r//6r/roo480cOBATZs2TT/84Q81duzY+Gc6Ojp0//33q6amRp2dnZo5c6aee+45FRYWpr3yYXG+peY0sTpXmzR/BlbhZBYhbsYQLItztUn9aGKq71de/baxbKk8piiVuc6mlJOuaTsNxzWmPpXM6VoNaTt9S+Uql2vRr9SbbpLV2TnVZ1XgafAe7Nq1S+Xl5dq7d6927NihU6dO6bbbblN7+59+UCtWrNBrr72mzZs3a9euXTp69KjuuuuutFccAAA5EfslBFIaWW/fvj3h740bN2r48OGqq6vT1772NbW2tuqll17Spk2bdMstt0iSNmzYoPHjx2vv3r2aOtXjiwYAAMhiVvesW1tbJUlDhw6VJNXV1enUqVOaMWNG/DPjxo3TqFGjtGdPz9OjnZ2damtrS1gAAOiNsw+Y2Sxh4LmzjsViWr58uW688UZdddVVkqSmpibl5uZqyJAhCZ8tLCxUU1NTj/upqqpSNBqNLyUlJV6rBADINk4alhDw3FmXl5fr/fffV01NjVUFKisr1draGl8aGxut9gcAwPnGU5z1kiVL9Prrr2v37t0aOXJkfH1RUZG6urrU0tKSMLpubm5WUVFRj/vKy8tTXl6el2oAALJctjwNnlJn7TiOli5dqi1btmjnzp0aPXp0wvZJkyZpwIABqq2t1bx58yRJhw4dUkNDg8rKvIeUoGfG0CAfs/943rfNfi3KumWx8rpfm3AZY9jRBPNvxXQ+rufq8buzCXVyaydTBqxhhurafK9uIW6mEMQhLt9PRlhcp5nICph2IZnKtpFSZ11eXq5NmzZp27ZtGjx4cPw+dDQa1cCBAxWNRrV48WJVVFRo6NChKigo0NKlS1VWVsaT4AAAeJRSZ71u3ZmcttOnT09Yv2HDBn3zm9+UJK1Zs0Y5OTmaN29ewktRAABIN6bBe+D04hn3/Px8VVdXq7q62nOlAADoFbJuAQAQdJE/Ljblg49EHgAABBwjawBAeDENDgBAwNFZo7dsYjZ941estA23Y/qV5s8mVtqw3S1FozEe2nRNuNSpdWXyNIvR1cnTN0rma9UU71y6whxobWoLt+u/1GsaUovvzi1G22uqUdf9mtrJ1A4ujLHSLnWyuRYz8u9IlqKzBgCEl22ay/MxdAsAgCCxzZx13mfdAgAAfYORNQAgvHjADACAgMuSe9ZMgwMAEHCMrNPALTTFrxR0xjAQt/SBNukdLcJEPMtU2JeBTTpKmzZ0C88yGTD306TbSu84nHSbWyiTMc2lW6W8hmdZhA3ZnM+QD5OPxFzD+WzCOE2haIb9uoahGVJ+hiFFZsQ5s9iUDwM6awBAeHHPGgCAgOOeNQAACAJG1gCA8GIaHACAgMuSzpppcAAAAo6RNQAgvLJkZB26zjqQ6Shd+FUn31LbuZQ1xgFb7NcmfaBbfGvS/drEhfuVHtBiv4fXTDVuL11t2Gj4fkxxx5I59tiN5/cQWFzj7t+795SfRjZx4x6vC9f6ZuAdBWnF0+AAACAIQjeyBgDgLN5gBgBA0GXJPWumwQEACDg6awAAAo5pcABAaEVkec86bTXxV+g66yCGZvnJa1iLa6q+vYaNGQpJMtXJNWTPEIpjbAu/ztWGRShN6QrTF+t936Wv/q/5A361o08pMt3Ktsw3hMBlIkWsnyzCLQPx+yF0CwAABEHoRtYAAMRlydPgdNYAgPDKks6aaXAAAAKOkTUAILR4gxkAAEGXJdPgdNZf4Dn7j4+8HtfX+hrCOQ7PH5R0m1tYkV/tb8zYlaHQFL+yfWUqK92ff/A/SbdtfPoOY1lzRi9D9iuXKDVTBjK3UDTT9oyFAhquVWPGOousW1YhoEgrOmsAQHgxsgYAINiy5Z41T4MDABBwjKwBAOGVJa8bpbMGAIQX96wBAAg27lkDAIBAYGT9BdmWftMzQyxpqSmm2YIpZlZyieG2iUu2OR/TcU1xsRYx2K5lPabIdIu3ff3K5L+dYVPNdRpmCE1uXWmIh/6N+VyGfJj8XqTpfQCSOc7a2MYW6U3dYrTN7zDw59+uUPybyDQ4AAABZzkNHpbOmmlwAAACjpE1ACC8mAYHACDgsqSzZhocAIAU7N69W3PmzFFxcbEikYi2bt3qWmbnzp267rrrlJeXp9LSUm3cuDGlY9JZAwBC62yctc2Sqvb2dk2cOFHV1dW9+vyRI0c0e/Zs3Xzzzaqvr9fy5cv1rW99S2+++Wavj8k0eLaySItnYhV2ZAgTGeZ5rzKeq036QL/SIbqmJTSEk7mVNaejNDGXG+ZTyFJ0tT+hgG7pWk1sUm+aztU1PNFt30m4hifapGvNUrNmzdKsWbN6/fn169dr9OjRevLJJyVJ48eP19tvv601a9Zo5syZvdoHI2sAQNZra2tLWDo7O9O27z179mjGjBkJ62bOnKk9e3ofx05nDQAILycNi6SSkhJFo9H4UlVVlbYqNjU1qbCwMGFdYWGh2tra9Pnnn/dqH0yDAwBCK13vBm9sbFRBQUF8fV5enmXN0ovOGgAQbmkIvyooKEjorNOpqKhIzc3NCeuam5tVUFCggQMH9mofTIMDAOCjsrIy1dbWJqzbsWOHysrMD/99EZ01ACC80nTPOhUnT55UfX296uvrJZ0Jzaqvr1dDQ4MkqbKyUgsWLIh//tvf/rb+8z//Uw8++KA++ugjPffcc/rZz36mFStW9PqYTIMDAEIrE/ms9+/fr5tvvjn+d0VFhSRp4cKF2rhxo44dOxbvuCVp9OjR+rd/+zetWLFCTz/9tEaOHKkXX3yx12FbEp31ecsmttI17tIUX2yKt/UxfaCJMR2lTyk9JfN3YGpjt/SNLeO9p0o0xUOb2sktPtsqNt+n78AqftjQTqYYbddYacO5usZ+e/xtubWD1e/Dp3cNBN306dPlOMl/Ez29nWz69On69a9/7fmYdNYAgPDKkneD01kDAEIrE9PgmcADZgAABBwjawBAeDENDgBAwGVJZ800OAAAAZfSyHrdunVat26d/uu//kuSdOWVV2rlypXxVGEdHR26//77VVNTo87OTs2cOVPPPffcOS8w91PrG6VJt0XvONxn9egLXkODJGUm5MLPY1qErhhZ1NnrcW1CeNxC9kyM4XwvmNvBr7SRmUhRarNvt3M1huXNN4d9DfkwknSba6pXA2N4lls7JPt+TndI+7Z5rlMqeMCsByNHjtTjjz+uuro67d+/X7fccovuvPNOffDBB5KkFStW6LXXXtPmzZu1a9cuHT16VHfddZcvFQcAIBNvMMuElEbWc+bMSfj7+9//vtatW6e9e/dq5MiReumll7Rp0ybdcsstkqQNGzZo/Pjx2rt3r6ZONf9XIwAAKeOetVl3d7dqamrU3t6usrIy1dXV6dSpUwkJtseNG6dRo0YZE2x3dnaek/QbAAD8Scqd9cGDB3XhhRcqLy9P3/72t7VlyxZdccUVampqUm5uroYMGZLw+cLCQjU1NSXdX1VVVULC75KSkpRPAgCQnc7es7ZZwiDlznrs2LGqr6/Xu+++q+985ztauHChfvvb33quQGVlpVpbW+NLY2Oj530BALIM96x7lpubq9LSM09cT5o0Sfv27dPTTz+t+fPnq6urSy0tLQmj6+bmZhUVFSXdX15envLy8lKvOQAAWcL6pSixWEydnZ2aNGmSBgwYoNraWs2bN0+SdOjQITU0NKSUYNvW+RaeZWIM13DLcOVXBizDNtdMYDbn41MYjzEkyS3EymNIkmvmJkN4kFsIjym08fD4YckLWoQVufJ6LbqUs8oEZhPuZ+B6zfjBpp0mePzNOqfcapU22RK6lVJnXVlZqVmzZmnUqFE6ceKENm3apJ07d+rNN99UNBrV4sWLVVFRoaFDh6qgoEBLly5VWVkZT4IDAPyRJU+Dp9RZHz9+XAsWLNCxY8cUjUY1YcIEvfnmm/r6178uSVqzZo1ycnI0b968hJeiAAAA71LqrF966SXj9vz8fFVXV6u6utqqUgAA9AojawAAgi3yx8WmfBiQyAMAgIBjZA0ACC+mwQEACDZCtzLsD4smq19u/jnrbVLBIT28puZ0jV21SIfoOV2oy36N8cM+xX5bxeK61OnU1uQxtaWG35YpPluSoiuSv9/Areynh5LHd5tiyt3iqG3i9k2pLI1pML2mlHQ5pvVxDazebxAEWTKy5p41AAABF9iRNQAAvRKS0bENOmsAQGhlyz1rpsEBAAg4RtYAgPDKkgfM6KwBAKGVLdPgge2sL97wnvpHBqRczhQmEsb0mcaQJFPIhUUIiWtZl7R5nvmV5tIUhiNzaIpVqkS/wl4svrthpvM17NcU8nWmbPI2jq42F9XKz5Ju+mzCl5JuswnjbF1pCIOSVHqHP+lNTaGApa+6XGs+/T4QDoHtrAEAcMU0OAAAwZYt0+A8DQ4AQMAxsgYAhBfT4AAABBydNQAAwcY9awAAEAjn3cg6iLHUnmOl3bZnKH2dVeyxgU07mdikFjTFzRrTZ8p7O5nawXW/LteEW1rJpMf0MTVtdLWpzv6c67DV5u/G67XoFrdvZHEtGtOq2rxzwU2ysqc7pH3bvO83FUyDAwAQbBHHUcTx3uPalO1LTIMDABBwjKwBAOHFNDgAAMHG0+AAACAQGFkDAMKLaXCki19hL17DcCTLOplCQQxhIDb1dQsvMZ2PMQzHJdTGr7Av037dvhvj+biEi/l2LdqE3Xm8nlzTgRqimdwMWJk81a5+Y5Fe1nQ+Lte4W6ig1/2a6mwVRthHmAYHAACBwMgaABBeTIMDABBs2TINTmcNAAivLBlZc88aAICAY2QNAAi1sExl26CzDjFj2IRbCImJTRYew3HdwqQ8h/C4bLcKV7JoC2MmJAO3cBkjn773w/MHGYuWvmoRwmMRnuWX6Grz+Sblcr2Y2tEYJiiLDHwWbej5t+Oc8nzM1I/lnFlsyocA0+AAAAQcI2sAQGjxNDgAAEHH0+AAACAIGFkDAEIrEjuz2JQPAzprAEB4MQ0OAACCgJF1iJlSTrqmB7SJy7SJw/aLT/G4xja2iRu3YIp9NaXllMyx36ZzdYuj9q2dLFKuWr2HwOt35/LbMMbeW7xLwOodBRmKZU8XngYHACDosuSlKHTWAIDQypaRNfesAQAIOEbWAIDwypKnwemsAQChxTQ4AAAIhMCOrP+waLL65eafs94q3eF5xnPKPCkzoR42ISIuZU0hS8bUgy77tWpjj1yvccNx3dJytr5Rmvy4q5OHOrmGSRnq7Jby0zXM0MMxXbl8d8ZQNFNImAtT+0fv8P77MLWxryligxD2xdPgAAAEG9PgAAAgEBhZAwDCi6fBAQAINqbBAQBAj6qrq/XlL39Z+fn5mjJlit57772kn924caMikUjCkp9/7gPUJnTWAIDwijn2S4peffVVVVRUaNWqVTpw4IAmTpyomTNn6vjx40nLFBQU6NixY/Hl97//fUrHpLMGAISXk4YlRU899ZTuu+8+LVq0SFdccYXWr1+vQYMG6eWXX05aJhKJqKioKL4UFhamdMzA3rO+eMN76h8ZkOlqhJdNmku/ytrEJbvUySaW2sQYv+oWb+tTHLZNOsroasNGQzu5prk08K2dghgD7HLM6B3+HNYm9tuXfwtOd0j7tnmvUwoisrxn/cf/bWtrS1ifl5envLy8cz7f1dWluro6VVZWxtfl5ORoxowZ2rMneUz7yZMndemllyoWi+m6667TD37wA1155ZW9ricjawBA1ispKVE0Go0vVVVVPX7us88+U3d39zkj48LCQjU1NfVYZuzYsXr55Ze1bds2/dM//ZNisZimTZum//7v/+51/QI7sgYAwFWa3mDW2NiogoKC+OqeRtVelZWVqazsT7N006ZN0/jx4/X888/rscce69U+6KwBAKGVrtCtgoKChM46mWHDhqlfv35qbm5OWN/c3KyioqJeHXPAgAG69tprdfjw4V7Xk2lwAAB6KTc3V5MmTVJtbW18XSwWU21tbcLo2aS7u1sHDx7UiBEjen1cRtYAgPDKwBvMKioqtHDhQl1//fWaPHmy1q5dq/b2di1atEiStGDBAl1yySXx+96rV6/W1KlTVVpaqpaWFj3xxBP6/e9/r29961u9PiadNQAgtCKOo4jFPWsvZefPn69PP/1UK1euVFNTk6655hpt3749/tBZQ0ODcnL+NHH9P//zP7rvvvvU1NSkiy66SJMmTdI777yjK664IpV6ej/Lxx9/XJWVlVq2bJnWrl0rSero6ND999+vmpoadXZ2aubMmXruued6HVPW1tamaDSq6bqT0C03PoUGWclUyjyb8BM/juknw/mYUoVKFiFuNiklLVJ+murkmnozAylkD88fZNytWwpTk4ykwfT42zntnNJObVNra2uv7gN7cbav+LPpq9S/f2pvA/ui06c79P92PuprXdPB8z3rffv26fnnn9eECYlf9IoVK/Taa69p8+bN2rVrl44ePaq77rrLuqIAAJwjloYlBDx11idPntS9996rn/zkJ7rooovi61tbW/XSSy/pqaee0i233KJJkyZpw4YNeuedd7R3r/f/ogQAoCdnp8FtljDw1FmXl5dr9uzZmjFjRsL6uro6nTp1KmH9uHHjNGrUqKRvduns7FRbW1vCAgAA/iTlB8xqamp04MAB7du375xtTU1Nys3N1ZAhQxLWm97sUlVVpUcffTTVagAAkDX5rFMaWTc2NmrZsmX653/+55TTeyVTWVmp1tbW+NLY2JiW/QIAssDZN5jZLCGQ0si6rq5Ox48f13XXXRdf193drd27d+vZZ5/Vm2++qa6uLrW0tCSMrk1vdkn2snQAANyk6w1mQZdSZ33rrbfq4MGDCesWLVqkcePG6bvf/a5KSko0YMAA1dbWat68eZKkQ4cOqaGhoddvdskk30Ij/GITkuRTqFPo2tBHxnAmQ5Yk1/Afw7OaXxr7mblSe3v/esPEcuZrYpjp+VGfwqTcMk2ZwthsQqhMXPdrOB/T9SJJLeOT9yrDPB7T7bjG79W07z7MupUtUuqsBw8erKuuuiph3QUXXKCLL744vn7x4sWqqKjQ0KFDVVBQoKVLl6qsrExTp5rjPwEASFmaEnkEXdrfYLZmzRrl5ORo3rx5CS9FAQAg3SKxM4tN+TCw7qx37tyZ8Hd+fr6qq6tVXV1tu2sAACDeDQ4ACDOmwQEACDjirAEAQBAwsgYAhFYmUmRmAp31F3iNA7ZK1WfDIkWjMbZS3vfrW6o+N15jw308prEdDWWHTHB5J4GhztHV5qImxrhkU2pNydwWbt+NT6lGh3wY8V7Y43sIXFOUGuKwXWOaXf6dScqt/Q3Xm+u/bS6x7n0iS+5ZMw0OAEDAMbIGAISXI7uc1OEYWNNZAwDCi3vWAAAEnSPLe9Zpq4mvuGcNAEDAMbIGAIRXljwNTmedBhlL/egxtaAUzHAym/PxizFtoUuIlSmsxZhK1CYcxiK9qU2oU+sbpUm3Re/wmJZTsrvG/WpHw3FtUmS6hgKafrN+7ddFsuu4u6uftM/zblMTk2QRpWf1cFofYhocAICAY2QNAAgtngYHACDosuSeNdPgAAAEHCNrAEB4ZcnIms4aABBeWdJZMw0OAEDAMbLuJWNcrI9x1p7jcW3SEvpV1qeUnpLUMj95akJT/LDbd2far1VMrYmPKSXdUh56ZRVL7ZVLOxl/O27pKC2O6xvD9271PgCLf7+SlT3tnPK8z5RlSZw1nTUAILQI3QIAIOi4Zw0AAIKAkTUAILxijhSxGB3HwjGyprMGAIQX0+AAACAIGFn3UqbSYPp2XL/CT7yGl0gapuRl3dphmMc6uTGGZ7nt19DGpnP1k6kdjWFdNmFSLt+d6bowhVi5haFZpcg0yFQYp0nGUt5mKowtgeXIWuEYWdNZAwDCi2lwAAAQBIysAQDhFXNkNZXN0+AAAPjMiZ1ZbMqHANPgAAAEHCNrAEB4ZckDZnTWSJ3HzFquoTQ2YSA2WcQysV+fyrqGMxlCfIzhPy4hPJ5DwmS+LjxnnZN8ywBnFSZlUSdjiJtNKKCpToEIzXLBPWsAAAIuS0bW3LMGACDgGFkDAMLLkeXIOm018RWdNQAgvJgGBwAAQcDIGgAQXrGYJIsXm8TC8VIUOmsAQHhlyTQ4nfX5ys/Udn7FHgeRTaxuBtrCNQY4A/HoNqkqreK3M5SG1DOX9jelC82q32SWorMGAIQXI2sAAAIuS95gxtPgAAAEHCNrAEBoOU5MjkWaS5uyfYnOGgAQXo5jN5XNPWsAAHzmWN6zprNGb9ikNDRyC9ewCfUIYChI6NIH2rR/ENOBZiIkzCZMzS82xwxgG3uu0+kOad8273XCOeisAQDhFYtJEYv7ztyzBgDAZ1kyDU7oFgAAAcfIGgAQWk4sJsdiGpzQLQAA/MY0OAAACAJG1gCA8Io5UuT8H1nTWWeY5zhqmWO0XdMS+hQrbayTxbm6MZ6vTZpLC761heF8rOL2fWoL1zp5/O6s9usmQ9eMkcc6md5BIPnUTs4p7/tMleNIsgndCkdnzTQ4AAABx8gaABBaTsyRYzEN7jCyBgDAZ07MfvGgurpaX/7yl5Wfn68pU6bovffeM35+8+bNGjdunPLz83X11VfrjTfeSOl4dNYAgNByYo71kqpXX31VFRUVWrVqlQ4cOKCJEydq5syZOn78eI+ff+edd3TPPfdo8eLF+vWvf625c+dq7ty5ev/993t9TDprAABS8NRTT+m+++7TokWLdMUVV2j9+vUaNGiQXn755R4///TTT+v222/XAw88oPHjx+uxxx7Tddddp2effbbXxwzcPeuz9w9O65RVnHs26O7qSLrt9Onk2yT59rSmsU5+PiFqOl/TcX1sp0y0hemYrsf12oYuZbu7+pmLejyu67m6fbcmNteMH8d0k4H2N+5TZ8r1xf3g006nVTKOs3Vta2tLWJ+Xl6e8vLxzPt/V1aW6ujpVVlbG1+Xk5GjGjBnas6fnaIs9e/aooqIiYd3MmTO1devW3lfUCZjGxsazr6NhYWFhYQnx0tjY6Ftf8fnnnztFRUVpqeeFF154zrpVq1b1eNxPPvnEkeS88847CesfeOABZ/LkyT2WGTBggLNp06aEddXV1c7w4cN7fb6BG1kXFxersbFRgwcPViQSUVtbm0pKStTY2KiCgoJMVy+waKfeoZ16h3bqHdqpZ47j6MSJEyouLvbtGPn5+Tpy5Ii6urqs9+U4jiKRSMK6nkbVmRS4zjonJ0cjR448Z31BQQE/hl6gnXqHduod2ql3aKdzRaNR34+Rn5+v/Px834/zRcOGDVO/fv3U3NycsL65uVlFRUU9likqKkrp8z3hATMAAHopNzdXkyZNUm1tbXxdLBZTbW2tysp6fqteWVlZwuclaceOHUk/35PAjawBAAiyiooKLVy4UNdff70mT56stWvXqr29XYsWLZIkLViwQJdccomqqqokScuWLdNNN92kJ598UrNnz1ZNTY3279+vF154odfHDHxnnZeXp1WrVgXu/kHQ0E69Qzv1Du3UO7RTdpo/f74+/fRTrVy5Uk1NTbrmmmu0fft2FRYWSpIaGhqUk/Onietp06Zp06ZNevjhh/XQQw9pzJgx2rp1q6666qpeHzPiOCF51xoAAFmKe9YAAAQcnTUAAAFHZw0AQMDRWQMAEHCB76xTTUN2vtu9e7fmzJmj4uJiRSKRc94t6ziOVq5cqREjRmjgwIGaMWOGPv7448xUNkOqqqp0ww03aPDgwRo+fLjmzp2rQ4cOJXymo6ND5eXluvjii3XhhRdq3rx557y04Hy3bt06TZgwIf5Cj7KyMv3iF7+Ib6eNevb4448rEolo+fLl8XW0FfwW6M461TRk2aC9vV0TJ05UdXV1j9t/9KMf6ZlnntH69ev17rvv6oILLtDMmTPV0eFT8oEA2rVrl8rLy7V3717t2LFDp06d0m233ab29vb4Z1asWKHXXntNmzdv1q5du3T06FHdddddGax13xs5cqQef/xx1dXVaf/+/brlllt055136oMPPpBEG/Vk3759ev755zVhwoSE9bQVfNfrt4hnwOTJk53y8vL4393d3U5xcbFTVVWVwVoFhyRny5Yt8b9jsZhTVFTkPPHEE/F1LS0tTl5envPTn/40AzUMhuPHjzuSnF27djmOc6ZNBgwY4GzevDn+mQ8//NCR5OzZsydT1QyEiy66yHnxxRdpox6cOHHCGTNmjLNjxw7npptucpYtW+Y4DtcT+kZgR9Zn05DNmDEjvs4tDVm2O3LkiJqamhLaLBqNasqUKVndZq2trZKkoUOHSpLq6up06tSphHYaN26cRo0albXt1N3drZqaGrW3t6usrIw26kF5eblmz56d0CYS1xP6RmDfYPbZZ5+pu7s7/kaYswoLC/XRRx9lqFbB1tTUJEk9ttnZbdkmFotp+fLluvHGG+NvC2pqalJubq6GDBmS8NlsbKeDBw+qrKxMHR0duvDCC7VlyxZdccUVqq+vp42+oKamRgcOHNC+ffvO2cb1hL4Q2M4aSIfy8nK9//77evvttzNdlUAaO3as6uvr1draqn/5l3/RwoULtWvXrkxXK1AaGxu1bNky7dixo88zPAFnBXYa3Esasmx3tl1oszOWLFmi119/Xb/61a8S0q4WFRWpq6tLLS0tCZ/PxnbKzc1VaWmpJk2apKqqKk2cOFFPP/00bfQFdXV1On78uK677jr1799f/fv3165du/TMM8+of//+KiwspK3gu8B21l7SkGW70aNHq6ioKKHN2tra9O6772ZVmzmOoyVLlmjLli166623NHr06ITtkyZN0oABAxLa6dChQ2poaMiqdupJLBZTZ2cnbfQFt956qw4ePKj6+vr4cv311+vee++N/3/aCn4L9DS4WxqybHTy5EkdPnw4/veRI0dUX1+voUOHatSoUVq+fLm+973vacyYMRo9erQeeeQRFRcXa+7cuZmrdB8rLy/Xpk2btG3bNg0ePDh+3zAajWrgwIGKRqNavHixKioqNHToUBUUFGjp0qUqKyvT1KlTM1z7vlNZWalZs2Zp1KhROnHihDZt2qSdO3fqzTffpI2+YPDgwedkR7rgggt08cUXx9fTVvBdph9Hd/PjH//YGTVqlJObm+tMnjzZ2bt3b6arlFG/+tWvHEnnLAsXLnQc50z41iOPPOIUFhY6eXl5zq233uocOnQos5XuYz21jyRnw4YN8c98/vnnzt/8zd84F110kTNo0CDnL/7iL5xjx45lrtIZ8Fd/9VfOpZde6uTm5jpf+tKXnFtvvdX55S9/Gd9OGyX3xdAtx6Gt4D9SZAIAEHCBvWcNAADOoLMGACDg6KwBAAg4OmsAAAKOzhoAgICjswYAIODorAEACDg6awAAAo7OGgCAgKOzBgAg4OisAQAIODprAAAC7v8DBXzdeWgqGvsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Setup scene\n",
    "\n",
    "scene_file = \"../scenes/floor_wall_backscattering/floor_wall_backscattering.xml\"\n",
    "tx_position = np.array([2., 0., 2.0])\n",
    "tx_orientation = np.array([0.0, 0.0, 0.0])\n",
    "tx_pattern = 'iso'\n",
    "fc = 3.5e9\n",
    "tx_slant_angle = 0.0\n",
    "\n",
    "cm_center = np.array([0., 0, 0.1])\n",
    "cm_orientation = np.array([0., 0., 0.])\n",
    "cm_size = np.array([5., 5])\n",
    "cm_cell_size = np.array([0.1, 0.1])\n",
    "num_samples = int(1e7)\n",
    "max_depth = 3\n",
    "\n",
    "scene = mi.load_file(scene_file)\n",
    "params = mi.traverse(scene)\n",
    "\n",
    "# Define Loss and set trainable variables\n",
    "\n",
    "def func_loss(x):\n",
    "    return -1e3*dr.sum(x)\n",
    "\n",
    "variables =[\n",
    "            params['mat-itu_marble.eta_r'],\n",
    "            params['mat-itu_marble.eta_i'],\n",
    "            params['mat-itu_marble.s'],\n",
    "            params['mat-itu_marble.lambda'],\n",
    "            params['mat-itu_concrete.eta_r'],\n",
    "            params['mat-itu_concrete.eta_i'],\n",
    "            params['mat-itu_concrete.s'],\n",
    "           ]\n",
    "for v in variables:\n",
    "    dr.enable_grad(v)\n",
    "    \n",
    "# Compute gradients with RBP\n",
    "    \n",
    "# Reset grads\n",
    "for v in variables:\n",
    "    dr.set_grad(v, 0.0)\n",
    "\n",
    "rbp_tracer = PathlossMapRBPTracer(scene, fc, tx_pattern, tx_slant_angle,\n",
    "                                  cm_center, cm_orientation, cm_size, cm_cell_size,\n",
    "                                  num_samples, max_depth)\n",
    "\n",
    "cm_test,_ = rbp_tracer(tx_position, tx_orientation, func_loss, loop_record=True)\n",
    "g_test = np.array([dr.grad(v).numpy()[0] for v in variables])\n",
    "print(\"RBP Grad: \", g_test)\n",
    "\n",
    "# Compute gradients with vanilla AD\n",
    "\n",
    "# Reset grads\n",
    "for v in variables:\n",
    "    dr.set_grad(v, 0.0)\n",
    "\n",
    "ref_tracer = MapTracer(scene, fc, tx_pattern, tx_slant_angle,\n",
    "                       cm_center, cm_orientation, cm_size, cm_cell_size,\n",
    "                       num_samples, max_depth)\n",
    "cm_ref,_,_,_ = ref_tracer(tx_position, tx_orientation, loop_record=False)\n",
    "l = func_loss(cm_ref)\n",
    "dr.backward(l)\n",
    "g_ref = np.array([dr.grad(v).numpy()[0] for v in variables])\n",
    "print(\"Ref. Grad: \", g_ref)\n",
    "\n",
    "print(\"Relative Error: \", np.linalg.norm(g_ref-g_test)/np.linalg.norm(g_ref))\n",
    "\n",
    "cm_ref_db = 10.0*np.log10(cm_ref.numpy()).T\n",
    "cm_test_db = 10.0*np.log10(cm_test.numpy()).T\n",
    "ale = np.abs(cm_ref_db - cm_test_db)\n",
    "plt.figure()\n",
    "plt.imshow(ale)\n",
    "plt.colorbar()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "4cd7ab41f5fca4b9b44701077e38c5ffd31fe66a6cab21e0214b68d958d0e462"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.0rc1"
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}