isuruf/m2m.ipynb

## m2m.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "import sumpy.symbolic as sym\n",
    "from sumpy.cse import cse\n",
    "from sumpy.kernel import *\n",
    "from sumpy.expansion.multipole import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 287,
   "metadata": {},
   "outputs": [],
   "source": [
    "knl = HelmholtzKernel(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "metadata": {},
   "outputs": [],
   "source": [
    "def translate_from2(self, src_expansion, src_coeff_exprs, src_rscale,\n",
    "        dvec, tgt_rscale):\n",
    "    if not isinstance(src_expansion, type(self)):\n",
    "        raise RuntimeError(\"do not know how to translate %s to \"\n",
    "                \"Taylor multipole expansion\"\n",
    "                           % type(src_expansion).__name__)\n",
    "\n",
    "    if not self.use_rscale:\n",
    "        src_rscale = 1\n",
    "        tgt_rscale = 1\n",
    "\n",
    "    logger.info(\"building translation operator: %s(%d) -> %s(%d): start\"\n",
    "            % (type(src_expansion).__name__,\n",
    "                src_expansion.order,\n",
    "                type(self).__name__,\n",
    "                self.order))\n",
    "\n",
    "    from sumpy.tools import mi_factorial\n",
    "\n",
    "    src_mi_to_index = dict((mi, i) for i, mi in enumerate(\n",
    "        src_expansion.get_coefficient_identifiers()))\n",
    "\n",
    "    for i, mi in enumerate(src_expansion.get_coefficient_identifiers()):\n",
    "        src_coeff_exprs[i] *= mi_factorial(mi)\n",
    "\n",
    "    result = [0] * len(self.get_full_coefficient_identifiers())\n",
    "    from pytools import generate_nonnegative_integer_tuples_below as gnitb\n",
    "        \n",
    "    for d in range(self.dim):\n",
    "        result = [0] * len(self.get_full_coefficient_identifiers())\n",
    "        for i, tgt_mi in enumerate(\n",
    "                self.get_full_coefficient_identifiers()):\n",
    "            \n",
    "            new_mis = []\n",
    "            for mi_i in range(tgt_mi[d]+1):\n",
    "                new_mi = list(tgt_mi)\n",
    "                new_mi[d] = mi_i\n",
    "                new_mis.append(tuple(new_mi))\n",
    "\n",
    "            for src_mi in new_mis:\n",
    "                try:\n",
    "                    src_index = src_mi_to_index[src_mi]\n",
    "                except KeyError:\n",
    "                    # Omitted coefficients: not life-threatening\n",
    "                    continue\n",
    "\n",
    "                contrib = src_coeff_exprs[src_index]\n",
    "\n",
    "                for idim in range(self.dim):\n",
    "                    n = tgt_mi[idim]\n",
    "                    k = src_mi[idim]\n",
    "                    assert n >= k\n",
    "                    from sympy import binomial\n",
    "                    contrib *= (binomial(n, k)\n",
    "                            * sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k))\n",
    "\n",
    "                result[i] += (\n",
    "                        contrib\n",
    "                        * sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(src_mi))\n",
    "\n",
    "            if d == self.dim-1:\n",
    "                result[i] /= mi_factorial(tgt_mi)\n",
    "        \n",
    "        src_coeff_exprs = result[:]\n",
    "\n",
    "    logger.info(\"building translation operator: done\")\n",
    "    return (\n",
    "        self.expansion_terms_wrangler.get_stored_mpole_coefficients_from_full(\n",
    "            result, tgt_rscale))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 289,
   "metadata": {},
   "outputs": [],
   "source": [
    "orders = list(range(4, 10))\n",
    "ops = [[],[],[],[]]\n",
    "for i, cls in enumerate([LaplaceConformingVolumeTaylorMultipoleExpansion, VolumeTaylorMultipoleExpansion]):\n",
    "    for order in orders:\n",
    "        expansion = cls(knl, order)\n",
    "        dvec = sym.make_sym_vector(\"d\", knl.dim)\n",
    "        coeffs = sym.make_sym_vector(\"c\", len(expansion.get_coefficient_identifiers()))\n",
    "        exprs = expansion.translate_from(expansion, coeffs, 1, dvec, 1)\n",
    "        ops[i].append(count_ops(cse(exprs)))\n",
    "        coeffs = sym.make_sym_vector(\"c\", len(expansion.get_coefficient_identifiers()))\n",
    "        exprs = translate_from2(expansion, expansion, coeffs, 1, dvec, 1)\n",
    "        ops[2+i].append(count_ops(cse(exprs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 290,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([115, 194, 305, 448, 644, 866], [82, 116, 154, 183, 235, 284])"
      ]
     },
     "execution_count": 290,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ops[0], ops[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 295,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEJCAYAAACOr7BbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydeVzVdfb/n597uXDZ90VkVURFQEVBURGXxAVzTZyWmZo2bV8mv9nMlDUzNU02v7KasmXSmWoKLcm1tFwyzUQ2V1xAEUFFZF8v3Hvfvz8uUO4oXC7g+/l48NDP+7Odey/c8znvc96vowghkEgkEonkYlSWNkAikUgknRPpICQSiURyWaSDkEgkEsllkQ5CIpFIJJdFOgiJRCKRXBbpICQSiURyWawsbUBb8PDwEEFBQZY2QyKRSLoU6enp54UQntc6rks7iKCgINLS0ixthkQikXQpFEU52Zrj5BSTRCKRSC6LdBASiUQiuSzSQUgkEonksnTpHMTlaGxspKCggPr6ekubIjETWq0WPz8/NBqNpU2RSLo13c5BFBQU4OjoSFBQEIqiWNocSTsjhKCkpISCggKCg4MtbY5E0q3pdlNM9fX1uLu7S+fQTVEUBXd3dxkhSiQdQLdzEIB0Dt0c+flKJB1Dl3QQiqLcqijKBxUVFZY25bKcPXuW3/zmN/Tu3ZuwsDCmTJnC0aNHLW1WuxMUFMT58+ctbYZEclOhO1lJQ2F1h9yrSzoIIcRaIcSDzs7OljblEoQQzJw5kzFjxpCbm8uhQ4d45ZVXKCoqMvu99Xq92e8hkUgsg6GygdLkIxS/t5fK71u1zq3NdEkH0ZnZunUrGo2G+fPnt4wNGjSIUaNGsWDBAsLDw4mIiCA5ORmAbdu2ER8fT1JSEqGhoSxcuJDPPvuMmJgYIiIiyM3NBeCee+5h/vz5xMXFERoayrp16wBYvnw5c+bM4dZbbyUhIQGAxYsXEx0dTWRkJIsWLQKgpqaGxMREBg4cSHh4eMv9Fy5cSFhYGJGRkTzzzDMAFBcXM3v2bKKjo4mOjmbnzp0AlJSUkJCQwODBg5k3bx6yG6FE0nHU7iumdl8xjmP9cftNvw65Z7erYrI0Bw4cYMiQIZeMr1q1iqysLPbu3cv58+eJjo5m9OjRAOzdu5fs7Gzc3Nzo1asX999/P6mpqSxZsoS3336bN998E4C8vDx++OEHcnNzGTt2LDk5OQDs2rWLffv24ebmxqZNmzh27BipqakIIZg2bRrbt2+nuLgYX19f1q9fD0BFRQWlpaWkpKRw+PBhFEWhvLwcgCeeeIKnnnqKUaNGkZ+fz8SJE8nOzuall15i1KhRvPDCC6xfv54PPvigI95SieSmpf5IKcIgsA1zxyG2B7b93bByt+2w+3drB/HS2oMcOl3ZrtcM83Vi0a0Drvu8HTt2cPvtt6NWq/H29iY+Pp49e/bg5OREdHQ0PXr0AKB3794tkUBERARbt25tuUZSUhIqlYo+ffrQq1cvDh8+DMCECRNwc3MDYNOmTWzatInBgwcDUF1dzbFjx4iLi+OZZ57h2WefZerUqcTFxaHX69Fqtdx///0kJiYydepUAL7//nsOHTrUct/KykqqqqrYvn07q1atAiAxMRFXV9frfh8kEsm10ZfUUb7uOPXZpVgHO2Eb5o6iVnWoc4Bu7iAswYABA/jyyy8vGb/adIyNjU3L/1UqVcu2SqW6IK9wcfVO87a9vf0F93nuueeYN2/eJfdJT09nw4YNPPfccyQkJPDCCy+QmprK5s2b+eKLL3jnnXfYsmULRqORXbt2YWt76S+jrCCSSMyHscFA1dZTVG0vQFGrcJ4chMPInhazp1s7iBt50m8r48aN449//CMffvghDzzwAAB79uzB1dWV5ORk7r77bkpLS9m+fTuLFy9uiQJaw8qVK7n77rs5ceIEx48fp2/fvmRmZl5wzMSJE3n++ee58847cXBwoLCwEI1Gg16vx83NjbvuugsHBweWL19OdXU1tbW1TJkyheHDhxMSEgJAQkIC77zzDgsWLAAgKyuLQYMGMXr0aD777DP+/Oc/880331BWVtZO75pEIgHQ5ZRTtfUUdoO9cJ4chNrJ5tonmZFu7SAsgaIopKSk8OSTT/Lqq6+i1WoJCgrizTffpLq6moEDB6IoCq+99ho+Pj7X5SD69u1LfHw8RUVFLF26FK1We8kxCQkJZGdnExsbC4CDgwOffvopOTk5LFiwAJVKhUaj4b333qOqqorp06dTX1+PEII33ngDgLfeeotHHnmEyMhI9Ho9o0ePZunSpSxatIjbb7+dqKgo4uPjCQgIaJ83TSK5iWk4U0Pj2RrsB3uh7e+G1xNRWPewv/aJHYDSlStRhg4dKi7uB5GdnU3//v0tZJH5uOeee5g6dSq33XabpU3pFHTXz1ly82CsbaTiu5PU/HwGtZMNPguGolh1TGGpoijpQoih1zpORhASiUTSgQijoGbPWSo35mGs02M/vAfOEwI7zDlcD9JBdBGWL19uaRMkEkk7oD9XS/nXOVgHOeEyLaTTTCddDukgJBKJxMwYKnXUHy7DPsYHjY89Xg8PQuPn0OmrAqWDkEgkEjMh9EaqdxZSufkUwijQ9nND7WSNtb+jpU1rFdJBSCQSiRmoP1JK+drj6M/Xoe3vhsvUXqidrC1t1nUhHYREIpG0M4aaRko+y0btZIP77wdg29fN0ibdEJ0vbd4NUKvVDBo0qOUnLy/vqsf/WjbbwcGhAyyUSCTtjbHBQPXuMwghUNtr8Lg/Au8no7qscwAZQZgFW1tbsrKyLG2GRCLpAIQQ1O0rpmL9CQyVDWh87LEJdMImwMnSprUZGUF0EMuXL+fRRx9t2Z46dSrbtm2znEESiaTNNJyupviDfZR+fgSVgwbP+ZHYBHZ9x9CMjCDMQF1dHYMGDQIgODiYlJQUC1skkUjaG2EQlHxyCKEz4DIzBPtoHxRV5y5bvV66t4P4ZiGc3d++1/SJgMmvXvUQOcUkkXRPhFFQm3kOu4GeKFYq3O/sj5WbFpWdxtKmmYXu7SA6EVZWVhiNxpbt+vp6C1ojkUiuF11eBeVrcmk8XQOA/RBvrP26xnqGG6V7O4hrPOl3JEFBQbz77rsYjUYKCwtJTU21tEkSiaQVGCp1VHyTR23mOdTO1rjd3g/bSA9Lm9UhdEkHoSjKrcCtzf0LugIjR44kODiYiIgIwsPDiYqKsrRJEomkFZQmH0GXV4njWH8cx/qjslZb2qQOQ8p9S7ok8nOWmJP6I6VoejqgdrCmsagGxarj232aEyn3LZFIJNfJr3tBO471x3liEBrvzqu2am6kg5BIJDc9na0XdGdBOgiJRHLTU7H+ODW7z3aaXtCdBekgJBLJTUnD6WoUazUaD1scx/pjN9gLmyBnS5vVqZAOQiKR3FQYaxup2HSSmt1nsA33MC12c9Fi5aK1tGmdDukgJBLJTcGVekFLrowU6zMDZ8+e5Te/+Q29e/cmLCyMKVOmcPToUUub1e78WqZcIunsVO8spDwlBytve7wej8J1eki3lchoL2QE0c4IIZg5cyZ33303X3zxBQBZWVkUFRURGhpq1nvr9XqsrORHKpE0Y6jUYahuxNrXAfsYH9TONthGeHT6XtCdBRlBtDNbt25Fo9Ewf/78lrFBgwYxatQoFixYQHh4OBERESQnJwOwbds24uPjSUpKIjQ0lIULF/LZZ58RExNDREQEubm5ANxzzz3Mnz+fuLg4QkNDWbduHWCSEZ8zZw633norCQkJACxevJjo6GgiIyNZtGgRADU1NSQmJjJw4EDCw8Nb7r9w4ULCwsKIjIzkmWeeAaC4uJjZs2cTHR1NdHQ0O3fuBKCkpISEhAQGDx7MvHnzuNIiyxdffJF7772XMWPG0KtXL956662WfZ9++ikxMTEMGjSIefPmYTAYWLFiBU8//TQAS5YsoVevXgDk5uYyatSodvhUJDcbQm+kctspzr6eRtnKowghUNlYYRfpKZ3DdSAfN9uZAwcOMGTIkEvGV61aRVZWFnv37uX8+fNER0czevRoAPbu3Ut2djZubm706tWL+++/n9TUVJYsWcLbb7/Nm2++CUBeXh4//PADubm5jB07lpycHAB27drFvn37cHNzY9OmTRw7dozU1FSEEEybNo3t27dTXFyMr68v69evB6CiooLS0lJSUlI4fPgwiqJQXl4OwBNPPMFTTz3FqFGjyM/PZ+LEiWRnZ/PSSy8xatQoXnjhBdavX88HH3xwxffh8OHDbN26laqqKvr27ctDDz1ETk4OycnJ7Ny5E41Gw8MPP8xnn31GQkICixcvBuDHH3/E3d2dwsJCduzYQVxcXPt9OJKbgrrDpVSsu7AXtHQKN0a3dhD/SP0Hh0sPt+s1+7n149mYZ6/7vB07dnD77bejVqvx9vYmPj6ePXv24OTkRHR0ND169ACgd+/eLZFAREQEW7dubblGUlISKpWKPn360KtXLw4fNr22CRMm4OZmamu4adMmNm3axODBgwGorq7m2LFjxMXF8cwzz/Dss88ydepU4uLi0Ov1aLVa7r//fhITE5k6dSoA33//PYcOHWq5b2VlJVVVVWzfvp1Vq1YBkJiYiKur6xVfb2JiIjY2NtjY2ODl5UVRURGbN28mPT2d6OhowNQ3w8vLCx8fH6qrq6mqquLUqVPccccdbN++nR9//JFZs2Zd93stuXmpyy6h5D+HsPKwxeP3A9B24XafnYFu7SAswYABA/jyyy8vGb+a5pWNzS+LclQqVcu2SqVCr9e37Lv4Kah5297+FykAIQTPPfcc8+bNu+Q+6enpbNiwgeeee46EhAReeOEFUlNT2bx5M1988QXvvPMOW7ZswWg0smvXLmxtL9WeudyT2L/+9S8+/PBDADZs2HDJa1Kr1ej1eoQQ3H333fz973+/5BqxsbEsW7aMvn37EhcXx8cff8yuXbv45z//eZl3TCL5BaPOgP5cLdb+jmj7uuE6uw92g71QrOQMepsRQnTZnyFDhoiLOXTo0CVjHYnRaBQxMTHigw8+aBlLTU0VL774okhISBB6vV6cO3dOBAQEiDNnzoitW7eKxMTElmPj4+PFnj17hBDign133323mDx5sjAYDCInJ0f07NlT1NXViWXLlolHHnmk5fyNGzeKmJgYUVVVJYQQoqCgQBQVFYnCwkJRV1cnhBAiJSVFTJ8+XVRVVYmioiIhhBAlJSXC1dVVCCHE7bffLl577bWWa2ZmZgohhHjsscfEX//6VyGEEBs2bBCAKC4uvuQ9WLRokVi8eHHL9oABA8SJEyfEwYMHRUhIyAX3zMvLE0IIsWzZMuHv7y8+/PBDodfrRb9+/cTgwYOv+D5b+nOWWB6j0ShqMovE6Zd/FoV/+UkYdHpLm9RlANJEK75jZQTRziiKQkpKCk8++SSvvvoqWq2WoKAg3nzzTaqrqxk4cCCKovDaa6/h4+PTMk3UGvr27Ut8fDxFRUUsXboUrfbShT0JCQlkZ2cTGxsLgIODA59++ik5OTksWLAAlUqFRqPhvffeo6qqiunTp1NfX48QgjfeeAOAt956i0ceeYTIyEj0ej2jR49m6dKlLFq0iNtvv52oqCji4+MJCAi4rvcmLCyMv/3tbyQkJGA0GtFoNPzrX/8iMDCQuLg4Tp06xejRo1Gr1fj7+9OvX7/rur7k5qHhdDXla3JpyKtE09MBl2m9byoZ7o5Cyn13Ee655x6mTp3KbbfdZmlTOgXd9XOWXJvGc7UUvZGOytYKp4lB3bIXtLmRct8SiaTbIIyCxsJqrP0d0XjZ4TIjBLsIj5t3oZsQ0AGVWdJBdBGWL19uaRMkEougy6ugfHUujedq8XlmKFauWhyG9bC0WZbB0AgHVsGut2HOf8C9t1lvJx2ERCLplBgqdVRsOEFtVrGpF3RSX9QuN6kMd0MNZHwCu96BilPg2Q9qzksHIZFIbj6M9XrO/r8MhN6A4zh/HMfcXL2gW6gpgdQPTD91peA/HKYshj4TQWX+Ml7pICQSSaehobAa654OqLRWuEwNxibYuVv1gm41ZSdh178g47+gr4PQyTDqSQgY3qFmSAchkUgsjv58Uy/ow6V4zo/EJsgZ+6E+ljar4zl7AHa+acozKApEJMHIx8HLMhV70kGYAbVaTURERMv2119/TVBQ0BWPDwoKIi0tDQ8PDxwcHKiuru4AKyUSy2PUNfWC/rGpF/SUYKz9HC1tVsciBOTtMDmGnO9BYw/D5kPsw+Dsd8GhZWdPY2VtjaObR4eYdk0HoSiKPVAnhDAqihIK9AO+EUI0mt26LoqtrS1ZWVmWNkMi6dQIo+Dcu1noi2qbekEHo3aytrRZHYfRAIfXmxxDYTrYecC4P8PQ+8DuFw0pIQQFh/aTvmE1uempDEpIZPy9869y4fajNVmO7YBWUZSewGbg98BycxrVHVm+fDmPPvpoy/bUqVPZtm1bq8+/5557ePzxxxkxYgS9evW6QO/pcvLer732WovM9lNPPcW4ceMA2Lx5M3fddVc7vCKJ5MZoPFeLMAoUlYLT+AA8HxqI29y+N49z0Osg/T/wrxhY8VuoLYHEf8JTB2D0ggucw+GdP/DpwidZ8Zc/Ungkm+Ezkxg2M6nDTG3NFJMihKhVFOU+4G0hxGuKomSa27CuTF1dHYMGDQIgODiYlJSUdrnumTNn2LFjB4cPH2batGncdtttV5T3Hj16NP/85z95/PHHSUtLQ6fT0djYKCW0JRbDUNNI5XemXtCut4ViP8Qbu0hPS5vVcdRXQNrH8PN7UF0EPpFw28fQfzqof/kqrq+uRuvgAMDJ/Xsx6BuZ8OCj9I8bi8a6Y8t8W+UgFEWJBe4E7ruO8yzO2VdeQZfdvnLfNv374fPHP171GHNNMc2YMQOVSkVYWBhFRUXAleW9f/e735Genk5VVRU2NjZERUWRlpbGjz/+eEEDH4nE3AijoCb1DJWbTmKs1+MQ64tt/5tIhrvqLPz8LqQtA10lBMfDzKXQa+wFq6FLCvJJ37Ca7O1bSVr0d3r06cvYex5AY6O1WD+L1nzRPwE8B6QIIQ4qitIL2HqNcyQXYWVlhdFobNmur6+/6vF/+tOfWpr7NDubX0toN2toiavIewcFBbFs2TJGjBhBZGQkW7duJTc3V2oYSTqU0s+yqTtYgk0vZ1ym9UbjY3/tk7oD53PgpyWw9wsw6qH/NFOpqu/glkOEEJzcm0H6htXk7c3ASmNN2Ohx2Dk7A2CttWyJ7zUdhBBiO6Y8RPP2ceBxcxrVXlzrSb8jCQoK4t1338VoNFJYWEhqaupVj3/55Zd5+eWXr3ndiRMn8vzzz3PnnXfi4OBAYWEhGo0GLy8vRo8ezeuvv87HH39MREQETz/9NEOGDJHdtSRmx1CpQ9FaobJWYz+8B7YDPW+eXtAF6bDzDcheB2prGHwXxD56wapnIQSKoqDX6Vj31mtYWdswcu5vibxlEnZOzhY0/kJaU8UUCjwDBP36eCHEOPOZ1f0YOXIkwcHBREREEB4eTlRUVLtc90ry3l5eXsTFxfHyyy8TGxuLvb09Wq1W5h8kZkXojVTtKKRqSz4OcX44TwhE2+fKnQe7DUJAzmZTRVLej6B1hrinTeWqDl4th1WXlZK1cT0F2QeYu+jvaLRa5jz/Ch7+AaitOp/w4DXlvhVF2QssBdIBQ/O4ECLdvKZdm5tJ7ltyIfJz7nxc0As6zB2XxODuvwraoIeDKbBzCRTtB0df0/qFIfeAzS/rOYpO5JKx/msO//QjRqOB3kOGMemhJ1uS0R1Ne8p964UQ77WDTRKJpJtSsSmPqi2nsPK8SXpBN9RC5ifw0ztQkQ8eoTD9X6aVz1YXluvm7cvkq5efR6O1ZWDCZAZPuhVXH18LGX59tMZBrFUU5WEgBdA1DwohSs1mlUQi6fQYdQaE3ojaXoNtuAcqrRUOI3y7dy/o2lKTcN7u95vE84bB5H9A6KQW8byG+joObP0eK2sNkeMn4R8Wwbjfz6N/3Fi09u0TMRjKy1E5O5s9p9MaB3F3078LfjUmgF7tb45EIunsCCGo21tMxYYT2PR2wW1uX6x9HbD2tcx0SYdQnv+LeF5jrckhjHwSAmNbDqksPkfGt2s5sGUTutoa+sSMIHL8JNRWVgyedGubTRBCUJeWRlnyCqo2biRg2cfYDb3mLFGbaE0VU7BZLWhCUZT+mEpqPYDNclpLIul8XNwL2n54N2/cU3TQlF/Y/2WTeN4cGPE4eIddcFjq6i/Z8cV/AQgdNpIhiTPo0advu5hgKC+nYvVqylaspCE3F5WDAy5z5mDl7d0u178arali0gAPAaObhrYB77dGi0lRlI+BqcA5IUT4r8YnAUsANfCREOJVIUQ2MF9RFBXw4fW+EIlEYl5qMoooW3kUlZ0VrrP6YDfUu3v2ghYCTv5kqkg6tqlJPG8eDH8YXPwBMBoMHN29E98+/XDy9KJHSChDEmcweNJUnDy8rnGD1pggqMvMpDw5mcpvNyJ0OrSRkfR4+W84TZ6Mys6uzfdoDa2ZYnoP0ADvNm3/tmns/lacuxx4B/hv84CiKGrgX8AEoADYoyjKGiHEIUVRpgELm86RSCQWRhgFxtpG1A7WaENccRjZE6dx/t2zF7TRCEc2mBxDwR6TeN7YP0P0L+J59TXV7N+8kcxv11FVUsyIOXcSe9vt+A+IxH9AZJtNMFRWUrFmLeXJyeiOHUNlb4/zzBm4zp2L1gJVe61xENFCiIG/2t7SVPp6TYQQ2xVFCbpoOAbIaVpwh6IoXwDTgUNCiDXAGkVR1gP/a809ugrNC2NefPFFXnzxxZaV0BePKYrCnXfeSVpaGhqNhpiYGN5//300mm74Bynp1OhOVFC+JhfFRo3nvEjUTta4TO2GqUe9DvYlw863oOQYuATClNdh0J1gbXpSF0Lwwyf/Zt/339Koqzclnu+dT6+otucAhBDU79tHWfIKKjdsQNTXox0wAJ+/vIRzYiIqe8utPG+NgzAoitJbCJEL0CS1YbjGOVejJ3DqV9sFwDBFUcYAswAbYMOVTlYU5UHgQYCAgIA2mNGxZGVlsWzZMsDUHyI1NZU5c+ZcMvbKK69w55138umnnwJwxx138NFHH/HQQw9ZzHbJzYWhQkf5NyeoyypG7WyD81h/S5tkHnRVkL7clHyuOgM+ETD73xA2A9RWCCE4dzwH714hKIpCQ10tocNHEjVlOl5BbXeUhupqKteupSx5BbrDh1Hs7HC+9VZc5s7FNnxA219fO9AaB7EA2KooynFAAQIxSX7fKJebtBRCiG2Y8htXRQjxAfABmBbKtcEOs5GXl8ekSZMYNmwYmZmZhIaG8t///peHH36Y2NhYGhsbee89Uw7e1tb2krEpU6a0XCsmJoaCggKLvA7JzYcur4LzHx9AGEX37QVdXQy7l8KeD00Kq0FxMP0d6D0eFAV9YyNHdmwmfcNqivOOc9erS/AO7s2EBx9rl7LSuv0HKF+RTMX6DYjaWmz698fnxUU4TZ2K2kIL565Ea6qYNiuK0gfoi+nL/bAQQneN065GAfDrRxI/4HQbrtcpOXLkCP/+978ZOXIk9957L++++y75+fncddddjB8/nj//+c/cdtttfPzxxxeM/e1vf2u5RmNjI5988glLliyx4CuR3AwYahpR22vQ+DpgN8gLx3i/7rcKuizPtLAt8xPTtFL/qTDyKfAbAoCutpaMb1azd9MGasrLcPcLYMKDj+HW09TVrS3OwVBdQ+X69ZQnJ1N/6BCKrS1OUyabcgsREZ1Wo+qKDkJRlHFCiC2Kosy6aFdvRVEQQqy6wXvuAfooihIMFAK/Ae64wWtdk5R/ZlwyFjLEi4gxfjQ2GFj39qXplH6xPeg/ogd11Q18+/6BC/bN/EPrNJT8/f0ZOXIkAHfddRdLlixh9erVvPjii8yYMYPp06cD8NZbb10y1szDDz/M6NGjpX6SxGw0nq+jYt1xGs/V4vPUEFTWalxn9bG0We3L2QOmUtUDX4GigoFzYcQT4BkKQGN9PRqtFoC0tSn07NufqCnTCYwc3OYv7vpDh0y5hbVrMdbWYhMaivfzf8Z52jTUjp2/terVIoh4YAtwuRUeArimg1AU5XNgDOChKEoBsEgI8W9FUR4FNmIqc/1YCHHweg3v7Fz8i6VqWmX54osvXrL/cmMvvfQSxcXFvP/+++Y1VHJTYuoFnU/Vj4UoViqcxge0rr9kV0EIyN8FO974pVR1+EOmUlXnngghyMtKJ33919SWl/Hb197Gxs6O+9/+CFtHpzbd2lhbS+WGDZQlr6B+/34UGxucJk/GZW4StoMGddpo4XJc0UEIIRY1/fcvQogTv97X9PR/TYQQt19hfANXSUS3J1d74tdYq6+639bButURw8Xk5+eza9cuYmNj+fzzzxk1alSrz/3oo4/YuHEjmzdvbnEsEkl7oS+t59zSvRgrG7CL8sJ5UjfqBW00wrGNJsdwajfYuV9Qqtqoq+fQd9+Q8c0aSgtPYe/qxuCJUzEa9KitNG1yDvVHjlCenEzFmrUYq6uxDumN9x//iPP0aaidO4+E9/XQmiT1V8DF35JfAkPa35zWoSjKrcCtISEhljLhmvTv35///Oc/zJs3jz59+lxXFdL8+fMJDAxskfCeNWsWL7zwgrlMldwkGOv0qGytULvYYNvPDbsh3tgEtu1pudNgaDStdt75JhQfBucAmLzY1IvB+pdFZTmpu/j+o3/hFdybyY/+gb6xo9oks22sq6Pym28pT06mbu9eFGtrHCdNxHXuXGyjoswWLeScqybI3Q4rtXkfIK+Wg+gHDACcL8pDOAFas1p1DYQQa4G1Q4cOfcCSdlwNlUrF0qVLb+hcvV7fztZIbmaae0HX7SvG++khqB2su0+eoaEGMj6BXe9AxSnwCoNZH8KAmaDWUHQ8h/QNq/EODmFI4nRCY0fh6O5Jz/4D2vTlrcvJoSx5BRWrV2OsrMQ6OBivhc/iPH06Vq7m6X9RVtPA2n2n+SqjkL2nyln++2jG9G37qu2rcbUIoi8mmQwXLsxDVAGd9otZIpGYuKQX9HBfFDM/cXYYtaWQ+qGpXLWuFNAKFq8AACAASURBVAJGQOL/gz4TMAojuempZKxfTUH2ATRaW9x7mgon1VYa/MLCr3Hxy2PU6ajauJGy5BXUpaejaDQ4JiTgMjcJu+hos0QLDXoj246cY1VGIZsPF9FoEPTzceTPif0J72n+aaur5SBWA6sVRYkVQuwyuyXdiKCgIA4cOHDtAyUSM2FsMFD83l4az9R0r17QFQWmhW3py5tUVSeb+jwHDG85ZNN7b3Hwh804engSf9e9hI9LaJPMtu74ccqTV1Dx9dcYKiqwDgzEa8ECnGfOwMqt/fteCCE4UFjJVxkFrNl7mtKaBjwcrPldbBCzo/wI8+24acHW5CAyFUV5BNN0U8vUkhDiXrNZJZFIbgijzoDKRo3KWo1Nbxccx/p3j17QxUdMpar7kk0VShFzYOQT4B1GZfE5Mj/9mMGTbsXJw5OI8ZMIHhxNn5hYVOobW+RnbGigatN3lCcnU7tnD2g0ON4yHte5c7GLiUExQ/FIUWU9KZmFrMoo4GhRNdZqFRPCvJk9pCdxfTzRWCD6a42D+AQ4DEwE/gLcCWSb0yiJRHJ9tPSC3nYKr/kD0fjYdw/dpFN7TBVJR9aDlS1E3w+xj4BLAKePZpP+v1c5lvoTAF6BwTjFjaVn3xsXtWvIy6NsxUoqUlIwlJWh8ffH8+mncZk1EysPj/Z6VS3UNRjYdOgsX2UUsuNYMUYBUQEuvDwznKkRvjhbWBSxNQ4iRAgxR1GU6UKI/yiK8j9MaxgsRleoYpJIOoq6w6VUrM1FX1KPNswdxaaLS2MIATnfw4434eQO0LpA/LMQMw/s3TEaDCS/8H+cPnIIGzv7Nstsi4YGqjZvpix5BbU//wxqNY7jxuEydy72I2LbPVowGgV78kpZlVHI+v1nqNbp6eliyyNjQ5gV5UewR+eZCmyNg2ju+1CuKEo4cBYIMptFraArVDFJJOZGCEHJp9nUHywx9YK+NxxtqHkqaDoEgx4OfW1yDEX7waknTPw7RP2O+kbIy0qn38h4VGo1Pfv2p9+IOAaMuQVr7Y1JgjScPEn5l19SvioFQ0kJGl9fPJ98AudZs9B4tX910MmSGlZlFLIqs4BTpXXYW6uZHNGD2VF+DAt2Q9UJe2u0xkF8oCiKK/BnYA3gADxvVqu6Idcj933fffeRlpaGEILQ0FCWL1+OQycT8ZJYDtFoQNGoURQF654O2AQ6de1e0I11kPWZSW67/CR4hML0dyFiDmXFxWR89gkHt22mUVePT0hfXLx9GH3njemFioYGqrZsoXzFCmp+2gVqNQ5jx+CalIT9yJEoN5izuBKV9Y2s33eGVRkF7MkrQ1FgZG8Pnp4QysQBPthZt+Yr2HJc1bqm7m6VQogyYDuyD/UNcz1y32+88QZOTqZKhaeffpp33nmHhQsXWsx2SedACEFdVjHl35zAdXYfbPu64TSu60jeX0JdOez5yFSqWlMMPYfCxFeg7xQqzhez9Y1/kJueilqtpt/IeKKmTMfF2+eGbnVxtGDl2wPPJx7HedZsNN7tGy3oDUZ+zDnPqoxCNh08i05vpLenPf83qS8zBvXE16XriCBe1UEIIYxNukkrOsiebkFb5b6bnYMQgrq6uq5fgSJpMxf3glbbd+EGUpVn4Od3IW0ZNFRByC0w6in0vjFUl5XiolJhY2dP8ck8hs+ay6CEROxdrn/qrKOjhcNnK1mVUUhKZiHFVTpc7DTMjfZndpQfkX7OXfLvWGme6rjiAYryPFAHJAM1zeNCiFLzmnZthg4dKtLS0i4Yy87Opr8FWvP9mry8PIKDg9mxY0eL3HdYWBj5+fkIIRg/fjxpaWktct+/HmuW+/7973/Phg0bCAsLY/369dh1UA/arkJn+Jw7ioqNeVRtO4XKzgrnicFdtxf0+Rz4aQns/QKMetNq55FPUmsfyN7vNrB30wbsXd246+9voigKwmi8oQTx5aIF1zlzzBItnK/WsSbrNF9lFHDwdCVWKoWx/byYHdWTsf28sLHqnAUDiqKkCyGu2Q6vNQ7ixGWGhRDC4tNNrXEQyS9dOjXTd3gcgyYm0qirZ9WrL16yf0D8LYSPuYXaygrWvvH3C/bNXfTqNe3Ky8tj9OjR5OfnA7Bly5YL5L6vlYNoxmAw8NhjjxEdHc3vf9+WHk3dj+7uIISx6fdDpVCz5yyNZ2pwuiWga/aCLswwaSQdWgNqa5M+0ohHKanTkL4uhewft6FvbCB40JAbltnuyGhBpzewOfscqzIK2HakGL1RENHTmdlRPbl1oC/uDjbtdi9z0VoH0ZqGQa1Sbu1IukKZa1vlvgHUajVz585l8eLF0kHcROiOV1C+Nhf7YT44DPfFPvrG5t0tihBwfJvJMRzfBjZOMOopxLD5GG3dUVtZcXrrJrJ3/EBY/DiiJk/H3e/6W5t2VG5BCEHmqXJWZRSwdu8ZKuoa8Xay4b64YGZH+RHq3fl7O9wIrYkg7ICngQAhxIPN3eWEEOs6wsCr0dmnmH766SdiY2N54IEH6NevH3/4wx+uea4QgtzcXEJCQhBCsGDBAgBef/11c5vdpegMn3N7o6/QUbHhBHV7i1G72OAyrTe2Ye6WNuv6MBoge61pcduZLHDwhuEP0xhxB4dS08n4Zg2R4ycxJHE6+sZGGuvrrltiuyOjhdPldaRkFvJVRgHHi2vQalRMHODDrCg/RoV4oO6KU320YwQBLAPSgRFN2wXASsDiDqIzc6Ny30II7r77biorKxFCMHDgwJbktaT7UpN2lvI1uaZe0OMDcIz361q9oPU6U25h5xIozQW3XnDrEqr9E8ja8j17lz9JfVUlXsG9cfExRURWGg1WmtZPmTXk51O+cqXZo4UanZ5vD5zlq4wCdh0vQQiICXZj3uheTInogaO2C07z3SCtcRC9hRBzFUW5HUAIUad0xXR8B3Ojct8qlYqdO3eawSJJZ0QYBIpaQe1sg02IKy5Te2HlZlE1/eujvhLSl8Gud6H6LPQYCHOWQ/9poFLzzV//RP7BfYQMHcaQKTOuW2a7o6IFo1Hw8/ESvswo4NsDZ6ltMBDgZscT4/swa7AfAe43Z5FIaxxEg6IotpjajKIoSm9AZ1arJJJuTmNxLRXrjmPlbY/LlGC0fVzR9ulCq6ArCkxy22nLQFcBwfEYp79LbpktWavWM/nReBxc3Yj/7X1Ya21x8elxXZe/crQwC423d7u9jNPldXyx5xRfpp3idEU9jjZWTBvoy+whfgwNdO2SpantSWscxCLgW8BfUZTPgJHAPeY0qqsj5b4lV8Ko01O15RRVO0y9oG26mjTGqT2mNQyHVgMC+k2lYehDHDh6noy3P6ei6CxOnl5UFJ3FwdUNr6DWFzt2VLSgNxjZeqSYz1Pz2XbkHAKI6+PJwin9SQjzRqvpQlN7ZqY1VUzfKYqSAQwHFOAJIcR5s1smkXQz6nPKKV1xxNQLeog3zpOCUDt2gV7QhkaTQ/j5PShMAxtnGP4QxDxIg9aTDx+9j/rqKnxD+zP6jnsIib4+me2OihYKy+tI3nOKFXtOcbayHk9HGx4eE8LcaH/83W7OKaRr0VohkHhgFKZpJg2QYjaLJJJuhjAKFJWC2skaKxcbnO/qj01AF+gFXVtqasyT+iFUnQa3XohJr3HaYSiFuSeIcQ3EGhg+6zf49u1Hj5C+rb50R0cL/9t9km1HiwEY3ceTl6YPYFw/L4v0WOhKXNNBKIryLhACfN40NE9RlFuEEI+Y1TKJpItjqGmkclMexjo97nf0R+Nlh9fDgyxt1rUpPmLSR8r6HPR1EByPYfLrHCu1J33tGs7mrEXr4EjkLZPQ2jswJHF6qy/dodFCaj7JaacoqtTh5WjDo2NDSBoqo4XroTURRDwQLpoWTCiK8h9gv1mtugZdYaGc5Oblcr2gm6OITosQkLvZNI2U8z2obSAyCYY/xOkKFWvf/AfVJedx7eHL+HsfYkD8eDTa1lVbdWS0sOXwOVNuoSlaiA/15C/TAxjfzwsrGS1cN61xEEeAAOBk07Y/sM9sFrWCrtgP4nrkvpt57LHHWLZsGdXV1ZYyW3KdNJ6rpfTzw12nF3RDLez7An5eCuePmBa2jf0zpX6TacAGH+8+uNhW4OEfyC33PUyvwUNbrY9kiha+pHzVKrNGCwVltazYc+qSaGFutD9+rt0jWqivaST/UAmlp2sYPr13h923NQ7CHchWFCW1aTsa2KUoyhoAIcQ0cxnXnbgeuW+AtLQ0ysvLLWav5Ppodu6qJpVVtzv6de5e0BWFsKepTLW+HHoMRMx4n1NKH9K/XcfxzOfw6zeAuS++ip2TM7Ofe6lVlzVFC1spX5H8S7QwZgyuc80TLfwvNZ8ffhUt/HV6AOO6SbRQUVxHbsY58vaf5+zxSoRRYOuoIWpiINbajukj0Zq7vGB2K7oZbZX7NhgMLFiwgP/973+kpMh6gM6M0Bup+rGQ+qNleD4Qgdpeg9fj1y8212EUpJnKVA9+TXOZKsMfJuecwk+ffUpx/qfYOjkzfNZvGJQwpdWXvVy04PH4Y7jMnt3u0ULynlOsaIoWvJ1seGxsCEndIFrQNxgoOFKGd5ATto7WnMouZVdKLh7+DkRNDCAowgOvIKcO7Tx3TS0mAEVRvDFFDgCpQohzZrWqlXR2LaYblftesmQJRqORp556CgcHBznFdBk6w+dcl11C+brjGErq0Q5wx21OKKoOerK7LgyNkL3GlF8o2GMSzov6HbX9bkfjE4LGRkvWpg1kbVzHkMQZ9B81Bivra5ffdmS0sLkpt9AcLYwJ9eT2mK4fLVSV1nPyQAkn95+n4HAZ+kYjY+7sy4C4nujq9DTWG3BwbX912PaU+04CFgPbMK2DiAMWCCG+bAc720RrHMS59y9Nl9hFeuAQ64uxwcD5ZQcv2W8/xBv7od4Yahop+TT7gn1e8yKvaVdb5L7PnDlDUlIS27Ztw8rKSjqIK2BJB2GobqBs5VHqj5Rh5WmLy7TenXMVdG0pZPzHVKZaWWjSRxr2EOc9RpH+3Xdk79jGmN/ez6CJiRgNBhSVqlWRz+UqkVxuu63do4VTpbWsSDtF8p5TnKsyRQtzh/p36WjBaBToahqxdbSmtrKBZf+3AwAnDy2B4R4ERrjTM9QFKzMv1mtPsb4/AdHNUYOiKJ7A94DFHURn5kblvjMzM8nJyaG5Qqu2tpaQkBBycnLMb7TkqrTkGbRWGKobcU4MNvWC7mxPsMVHYfd7vypTHY2Y8jp5tZ6kf7OGk/v+DytrGwbEjycgYiDANRe2XbESac4c7EeNardoobE5t7A7n+3HfokWXh4WyNi+nl0yWmhOMJ/cX0L+wVJ8ejuT+HAkdk7WjP1tP3yCnXHtYdcppyVb4yBUF00plQBd5lO62hO/ylp91f1qe02rIobLkZ+fz65du4iNjeXzzz9n1KhRrTovMTGRs2fPtmw7ODhI52BhmntBV/90Go8HIky/N48M6lxlq5ctU52DMXoeKt9IFGDX889QVXyOUb/5HZG3TGqVzHbDyZO/RAulpWarRDpV+ktu4VyVDh8nLY+N68PcaH96dqEezhez7bPDHNpxGiFA66AhMNydXoM9W/aHjfS1oHXXpjUO4ltFUTbyy0K5ucAG85nUPbhRuW9J56KhsKkX9MlKNH4OGKsbUbmpO49zaKiFfckmx9BSpvonqnrNIGvHLg799Z/89h9vYefkzNQn/g97F1fUVleXqxYNDVRt3kzZihXU7voZ1Gocx43FJSkJ+xEj2jVa2Jxtyi00Rwtj+3pxe0xAl4sW9A0GCo+Wk7f/PIVHypjzx2g01mo8AxyJmhRokQRze9AaLaYFiqLMwiS1oQAfCCFkac01uFG574uR+QfLIPRGytfmUpN6FpWdBtfZfbAb0ol6QVcUwp6PTFLbdWXgEwkz36fILpL0bzdw5P0FGI1GQoYOp7G+DpyccfK4es+Ehrw8ylaupCLlawylpWh8ffF88glTtODVfv0WulO0cO5kJXvWnWhJMFtZq/Dr50Z9dSMaNzUD4npa2sQ20aqSCyHEKmCVmW1pNXIltcTsqBX0JfU4jPDF6ZZAVLadpDqpIL1JTfVrEEbolwjDH4aAWCpLivn00fvQ2GgZlJDI4MnTcPG+ertSY0MD1d9/T1nyCmp3726KFsaZooWRI1q9KO5aNEcL/0vN58djxSjAmL5e3BETwJguEi0YjYKi4xXkHSghcIAbvk2FCaVnaug/ypegcHd8OyDB3JG0qsy1s9JZy1wl5sccn7PueAUVG/Nwu6MfVs42nUcew6D/VZlqakuZakPkbzmQeYTyc2cZd888AI7t/omAiIHY2F199bbuxAnKV35JRUoKhrIyNH5+uMyZg/PMGe0eLXyxJ58VaQUUN0ULc6NNlUhdIVowGozkZJxrSTDX1zSiqBRiZ/Zm8ISAC6oRuxLtWcUkkXRrLu4FbSjXYeVsY3nnUFcG6c1lqgXgGgyTX6Oi5wQyt2xh/3Mv0FBXi1//cAz6RtRWGvoMG3HFyxkbGqj67jvKk1dQm5oKVla/RAsjYts5Wijif6mnWqKF5txCZ48WhBCUnq6hulxH4AB3FJXCzi9zMBoEgRHuBIa7ExDmho2dKY/T1RzD9XJFB6EoymYhxHhFUf4hhHi2I41qKxdrGkm6F+0V9QohqPqhgKrNpgWMnaYXdPFRk5rq3s+hsRaCR0Pi69AngcO7drDhmSdRVAqhw0cxZMp0fEJCr3o53fETlK9caYoWysvR+Pvj+fTTuMycgZWn51XPvR4uFy083pRb8O3E0ULzCuaTB0rI23+e6lId9i423P33ESiKwuwFQ3Bw03a5BHN7cLUIooeiKPHANEVRvsCUoG5BCJFhVstuEK1WS0lJCe7u7tJJdEOEEJSUlKBtpZLo1VAUBX1xHdpQV5wTLdwLWgjI3dJUpvpdS5mqYeiDHMsrx17vir9Kjf+ASIZOm8XgiVNxdPe44uWMOh1V331PeXIytXv2mKKF8eNxSZqDfWz7Rwuf7c5nR875lmjhjmEBxId23mihqrQeBxdTlLjjyxwObi/EylqFf383hk4OIjD8Fx0tJ4/O69zMzRVzEIqi3Abch6l6Ke2i3UIIMc7Mtl2Ty+UgGhsbKSgooL6+3kJWScyNVqvFz88Pjebq5ZqXo7G4lor1J3C6JQBrP0eEwWjZhW7NZaq7l0LxYbD3guj7qe8/l30/p5G5cR3VJefpP2oMUx575pqX0x0/TvmKlVR8/XVLtOCSNAeXmTOx8riyQ7lejpytYlVmAasyCimu0tHDuSm3MLRzRgu/TjCf3H+eksIakv4YjWeAI+cLqqmt0HW7BPPVaHMOoklK40tFUZ4XQvy1Xa0zIxqNhuDgYEubIelkGHV6KrecorqpF7S+tB5rP0fLOYfK06bcwq/LVGcshfBZ7F67mp8X/AG9TkdA+MAWme0rYdTpqNr0nSlaSEszRQu33IJr0hzshg9vt2ihuErHmr2nWZVRwMHTlVipFMb09WzKLXih7mRTMM1TzcX5VaxekomuRo+iUujR25kRs0KwczbpTXn4OYCfg4Wt7Zy0Zh3EXxVFmQaMbhraJoRYZ16zJJL2o3ZfMeVrj2Os6gS9oAvSTTIYB1PAaIB+iYhhD3GqxpEeof3QWNmgdXCkb2wcQ6ZMxzPwyg87utzcX6KFigo0AQF4PfMHnGfMaLdoob7RwPfZRazKKOSHo8UYjIKIns4sujWMWwf64uHQ/kJyN0pzgjlv/3lOHighIMyNoVOCcfGxIzjCg4CLEsySa9OalqN/B2KAz5qGnlAUZaQQ4jmzWiaRtBP64jrUzta4/9ZCvaANeji81pRfOLUbrB0hZh76wfdw+NBJMt79guL8PBLmP07E2AQGTpjMwAmTL3spo05H1caNlK1YQV1aOmg0ON4yHtekJOyGDWuXaEEIQdrJMlZlFLBu3xmq6vX4OGl5IK4Xs6J6Eurt2OZ7tDc7vzxGTsY5qkt1AHj4O2Db9BCgsVYz/p4wS5rXZWmNmus+YJAQwti0rQYyhRA3JlLUjlwuByGRNPeC1oa6YjvAA2EwgqJ0fNlqXRlk/Bd2f9BUphoEwx7CEJ7E7m++Ze93G6itKMfDP5CoxOn0H3llmW1dTo5JE+nr1RgrKtAEBuCalGSKFtzd28XcvPM1rMosJCWzgFOlddhZq5kU7sPsKD+G93LvNFNIzRLZZWdriEsyVXBt/OgAhkYjgeHuBIZ7mEUiuzvR3usgXIDSpv8737BVEokZEQZTL+iKTScROj1qVy220PF5hvPHTEnnrP+ZylSD4mDKYmq8hmHv5o5KCE5k7sG7VwhDpswgIGLgZSvujPX1TdHCSurSTdGC04RbcElKwi4mpl2ihYraRtbtP82qjELST5ahKDCytwdP3RLKxAE+2Nt0jqVSJaerOZpaxMn9JZQUmuRnnDxtGT7DgMZazcT7wy1sYfekNZ/+34FMRVG2Yip1HQ3I6SVJp0J3spLylBwaz1qoF7TRaFJT3b20SU3VGiLmIGLmcaKogfQvV3P66L958N3l2Do4kvTiq2isL/+Uqzt2zKSJtHoNxooKrAMD8VqwAOeZM7Byc2uzqY0GIz8cKWZVZgHfHzpHg8FIHy8Hnp3UjxmDfenhbPkqpGaJbL++btg5WXMmp4KsTfn0CDElmIMi3XHx7pwS2d2J1iSpP1cUZRumjnIK8KwQ4uzVzzIvUotJcjH60nqM9Xrc7uyHbXgH9oLWVZsWtO1+H0qOtaipNobfzqGM/aQvfo+y0wU4uLoxfOZcVCpTGeXFzsFYX0/lt99SvmIldRkZKBoNjhMm4DJ3LnYx0W1+PUII9hdWsCqjkDV7T1Na04C7vTV3Dg9g1mA/wns6WfTLtjnB3LxY7WxuBULA2N/2I2ykL6Ex3vQZ6iUTzB1Mt9NiktwciEYjVTsKUGmtcIj1Na2u1htROqqOvfSESU014xPQVYBvFAx/CNF/OorGhpKCfJb/4eGmaaTphMaOuqzMdv3RoyZNpNWrMVZWYh0UhEtSEs4zprdLtHC6vI6vswpZlVFIzrlqrNUqJoR5MyuqJ6NDPdFYcA2IvsFAfU0jDq5aasp1LF+4EzAlmIMiTN3VvAOdLC950g1pt5ajnRnpIG4+hBDUZ5dSvt7UC9puiDduc64uNdGON4e8H+HnpXBkA6jUEDYdhs2nqMGV9A2rQVGY8ugfACjOz8PDP/CSJ3NjXR2V326kfMUK6jIzTdHCxIm4JM3BLrrt0UKNTs+3B86yKrOAn3JLEAKGBroyK8qPxIgeOFvwKbw5wZy3/zyFh8vw6+9G4sOmepec9HP49HKWCeYOQIr1Sbod+vN1lK/NNfWC9rLF477wjukF3VgH+1eappGKDoCtG8Q9jXHIveQeySf9o68oPHwQa1tbIm+Z3LJAyzMg6ILL1B89alq3sGaNKVoIDsbr2WdN0YJr216HwSjYlVvCqowCvjlwlrpGAwFudjw+rg+zonoS6N6B+Zhf8WtdtE3/PsixPUWAqQdz/1G+9Br4y3qNkCHtpyIraR+u6iAURVEB+4QQskRAYnEM1Q3oTlbinNgLhxE9zF+d1NKUZznUlYJ3OEx7ByJuA40te1JWsOOL/+Lk6c2Y391P+NgEbOzsLriEsa6Oym++NUULWVko1tY4TpyIa9IcbIcObXO0cLSoiq8yCvg6s5CiSh2OWitmDO7J7KieDAl0tUheob6mkfyDJeTtL+FMTjl3vjQcK2s1AQPc8AxwJChCJpi7Cld1EEIIo6IoexVFCRBC5HeUURIJmJ4+a7OKMZTU4XRLIDZBzvRYGINKa8bAVwg4lWpa7XxoDSCg7xQYNp8Kuz5kblxLkDqboIFRDIgfj6tvT0Kih7ckn5upP3KU8uRkKtauxVhVhXWvXngtfBbn6W2PFs5X61iTdZpVmQUcKKxErVIYE+rJC1P9GN/fC62F9ITO5JSz6+vclgSzraOGgAHu6Or0WFmr6Te8h0Xsktw4rflL6wEcVBQlFahpHhRCTDObVZKbnl/3grYOcMRxrD+KWmU+56BvMMlf7H4PTmeCjbMp6Rx9P4XnaslIWU1O6mIUlYKtkwtBA6NwcHMndNjIlksYa2qo/OYbylaupH7vPlO0MGkirklJ2A4Z0qYn5vpGA5uzz7Eqo4BtTZIX4T2deGFqGNMGdbzkRYtE9v4Sekd54tfPDbVGRaPOwJDJQQRGuOMV2PV6MEsupDV/bS+Z3QqJpAlDTSOVG/Oo2dNBvaCrz0Hax6af6iJw7wNTXoeBt4ONA2v/3ysc2/0TWnuHK8ps1x04SPnKlVSuW4expgbrkN54P7cQp2nT2hQtdDbJC4PeSPbO0+QdKKGwuQezjRr3nvb49XPDK9CJuX+K6VCbJOalNesgflAUJRDoI4T4XlEUO+Dm0MSVdDiiXk9tVrH5e0GfzjItajvwFRgaIGQCDJ9PnVc0B7dvYZDKBisgZOhwAgYMZED8eDS/6kFhqKqict06ylauRHcoG0WrxWnSJFySkrAdPKhN0cLJkhpWZRSSkllIfmkttho1k8N9mBXlR2zvjpO8MBqMnD1RSX1VI70Ge6JSK+zZkIeVRmXqwRzhTs8+rqg1nbPng6TttEas7wHgQcAN6A30BJYC481rmuRmQXe8nLpDpbhM7YWVuy09FkajMkcpZrNo3u73IX8XaOxhyD0Q8yClDXZkbFjNwe3vodfpcOvpR6/B0YSN/qXtiRCCuqwsyld+SeU33yDq6rDp3x/vF57HeepU1E43LgTYLHmRklFIWpPkxYje7jwxvg+TwjtO8uLXCeb8QyXoavQ4umsJHmRafDj3TzHYOmpkgvkmoTW/dY9gUnPdDSCEOKYoiqxHk7SZi3tBO47xQ+1g3f7OobYUMv4DqR+ZRPNcAmHiKzD4LnRGKza8/TrHM/agtrKi36gxl8hsG8rLqVizlvKVt/v9vwAAIABJREFUK9EdO4bKzg7nqVNxSUpCGz7ghr8sLyd5EdLBkhdCCErP1ODmY4+iUtj1dS6HfjyNraOG4AgPAiM88A9za3mNdk4WkkmXWITWOAidEKKh+RdEURQroOuurpNYnOZV0FVbTiEE5usFXXQIUv9/e+cdH+WV3f3vndFo1GfURl2ogUAINQQCYcCAsQEDBgzYa29xicu7726cxG+SzcbJtmySbUl2s/tu1q/XZat3MRgb27hiY6pt1JBEE6iAJNQ1oy5Nue8fjzQILEAICbX7/Xz4oHnmmec592o0Z84995zfr6DoT+Do1rSd1/0QR9xKGqvPE+FlwlNKnA4Hi7d+gfTV6/A1azkDKSVdn32GdccrtL/zDrKvD6958wj/3ncJWLsOvd/I6gqklJTUtLEzv5o9RbU0d/YR5OvJAzmx3Jt1a1peDE4wV5ZoGswD6mrpK2KYkxuhKpgVwPAcxH4hxDcBbyHEauCrwJ6xNUsxlZEOFx2HasdGC9rlhLJ3Ne2Fiv3g4QVp2yHnSTqNkRS9t5ei/3gMe28vT/zyJYw+Pmz9x0uCiY6WFmyv7sb6yiv0VVSg8/fHvHUr5u3b8Jo9e8RmXbR1s7tAU2Mr6295cUeKhS2Z0SxPHvuWF9IlETpBXYWN1/6jwJ1gjpkdyIJ18fgHa7+DoMjxKahTTEyGowehQ9OmvhOtWd87wHNyAvToUK02Jg/2xi46DtdiXp+I0Auc7X2jq+rW0wYFv9MihtZKCIiCBX8B8x/C2t7H0V0vc+rgRzgdDhKyFpC17h5iU7U229LlouvoUVp37KD9/Q/Absc7Kwvztm0ErLkLnffIlnqu1vJic1YU6+dFjmnLi4EEc1VxM1UlTSRkWli4Ph57n5Ojr55jhkowT2tGrdVGf7HcS2g5CAmcngjOQTE5uFIL2ndBOJ6RfqPnHJrPaUnnwt9DXwfE5MCqbyGT76bP7sDo40tfQzmnDx8gdcVqMtduJDgqBgB7QwO2Xa9i3bkT+4UL6E0mgh54APO2rRhH2CnY6ZIcPtfEq/k17pYXMUHe/OXKmWzOjCIuZOy/oX/425OcK2ykt9OBTieImGnCHKY5OYOnnqX33aLeVYpJz3B2Md2NtmvpHFoEES+EeEJKuXesjVNMXgaqoG1vVYy+FrSUcG6ftk217F3QGSD1Xsh5AntwCqUf7yP/108RnjiTdV97GktcAk/+6rcYfXyQTicd+/fTumMHHR9+BE4nPjk5hD71FP6r70BnHFnB2em6dnblV7O7cHDLi0i2ZEWTPUYtLwZrMLc1drPiS3MAcLnkZQlm41htFVZMeYbzzvkJsEJKeRZACJEIvAmMm4NQehCTABe0f3RhdLWg+zqh6GUtYmg6Db6hsPwbkP0I7XY9hW+/wfH3f0hPZwdhCTNJnH+paEtns9H4/AtYd+7EcfEi+uBggh9+CPPWrXjGxY3InIb2Hq3lRX4NJy624aETLL8FLS/qK9s4dfiiO8EMEBrrj8PuxMOgZ9VXlP6yYnQYTg7iYynlskGPBbB/8LHxQuUgJhbOTjvtH14g4I5YdF4eONt60fl53vxuGOt5+PRZTd+5xwYRGbDof8HczeChfeP/6Le/Jv/N10hasIisu+8hKjkFHA53tNB54CBIiW9uLubt2/FfcTviKvrP16K7z8m7J+rYlV/DgbJGXBLSok1szoxiQ/rYtLxob+mhqriJ+IxQfE1GSg/UcPCVs8TMDtR0E1KD8TWrFtmK4XPTehBCiC39P64GZgB/RstBbEPLQzw9SraOGOUgJgbSKen8pF8Lus9B8IMpeM8NvsmLSqg6rPVGOvUmICBlI+Q8iSsqm3PHPiXvrd0s2nI/celZdFpbcfT1YrKE03fhAtZXdmLbtQtHYyMeFgume7dgvncrntFRN2yKyyU5WtHszit09DqINHmxKTOKLVlRJFlGt+XF4ARzZXETLbVaC7Q7Hk4hOSccR58TIYRKME8zpJRUd1STX59PQUMBj857lBj/mBFdazSS1BsG/VwPLO//uRG4BU34FZOB3nIr1tfLNS3oJDPmDQkYwm4iEWvvgZJXtPxCXTF4B8KSp2DBX9BrCKLkw/co+NET2BrqMVnCcPT1AeDj40v7kaOc37GDzsNHQKfDb9kyzNu347dsKcLjxtfhzza0syu/ht0FNdTaevAzerhbXuTEB41qI7qeDjs9XXbMFh+62uy8+uN8LcGcZCL33iR3i2wAj9GuF1FMSJwuJ2XWMvLq8yhoKKCgvoCG7gYA/D39WT1j9YgdxHBRinKKm6LphRLs9V2Y7k7AOzV45MnYtotw7Ndw7AXoagJLCuQ8AfO2g6cPUkpefPqrtNRcIGp2CvPXbSJxQQ72qvOaZOfu3ThbWvCIjNDqFrZswRAefuPj6ehlT1EtrxbUcLzahk7AslmhbM6M4s6UcLxH6cNZSklzTSdVJU1UFTdTV25jxrwQt7paVWkz4QkmlWCeRvQ6eyluLKagoYC8hjyKGorosHcAEOYTRlZYFvMt88kMyyTJnIROjDyCHDXJUSFEPPB1II5BEcdEaPetHMStZ6AK2ifdgkeQF872PoRRP/Iq6OpjWlHbid1akVvyWsh5Ehm3lJozJzmx/wNWPfpV9B4enMv7BF9zEJboWNrffRfrn3fQ9dln4OGB/4oVmLdvwzc3F6G/MVsGt9Lef6YRh0syNzKAzZlRbMyIxOI/OoV8LqcLXX9B3Ju/KKKyuBm4pMEcnx6CZcYoJPMVkwJbr43ChkLyG/LJr8+ntLkUu8sOQJI5iUxLJllhWWRZsoj0ixzVe4+m5Ohu4Ndo1dOumzVMMTlxa0G/UY6zpQeh1+G/LHpk21YdfXDydc0x1BwDYwAsfAIW/gXOgFjOHD1I3v97mvryMrx8/chcs4HQGfFE+QdifXkHZa+/jstmwxAbS+jf/A3mzZvwCA29IRNcrkuttN8s1lpphwUYeXRpPFsyo0kOH528wkCCubK4mboKGw/92xI8PPXMygknPj1UJZinEXWdde7lorz6PM5azwLgofMgJTiFL875IpmWTDItmZi9zONsrcZwHESPlPJnY26JYsJib+zCuqec3jM3qQXd2aQtIX32HHTUQVAirP0RZHwBjP60NTbwx689QkdrC4ERUax69KvMWbCI7n0fUfn339QkOw0G/Fevxrx9Gz4LFyJ0NxZmVzR18mp+NbsKaqhu7cbHU8+auaPfSvvCqRYO/rnMnWAOCPEieWE49j4nHp56ZmaHjcp9FBMTl3RRbi3XooOGfArqC6jtrAXAx8OHDEsGa+LWkBWWRWpIKt4eY9+YcSQMx0H8VAjxLeBdoHfgoJQyf8ysUkwoOg7X0lfVhml9An6LR6AFffG4VrtQvAOcvZC4Cu75OSSuoqWuluaiYmYuzMU/JJS4jGxmLlxMhNEb6ys7qfzmt3B1dGiSnX//95g23bhkZ2tnH28cr2Vnfg2FF6zoBCxJCuHpO2dxZ8rNt9Ie3CJ7zuIIYlKC8PIx4O1vuCzBrFpkT13sTjulzaUUNBRou4waC7D12gAI9gomKyyLL8/9MpmWTGYFzsJDNzlyS8Oxch7wJWAll5aYZP9jxRRkoAraI9gLY2wAptUzCFgZe2PLSU4HnH5L241UdQgMPpD5Rch5Ahkyi/PFReT98LtUFBzDx2QmIWsh9PSw0ByG9fs/oLK0FGE0ErDmLk2EJyvrhj5gex1OPjzVwK78Gj483YDdKZkd7s8/rJ3NPRlRhJtuLq/gtLso/OC8O8Hs1mBOCQK0wrVNf511U/dQTFw6+jooaixy5w+Km4rpdWrfn+MC4lgZs9KdP4jxj5m0Xw6G4yA2AwlSyr6xNkYx/gzWgvbJDsMYG3Bj+gzdrVpB26fPge08mGPhzn/RnIN3IBdOFLPvh1+n6XwlPiYzi7c+QHLkDOq//W3a3tqL7OrCOGsWYc88g2nDevQm07BvLaUk/3wru/JreOP4RWzddkL9jXxlcRxbsqJJiRx5Atje56TmdCu9XQ6Sc8LReQiKP6rBJ8DTrcGsWmRPXZq6m8ivz3c7hNOtp3FJFzqhY3bQbLbN2sb8sPlkWDII8Q65/gUnCcNxEEWAGWgYY1sU48hlWtC+BgK3zsQn6wbWyRtPa9FC0ctg74K4pbDmXyF5HZ1tbbi6nPh7g8FTS8iu/srjhDda6fjtn2g4cwbh40PAurUEbtuGV1raDX3jOt/cxa6Cal4tqKGquQsvg4675oazOTOK25JC8BhhK+3BCebq06047S7MYT4k54QjhOCBb+fg6TU5lgoUw0dKSVVblTuZXNBQwPn28wB46b1IC03j8bTHybRkkh6ajq9h6rZIH867Oww4JYT4jMtzEOO+zVUxenTlN9B5rA6/JVHuVhnXxWnXlpGOvQDlH4LeCGnbIOdJCJ9H4/lK8n7135w6+BGzb7udu558igBrG3cIH9r/4Z9p7u3Fa+5cwr/9bQLW343ez2/Y9g4l0bk4IZivrUhi7bwI/EaQV3A5XdRXtBGeaEIIwbG9lZw4UEtAqDdzl0YSlxpC5MxLu0uUc5gaOF1OTrWcckcH+Q35tPS0AGA2msm0ZLI9eTuZlkzmBM/BoBu7Nu0TjeHUQSwf6riUcv+YWHQDqDqIm6O33IrL7sI7OQjpcOFo6cFg8bn+C63nIe8lKPgtdNRDQDRkPwTzHwbfECqL8vns9Z2cLynCw2gkZdFtJLo8kG+9TV95OTpfXwI2btCihZThN5brc7jYf6aRXfnVfHDykkTnlqwoNmVEEWm+8Z0gPR12qkqbqSpp5nxpM71dDre6mrWhC+mSKsE8BanrrONQzSEO1x7m6MWjtPW1ARDlF0WWJcudP4g3xU/J3/1o6kGMuyNQjC4Oay+2t8rpPt6E54wAvGYFIjx013YOTofWWvvY83D2fRACZt4J2Y9A0h3Y7XY8PI0IoDz/M1pqq8lZsoLI8vP0/fIFeu12vDMyiPj+9wlYuwadzzAcEVq4X1Rt49X8al4vqqW1y06wrycPLoplS+aNS3RKKXE5JHqDjtqzVnb/JN+dYI5P01pkmyyaozEPx1kqJgVd9i6O1R/jSO0RDtUeosJWAYDF28LK2JUsjljM/LD5hPmq7ceDGY4eRDuXNKg9AQPQKaVUJZ+TDGl30X6gmvYPL9eCvuYHrK1GSzrn/wbaa8E/Apb/HWR+CcwxtDc3UfDybyl+/202Pv1NIiwRzO52En2uDueHn2I3mTDffz/mbVvxmjV8oZrq1i52F9SwK7+G8qZOPD10rE4J496sKJbOvDGJTnufk5pTrVSWNFNV3MSc3AgWbkggNNaf+eviiEsNwTLDXyWYpxAu6eJM6xkO1x7mcM1h8hvysbvsGPVGssOy2TpzK7mRuSSaE6dkhDBaDCeCuKykVAixCVh4ldMVE5ieslba3q3COzUY07praEG7nHD2A8h7Ac68rXVWTVoF634Es9aA3oO6s2fI+82POHP0INIliUtIwvbzX9Bx4DA4HPgsWID561/Df/VqdF7D21La1mNnb/FFduXX8EmFtga8MD6IJ5YnsHZeBAFeN7b2K6Xknf9XQmVxM85BGsyh/e0sDJ56cjYk3NA1FROXpu4mjtQe4XDtYY7UHqG5R2tlMjNwJg/OeZDFkVqUYNSryvXhMqJmfUKIo1LKRWNgzw2hchDXx97Yhb22E5/0UKSU2Ks78Iy5ShuJ9jotr5D3G22Lqq9F2546/ysQGOc+zemw8+xXH8be00NScBiRxacwVteiDwzEtHkz5q1bMSbED88+p4sDZY3syq/hvRP19DpcJIT4sjkzik2ZUcQEDW+Zx+V0UVfeRlVJE+0tvdz56FwA9v/hNDq9UBrMU5A+Zx/5Dfluh3Cq5RQAgcZAFkcuJjcyl8WRi7H4WMbZ0onHqOUgBulCAOiAbC4tOSkmKK4eB237ztNxsBadrwHvlGCEQfd55+ByQcVHWm7h9F5wOSB+Odz5XUi+Gzw86e3qouTN3Zw79glbnv5HOvftI6db4pF3AoM8ge+SJZif/j/4r1o1LBEeKSWltW3szK9mT1EtTR19BPoYuG9BDJszo8iIMQ877K8500rpgVp3glnTYDbjtLvQG3QsfyB5BLOnmIhIKaloq+BwzWEO1R4irz6Pbkc3HjoPMi2ZPJX1FLmRucwOmn1TnU4VlxjOPr3BuhAOoBK4Z0ysUdw00iXpKmzAtrcCV7sdn2xNC1pc+c25oxEKfwd5L0JrJfgEw6KvwvyHIDgRAFtDHfl791Dy4bv0dXcT6u1H6apVeNraCIiJwfz1r2PatAlDRMSwbLto62Z3QS278qspa+jAU69j1RwLmzOjuD3ZgqfHtf+oB2swz14cga/JSGtdF9WnWtwJZqXBPLWw9do4evGoO7lc11kHaNXKm5I2sSRyCdnh2VO6FmE8GU4O4uFbYYhidHA0dNG64wyGaH9Cvjz38ohBSqg8oNUtnNwDLjvMuA1W/hPM2eCW7wSoO3uGPzyjiQZG2yUx5dUESkHAXXdi2nIvPguyh9Uor6PXwdsldezKr+ZIeTNSQvaMQL6/OZX18yIxXadK22l3ceFUC5XFWoK5o1UrxTGF+pA038Kc3Ajm3hapEsxTBIfLQXFTMYdqDnGk9gglzSW4pAt/gz85ETk8nvY4uZG5RPnduDKg4sa5qoMQQvzzNV4npZTfGwN7FCPA2dFHz+lWfOeHYQj3JfSJNDxjB7V96GqBwj9oSefms+BlhoWPadFCqLYE43Q4OHPwI/q6ukgweGPfuZNZdS1ENdkInJuK6R++ScDatej9r98G2+mSHDzbxKv51bxdWkeP3UVskA9PrZrJ5swoZgRf+9teW3M39l4nwZF+dHfYefMXx/Ew6omdE8SC9cHMmHupRbb+OlGHYuJT3V6t7TaqPcwnFz+hw96BTuhIDUnlibQnyI3MJTUkddI0uJtKXGvGO4c45gs8CgQDykGMM9Ip6Txai+2980i7E68kM3qTEWOcqV/T+YiWWzjxmtZFNSYHlv0tpNwDBm2vf3dHO8ff20vBW6/R2WYjsNeB4VQVHsHBLLhnE+YtmzEmJQ3LnhO1bbxaUM1rhbU0tPcS4OXBlqxo7s2KIis28Kp5hcEJ5sriZlpqO4lL09TV/AKNbPk/WVhmBKgE8xSh097Jpxc/dTuFgTYWEb4R3BV3F7mRueRE5GAyDr8Pl2JsGNYuJiGEP/AUmnP4M/ATKeW49WYSQmwANiQlJT1WVlY2XmaMKz3nrFhfP4ejvutyLejuVij6kxYtNJ7SxHjS79eihbC5l12jaO8ePvztczidToLbu4hvbic+awGBW+/Fb9kyhOH620prrN280S/ReaquHYNecHuyhXuzolgx24LRY2h1N3uvE4NRe273f+ZTc9rqTjDHzQsmbl6IW4NZMblxSRcnm09yuFZLLhc1FOGQDrw9vFkQvoDcyFxyI3OJC4hTNQm3iFGRHBVCBAF/AzwIvAT8VErZOmpW3iTTdZurq8vOxX//FJ2vAfPdCXilBCFq8jSnULILHN0QNV9rfZG6BTy1JR0pJVXFhRjrG5H7PqJq/z7O+3gy08ufGZu3YNq4cVjKbI3tvbxVfJE9RbUcq9LeDhkxZrZkRbE+LZIg38/vZBqswVx5vJmmmg4e+dFtGDz1lBc24nJKlWCeQtR31nPk4hEO1xzmyMUjWHutAMwJmsPiyMUsiVxChiUDT/0IFAkVN81Nb3MVQvwI2AI8C8yTUnaMon2KG0TaXXQdb8Qny4LOx0DII6l4BkvEqVfgVy9AfQl4+mnRQvbDEJHufq2jr4+St14n//WdtHa2E9doJdXWQ9y6daRv2Yx3RsZ1v7lZu/p4p7SO14tqOXKuGZeE2eH+/O1dyaxPi7hmXqGyuIn9fzjtTjCHxvqTcUcMLocLPPUkZNyYXKhi4tHj6CG/Pp9DtVp/owE5zWCvYJZGLSU3KpdFEYumVCvs6cC1vq49jda99RngHwd9gAi0JLVqtXELkFLSc6IF65uaFrRHsBdGz3KMx5+H4p1g74TwNFj/nzBvGxgvJZFlXx8Hf/oTivKO0Ctd+Hf3ku0XROr/fpjAtWuv2w+po9fB+yfq2VNUy8dljdidkrhgH/73iiQ2pEcyK+zzCeu25m6qijV1tfSV0cT2J5QtMwK0BHNqML4mVck62bH12jhnPUdxUzGHaw+TV59Hr7MXg87A/LD5bEzcSG5kLrMCZ6llo0nMVR2ElFJlBMcZe0MX1jf6taBDvQhZVovxnWfgYqGm0JZ6rxYtRGZpzfP6uXjoIGL/x7TteYM6Hz1mP1/mZS1k9lcewRgXd8179tidfHS6gT1FF/ngVD09dhcRJi8eXhLPhrTIIZvj2fucHHuzwp1gBggI9aavxwlAaIw/a5+cN7qTo7gltPS0cM56jnJrOedsl/5v6m5yn5NgSmDbrG3kRuaSHZ49YfWVFTeOWvCdoEiXpOnFUlwdPZjij+PX9EPEp1awzIV1P4a07eB1aZeHo7WV0heeo+joARqFi0UVdczIXcqaLZvxW7IEoR86WQxau4uDZU3sKarl3RP1dPQ6CPHzZHt2DBvSI5kfG4huUJ3BQItsp91Fym2ReBh0lH3WQECoN0u2RjAjVWkwTyaklDT3NHPOek5zBrZy98+tvZdSjj4ePiSaE1kSuYREcyKJ5kSSA5NVB9QpjHIQEwjpknSXNOGd5I04/RpBPu/i0fkh+sZemLtZSzrHLHRHC9Llou3gQQp+/xKn6qvpNBrwAubPTiP9J/+DX3T0Ve/ldEk+qWhmT9FF9pZcxNplJ8DLg7vnRbAhPZJFCUGXKbG11HZScbyRyuPN1FdoGswhMX6k3BaJEIIHv7cI/QiV2xS3BiklDV0N7kjgrPWs2xkM6CEA+Bv8STAnsCJ2BYmmRLczCPMJU05/mqEcxAShr7od684S+i46CPT+Fb5yD8aQWbDmG1ri2Sdo0Lk1tOzcScdrr9Fz8SLH58YRYDaTu3Yjc++9D73H0L9WKSUFF6y8XljLm8UXaWzvxcdTz+qUMDakRbJ0Voh7W6q9z0nVyWZi5wYhhKDowwucOFD7uRbZAyjnMHGQUlLXWcc52zl3JDDgFDrsl/aamIwmEk2J3BV3F4nmRBJMCSSaEwn1DlWOQAGMsJvrRGEqbHN1Wttp23GIznM+6LBi8vw9PvP8EAsegRm57mjB1dND+3vvU77jZU7WXaDd28ia4BgCt94LmRn4h0cM+UctpeTExTb2FGnbUmus3Xh66FiRHMrG9ChWzrbg7ak5hbambqpKtARzzRlNg3lAXa2tuRu9XueuYFaMPy7poraj1u0ABnIF5bZyuhxd7vOCvIIucwCJpkQSzAkEewUrRzBNGbVurooxoqkM8l6k5UA0vY7Z+Pl+RMDSIHQL/i/4BgP9O5iKi2nZuZMzH31AhZ8nrb7eGEICSV26gsiHn8BwFa2Fc40d7CmqZU9RLecaO9HrBEtnhvA3q2exem4YAV4GXE4XDrsLgOpTLbz2X4UAmAZpMAdFaNtXA4JV4nG8cLqcVHdUfy4/UGGroMfZ4z4v1DuUBHMCm5I2uZeFEkwJBHoFjqP1ismMchC3EkcfnNpDz/73MTS8gV7fhSnpy4h5czBkPgP9ze8cLS3YXn8d285d9JaVUR9iJj8qGP8AM7dv2sq8lXfi6f35LaoXWrp447gWKZy42IYQkBMfxKO3JbAmNZwgX08twVzUzNF+Dea0FdEs3JBAWIKJJVuTVAXzOGJ32bnQfkHbKTRoWajCVkGfq899XphPGEnmJLLDs905gnhTvGpNoRh1lIO4FbSUQ95LOI69ha19M92uB/BPvB3T/bfj6a/tAJEOBx0ff4xt5y7qDx2g0uyLX6iFhd/5Dkl33UnMudPEZ2aj012+G6mhrYc3+6ua889r1aqZsWb+eX0Kd6dFEBagRRhSSl77rwJqTrde0mBODyFyphnQ1NUy7oi9hZMyfbE77VS1VV22bfSc9RyVbZU4XA73eVF+USSYElgcudi9PJRgSsDP028crVdMJ5SDGCucdk2A59jzyHMHaXfeS7vzB0idBwGrYvC/PRcMenorKrDtehXr7t00drVTGRlK3cwodDodGWvWEnjfdgAS5+e4L93a2cfekjr2FNVytEJroT0nIoC/W5PMhrRIwv2M1Jxq5eSeSo7Z+rj7q2kIIQiJ9iM80UTcvBAssUqDeazpdfZSaauk3Na/Y6jfGZxvO49TajUiAkG0fzSJpkSWRS9z5wjiTfH4GFQkpxhflIMYbaznIe8lTbqzox4CorGG/ZLO8+GaFvTdCeg8ndhefw3rrl105+WBXs/Z7Hmc6WnHy8+PnNXrSL9zHf5Bl9oStPfYea+/qvlAWRMOlyQhxJe/XDmTDekRJFn8qSpppvD3Ze4Es8GoJyYlCKfThV6vY8nWmeM4MVMbp8vJWetZChoKKGwspKSphAvtF3BJLcejEzpi/WNJMCVwR+wd7hxBXEAcXh7D0+xWKG41ykGMBk4HlL2rNcsrew+EwB6zHbF8Cx7z78DfaseruRvZXkHDj79H29636e3toTZpBrO/+iTR938Bs62V2MpyUpatwGDUPjB67E72nWrg9cJa9p1uoM/hIsrszaNL41mfGkFQt6SqpJmIfpnP9pYebA1dpC6NYsa8YCKTzKpF9hjR0dfB8abjFDYUUthQyPGm43TatSryUO9Q0kLTWBu/1r1jKC4gTjWmU0w6lIO4GWw1WqSQ/xtoqwG/cFy536CtYy0dx9rx9gwiILYZ22uvYdu1i77KSrpM/lRnp1LZ2YbDYSd8VgLxFgsRFgsRM5Ppc7j44KQWKbx3op7OPieh/kYeWBjLujlhBDTbOV/SzCfvHb+kwZxkJt5sJOW2SFKXKaWt0UZKSU1HDYWNmjMoaCigrLUMiUQndMw0z2R9wnoyLBk4FwsUAAARj0lEQVRkWjKJ9I1U20cVUwLlIG4UlxPO7dOEeM68rQnzJK5ErvkBXV3zsb1zHldHG4ZwB10Hf0XDv74PLhfe8+dTOjeBquoq9F1tzFm6gvnrNhISG4fTJTla3syeolr2ltRh67Zj8jawMT2CVRFBzIsyETYjgI7WXl76ySF3gnlG6uUazDqVUxgV7E47J1tOatFBv1No7G4EwNfgS1pIGqvSV5FhySAtJE0ljRVTFlUoN1za67RoIe83YDsPvqGQ+SWY/xUIjKP9YA22N8oR+na6P3sR+/lihMVCx/IlpD36OJ5xcXz8+xcwGL1IX70WL38T+edb2VNUy5vFdTR19OLrqeeu2WHcbvbHt9lO9YkWOlp73epqAM01HQRF+KoE8yhi7bFeFh2UNpfS69Rak0f5RWmRQWgmGZYMksxJ6HVX72ulUEwGRkUwaKIz5g7C5YKKj+DYC3D6LXA5IH651kE1+W6cPWCva6X72IdYd+3B0eqJoy4fj+VLqZ4RwckzJ+hus/GlH/wMS1wCUkpKa9vYU1TLG8cvUmPtxuih487EENZlR7NitoU3/7OQunKblmCeE8SMeapF9mgipaSircKdOyhoKKCyrRIAD50HKUEppFvSybRkkh6ajsXHMr4GKxRjgKqkvhk6GqHw95D3IrRWgHcQLPpfWrO84ERcdietfzpK9/EenO0NdO37LsbkZHy2305RRzJnjh3F2VRFQtYC5q/bhNUrhN+9e5o9xy9S0dSJQQjWRgTymCkQj/peOvM7ueNBCwaDnux1ceh0gsiZKsE8GnQ7uiltKnVHCIWNhdh6bYDWiygjNIN7ku4h05LJ3OC5akeRQjEI5SAGkBIqD2q5hZN7wGWHGUtg5TMwZwN4GLHX1NDy4xfpPe+HzicMZ2sFhggrwc8/R+jiXOy9PVz46yeZt+ouwhffyUd1kl+9X8upugPoBCxODOahGAt9RxvpO9FNt66HiJlmUm+LRLq0SG5GavA4T8TkpqGrwR0ZFDUWcbL5JA6pFZ/Fm+JZGbNSiw4s6cQHxKtkskJxDZSD6GqBwj9o0UJzmaaxsPAxmP8QhCZrTfLefh/brp30nGvHZ/FfIoztGGa10WAJofD9TxE7zvLQosW09Ap09z3Dz0obKHrxFKEuwVJfP+4xmMndkED6wgiaqtsp6oG4eSHEzAnCU2kwjxiny0mZtewyh1DTUQOAUW8kNSSVh1IfIiM0g/TQdMxe5nG2WKGYXEzPTycp4fxRrW6hdDc4eyEmB5b+D8zdhPTwoqf0BLZffhfbW++A9EXv58K8bTMy0pszneco+nAvvZ2dhMTPxD7nNu5/9gifVlnxcMEWnS932v2gywntDoJjvTF5GQAIifZn1VdSxnkCJicdfR0cbzxOYaPmEI43Hnd3LQ31DiXDksGDcx4kIzSD2UGzMegN42yxQjG5mZ4O4o2/0iIGYwBkfVlLOofNxdHaStsfd2DduYve06fxiMnGd/k/IYy+WL6xAIOPFycP7eez/96F96wMykzp/K7em7hCHSF+ffzVqlncnRbOJ78oJSjOlxmpKsE8UgZqDwYig6FqDzYkblC1BwrFGDI9dzFVfAytlZB6L1JnpPPQIay7XqV93z6w2/GavwyvlK04271w+UNp12H8MmOxzlzOnoILVB2vIsxpItllwGzXLhk508zmp7MA7cNNfVjdGAO1B4MdwoDu8UDtwUDuQNUeKBQ3h9rFdC3il9EnYrH+/Flsu3fjaGhAHxRE0IMP4rdqA62vNmPvcnGq5xNKKj7G7mPibKsX731WgCXAyH0hMXhXdxM1M5AZ84I/1yJbOYfr09rT6nYEhQ2Fn6s9yInIUbUHCsU4My0dRP2//RstL/0GdDr8li0j4Jln8ErOxjgjEIfTxcnXn+f02XdpNphx+a8jWJ9IhkvP/ffM5PZF0fS092Hw1KsE8zCRUlJhq3DnDgobCj9Xe7A9eTuZlkwyQjMI9QkdX4MVCgUwTR2Ez6JF6ENCMG28B1efFw1/KsG27zh/SoCX6/uIshpZ7buZEI8wvMK8SZlvITE9lNAYrUW2yilcm2vVHpiNZlV7oFBMEqalg/BfsQJj5kJOPHuYoGYjPc5OjrceprIthkW3ZXP3+rlEt7tITA9VzmAYXK/2YFXsKjJCM8iwZBAXEKeW4BSKScK0dBC/33eanL1VBOo8OdNWTKnLRUD6ev72rpnEJQeNt3kTmoHag4GlosKGQmo7awFVe6BQTDWmpYNo13uQ56zDEeTH0iceYFWCEnW/GgO1BwWNmkMYqvbgiylfVLUHCsUUZFo6iCeWJSCWJ463GROO69UezAqcxYbEDVoy2ZKhag8UiinOtHQQ6kNN43q1B+mh6dwRewcZlgzmhcxTtQcKxTRjWjqI6cr1ag8WRSxyJ5NV7YFCoZgwDkIIsQm4G7AAv5BSvjvOJk1qXNJFpa3ymrUH9yXfR4YlQ9UeKBSKIRlTByGEeB5YDzRIKVMHHV8D/BTQA89JKf9dSrkb2C2ECAR+DCgHcQN0O7opaSpxRwhFjUWfqz3YlLSJDEuGqj1QKBTDYqwjiBeBnwO/GTgghNADvwBWA9XAZ0KI16WUJ/pPeab/ecU1aOhquGyr6amWU+7agwRTgqo9UCgUN82YOggp5cdCiLgrDi8EzkopywGEEC8D9wghTgL/DuyVUuaPpV2TjWvVHnjpvdy1B5mWTNJC0lTtgUKhGBXGIwcRBVwY9LgayAG+DtwBmIQQSVLK/xnqxUKIx4HHAWJjY8fY1PGhva+d4sbi69YeZFoySQ5MVrUHCoViTBgPBzHUWoeUUv4M+Nn1XiylfBZ4FrR236Ns2y1HSkl1R7U7MihsLFS1BwqFYkIwHg6iGogZ9DgaqB0HO8YFu9POiZYTlzmEz9UezLiDjFBVe6BQKMaX8XAQnwEzhRDxQA1wP/DAONhxS2jtaXU7gitrD6L9olkUsUgTwglNV7UHCoViQjHW21z/CNwOhAghqoFvSSl/LYT4GvAO2jbX56WUpWNpx61C1R4oFIqpxFjvYvrCVY6/Bbw1lve+FajaA4VCMZWZMJXUN4IQYgOwISkp6ZbeV9UeKBSK6YSQcvJuBMrOzpbHjh0bk2s7XA7OWs9es/Ygw5Khag8UCsWkQwiRJ6XMvt55kzKCGAuuVXtg8bZcXnsQlIxBp2oPFArF1GZaOojh1B5sTNzojhAifCPUcpFCoZh2TEsH8Z0j32Fn2U4A/Ax+pIWmuWsP0kLT8DX4jrOFCoVCMf5MSwdxZ9ydpASnqNoDhUKhuAaT0kHc7C6m3Mjc0TVIoVAopiC68TZgJEgp90gpHzeZTONtikKhUExZJqWDUCgUCsXYoxyEQqFQKIZEOQiFQqFQDIlyEAqFQqEYEuUgFAqFQjEkk9JBCCE2CCGetdls422KQqFQTFkmpYNQ21wVCoVi7JnU3VyFEI1A1QhfbgJUCDJ6qPm8nOkyH1NlnJNxHDdj8wwp5XUVyya1g7gZhBDPSikfH287pgpqPi9nuszHVBnnZBzHrbB5Ui4xjRJ7xtuAKYaaz8uZLvMxVcY5Gccx5jZP2whCoVAoFNdmOkcQCoVCobgGykEoFAqFYkiUg1AoFArFkExKPYjxQgjhC/xfoA/4SEr5+3E2aVKj5vPzTIc5mUpjnGxjuVF7J20EIYTQCyEKhBBv3MQ1nhdCNAghSoZ4bo0Q4rQQ4qwQ4hv9h7cAr0gpHwM2jvS+Ew0hRKUQolgIUSiEOHYT15ky8ymEMAshXhFCnBJCnBRCLB7hdSbknAghkvt/3wP/2oQQfzXCa437GIUQfy2EKBVClAgh/iiE8BrhdW7JWIQQT/XbWjrSeb8V9k5aBwE8BZwc6gkhhEUI4X/FsaHk514E1gzxej3wC2AtkAJ8QQiRAkQDF/pPc47Y8onJCillhpQy+8onpul8/hR4W0o5G0jnivfaZJ8TKeXp/t93BjAf6AJevcLGSTFGIUQU8JdAtpQyFdAD919xzoQZixAiFXgMWIj23lovhJg5Ee2dlA5CCBEN3A08d5VTlgOvDXyLEEI8BvzsypOklB8DLUO8fiFwVkpZLqXsA14G7gGq0SYYJuncjZBpNZ9CiABgGfBrAClln5TSesVpU2lOVgHnpJRXdiWYTGP0ALyFEB6AD1B7xfMTaSxzgKNSyi4ppQPYD2yeiPZOmD/KG+S/gL8DXEM9KaXcAbwNvCyEeBB4BNh+A9eP4pKXBW1So4BdwL1CiF8yOQtrroYE3hVC5AkhPleZOQ3nMwFoBF7oX8Z8rn/t1s0Um5P7gT9eeXCyjFFKWQP8GDgPXARsUsp3rzhnIo2lBFgmhAgWQvgA64CYiWjvpEtSCyHWAw1SyjwhxO1XO09K+UMhxMvAL4FEKWXHjdxm6EvKTuDhGzJ4crBESlkrhLAA7wkhTvV/M3EzzebTA8gCvi6l/EQI8VPgG8A/DT5pKsyJEMITbS36H4Z6fjKMUQgRiPbtOB6wAjuEEF+UUv7uihtPiLFIKU8KIX4AvAd0AEWAY4jzxt3eyRhBLAE2CiEq0cKmlUKI3115khBiKZCKtq76rRu8RzWXe/RoPh+yThmklLX9/zegzdfCK8+ZZvNZDVRLKT/pf/wKmsO4jCkyJ2uBfCll/VBPTpIx3gFUSCkbpZR2tG/JuVeeNJHGIqX8tZQyS0q5DG2JqOzKcyaEvVLKSfsPuB14Y4jjmcApIBHNCf4B+JerXCMOKLnimAdQjvaNxBPNw88d7/GO0Rz6Av6Dfj4MrJnu8wkcAJL7f/428KOpOCdoX7Ievspzk2KMQA5QipZ7EMBLaNHfhB0LYOn/P7bfrsCJaO8tfTOOwRvjdoZ2EEuAeYMeG4DHhjjvj2hrlnY0j/vooOfWAWeAc8A/jvdYx3AOE/rfPEX9f2SfG+t0nE8gAzgGHAd2D/EHPOnnpP8DtRkwXeX5STNG4Dv9H6glwG8B40QeC9oXkBP9f3erJurcq2Z9CoVCoRiSyZiDUCgUCsUtQDkIhUKhUAyJchAKhUKhGBLlIBQKhUIxJMpBKBQKhWJIlINQKBQKxZAoB6GY9gghooUQrwkhyoQQ54QQP+1vQTGc174ohNh6k/ffJIT45+ucM08I8eLN3EehuFGUg1BMa4QQAq01w24p5UxgFuAHfH+Ic2+6d9lVrvF3aCIuV0VKWQxECyFib9YGhWK4KAehmO6sBHqklC8ASCmdwF8DjwghfIQQDwkhdggh9qB1vBVCiJ8LIU4IId4ELAMXEkLMF0Ls7++K+44QIqL/+EdCiH8VQuxH0zFh0GtmAb1Syqb+xy8KIf5HCHFACHGmvznlAHu4QudAoRhLJl03V4VilJkL5A0+IKVsE0KcBwYEWhYDaVLKFiHEFiAZmAeEobVLeF4IYQD+G7hHStkohLgPLQp5pP8aZinl8iHuvwTIv+JYHJoeQCLwoRAiSUrZg9b64xvAD29mwArFcFEOQjHdEWh6GNc6/p6UckCUZRnwx/5Io1YIsa//eDJa5833tFUr9Gg9cgb401XuH4GmPTGYP0spXUCZEKIcmA0UAg1A5HAHplDcLMpBKKY7pcC9gw/0K8rFoDU5mw90XvGaqzmUUinl1bSrr7zGAN2A6TrXH3js1X++QnFLUDkIxXTnA8BHCPFlcGv5/gR4UUrZNcT5HwP3CyH0/TmGFf3HTwOhQojF/dcxCCHmDuP+J7m0lDXANiGETgiRiNZt93T/8Vlo3UoViluCchCKaY3U2hlvRvtQLkNrj9wDfPMqL3kVTdylGE3pa3//dfqArcAPhBBFaEtCnxOtGYKPgcz+3VQDnO6/7l7gyf78A2jO6M3hj06huDlUu2+FYpzplzTdI6V8v7/W4Q0p5StXnGNEcxq3SU3oXqEYc1QEoVCMP/+KJt5zLWKBbyjnoLiVqAhCoVAoFEOiIgiFQqFQDIlyEAqFQqEYEuUgFAqFQjEkykEoFAqFYkiUg1AoFArFkCgHoVAoFIoh+f+wyTnc8aDz9gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.loglog(orders, ops[0], label=\"Compressed\")\n",
    "plt.loglog(orders, ops[1], label=\"Full\")\n",
    "plt.loglog(orders, ops[2], label=\"Compressed-new\")\n",
    "plt.loglog(orders, ops[3], label=\"Full-new\")\n",
    "for q in [2,3,4]:\n",
    "    plt.loglog(orders, np.array(orders)**q*ops[0][0]/orders[0]**q, '--', label=\"p**{}\".format(q))\n",
    "plt.legend()\n",
    "plt.xlabel(\"Order (p)\")\n",
    "plt.ylabel(\"Number of operations\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "metadata": {},
   "outputs": [],
   "source": [
    "sym.UnevaluatedExpr = lambda x:x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1049\n",
      "1400\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "#expansion = VolumeTaylorMultipoleExpansion(knl, 4)\n",
    "coeffs = sym.make_sym_vector(\"c\", len(expansion.get_coefficient_identifiers()))\n",
    "new_m2m = translate_from2(expansion, expansion, coeffs, S(1), dvec, S(1))\n",
    "coeffs = sym.make_sym_vector(\"c\", len(expansion.get_coefficient_identifiers()))\n",
    "old_m2m = expansion.translate_from(expansion, coeffs, S(1), dvec, S(1))\n",
    "print(count_ops(cse(new_m2m)))\n",
    "print(count_ops(cse(old_m2m)))\n",
    "for x, y in zip(new_m2m, old_m2m):\n",
    "    print(expand(x-y))"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}