diff --git a/docs/source/api.rst b/docs/source/api.rst index 1a0bb63f1..64abe87ce 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -94,6 +94,16 @@ Continuous nodes predict_precision continuous_node_prediction +Dirichlet processes +------------------- + +.. currentmodule:: pyhgf.updates.prediction.dirichlet + +.. autosummary:: + :toctree: generated/pyhgf.updates.prediction.dirichlet + + dirichlet_node_prediction + Prediction error steps ====================== @@ -161,6 +171,21 @@ Continuous state nodes continuous_node_volatility_prediction_error continuous_node_prediction_error +Dirichlet processes +^^^^^^^^^^^^^^^^^^^ + +.. currentmodule:: pyhgf.updates.prediction_error.nodes.dirichlet + +.. autosummary:: + :toctree: generated/pyhgf.updates.prediction_error.nodes.dirichlet + + dirichlet_node_prediction_error + update_cluster + create_cluster + get_candidate + likely_cluster_proposal + clusters_likelihood + Distribution ************ @@ -238,6 +263,8 @@ Utilities for manipulating neural networks. list_branches fill_categorical_state_node get_update_sequence + concatenate_networks + add_edges Math **** diff --git a/docs/source/learn.md b/docs/source/learn.md index 86553d570..ea4173432 100644 --- a/docs/source/learn.md +++ b/docs/source/learn.md @@ -178,6 +178,17 @@ A generalisation of the binary Hierarchical Gaussian Filter to multiarmed bandit :::: +### Non-parametric predictive coding + +::::{grid} 1 1 2 3 + +:::{grid-item-card} Self-organizing neural network using Dirichlet Process nodes +:link: example_3 +:link-type: ref + +::: +:::: + ## Exercises Hand-on exercises to build intuition around the main components of the HGF and use an agent that optimizes its action under noisy observations. diff --git a/docs/source/notebooks/0.3-Generalised_filtering.ipynb b/docs/source/notebooks/0.3-Generalised_filtering.ipynb index d089871c8..16c0b94fd 100644 --- a/docs/source/notebooks/0.3-Generalised_filtering.ipynb +++ b/docs/source/notebooks/0.3-Generalised_filtering.ipynb @@ -28,7 +28,15 @@ }, "tags": [] }, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.\n" + ] + } + ], "source": [ "import jax.numpy as jnp\n", "import matplotlib.animation as animation\n", @@ -140,13 +148,6 @@ "tags": [] }, "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.\n" - ] - }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+MAAAESCAYAAACW1QRFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAADvC0lEQVR4nOzdeVyVZf7/8dfZ2PfVBRBQZBMEF1RM09xSK9u1Rcts+9avzMZpnWlqmqlpt2naN8tMrTHNEvdUUnEXV0QUBRUQQWST5Wy/PxhOHgFl5+b4eT4e56Hc5z73fX14cx24zn3f160ym81mhBBCCCGEEEII0W7UHd0AIYQQQgghhBDiaiODcSGEEEIIIYQQop3JYFwIIYQQQgghhGhnMhgXQgghhBBCCCHamQzGhRBCCCGEEEKIdiaDcSGEEEIIIYQQop3JYFwIIYQQQgghhGhnNjMYN5vNlJSUILdNF0IIIYQQQgihdDYzGC8tLcXd3Z3z5893dFPEJYxGI4cPH8ZoNHZ0U8QlJBtlk3yUS7JRLslG2SQf5ZJslEuyUbaW5GIzg/FacmRcecxmMzk5OZKNAkk2yib5KJdko1ySjbJJPsol2SiXZKNsLclFZbaRVEtKSnB3d6e4uBg3N7eObo4QQgghhBBCCNEgmzsybjKZOroJ4hImk4msrCzJRoEkG2WTfJRLslEuyUbZJB/lkmyUS7JRtpbkYnODcbmWQnmMRiOpqamSjQJJNsom+SiXZKNcko2yST7KJdkol2SjbHLNuBBCCCGEEEII0YloO7oBQgghhBBCiKuD0WhEr9d3dDM6Fb1ej1arpbKyUo6OdwCdTodGo2mTbdvcYFylUnV0E8QlVCoVvr6+ko0CSTbKJvkol2SjXJKNskk+ytXW2ZjNZvLy8uQ2xM1gNpvp3r07J0+elL7TQTw8POjSpUu93/+WZCKzqQshhBBCCCHaVG5uLufPn8fPzw8nJycZVIpOwWw2c+HCBfLz8/Hw8KBr166tun2bOzIup24oj9FoJCMjg7CwsDY7xUM0j2SjbJKPckk2yiXZKJvko0z6aj3ZR06TdeIEg0cMwsnFsVW3bzQaLQNxb2/vVt321cBsNlNZWYmDg4N8iNEBHB1r+kN+fj5+fn513ruMRmOz389sbgI3mfJfeUwmE+np6ZKNAkk2yib5KJdko1ySjbJJPspiMplY/vUaHhv+LH+5/XU+n72Ax4c/yw/v/4yh2tBq+6m9RtzJyanVtnk1MZvNVFVVYSMnNHdKtT+79c130JL3M5s7Mi6EEEIIIYS4PEO1gQ/+9AXb1+wBwMnVEaPJSGV5FUs+TiJtxxGe+eQJHF0cWm2fclRXdFZt9bNrc0fGhRBCCCGEEA0zm818/tJ3bF+zB52dlgdevpsPk//FzX8dxeNvPYCjiwOHdx5lzsxPMRrkElAh2orNDcbVapsrqdNTq9UEBQVJNgok2Sib5KNcko1ySTbKJvkow4aftpC8NAW1Rs3T//k/xky5Fp2djh49ejBkwgBe/HoW9o527Nt8iKWfrujo5l71VCoVdnZ2cmaBQrXk/UxmUxdCCCGEEOIqUZhXxOwJf6PyQhV3zb6Vmx4cV+96m37Zxod//gq1Rs1ri1+gR0Rgs/dZWVnJ8ePHCQkJwcGh9U57V6oRI0YQFxfHnDlzOropopW01c+wzX0sKbOpK4/RaGTPnj2SjQJJNsom+SiXZKNcko2yST4db8E7P1F5oYqwuFBumD7GsvzSbIbekEDC2H6YjCa+fe2Hq3bysPvvvx+VSoVKpUKn0xEaGsrs2bMpLy9v8DU//fQTr776aqu1ofb2WldrBkrXkvczmxuMy+ycymMymcjOzpZsFEiyUTbJR7kkG+WSbJRN8ulYxw9ms/mX7ahUKu7/yxTUmj+GApdmo1KpuPfZ29HZ6zi0/Qh7NuzvqGZ3uOuvv57c3FwyMzP5xz/+wUcffcTs2bPrrFc707aXlxeurq7N3p/RaLTqI2azmerqahmMK1RL3s9sbjAuhBBCCCGEqGvJJ8sBGDJhAKF9elxxfd/u3lw/dSQAP328vNUHg+Xl5Q0+KisrG71uRUVFo9ZtLnt7e7p06UJgYCB3330399xzD0uXLuXll18mLi6Or776itDQUOzt7TGbzYwYMYKnnnrK8vqioiKmTZuGp6cnTk5OjB8/noyMDMvzc+fOxcPDg19//ZWoqCjs7e3JyspqdntF5yGDcSGEEEIIIWzcqaM57FiTikql4pb/m9Do102cPgY7Bx3H9p3gQMrhVm2Ti4tLg4/bbrvNal0/P78G1x0/frzVusHBwfWu11ocHR0tR8GPHj3KDz/8wOLFi0lNTa13/fvvv5+dO3eybNkyUlJSMJvNTJgwweqe1RcuXOD111/niy++4ODBg/j5+bVae4Vy2dx9xmV2TuVRq9WEh4dLNgok2Sib5KNcko1ySTbKJvl0nJXz1gMwYHRfAnp1q/N8Q9m4e7sx8vZrWPXdelZ+u46YxMh2aa9Sbd++ne+//55Ro0YBUF1dzbx58/D19a13/YyMDJYtW8bmzZtJTEwEYP78+QQGBrJ06VLuuOMOoOYU948++oi+ffvW2YZKpcLe3l5mU1eolryfNfmVycnJ3HjjjXTr1g2VSsXSpUsvu/6GDRsskx5c/Dh82PqTtcWLF1tOy4iKimLJkiVNbRoAGo2mWa8TbUej0RARESHZKJBko2ySj3JJNsol2Sib5NMxyorL2bRsKwDXTx1V7zqXy2bcvTWnqu/ZeIAz2Wdbr11lZQ0+Fi9ebLVufn5+g+uuWGF9+7UTJ07Uu15z/frrr7i4uODg4MCQIUMYPnw4H3zwAQA9evRocCAOkJaWhlarZdCgQZZl3t7ehIeHk5aWZllmZ2dHbGxsvdtQqVQ4OjrKYFyhWvJ+1uTBeHl5OX379uU///lPk16Xnp5Obm6u5REWFmZ5LiUlhcmTJzN16lT27t3L1KlTufPOO9m2bVtTm4fBYGjya0TbMhgMbNmyRbJRIMlG2SQf5ZJslEuyUTbJp2OkJO2kqqKawLBuRA4Mq3edy2XTNdifvsOiMZvNbPhpc6u1y9nZucHHpbePuty6jo6OjVq3uUaOHElqairp6elUVlby008/WU4jv9J2G7rO3mw2Ww2uLzfYNpvNlJWVyQRuCtWS97Mmn6Y+fvz4OtdlNIafnx8eHh71PjdnzhzGjBnD888/D8Dzzz/Pxo0bmTNnDgsWLKj3NVVVVVRVVVm+LikpAWpOFam9/kKtVqPRaOrMSFi73GAwWP1QazQa1Gp1g8svvq4DQKut+fZdGkBDy3U6HSaTyWr6e5VKhVarbXB5Q23vTDXp9XrOnj2LXq+3mZouXd5ZazKbzZw9e9Zqhs7OXpMt5VTbd8xms83UdHHbO3NNRqPRqu/YQk22klNj+k1nq6mWLeRUXV1t1XdsoabOkNPGJVsAuGbSIAwGQ71tr+07JpOpzvdAp9Mx/OYh7P39IL//vJWb/288Go2m0TXp9XpLu8xms1Uba8+aberyS2ewrh3IXjpgbWi5Wq1u1D7NZjPOzs706tXLsrx23w3VVMtsNhMREYHBYCAlJYWhQ4eiUqkoKCjgyJEjhIeHYzKZLK9tqCaTyYTBYMBkMlna2JKa2nN5e+XU1ssvzunS/tGSD0na7Zrx+Ph4KisriYqK4i9/+QsjR460PJeSksKsWbOs1h83bhxz5sxpcHuvv/46r7zySp3l69evx8nJCYCgoCDi4+PZt28f2dnZlnXCw8OJiIhg+/btnD37x6k2cXFx9OjRg+TkZEpLSy3LhwwZgp+fH6tXr7Z60xo5ciSOjo4kJSVZtWHChAlUVFSwfv16yzKtVsvEiRMpKCggJSXFstzV1ZXrrruOkydPWk364OvrS2JiIhkZGaSnp1uWd+aadu/ezdChQ22qps6e07BhwwBYs2aNzdRkaznVsqWabCGnwsJC4I++Yws12VJOtWypJlvLac2aNTZXEygzJ+MFM8f2nUClggsOJSQlJTVYE9ScOv7777/XqSkoths6By2FuUV8+9F8evYNbnRNWq2WgIAAACoqKqiurrasb29vj6OjI+Xl5Vbfd0dHR+zt7SktLbUa0Dk7O6PT6SwH4i7+HqtUqjrL3dzcMJvNVt8XAA8PDwwGg9Us62q1Gjc3N6qrqy0ztOv1eku7KisrrQ4G1i6/tKba9paXl+Pv78+ECRN46KGH+PDDD/H29mb27Nl07dqVkSNHUlJSYlm/oZpq2177b0trgppMXFxc6tRkZ2eHk5NTp8uprWuqrq62/P/S94hhw4bh5eVFc6jMLRjKq1QqlixZws0339zgOunp6SQnJ9O/f3+qqqqYN28en3zyCRs2bGD48OFAzTdo7ty53H333ZbXff/990yfPt3qG3mx+o6MBwYGkpubi7e3NyCflCqlJr1ez5o1axg7diyOjo42UdOlyztrTWazmaSkJMaMGYNOp7OJmmwpp9q+M2HCBEs7O3tNF7e9M+dUVVXFypUrLX3HFmqylZwa0286W021bCGniooK1qxZY+k7tlCT0nNa/J9f+emj5cQOjWL2J4812PbavjN+/Pg6p0vX1vTZX+excfEWht8yhIdevbfRNVVWVnLy5Emr239dvL6Sj7hOnz6d4uJili5dWmf9V155hZ9//pk9e/ZYLb/uuuuIi4vjvffew2w2U1RUxFNPPcUvv/xCdXU1w4cP5/3337dctvvNN98wa9Yszp07V2/bTSYTpaWlloGs0o8iX7zcVo6MV1ZWcuLECUJDQ+tcI242m7Gzs6M52nwwXp8bb7wRlUrFsmXLgJrB+DfffMNdd91lWWf+/PnMmDGjzj0GG1JSUoK7uztFRUUNng4vOobJZOLkyZMEBgbK7KkKI9kom+SjXJKNckk2yib5tC+z2cxTY/5C/qkC/t9bMxh6Y0KD6zYmm7QdR/j71HdwdHHgk01vYefQuAFIZWUlx48fJyQkpM614OLKzGYz1dXV2NnZ1fmgRLSPy/0Mm0ymZr+fdci74ODBg61udN+lSxfy8vKs1snPz8ff37/J25Y3duVRq9X06NFDslEgyUbZJB/lkmyUS7JRNsmnfR3Zc4z8UwU4ONkzYHTcZddtTDbh/Xvh292birJKdv22t5VbKxqiUsmtzZSsJe9nHfJOuGfPHrp27Wr5esiQIVbXrELNufi19+JriktPFRIdz2Aw8Ntvv0k2CiTZKJvko1ySjXJJNsom+bSv7av3ADBgdBz2jpc/it2YbNRqNUMmDARgx9rUVmunuDyz2UxJSUmLJgoTbacl72dNnsCtrKyMo0ePWr4+fvw4qampeHl5ERQUxPPPP8/p06f59ttvgZqZ0oODg4mOjqa6uprvvvuOxYsXW907cObMmQwfPpw33niDSZMm8fPPP7N27Vo2bdrU5ILkh1R5aidjkGyUR7JRNslHuSQb5ZJslE3yaT9ms9ly9HrAdX0btX5jshkwqi/LPl9JavIBDNUGtHbtNh/0VctsNltm85aj48rTrrOp79y502om9KeffhqA++67j7lz55Kbm2s162N1dTWzZ8/m9OnTODo6Eh0dzfLly5kwYYJlncTERBYuXMhf/vIX/vrXv9KzZ08WLVrEoEGDml2YEEIIIYQQV6vTx3I5k30WrU5L7DXRrbbdnrHBePi6cf5sCYe2HyH2mqhW27YQV5smD8ZHjBhx2dH/3Llzrb5+5plneOaZZ6643dtvv53bb7+9qc0RQgghhBBCXGLnupqj4n2GRODo0nqTpqnVavqNiOW3Hzex87dUGYwL0QI2N3vGpVPNi46n0WgYMmSIZKNAko2yST7KJdkol2SjbJJP+6k9Rb1/I05Rh6Zl039U3P/2sU8uOWgHKpUKZ2dnOUVdoVryfmZzg3GZnVN51Go1fn5+ko0CSTbKJvkol2SjXJKNskk+7aMov5ije48D0G9kbKNe05Rs+gwOx97RjnN5RRw/mH3F9UXLqFQqdDqdDMYVqtPNpt6W9Hp9RzdBXEKv17N8+XLJRoEkG2WTfJRLslEuyUbZJJ/2kZq8H6i5vtvL36NRr2lKNnYOdpbT03ev39fsdorGMZlMnD9/HpPJ1NFNEfVoyfuZzQ3GhTLJLUyUS7JRNslHuSQb5ZJslE3yaXv7Nh8CIG54nya9rinZxA2PAWD/lrQm7UN0LJVKxdKlSwE4ceIEKpWK1NTUVt9PcHAwc+bMqXe/bb2vzkQG40IIIYQQQtgIk9HEgS2HAYhNbLvJ1WISIwE4uu84F0or2mw/HUWlUl32cf/993dY21pr8BkYGEhubi59+lz5Q5umDtx37NjBww8/3MIWWps7dy4eHh7tsq/2IjcGFEIIIYQQwkYcP5RNWXE5ji4O9IwNbrP9+Hb3pksPP/Ky8jm4LZ2Bo+PabF8dITc31/L/RYsW8dJLL5Genm5Z5ujo2KTtVVdXY2dn12rtaw0ajYYuXbq06jZr6/T19W3V7V5Oe+6rtdnckXGtVj5fUBqtVsvIkSMlGwWSbJRN8lEuyUa5JBtlk3zaXu1p49GDI9BoGz/Lc3OyiR0a9b99HmpaIzuBLl26WB7u7u6oVCrL1zqdjkcffZSAgACcnJyIiYlhwYIFVq8fMWIE/+///T+efvppfHx8GDNmDADLli0jLCwMR0dHRo4cyTfffINKpeL8+fOW127ZsoXhw4fj6OhIYGAgM2fORK1Wo1KpGDFiBFlZWcyaNctylL4hGRkZDB8+HAcHB6KiolizZo3V85ce7S4qKuKee+7B19cXR0dHwsLC+PrrrwEICQkBID4+3tIOgPvvv5+bb76Z119/nW7dutG7d2+g/qP3ubm5jB8/HkdHR0JCQvjxxx8tz23YsKHO9yE1NRWVSsWJEyfYsGED06dPp7i42FL3yy+/XO++srOzmTRpEi4uLri5uXHnnXdy5swZy/Mvv/wycXFxzJs3j+DgYNzd3ZkyZQqlpaUNfi8vpyXvZzY3GBfK1NRPD0X7kWyUTfJRLslGuSQbZZN82tb+/10vXnsaeVM0NZuYoZH/22fTrhs3m81UXqjqkEdr3IqtsrKS/v378+uvv3LgwAEefvhhpk6dyrZt26zW++abb9BqtWzevJlPP/2UEydOcPvtt3PzzTeTmprKI488wosvvmj1mv379zNu3DhuvfVW9u3bx6JFi9i8eTNPPvkkAD/99BMBAQH8/e9/Jzc31+oI/sVMJhO33norGo2GrVu38sknn/Dss89etq6//vWvHDp0iBUrVpCWlsbHH3+Mj48PANu3bwdg7dq15Obm8tNPP1let27dOtLS0lizZg2//vrrZbd/2223sXfvXu69917uuusu0tIa97OTmJjInDlzcHNzs9Q9e/bsOuuZzWZuvvlmzp07x8aNG1mzZg3Hjh1j8uTJVusdO3aMpUuX8uuvv/Lrr7+yceNG/vWvfzWqLa3J5j6WlElBlMdgMJCUlMSECRPQ6XQd3RxxEclG2SQf5ZJslEuyUTbJp21VXqgifc8xoOmD8eZkE5UQjlqjJi8rn/xTBfgF+DTqdVUV1Uzv92ST2tdavt79bxyc7Fu0je7du1sNBJ944glWrlzJjz/+yKBBgyzLe/XqxZtvvmn5+rnnniM8PJy33noLgPDwcA4cOMA///lPyzpvvfUWd999N0899RQAYWFhzJkzh5EjR/LRRx/h5eWFRqPB1dX1sqeYr127lrS0NE6cOEFAQAAAr732GuPHj2/wNdnZ2cTHxzNgwACg5ohzrdpTwb29vevs19nZmS+++OKKp+HfcccdPPjggwC8+uqrrFmzhg8++ICPPvrosq8DsLOzszpDoSFr165l3759HD9+nMDAQADmzZtHdHQ0O3bsYODAgUDNhxVz587F1dUVgKlTp7Ju3TqrLBrLYDA0+/1MjowLIYQQQghhA9J2HMGoN1qu525rTq6OhPWtOX259oj81cBoNPLPf/6T2NhYvL29cXFxYfXq1WRnW99zvXZQWys9Pd0yGKyVkJBg9fWuXbuYO3cuLi4ulsf48eMxmUwcP3680W1MS0sjKCjIMhAHGDJkyGVf83//938sXLiQuLg4nnnmGbZs2dKofcXExDTqevhL9z9kyJBGHxlvrLS0NAIDAy0DcYCoqCg8PDys9hUcHGwZiAN07dqV/Pz8Vm1LY9jckXEhhBBCCCGuRge31kww1icx8rLXEremPkMiSd99jIPb0hk1eXijXmPvaMfXu//dxi1reN8t9c477/Dee+8xZ84cYmJicHZ25qmnnqK6utpqPWdnZ6uvzWZznVwuPW3eZDLxyCOPWE5Lr11WVlZGz549G93G+k7Hv9LPxPjx48nKymL58uWsXbuWUaNG8fjjj/P2229f9nWX1tkUtW1Sq2uOEV/c7ubcv7u+73F9yy89kq1SqTrkPu4yGBdCCCGEEMIGpO04AkBUQu9GrW80Gjl79izl5eXodLpmDUaiEnqz+MOafTc0ELqUSqVq8aniHen3339n0qRJ3HvvvUDNYDkjI4PIyMtfGhAREUFSUpLVsp07d1p93a9fPw4ePEivXr0sy0wmEyUlJZajz3Z2dhiNxsvuKyoqiuzsbHJycujWrRsAKSkpV6zN19eX+++/n/vvv59hw4bx5z//mbffftuy7yvt93K2bt3KtGnTrL6Oj4+37BdqJnnz9PQEqHMbtabUffLkScvR8UOHDlFcXHzFfDqCzQ3GZXZO5dFqtUyYMEGyUSDJRtkkH+WSbJRLslE2yaftXCir4PihmtOkIwdefjBeXV3Npk2b2LVrFxcuXLAsd3BwYMOGDQwbNgwHB4dG7bdX3xC0Oi3nz5aQl5VP12D/5hfRSfTq1YvFixezZcsWPD09effdd8nLy7viYO+RRx7h3Xff5dlnn2XGjBmkpqYyd+5c4I8jxM8++yyDBw/m8ccf56GHHsLZ2ZlDhw5Zrq+GmlOsk5OTmTJlCvb29pZJ1i42evRowsPDmTZtGu+88w4lJSV1Jou71EsvvUT//v2Jjo6mqqqKX3/91VKTn58fjo6OrFy5koCAABwcHHB3d2/S9+3HH39kwIABXHPNNcyfP5/t27fz5ZdfWr6ngYGBvPzyy/zjH/8gIyODd955x+r1wcHBlJWVsW7dOvr27YuTkxNOTk516o6NjeWee+5hzpw5GAwGHnvsMa699to6lw20FplNXSheRUVFRzdBNECyUTbJR7kkG+WSbJRN8mkbR3Yfw2wy4xfog3cXzwbXy83N5eOPPyYlJYXY2FimT5/OzJkzmT59umWSq48//pjTp083ar92Dnb06hsMQNqOjNYoRfH++te/0q9fP8aNG8eIESPo0qULN9988xVfFxISwn//+19++uknYmNj+fjjjy0DZHv7mjMFYmNj2bhxIxkZGQwbNoz4+Hheeuklq0nL/v73v3PixAl69uzZ4D221Wo1S5YsoaqqioSEBB588MErTk5mZ2fH888/T2xsLMOHD0ej0bBw4UKgZsD573//m08//ZRu3boxadKkxnyrrLzyyissXLiQ2NhYvvnmG+bPn09UVM3t8XQ6HQsWLODw4cP07duXN954g3/84x9Wr09MTOTRRx9l8uTJ+Pr6Wk2OV0ulUrF06VI8PT0ZPnw4o0ePJjQ0lEWLFjW5ve1BZW6N+f0VoKSkBHd3dwoKCvD29u7o5oiL6PV6mTlVoSQbZZN8lEuyUS7JRtkkn7az4J2fWPb5Kq69NZFHX7uv3nVOnDjB999/j4+PD3fccYfldGD4I5uhQ4fy888/k5eXx5QpUxp1nfIP7//Mko+TuOamQTz+5gN1nq+srOT48eOEhIQ0+oj71eKf//wnn3zyCSdPnmxwndrT1N3c3CzXVov2dbmfYb1eL7OpCyGEEEIIcbVK215zvXhDp6jn5eXx/fffExAQwP333281EL+Yu7s706ZNIyQkhIULF3Lq1Kkr7rt2n7XXjYuGffTRR+zYsYPMzEzmzZvHW2+9xX331f/hibB9MhgXQgghhBCiE6u8UEXmwSwAogaG1Xm+oqKCBQsW4OPjw5QpU654GyqdTscdd9xB165dWbhwIaWlpZddPywuFI1WTWFuEWdPFza/kKtARkYGkyZNIioqildffZU//elPvPzyyx3dLNFBZDAu2oVM1KJcko2yST7KJdkol2SjbJJP68tIzcRoMOHTzQvfAOvJvMxmM7/++ivV1dVXHIhfnI1Op2Py5MmoVCp++umnyx7xdnCyJ7RPMPDHEXpRv/fee4+cnBwqKys5cuQIf/3rX6VPXMVsbjAu1x8pj06nY+LEiZKNAkk2yib5KJdko1ySjbJJPm3j8M6aidMiBtQ9Kr5v3z4OHTrEDTfcgJubW4PbqC8bZ2dnbr31Vk6cOMHu3bsv24bI/x2RT9t5dUzi1p7UajUeHh5yvbhCteT9zOYS7YibtYvLM5lM5OfnSzYKJNkom+SjXJKNckk2yib5tI0je44BEN6vl9Xy8vJyVqxYQWxsLNHR0ZfdRkPZhISEEBcXx9q1aykrK2vw9bX7zkjNbE4J4jLMZjN6vV6ux1eolryf2dxgvCU3ohdtw2g0kpKSItkokGSjbJKPckk2yiXZKJvk0/qMBiNH9x4HoHe89czn69evR6VSMW7cuCtv5zLZjBkzBpVKxapVqxp8fa+4UAByMvMoLap/0C4fwjSP2WymvLxcBuMd6HI/uy15P5MLFIQQQgghhOikTmbkUHmhCkdnBwJ6dbUsz8/PZ/fu3YwZMwYnJ6cW7cPJyYkxY8awbNkyBg0aREBAQJ113Dxd6BbiT87xM2TsPU6/ETGW5+zs7FCr1eTk5ODr64udnR0qlapFbbqamEwmqqurqayslFPV25nZbKa6upqzZ8+iVquvOPlhUzV5MJ6cnMxbb73Frl27yM3NZcmSJZe9yf1PP/3Exx9/TGpqKlVVVURHR/Pyyy9bfUI3d+5cpk+fXue1FRUVci9CIYQQQgghGlB7WnivviGoNX8M1NasWYOnpycJCQmtsp++ffuSkpLC+vXrmTp1ar3rhMWF1gzGU49ZDcbVajUhISHk5uaSk5PTKu25mpjNZioqKnB0dJQPMTqIk5MTQUFBrf5hSJMH4+Xl5fTt25fp06dz2223XXH95ORkxowZw2uvvYaHhwdff/01N954I9u2bSM+Pt6ynpubG+np6Vavbc5AXH5AlUelUuHq6irZKJBko2ySj3JJNsol2Sib5NP6aq8Xv/gU9ezsbI4ePcodd9yBRqNp1HaulI1arWbkyJH88MMPnDhxguDg4DrrhMX3ZOOSFI7sqXvduJ2dHUFBQRgMBrlMoYkMBgM7d+4kPDxcZl7vABqNBq1W22DfaMn7WZPTHD9+POPHj2/0+nPmzLH6+rXXXuPnn3/ml19+sRqMq1QqunTp0ujtVlVVUVVVZfm6pKQE+GOCA6h509BoNBiNRqvz/GuXGwwGq2svNBoNarW6weW1261V2xkMBkOjlut0Okwmk9UbkEqlQqvVNri8obZ3tpqGDRtmed5Warp4eWeuacSIERiNRss+bKEmW8pp+PDhNleTLeSkVqsZNmyY5XeOLdRkSzkNGzbssm3vjDWBbeRkNput+o4t1NTROdUOxsPiQi3Lk5OT8fHxITw8HKDRNY0cObLe9Wtr6tmzJ126dGHdunVMmzatTq2hfYIAOLb/BNVV1ajUfwxSatt+6ZHFqyWnltSk0WgYOnSoTdXU2XIyGAxXrKk52v2jFZPJRGlpKV5eXlbLy8rK6NGjB0ajkbi4OF599VWrwfqlXn/9dV555ZU6y1evXm25LiYoKIj4+Hj27dtHdna2ZZ3w8HAiIiLYvn07Z8+etSyPi4ujR48eJCcnU1paalk+ZMgQ/Pz8WL16tdUPyMiRI3F0dCQpKcmqDRMmTKCiooL169dblmm1WiZOnEhBQQEpKSmW5a6urlx33XWcPHmS1NRUy3JfX18SExPJyMiwOmNAapKaWrOmESNGkJqaysmTJ22mJlvLycHBgTFjxthUTbaQU35+Ptu2bbOpmmwpJ7VabXM12WJOUlPLa6osqyL/ZAEqVc1p6snJyZw5c4Zjx47Ro0cPCgsLm1RTTEwM3t7ebNiwocGaHB0dOX78OL/++iu33HKLVU1mkxk7Rx1VF6rYvG4rJfoiyUlqsvmaRowYgbu7O82hMrdgWj6VSnXFa8Yv9dZbb/Gvf/2LtLQ0/Pz8ANi6dStHjx4lJiaGkpIS3n//fZKSkti7dy9hYXXvlwj1HxkPDAwkNzcXb29vQD7ZUUpNer2eNWvWMHbsWBwdHW2ipkuXd9aazGYzSUlJjBkzxnKPxM5eky3lVNt3JkyYYGlnZ6/p4rZ35pyqqqpYuXKlpe/YQk22klNj+k1nq6mWLeRUUVHBmjVrLH3HFmrqyJx2/baP92d+RkCvrrz168sYDAYWL15MXl4ejz76qOVMnsbUVNt3xo8fX+e024trMpvNfP7553h5eXHXXXfVqemtRz9k/+Y07vvLZEZN/uPMyKs5p5bWdGm/sYWabCkns9nc7Ind2vXI+IIFC3j55Zf5+eefLQNxgMGDBzN48GDL10OHDqVfv3588MEH/Pvf/653W/b29tjb29dZXvvGfjGNRoNGU/d6mYZOKWhoeUM3dG/KcrVaXef0nMstb6jtnbGmi0+1sZWaanXWmmrfTJrSb5Re0+WWS01SU0NtbOry2nZf2nc6e022mJPUpMyaLu47tlLTxdqrpsz9J4A/rhcvKSnh8OHD3HDDDVZ/Jzel7SqV6oo1DRkyhF9++YXCwkK8vb2tauod35P9m9M4uvc41997XZNrauzyzpRTY5c3pd80tH5nq6kxy5Ve06WD86Zot7nxFy1axIwZM/jhhx8YPXr0ZddVq9UMHDiQjIyMdmqdEEIIIYQQnUvtRGlh/7vH944dO3B0dCQ2NrZN9xsbG4uzszNbt26t81ztBwMZ9UziJoSw1i6D8QULFnD//ffz/fffM3HixCuubzabSU1NpWvXrldc91IyO6fyqFQqfH19JRsFkmyUTfJRLslGuSQbZZN8Wo+h2kDmgSygZgCs1+tJTU0lPj6+waOOl9OUbLRaLQMHDiQ1NZWKigqr53r1DUGlUpF/qoDzZ4ub3A5Rl/QbZWtJLk0ejJeVlZGammq5WP748eOkpqZaLoZ//vnnmTZtmmX9BQsWMG3aNN555x0GDx5MXl4eeXl5FBf/0TlfeeUVVq1aRWZmJqmpqcyYMYPU1FQeffTRJhfUktnsRNvQarUkJiZKNgok2Sib5KNcko1ySTbKJvm0nhOHT6Kv0uPi7kzXEH/2799PZWUlAwYMaNb2mprNgAEDMJlM7Nu3z2q5k4sjAWHdgD/ugS5aRvqNsrUklyYPxnfu3El8fLxlpvOnn36a+Ph4XnrpJQByc3OtZqn79NNPMRgMPP7443Tt2tXymDlzpmWd8+fP8/DDDxMZGcnYsWM5ffo0ycnJJCQkNLmgiy/iF8pgNBo5fPiwZKNAko2yST7KJdkol2SjbJJP66kd6Naeor59+3Z69+6Np6dnk7aj1+upqKhAr9c3KRtnZ2ciIiLYtWuX1SRX8Mep6kdkMN4qpN8oW0tyafIwfsSIEXU63MXmzp1r9fWGi26N0JD33nuP9957r6lNqdfFM+QJZTCZTKSnp9OzZ896J00QHUeyUTbJR7kkG+WSbJRN8mk9lsF4fCinT5/mzJkzjBo1qlGvzc7OZvfu3WRmZlpu3aTVarG3t2fo0KH079+/UbND9+vXj++++45Tp04RGBhoWR4WF8q6RclyZLyVSL9RNpPJ1Oxc5FwHIYQQQgghOplj/5tJvWdMMKmpqbi5udGzZ8/LvqakpIRffvmFo0eP4uXlRUxMDP7+/mi1WoqKiti+fTtr1qxh06ZNjB49mri4uMteDxsaGoqHhwe7d++2Goz3ig0G4MShbIwGIxqtDCCFqI8MxoUQQgghhOhEys6Xk3+yAICg8O4kfbGMhISEem/zVCstLY1ly5ah0+m44447iIyMtBpo6/V6CgoKGDp0KL///jvLli0jIyODW265pcEJ4VQqFf369SM5OZlx48bh4OAAQNcQfxydHagor+TUsVx6hAe0YvVC2I52u7VZe7ncm5DoGGq1mqCgIMlGgSQbZZN8lEuyUS7JRtkkn9ZRO4u6f5AvJ3Ozqaqqom/fvg2uv2PHDn744QdCQkJ49NFHiYqKqnPEuzYbT09PbrnlFu68806OHj3KvHnz6syYfrG4uDiMRiMHDx602lZInx41bf3fEXzRfNJvlK0ludhconIdhfJoNBri4+MlGwWSbJRN8lEuyUa5JBtlk3xax7EDJwAI7VNzinpQUBDe3t71rrt161aSkpJISEjgjjvuwMnJqd71Ls0mMjKSadOmUVBQwNy5c6msrKz3da6uroSGhtaZVb32VPVj+040vUBhRfqNsrUkF5sbjMssg8pjNBrZs2ePZKNAko2yST7KJdkol2SjbJJP66g92tw9zJ/MzMwGj4ofPHiQVatWkZiYyPXXX3/Z67/ryyYgIIDp06dTUlLCDz/80GBuMTExZGdnc/78ecuy0D7BwB8fHIjmk36jbC3JxeYG4zKbuvKYTCays7MlGwWSbJRN8lEuyUa5JBtlk3xaR+1p6iZHAxqNhujo6Drr5ObmsnTpUmJiYhg9evRlB+LQcDa+vr5MmTKFrKwsfv3113rvqhQZGYlOp7M6Ol57ZDw7/TTVldVNLVFcRPqNsrUkF5sbjAshhBBCCGGrivKLOXfmPCq1ioLyM/Tu3Rt7e3urdSoqKli4cCF+fn7ceOONVxyIX0mPHj246aabSE1NJTU1tc7zdnZ2REREsH//fstg3auLJ+4+bpiMJk6knWzR/oWwVTIYF0IIIYQQopOoPe27Sw8/8gvO1HtUfPny5VRXVzN58uQGZ0Jvqr59+xIXF8eKFSsoLCys83xsbCwFBQXk5uYCNTOt9/zfJG7HZBI3Ieplc4NxmWVQedRqNeHh4ZKNAkk2yib5KJdko1ySjbJJPi13/H+nqLt3dUan0xEWFmb1/P79+zl48CATJ07Ezc2t0dttTDbjx4/H1dWVn376qc51sqGhoTg5OVnNqt4zNgSQwXhLSb9RNplN/SIyy6DyaDQaIiIiJBsFkmyUTfJRLslGuSQbZZN8Wq52YGuwryY8PBw7OzvLc2VlZSQlJdGnTx/69OnTpO02Jhs7OztuvfVWcnNz2bZtm9VzarWaiIgI0tLSLKeqh8b878i4zKjeItJvlE1mU7+IwWDo6CaISxgMBrZs2SLZKJBko2ySj3JJNsol2Sib5NMyZrPZMhjX6yrrnKK+Zs0aVCoV48ePb/K2G5tN9+7dGThwIBs2bKC4uNjquaioKIqKisjLywOg5/9mVM/LyqesuLzJbRI1pN8oW0tysbnBeH0zPIqOZTabOXv2rGSjQJKNskk+yiXZKJdko2yST8ucPV1I2fly1Bo1bv4u9OrVy/JcdnY2+/btY/To0Q3eS/xympLNyJEjsbe3Z+XKlVbLg4ODcXR05NChQwC4errgF+gD/DEDvGg66TfK1pJcbG4wLoQQQgghhC2qHdA6edsTFR2JVqsFam6tlJSURPfu3YmPj2/zdjg4OHD99ddz+PBhjh07Zlleezr1oUOHLAOUnjFy3bgQDZHBuBBCCCGEEJ1A7YBW66ayOkU9NTWVM2fOMH78+BbfxqyxoqKiCAwMZO3atVZHBqOiojh37hxnzpwBoOf/rhvPlMG4EHXY3GBcJjZQHo1GQ1xcnGSjQJKNskk+yiXZKJdko2yST8vUDmhd/JwIDQ0FQK/Xs2HDBvr06UP37t2bve2mZqNSqRgzZgx5eXns37/fsjwkJAQHBwfLqeq1M6oflUncmk36jbLJBG4XkSn/lUetVtOjRw/JRoEkG2WTfJRLslEuyUbZJJ/mM5lMHD+YDUBk/96WAcCOHTsoLy9n5MiRLdp+c7IJDAwkMjKS3377zTKJ1aWnqgdHBqJSqzh/tphzZ4pa1MarlfQbZZNbm11EZhlUHoPBYPUmLZRDslE2yUe5JBvlkmyUTfJpvtwT+VSUV6LSqEgY3h+AyspKfv/9d/r164eXl1eLtt/cbEaNGkVJSQm7du2yLIuKiqKwsJD8/HwcnOwJDKs5Yi+3OGse6TfKJrOpX0RmGVQes9lMaWmpZKNAko2yST7KJdkol2SjbJJP89Weou7opSOsdxgA27dvx2AwMHz48BZvv7nZeHt7ExMTw+bNmy2DktDQUOzt7UlLSwP+uG5cJnFrHuk3yiazqQshhBBCCGHDagey/qG+2NvbU11dzbZt24iPj8fV1bVD2zZ8+HDKysrYvXs3UHOqelhYGEeOHAEg9H/3G5fBuBDWZDAuhBBCCCGEwh3ddxyAqP7hAOzevZvKykoSExM7sllAzdHxPn36WB0d7927N7m5uZSUlNAzNhiouTWbyWTqwJYKoSw2NxiXWQaVR6PRMGTIEMlGgSQbZZN8lEuyUS7JRtkkn+Yx6I1kpZ0EYPCoBIxGI1u2bCEmJgYPD49W2UdLsxk+fDglJSXs27cPgLCwMNRqNenp6QSGdUdnp+VCaQVnss+2SnuvJtJvlK1dZ1NPTk7mxhtvpFu3bqhUKpYuXXrF12zcuJH+/fvj4OBAaGgon3zySZ11Fi9eTFRUFPb29kRFRbFkyZKmNg2Q2dSVSK1W4+fnJ9kokGSjbJKPckk2yiXZKJvk0zynjuZg0BvR2msIjezB3r17KS0tZejQoa22j5Zm4+PjQ3h4OCkpKZjNZhwcHOjRowdHjhxBq9PQIzIQgMz9Wa3W5quF9Btla9fZ1MvLy+nbty//+c9/GrX+8ePHmTBhAsOGDWPPnj288MILPPnkkyxevNiyTkpKCpMnT2bq1Kns3buXqVOncuedd7Jt27amNg+9Xt/k14i2pdfrWb58uWSjQJKNskk+yiXZKJdko2yST/Mc2XMUgC6hvgBs3ryZyMhIfH19W20frZHNkCFDKCgo4OjRmvaGh4dz/PhxqqqqCO3zv0ncDpxojeZeVaTfKFtLctE29QXjx49n/PjxjV7/k08+ISgoiDlz5gAQGRnJzp07efvtt7ntttsAmDNnDmPGjOH5558H4Pnnn2fjxo3MmTOHBQsWNLWJQoHkVgzKJdkom+SjXJKNckk2yib5NN3erQeAmvuLp6Wlce7cOcvf0a2ppdkEBQXRrVs3UlJSCAsLIzw8nJUrV3Ls2DF6xgQDNdeNi6aTfmObmjwYb6qUlBTGjh1rtWzcuHF8+eWX6PV6dDodKSkpzJo1q846tQP4+lRVVVFVVWX5uqSkBKj5ZKL20wm1Wo1Go8FoNFpNFlG73GAwWE1Fr9FoUKvVDS6/9FMPrbbm23dp52houU6nw2QyYTQaLctUKhVarbbB5Q21vTPVVLsPg8GATqeziZouXd6ZatLr9SQlJfHrr7+yf/9+8vLyCAoKYuDAgdxxxx0MHjxYclJITRdvz1Zqurjtnb0m+CMjW6nJFnJqTL/pbDXVsoWcatev/dcWamqPnDIPZAMQ2T+c7du3ExQUhK+vr6WtrVFT7Tpms7nO+k2padCgQSxZsoTTp0/j5+eHn58faWlpDOwzCIATh7KprKhEo9VY1WoLOV26vLVqurTf2EJNtpRTS25t1uaD8by8PPz9/a2W+fv7YzAYKCgooGvXrg2uk5eX1+B2X3/9dV555ZU6y9evX4+TkxNQ8+lcfHw8+/btIzs727JOeHg4ERERbN++nbNn/5hEIi4ujh49epCcnExpaall+ZAhQ/Dz82P16tVWPyAjR47E0dGRpKQkqzZMmDCBiooK1q9fb1mm1WqZOHEiBQUFpKSkWJa7urpy3XXXcfLkSVJTUy3LfX19SUxMJCMjg/T0dMvyzlzT7t27GTp0qE3V1Blz+vTTT1mxYoXV63Jycti6dSsffPABr776Kn/5y186VU21bCmnS9lSTbaQU2FhIQBr1qyxmZpsKadatlSTreW0Zs0am6sJ2ianzKOZnM+tOehUYjhHdnY2sbGxVu1pak21tyJbv349RqMRlUqFTqfD3t6esrIyfv/992bX5O3tjbu7O2vWrMHT0xOVSkVaWhpBQUE4ONlTeaGKRd/8iEdXN5vKqS1/9mp/19T+aws12VJOw4YNw8vLi+ZQmVswlFepVCxZsoSbb765wXV69+7N9OnTLaegQ811Ltdccw25ubl06dIFOzs7vvnmG+666y7LOvPnz2fGjBlUVlbWu936jowHBgZy9uxZ3N3dAflkRyk1mc1mysrKcHV1lSOuCqgpLS2NcePGMXnyZIYPH469vT2FhYWsXr2a3bt3s3PnTpycnDpVTbVsKSfA0nc8PT0xm802UdPFbe/MORmNRoqLi3FxcUGlUtlETbaSU2P6TWerqZYt5KTX6ykrK7P0HVuoqa1z2rJ2Ox/NnIuzuxMTnh1GVlYWjz/+OCqVqkk1VVVVkZ6ezoEDBzh16lSddtfy9fWlZ8+exMXFWQYYTa1px44drF27lscee4yysjLmzp3L1KlT+e5vS0jbcYQH/34Pw28ZYlM5tfWR8Yv7jS3UZEs51S5vjjY/Mt6lS5c6R7jz8/PRarV4e3tfdp1Lj5ZfzN7eHnt7+zrLdTodOp3OaplGo0GjqTvlfG3ojV1+6Xabs1ytVtcbVkPLG2p7Z6rJbDbj6upqaYMt1HQpJddUVlbGmjVrmDBhAmq1mtjYWLKzs9HpdJjNZgwGA1qtlmnTpln+X2v+/PlMmzbNqmYl1FTLlnKqdXHba/vOpcsb03al1tSY5Z2lptr3tYv/IO7sNdlCTo3pN52tpovZQk2X9h1bqOlSrVnTvpT9AIT0CeLQoUOMGDECOzu7ettSX01VVVVs3bqVlJQUqqqqCAkJ4brrriMgIABPT0/s7e0xGo2UlJSQk5NDdnY2e/fuZevWrYSFhTFq1Cj8/f2bVFO/fv3YsGEDe/fuZeTIkbi4uHD06FFC+/QgbccRTqSdYtSd1t+3zp5TW//s1fc7p7PX1NjlSq+pJaept/n8+EOGDLGcUlFr9erVDBgwwFJIQ+skJiY2eX8NfconOo7BYCApKUmy6QDl5eVMnDiRG264gXnz5lmW1/a9S7O5+E1n5syZTJ8+nf/7v/9r0ZuMaD7pO8ol2SiXZKNskk/THd13AgAHLx0qlYp+/fo1+rWHDh3igw8+4Pfffyc+Pp6nnnqKadOmMWTIEAIDA3FxcUGn0+Hg4ICnpydZWVmMHz+eWbNmMWnSJAoLC/n000/59ddfrc5IvRJ7e3v69u3L7t27MZlMhIWFkZGRQc+YmhnVM/efaMq34Kon/UbZWpJLkwfjZWVlpKamWs7PP378OKmpqZbz759//nmmTZtmWf/RRx8lKyuLp59+mrS0NL766iu+/PJLZs+ebVln5syZrF69mjfeeIPDhw/zxhtvsHbtWp566qlmFybE1c5gMHDLLbeQnJyMq6srERERTXr9oEGDUKvVfPbZZ7z++utt1EohhBBCNKS8vJxzp4oBKNafIyYmxjI30uVUVVWxePFifvzxRwIDA3nyyScZN26c5VLOK9HpdMTFxfHYY48xduxY9u3bx6effsqpU6ca3fYBAwZQXl7O4cOHCQsLo6CgAJ8gTwCy0k+hr5bbdAnR5MH4zp07iY+PJz4+HoCnn36a+Ph4XnrpJQByc3OtLowPCQkhKSmJDRs2EBcXx6uvvsq///1vq9sxJCYmsnDhQr7++mtiY2OZO3cuixYtYtCgQS2tT4ir1vPPP8+aNWtwdnZm5cqVDBw4sEmvnzp1Kv/+978BePHFF+U2g0IIIUQ7O7j/ENUl/5sPwcFAQkLCFV9z7tw5vvzySzIyMrj11lu58847cXNza9b+NRoNgwcP5tFHH8XJyYmvv/6a3bt3N+q1fn5+BAUFsWPHDkJDQ1Gr1RRdOIeLuzNGvZHs9NPNapMQtqTJ14yPGDHisqeszp07t86ya6+99ood9/bbb+f2229vanOEEPX48ccfefvtt4GaPjlkyJBmbefxxx/nxIkTvP322zz00EP079+f3r17t2ZThRBCCNGAXb/vAcDBzY6Q3sF06dLlsuufPHmS77//HmdnZx588EF8fHxapR1eXl5Mnz6dFStW8Msvv1BYWMjo0aOtrl+uz8CBA1m8eDElJSUEBQVZrhvft/kQmQeyLPceF+Jq1aLZ1JWkpKQEd3d3zp8/3+hTcET7uHiSsCu9aYuWy8rKok+fPpSVlfHnP/+ZN998s8F1G5ON0WhkzJgxrF+/nn79+pGSktLgxDGidUnfUS7JRrkkG2WTfBrPaDTy5/v+Qu7Oc7gEOPD429OJi4trcP3MzEwWLlxIt27dmDJlCg4ODk3aX2OyMZvNbN26ldWrV9O/f38mTpx42RyNRiPvvfce0dHRuLu7s379ekIdIln2+SpG3DaUR/45rcHXij9Iv1E2s9nc7FzafAI3IQAqKio6uglXjfnz51NWVsbQoUN57bXXrrj+lbLRaDTMmzcPLy8v9u/fz9atW1urqaIRpO8ol2SjXJKNskk+jXPy5EnK8i8A4OLnSHR0dIPrZmZm8v3339OjRw/uueeeJg/Ea10pG5VKxZAhQ7jpppvYtWsXv/zyy2XPmNVoNMTHx7N3716Cg4MxGAw4+zkCcEwmcWsS6Te2yeYG4zLLoPIYDAbWr18v2bST559/ngULFjB37twGb8lQq7HZdO/enQULFrBt2zaGDx/ems0VlyF9R7kkG+WSbJRN8mm8jIwMqopq7qEcPTCiwVtBnT59moULFxISEsLkyZMbXO9KmpJNfHw8N998M3v27GH16tWXXbd///5UVVWRm5uLh4cHem3NoPLU0RwqLzR+hvarmfQbZWtJLm1+n3EhRPtSqVRMmTKl1bc7duzYVt+mEEIIIep3aF8a1aU1M46PuuHaetcpLCxk/vz5dOnShTvvvPOKH8K3pr59+1JVVcWKFStwdXVt8JbEHh4e9OrVi9TUVHr16sWxY8fw8HXn/Nlisg6fJLxfr3ZrsxBKY3NHxoW4Wn355ZcUFBS0y7727NnDokWL2mVfQgghxNXm/PnznMrIBcDR056QsJA661RWVrJgwQKcnZ256667mn1EvCUSEhK45pprWLNmDYcOHWpwvfj4eE6dOoWfnx9FRUUE9O4KyKnqQshgXLSL9vyk9mq0ceNGHnzwQSIjIykuLm7Sa5uazdatWxk4cCAPPfQQZ86cadJrRdNJ31EuyUa5JBtlk3yuLCMjg6pzNae+BkcG1nneZDLx008/UV5ezpQpU3B0dGyV/TYnm+uuu47o6GiWLl3a4N8F4eHhODk5UVBQgEajweV/141n7s9qUXuvJtJvbJPNDcY74lNBcXk6nY6JEydKNm3EYDDwxBNPAHDHHXc06W4CzckmISGB+Ph4SktL+dvf/tbk9orGk76jXJKNckk2yib5NE5GRgaGkpqJ0eISY+o8v379eo4ePcptt92Gt7d3q+yzudmoVCpuuukmvLy8WLhwYb0TjWk0GmJiYjh48CDBwcEY7asByDwgg/HGkH6jbC3JxeYG4yaTqaObIC5hMpnIz8+XbNrIt99+y/79+/Hy8uLVV19t0mubk41arebdd98F4PPPP+fAgQNN2qdoPOk7yiXZKJdko2ySz5Xp9XqOHz9umUm9d7z1NdVpaWls2rSJ0aNH06tX611v3ZJs7OzsmDJlClVVVfz888/1zrAeHx9PeXk5bm5ulBiKAMg9cYbykgstbrutk36jbC3JxeYG40ajsaObIC5hNBpJSUmRbNpAZWWl5ej0iy++2ORPx5ubzbBhw7jtttswmUzMnj27Sa8VjSd9R7kkG+WSbJRN8rmyEydOUFFShf6CEZVaZXWaenFxMcuWLSMyMpIhQ4a06n5bmo2Hhwc33XQT6enp7Ny5s87z/v7+dOvWjaKiIlQ68PBzA+D4QTk6fiXSb5StJbnY3GBciKvJhx9+yKlTpwgICOCxxx5r132/+eab6HQ6Vq1axcaNG9t130IIIYStysjIQF9cc6QtMKwbDk72QM3Rt6VLl2JnZ8eNN96ISqXqyGbWKyIigoEDB7Jq1ap6rx+Pi4sjKysLT09P3Lo4A3Kquri6yWBciE6quLiY1157DYBXXnkFBweHdt1/aGgoDz74IAB//etf6z0lTQghhBCNZzabOXLkCOVna6677hkTbHlu8+bNnDhxgltuuaXVJmxrC2PHjsXHx4f//ve/6PV6q+diYmJqJnBzccHsWDNB3TGZxE1cxWxuMK7ETwmvdiqVCldXV8mmlRmNRiZPnkxsbCzTpk1r1jZams2LL75IUFAQt956q1zH1Aak7yiXZKNcko2yST6XV1hYSHFxMRWFNYPY0D7BAJw+fZoNGzZwzTXXEBwc3Cb7bq1stFott912G+fPn2ft2rVWzzk4OBAZGcn58+fBuebvhswDJ1q0v6uB9Btla0kuKrONHM4qKSnB3d2d4uJi3NzcOro5QrSb6upq7OzsOmz/RqMRjUbTYfsXQgghbEVKSgqrVq3i2NJ8DJVGXvvpRQJ7d+Ozzz5Dq9UyY8aMTvM7d+vWraxatYr777+fHj16WJYfP36cb7/9FpVJTdqiUwB8svkt3L3l73dx9bG5I+NydE55TCYTWVlZkk0baclAvDWy6Sx/FHRG0neUS7JRLslG2SSfy0tLS0NfZsRQaURnryMwrDubNm2isLCQSZMmtenv3NbOJiEhgcDAQJYtW2Z1unpwcDAeHh64erjg7FVziZ1cN3550m+UTWZTv4jMMqg8RqOR1NRUyaaVlJSUcN9997Fr164Wb6u1sjGZTPz3v/9l8uTJ8ouiFUnfUS7JRrkkG2WTfBpWXV3NqVOnqCqquZY6OCKAwnMF/P7771xzzTX4+/u36f5bOxu1Ws1NN91EcXEx69evtyxXqVTExsZy4cIFNG41p/fKYPzypN8om8ymLsRV5NNPP+Xbb79l6tSpipk0rbi4mAceeIAffviB5cuXd3RzhBBCiE4nMzMTs9mMvdEJgJA+PVi2bBne3t4MGzasg1vXPD4+PowcOZKtW7dy6tQpy/LY2FgMBgMOXjoAMvef6KAWCtGxZDAuRCdSVVXFe++9B8AzzzyjmIk8PD09LbdWe+211xTzIYEQQgjRWezduxcAQ8n/foc6GcnJyeGmm25Cq9V2YMtaZsiQIXTt2pVff/3Vcvact7c3AQEBuPrXfPBw7ECW/O0grko2NxhXyuBE/EGlUuHr6yvZtIJ58+aRm5tLQEAAd999d4u315rZzJo1CwcHB7Zu3Sr3HW8l0neUS7JRLslG2SSf+pnNZjIzM9FpdeRl5gNw4mwmAwcOJCAgoF3a0FbZqNVqJk6cyJkzZ9i+fbtleWxsLBoXQAXFBSWcO3O+VfdrS6TfKFtLcrG5wXhn/uTQVmm1WhITEyWbFjIajbz55psAPP30060yg3prZuPv788DDzwAYLn/uWgZ6TvKJdkol2SjbJJP/c6cOUN1dTV+bl2prtSjtdfi4u3Idddd125taMtsunXrxoABA1i/fj2lpaUAREdHo9FpsHev2Z+cqt4w6TfK1pJcbG4wLhMbKI/RaOTw4cOSTQstWbKEjIwMPD09eeihh1plm62dzZ///Gc0Gg1r1qxhx44drbLNq5n0HeWSbJRLslE2yad+W7duBcDLwRcAOw81Y8aMwcHBod3a0NbZXHfddeh0OlavXg2Ak5MTvXv3xsG75uDCMRmMN0j6jbLJBG4XkZmclcdkMpGeni7ZtIDZbOZf//oXAE888QQuLi6tst3WziY4OJh77rkHgNdff71Vtnk1k76jXJKNckk2yib51C8jIwN7e3sKss4B4BfsTd++fdu1DW2djaOjI2PGjOHAgQNkZmYC0Ldv3z8mcZMZ1Rsk/UbZ5NZmQtg4o9HIPffcQ+/evXniiSc6ujmX9eyzz5KYmGg5ZV0IIYQQDSsqKuLChQuEhoayf3saACMmDLPJ64NjY2MJCgpixYoVmEwmwsLCcPF1BGoG4zKJm7jaNGsw/tFHHxESEoKDgwP9+/fn999/b3Dd+++/H5VKVecRHR1tWWfu3Ln1rlNZWdmc5glhc7RaLbNmzeLw4cP4+Ph0dHMuKyoqis2bN3PDDTd0dFOEEEIIxdu0aRMAUZHRFJysOTI+4Np+HdmkNqNSqbj++uspKChg586daLVa4obEolJDeckFzmSf7egmCtGumny1+aJFi3jqqaf46KOPGDp0KJ9++injx4/n0KFDBAUF1Vn//ffft5xeC2AwGOjbty933HGH1Xpubm6kp6dbLWvOdTJqtRzsVxq1Wk1QUJBk0wraYoZTyUa5JB/lkmyUS7JRNsmnrsOHD6PT6UhZtw3M4O7jhpe/R6vvx2g0cu7cOc6ePUtZWRmVlZVUVlaiUqmws7NDp9OhVqvJz8/Hz88PnU7X6m0A6Nq1K3FxcWzYsIGYmBgGDU5gued6Kgv1ZB7IoksPvzbZb2cm/UbZWpJLkwfj7777LjNmzODBBx8EYM6cOaxatYqPP/643mtE3d3dcXd3t3y9dOlSioqKmD59utV6KpWKLl26NLU5dWg0mhZvQ7QujUZDfHx8Rzej03rhhReIiYnhjjvuaPVZNNsym3PnzvHhhx9iMpn429/+1ib7sHXSd5RLslEuyUbZJB9r+fn5XLhwgYCAAPatqzlFvWdMcKt8+K7X68nOzubYsWMcP36c/Px8y7WtGo0Ge3t7y4Gv6upqqqqq0Ov17N69GwAvLy9CQkIICQkhNDQUR0fHFrep1nXXXcehQ4dITk5m7NixuPg6Ulmo59j+EyROHNhq+7EV0m+UrSXjzyb9ZV9dXc2uXbt47rnnrJaPHTuWLVu2NGobX375JaNHj6ZHjx5Wy8vKyujRowdGo5G4uDheffXVy/7QVVVVUVVVZfm6pKQEgMrKSsubhVqtRqPRYDQarS6sr11uMBisrk3RaDSo1eoGl+v1eqs21A6MDAZDo5brdDpMJpPVjHsqlQqtVtvg8oba3plqMhqNHDx4kD59+mBvb28TNV26vK1qOnbsGP/6178wm81ERkZaLu9orZpUKhV79+4lKirK8kbSWjXt3LmTl156CScnJx555BG8vb2tarWlnNqqptq+ExcXh0qlsomaLm57Z85Jr9ezb9++mlvzaDQ2UZOt5NSYftPZaqplCzlVVVVx8OBBS9+xhZpaklNKSgpQ8/eruazm6FpIdBBGo7FZNen1erKysti3bx+HDx9Gr9fj6upKcHAwffv2xcfHBx8fH8uBsovbbjQa2bt3L126dKGgoICcnByOHz/Orl27UKvV9OrViz59+tCrVy+0Wm2LcnJwcCAxMZHk5GTi4+MJjQmm4PA+Du44bPn+Kymnjv7Zu7Tf2EJNtpRT7X6bo0mvKigowGg04u/vb7Xc39+fvLy8K74+NzeXFStW8P3331stj4iIYO7cucTExFBSUsL777/P0KFD2bt3L2FhYfVu6/XXX+eVV16ps3zdunU4OTkBEBQURHx8PPv27SM7O9uyTnh4OBEREWzfvp2zZ/+4NiUuLo4ePXqQnJxsuQciwJAhQ/Dz82P16tVWPyAjR47E0dGRpKQkqzZMmDCBiooK1q9fb1mm1WqZOHEiBQUFljdeAFdXV6677jpOnjxJamqqZbmvry+JiYlkZGRYnb7fmWuqrKxk6NChNlVTW+f0448/YjabSUhIICsri6ysrFatadiwYZw8eZKTJ0+2ek19+/alX79+7N69m1mzZjFlyhSbzakta4KaCW9ycnJspiZbyOns2bOcOnWKU6dO2UxNtpQT1PSboqIim6nJ1nI6deqUzdUETcvJbDZz4MABVCoVBQUFlBfUzJVUcOEM+/bta1JNffr0obCwkPXr11NZWYmdnR3e3t5ce+21hIeHk5SUxJkzZzhz5sxlawIICwuztDEoKIiAgAC6d+/Orl27+Omnn9BoNPj4+BASEsK4ceOanZPJZEKj0fDzzz8z+sYRbP9xHyfTT/Prr8tRq1WKyakpNV0pp+bWtGrVKgDL7xxbqMmWcho2bBheXl40h8rchGkLc3Jy6N69O1u2bGHIkCGW5f/85z+ZN28ehw8fvuzrX3/9dd555x1ycnKws7NrcD2TyUS/fv0YPnw4//73v+tdp74j44GBgeTm5lqOwMknO8qoSa/Xs2bNGsaOHYujo6NN1HTp8raoKScnh969e1NdXc2GDRtITExs9ZrMZjNJSUmMGTPGcm1Ya9b03//+l8mTJ+Pl5cWxY8dwdna2uZwuXt7aNdX2nQkTJlja2dlrurjtnTmnqqoqVq5caek7tlCTreTUmH7T2WqqZQs5VVRUsGbNGkvfsYWampvT8ePHWbBgARqNhi6+XVn33nYAPvz9X7h7uTWqptqzVrdv386FCxeIioqiX79+BAQEoFKpmlRTbd8ZP358ndPka2s6c+YMu3fvJjU1FZVKRXx8PEOHDrU6hb0pOaWlpbFkyRLuvutuXr71HUwGM/9c/DyBvbsrJqem1lRfTpe2vak1XdpvbKEmW8rJbDZfdmx7OU06Mu7j44NGo6lzFDw/P7/O0fJLmc1mvvrqK6ZOnXrFxqrVagYOHEhGRkaD69jb22Nvb19nee0b+8U0Gg0aTd1z+Rs6naCh5Q1NZNGU5Wq1GrW67kX+DS1vqO2dsabaNthSTbXaoqYPP/yQ6upqrrnmGq699tpGt72h5fW1vfbNpCn9pik13XbbbfTs2ZNjx47x7bff8uSTTzar7Q0tV0JOV2pjU5dLTZ2nptrnLn6+s9dkizlJTcqs6eK+Yys1XawxNe3fvx+oOT28Z5ferGM7XXr44eXraVn/crXu37+ftWvXcuHCBeLi4khMTGzw6FxT2q5SqRqsqWvXrkycOJGRI0eybds2tm/fTmpqKsOGDWPw4MFW39PG5BQTE8OuXbtY99s6fII8yc88x/5thwiNDras39E5NdT2xixv7Z+9S3/n2EJNjVmu9JouHZw3RZOmfrOzs6N///6sWbPGavmaNWusjtrVZ+PGjRw9epQZM2ZccT9ms5nU1FS6du3alOYBMpu6EqnVasLDwyWbJigqKuLjjz8GqDNHQ2tq62w0Gg2zZ88G4J133mnRm9XVSPqOckk2yiXZKJvkU6O6upq0tJoJ2yIjIynIKgIgLC70iq/Ny8vjyy+/ZMmSJQQEBPD4449zww03NPs02VpNycbJyYmRI0fyxBNPEBcXx/r16/nPf/7DkSNHmrRPlUrFmDFjOHPmDAHhNX/3b123s1ntt2XSb5StJbk0+ZVPP/00X3zxBV999RVpaWnMmjWL7OxsHn30UQCef/55pk2bVud1X375JYMGDaJPnz51nnvllVdYtWoVmZmZpKamMmPGDFJTUy3bbIr6Pu0QHUuj0RARESHZNMHHH39MWVkZffr0YcKECW22n/bI5v7778ff35/s7GwWLVrUZvuxRdJ3lEuyUS7JRtkknxqHDx+2nF47fvx4ju49Dlx+MG4wGPjtt9/4/PPP0ev13Hfffdx55514eno2+JqmaE42Tk5OjB8/nsceewwfHx8WLFjAkiVLuHDhQqO3ERAQQGRkJGXmYgDyjuU3ue22TvqNsrUklyYPxidPnsycOXP4+9//TlxcHMnJySQlJVlmR8/NzbW6MB6guLiYxYsXN3hU/Pz58zz88MNERkYyduxYTp8+TXJyMgkJCU0u6NLrDETHMxgMbNmyRbJpgkGDBjFs2DCeffbZVr+3+MXaIxsHBwf+/Oc/c//99zNwoNyupCmk7yiXZKNcko2yST41du6sOfrbvXt3nJ2dObrv8oPx3NxcPvvsMzZv3sywYcN4+OGHCQ4ObtU2tSQbb29v7rnnHiZNmsSRI0f46KOPmnSU/LrrrsPsVLPfyvN6MjMym9wGWyb9Rtlakkuz5mB/7LHHeOyxx+p9bu7cuXWWubu7X/YTsvfee4/33nuvOU2pownz0Yl2YjabOXv2rGTTBKNGjWLUqFFt/j1rr2z+9Kc/ten2bZX0HeWSbJRLslE2yadm0uHau5gkJiaSk5nHhdIK7B3tCAzrZrWu2Wxmx44drF69Gl9fXx5++OErztPUXC3NRqWqmQG9V69e/PLLLyxYsIDBgwczevToKx459PHxISFxIFmrV1BdZmDF4jU8/twjzWqHLZJ+o2wtyUUuPBBCwdryqLgQQggh2t/WrVuBmutMw8LCyEitOQoc2icYjfaPQWtlZSU//vgjK1asoH///syYMaPNBuKtycXFhSlTpjBu3Di2b9/Ol19+SVFR0RVfd+211+LoWzPJ85E9x6xmvRbCVslgXAgFWbp0Ka+88gqFhYUd3ZQ2sX//fqZOncqWLVs6uilCCCFEuzObzezevRuVSkVISAg6nY6Meq4XLyws5IsvviAzM5M777yT8ePHNzizsxKpVCoGDx7MjBkzqKys5PPPPycz8/Knnru6uhI9MAKA8vzKK64vhC2wucG4TGygPBqNhri4OMnmCsxmMy+//DIvv/wyn3zySbvss72zef/99/nuu+9444032mV/nZ30HeWSbJRLslG2qz2fffv2UVVVBUDv3r0BLEfGw+JCAMjMzOSLL74A4KGHHiIyMrJd2tYW2XTr1o2HHnqIrl278t1337Ft27bLntI75ubrAKgorGbXzl2t1o7O7mrvN0rXklxUZhu5+KCkpAR3d3eKi4txc3Pr6OYI0WQrV65k/PjxODs7k52d3eJblChReno6kZGRmM1mDh48SFRUVEc3SQghhGg37733HqWlpZjNZmbOnImdxp4HE2ZhNpv5eNNbZBxPZ8WKFYSGhnL77bfj4ODQ0U1uFSaTibVr15KSkkK/fv2YOHFivbeDMuiN3N/vCYx6Iz1v9Oev/3gRe3v7DmixEO3D5o6MyyyDylN7Kw7J5vL+9a9/AfDII4+020C8vbMJDw/n5ptvBuDNN99sl312ZtJ3lEuyUS7JRtmu5nyOHj1KSUkJrq6u+Pn54eHhwdF9xzGbzfgG+JCyYzNJSUkkJCRw9913t/tAvC2zUavVjB07lptuuok9e/awaNEi9Hp9nfW0Og09Y4IBKM+v4tChQ63els7oau43nUFLcrG5wbiNHOi3KWaz2fIpsKjf1q1b2bhxIzqdjlmzZrXbfjsim2effRaA+fPnW2aTFfWTvqNcko1ySTbKdrXmYzabSUpKAqC6utpyinrt/cUdvLTs2LGDiRMncv3119d71Lg92tjW2cTHx3P33Xdz/Phxvv3223rvttQ7vicAlYXVbNu2rc3a0plcrf2ms5DZ1IXo5GqPit97770EBAR0cGva1qBBg7j22msxGAytdktDIYQQQsnS0tIoKirC1dWVyspKy2D8yJ6jAOjtKrn99tsZMGBARzazXfTq1Yv77ruPc+fO8dVXX9WZaT2sb8218xUF1Zw5c4Zz5851RDOFaBcyGBeigx04cICff/4ZlUrFM88809HNaRfPPfccAJ999pn8khVCCGHTaq+Xhpr7aTs5OdG9e3fKy8s5uCMdgFvvnXRVzaPSvXt3ZsyYgclk4uuvv6agoMDyXK++NbPKVxUbMOpNbNiwoYNaKUTbs7nBuMwyqDwajYYhQ4ZINg1wdXXlgQceYMqUKURERLTrvjsqm3HjxjFu3DhefPFFdDpdu+67M5G+o1ySjXJJNsp2NeaTmppKUVERKpWKkpISwsLCqKio4JM5n2GoNKK10zLkuoSObma7Z+Pl5cUDDzyAg4MDc+fOJT8/v2a5vwc+3bzADFXn9Bw8ePCqv+f41dhvOhOZTR2ZTV10fmazGZVK1dHNEEIIIUQr0ev1fPDBBxgMBvz9/Tlx4gQ33XQTKSkpZO/JI+v3M0QMCONv383u6KZ2mPLycubNm0dpaSlTp06lS5cufPCnL9iyfAfefVzwjXFj+PDhjBw5sqObKkSrs7kj4/XNzCg6ll6vZ/ny5ZLNFXTEQFyyUTbJR7kkG+WSbJTtastnx44dlJWVUVFRgbu7O2q1ms2bN1NRUYG/czcAIgb06uBW1uiobJydnbnvvvtwd3fnm2++4fTp00QMCAOgssCAWq1m69atGI3Gdm2Xklxt/aazaUkuNjcYF8okt2KoKycnh6lTp7Jv374ObUdHZmMymViyZAl33HGH/Iw0QL4vyiXZKJdko2xXSz6VlZVs2rQJf39/HBwcKCoqQqPRUFVVxX333ceJAzV3FInoH9bBLf1DR2Xj6OjItGnT8PHxYd68eXj3cAdqZlR3sHegurqa7du3d0jblOJq6TdXGxmMC9FB3n33Xb777jsef/zxjm5Kh6msrOThhx/mv//9Lz/++GNHN0cIIYRoNVu2bKG6upry8nJ69epFdnY2KpWK++67D7VRQ/6pAlRqFWHxoR3dVEVwcHDg3nvvxdvbmzW/r8LJ1RGj3kThqfMAbNy4kerq6o5tpBCtTAbjQnSAwsJCPvnkEwBeeOGFDm5Nx3FycuLJJ58E4I033pD7ZwohhLAJZWVlbN26lYiICEpLS8nMzARg8uTJ+Pj4cHhnzS3NgiMCcXJx7MimKoq9vT333nsv7u7u6DxrhikX8qtxc3OjqqqKlJSUDm6hEK1L29ENaG1arc2V1OlptVpGjhwp2Vzkgw8+oLy8nLi4OK6//voOa4cSsnn88cd544032Lt3L6tWrerQ74fSKCEfUT/JRrkkG2W7WvJJTk5Go9FgNpstp6Z7eXkRGlpzFDx9V81gPLwNrhe/cOECxcXFlJSUUFJSQmlpKZWVlej1eqqrqy3Xt6rVajQaDRqNBjs7O5ydnfH09OTIkSO4ubnh7e2No6Nju89p4+joyNSpU8na8SbF2eWYSlU4OjpSUlLC5s2bGTBgAM7Ozu3apo52tfSbzqoluUiiol04OsqnvrXOnz/PnDlzgJqj4h09g3pHZ+Pl5cXDDz/Me++9x6uvvsq4ceM6/HuiJB2dj2iYZKNcko2y2Xo+RUVF7Nq1i6FDh7Jp0yZUKhV2dnZERkZa1jn8v8F4RL/mD8bNZjP5+fmcOnWK/Px8y+PChQuWddRqNS4uLjg4OGBnZ4ednR06nQ6VSoXRaESv12M0GqmqqqK8vJzy8nKr24g5ODjg4+ODj48P3bp1o3v37vj7+7f5LbacnZ2Z/MDtvLH5P5w/VUpBQQFqtRqTycTvv/9+VX5wb+v95mplc4NxmdxAeQwGA0lJSUyYMEHuKQ289957FBcXEx0dzW233dahbVFKNn/+85/5+OOP2bJlC6tXr2bcuHEd1hYlUUo+oi7JRrkkG2W7GvJZv349zs7OHD58GLPZzLXXXsuGDRuIiIgAoLzkAiePnAYgfEDjJ28zm83k5OSQmZlJdnY2p06dorKyEpVKhbe3N35+fgwcOBA/Pz88PDxwc3PD2dm50R9w187YPWrUKC5cuEBhYaHlcebMGfbt24fJZEKj0dCtWzdCQkLo2bMn3bt3b5PBeZ+Bkejsdeir9Fw4V0lIZA9KSkrYuXMngwcPxsPDo9X3qVRXQ7/pzAwGQ7NzsbnBuBBKdu7cOd577z0AXnnlFdRqmbYBoGvXrjz22GO8++67vPTSS4wdO1aOjgshhOh08vLy2L9/P926dSMnJ4cuXbpQUVGBq6sr3bt3ByB99zHMZjNdevjh4eN22e1VV1dz9OhRMjIyyMjIoLy8HDs7OwIDAxkyZAhBQUF079691QZoKpUKBwcHXF1d8ff3t3pOr9eTl5fH6dOnOXnyJDt27CA5ORk7OztCQkKIiIggPDy81Y7gau20hPUN4dD2I1Sc1VPgW0B1dTUODg5s2LCBm2++uVX2I0RHksG4EO3IwcGBl156iVWrVnHLLbd0dHMU5ZlnnmHr1q0888wzHd0UIYQQolnWrVuHvb09OTk5ACQmJrJu3ToiIiIsHzKn78oAsNxL+1IGg4GjR49y4MAB0tPTMRgM+Pr60rdvX8LCwggMDGzz08Tro9PpCAwMJDAwkMGDB2MymcjNzSUzM5OMjAx+/vln1Go1oaGhREZGEhUVhYODQ4v2GTEgjEPbj2Cnd6S6uhq1Wo2vry979+4lMTERPz+/VqpOiI4hg3Eh2pGTkxOzZ89m9uzZHd0UxfH392fz5s0d3QwhhBCiWbKysjh6tOZa8J49e5KTk4OHhwfFxcWWU9Thj8nbIvpbXy+ek5PDrl27OHjwIFVVVfj7+3PttdcSHR2Np6dn+xXSSGq1mu7du9O9e3eGDRtGaWkpaWlppKWl8euvv7JixQoiIyOJj48nODi4WWe8hf/ve1SScwGnSGdMJhOnT5/Gw8OD3377jSlTprR2WUK0K5XZRu4lVFJSgru7O+fPn8fd3b2jmyMuYjabMRgMaLVaOfVYYSQbZZN8lEuyUS7JRtlsNR+z2cwHH3xAUVERI0aMYMeOHURHR2NnZ8euXbv405/+hEajobpKz4wBT2HQG5iz+h94+Ltx4MABdu3aRW5uLm5ubvTt25eYmBh8fX3bvYbWyqa0tJR9+/axZ88eCgsL8fDwoF+/fvTv3x8nJ6dGb6eirJIHB83CZDQRMyWYgNDuHD9+HB8fHwoKCpg+fTpBQUEtamtnYKv9xlaYzeZm5yIXrIp2UVFR0dFN6FD5+fkkJiby008/Ke5e2krLpry8nDfeeIPRo0dbzeh6tVJaPuIPko1ySTbKZov5/PbbbxQVFREeHo6fnx/l5eX069ePw4cP07t3b8tp5Uf3HsegN+Dm7crew3t47733WL58Oa6urtx1113MnDmT6667rt0H4rVaKxtXV1eGDh3K448/zgMPPEBwcDDJycm89957LFu2jDNnzjRqO44uDgRHBtZsU+VJUVER3t7eFBQU4OLiwrp16xT3d1VbscV+I5o5GP/oo48ICQnBwcGB/v378/vvvze47oYNG1CpVHUehw8ftlpv8eLFREVFYW9vT1RUFEuWLGlO02Q2dQUyGAysX7/+qs7mzTffJCUlhX/9618d3RQrSsymsrKSf/7zn6xbt46lS5d2dHM6lBLzETUkG+WSbJTNFvPJyspi06ZNODk5cccdd5Camkq3bt1Qq9UUFBRY3dJs29qdNf9xNbB792769+/PzJkzueuuu+jdu3eHTuzaFtmoVCoCAwOZNGkSs2bNYvjw4Rw9epRPPvmEefPmkZWVdcVtRA0KB6CywMD58+fp27cvAGVlZWRnZ5ORkdFq7VUqW+w3tqQluTS5xy9atIinnnqKF198kT179jBs2DDGjx9Pdnb2ZV+Xnp5Obm6u5REW9sekFSkpKUyePJmpU6eyd+9epk6dyp133sm2bduaXpEQCnPy5Ek+/PBDAF5++WU5vegKvL29mTlzJgB/+ctf5BePEEIIxSooKGD+/PkA3HnnnZSXl5ORkUF8fDyHDx9Gp9MRGhrKqVOnmD9/Pr8nbQFg4Mh+zJo1izFjxlw1l1c6OTkxbNgwZs6cyW233UZ5eTlz587l66+/5ujRow0e4e4zuOZ6+6wDp7G3t6eyshIHBwe6dOkCQFJSkpxJJzqtJg/G3333XWbMmMGDDz5IZGQkc+bMITAwkI8//viyr/Pz86NLly6Wx8WzQM6ZM4cxY8bw/PPPExERwfPPP8+oUaOYM2dOkwsSQmn+9re/UVlZafngSlzZ7Nmz8fb2Ji0tjblz53Z0c4QQQog6SktLmTdvHkajkd69e9OjRw/27t2LRqOhT58+HD58mICAAH788Ue+/PJLCs+eo6rICMCkeydib2/fwRV0jNrvzyOPPMKUKVMwGo3Mnz+fL774gvT09DqD8vD+vdDoNBTmniPAN4j09HRiYmIoLS0lJCSE4uJi1q9f30HVCNEyTZpNvbq6ml27dvHcc89ZLR87dixbtmy57Gvj4+OprKwkKiqKv/zlL4wcOdLyXEpKCrNmzbJaf9y4cZcdjFdVVVFVVWX5uqSkBKi5B6JerwdqZnnUaDQYjUarT8xqlxsMBqsOr9FoUKvVDS6v3W4trbbm23fpkbuGlut0OkwmE0aj0bJMpVKh1WobXN5Q2ztTTXq93tI2nU5nEzVduryhtqemploGk6+//jomk0lRNdXu4+LtK+Fnz8nJieeff57Zs2fz0ksvMWXKlDp/tFwN/am27wA2U9PFbe/sNV3cd2ylJlvIqTH9prPVVMsWcqrNp/Z1nbWmyspK5s+fT1VVFWazmREjRlBdXc2ePXuIjo4mKyvLcnszLy8vbr75ZgznIMV4AL8AHzz83NDr9YqqqTYbs9lcZ/22yik0NJSePXuSnZ3Nxo0bWbhwIQEBAYwYMYKQkBDUajVaOw29YoJJ330MSjUUFhcyatQoduzYQVxcHHl5eWzatImIiAi6detmk/3p0n5jCzXZUk4tmbegSYPxgoICjEYj/v7+Vsv9/f3Jy8ur9zVdu3bls88+o3///lRVVTFv3jxGjRrFhg0bGD58OAB5eXlN2ibUDGxeeeWVOsvXr19vmaUxKCiI+Ph49u3bZ3UafXh4OBEREWzfvp2zZ89alsfFxdGjRw+Sk5MpLS21LB8yZAh+fn6sXr3a6gdk5MiRODo6kpSUZNWGCRMmUFFRYfUpnVarZeLEiRQUFJCSkmJZ7urqynXXXcfJkydJTU21LPf19SUxMZGMjAzS09MtyztzTXv27LG5mq6U0xNPPIHZbCYxMZGCggIyMjIUV1NMTAxr1qxpdE3tlVOPHj3w9/cnNzeXN954gwEDBrRZTkr/2dPpdGRlZdlUTZ09p+LiYoxGo6Xv2EJNtpaTTqcjPz/fpmqypZzWrFnTaWsymUxkZmZa/iD39PRk+/btlJaWUlRUhKenJwsXLgQgICAAb29vCgsLOb2tAACXbn9sSyk1NZQTtN/PXm5uLjqdjpycHL777jsCAgK44YYbOHjwIDqvmpN5924+hFc/B86cOYOTkxPr1q0jICCAjIwMvv32Wx566KE6t0m1hf5U+7um9l9bqElJP3utUZOdnR3N0aRbm+Xk5NC9e3e2bNnCkCFDLMv/+c9/Mm/evDqTsjXkxhtvRKVSsWzZMgDs7Oz45ptvuOuuuyzrzJ8/nxkzZlBZWVnvNuo7Mh4YGEh+fj4eHh6AfLKjlJpMJhOFhYX4+PhgZ2dnEzVdury+tu/evZuEhAS0Wi179+4lLCxMcTWp1Wry8/Px9PS0TBqjpJ+9hQsXMm3aNDw9PTl+/LjV7VCuhv5kMpk4d+6c5cNKW6jp4rZ35pwMBgP5+fl4e3ujVqttoiZbyan2d07t9aS2UFMtW8ipurqawsJCS9/pbDWZzWZ+/vln0tPTiY6O5uDBgzz66KM4ODgwd+5cCgsLcXJyQqvV0rVrV2677TZL21+87XWy00/x2FvTGXx9f8XUVMtkMnH+/Hl8fHzqXIPdnjmZzWYOHz7Mxo0bOXfuHH369CHQM5j3HvsMFw9nxv0pkbMFZ+nfvz8rV67k8ccf57fffiMtLQ1PT0+mTZtm9feCLfSnS/uNLdRkS+97tdtpjiYdGffx8UGj0dQ5Yp2fn1/nyPblDB48mO+++87ydZcuXZq8TXt7+3qvtal9Y7+YRqOp9xtUG3pjl1+63eYsV6vVlkFPY5Y31PbOVJNer2fHjh1MmDDhsm3vTDVdqr62DxgwgF9++YVDhw4RFRXVqLa3d016vZ5t27YxYcKERveb9szpnnvuYd++fcyYMaPeCW5svT/p9Xq2b99uyccWamrM8s5Qk9lstryvXfx8Z67JVnK6+HdOQ/2ms9V0sc6ek1qtrtN3OlNNq1at4tChQ0ycOJFVq1YxaNAgTpw4wW+//UZFRQWhoaGMHj2azz77jHHjxlm2VVxYQnb6KQD6Do1W5O9cvV7P1q1b6/2bANo3p9jYWKKjo9mzZw8bNmwg7dBhtPZays6X4+3sy8FDBwkICECn07Fv3z7Gjh1Leno6paWl/Pe//2XatGl1aujomhqzvCn9pqH1O0tNtvS+d/HlUU3VpAnc7Ozs6N+/v9UprVBzykRiYmKjt7Nnzx66du1q+XrIkCF1trl69eombVMIpVGpVNxwww0888wzHd2UTkutVvPWW28RERHR0U0RQghxlUtJSWHr1q2MHz+eU6dOodFoSE9PZ/ny5bi5uaHVarnjjjvIyMhAp9NZ3Tno0LaaU2aDwgNw83LtqBI6FY1Gw4ABA3jiiScYNDgBB++agdHRPSewt7fnyJEjxMbGsnv3blxcXEhISABqLn9dvHhxnaP7QihRk2dTf/rpp/niiy/46quvSEtLY9asWWRnZ/Poo48C8PzzzzNt2jTL+nPmzGHp0qVkZGRw8OBBnn/+eRYvXsz/+3//z7LOzJkzWb16NW+88QaHDx/mjTfeYO3atTz11FMtr1CIdlZSUmJ13YloPfXNsiqEEEK0tQMHDrB69WqGDh2Kp6cne/fupaqqCmdnZ2bMmEFFRQUxMTE4ODhw8OBBIiIirI6eHUipuZSzz+Dwjiqh07K3t2fMmDGMubVm8udta3eh1WpJTU1lwIABlJWVcfjwYa655hrUajU9e/bkyJEjLF++XP5mEIrX5MH45MmTmTNnDn//+9+Ji4sjOTmZpKQkevToAUBubq7VhfHV1dXMnj2b2NhYhg0bxqZNm1i+fDm33nqrZZ3ExEQWLlzI119/TWxsLHPnzmXRokUMGjSoyQXJPZyVR6VS4erqetVk8/e//53evXvz/fffd3RTrqgzZfPUU08RGRnJkiVLOrop7aYz5XO1kWyUS7JRts6Yz/Hjx1m6dCkxMTHY2dmxYMECVCoVkyZN4v7776esrIySkhIGDhxIfn4++fn5REdHW23jwNaawXj0kMiOKKFRlJ5Nwqia6+z1RSbUKg3FxcUkJycTGBjIjh07cHJyIjExkaNHjzJ69Gh2797Nxo0bO7jVrUPp2VztWpJLkyZwU7KSkhLc3d0pLi7Gzc2to5sjrlKHDx8mJiYGg8FAUlKS3Fe8Fb300ku8+uqrBAcHc+jQIRwdHTu6SUIIIWzcmTNn+Prrr/H09KS6uppz584BcNttt9GnTx8AvvvuOyorK3nwwQdZv34927ZtY/bs2ZbrTc+eKuDJ0S+i0ar5fOt7OLo4tFr7zGYzJpMJlUpledgqk8nEo0P/TGlRGX/5dhbL1i3BYKi5bW5VVRWPPvoonp6efPDBB/To0YMuXbqwbt06Jk6caHVHFiGUpEkTuHUGcn2I8phMJk6ePElgYGC9kyzYCrPZzKxZszAYDNxwww2dYiDembJ59tln+frrrzlx4gRvv/02f/3rXzu6SW2uM+VztZFslEuyUbbOlE9xcTHfffcdKpWKvLw8AgMDMZlMeHp6Wo58nzt3jmPHjjFp0iTMZjMHDx4kMjLSauKnfVvSAOgZE1LvQLy6upqSkhLLo6ysjMrKSioqKqweVVVVGI1GDAaD5d9L/+6tHZCr1Wrs7Oywt7e3/Fv7fycnJ1xcXOo8nJycOH36tGKzUavVRA8OZ+uKXRzaeoS4uDgOHqyZyC09PZ3vv/+e++67j1GjRvHzzz8zcOBAEhISSEpKwsXFpVPPP9OZ+s3VyGQyNTsXmxuMXzy9vVAGo9FIamoq3bp1s+k3kMWLF7Ny5Up0Oh3vvvtuRzenUTpTNs7Ozrz99ttMmTKF1157jbvvvpuePXt2dLPaVGfK52oj2SiXZKNsnSWfCxcu8NVXX1FeXo69vT2TJk2iurqaFStWMHnyZMsR6O3bt+Po6Eh0dDR5eXkUFhZy/fXXW20rdeN+AHrEdmfXrl0UFhZSWFjI+fPnKSkpqXMbXwcHBxwdHS0PFxcXfHx8sLe3R6vVotFo6vxbe4TcZDJhNpsxGo3o9Xqqqqqorq62/HvhwgUKCgooKyujvLzcar+1t8zq3r07Xl5eeHp64uXlhZeXF97e3g3OLN2e+l4TzdYVu0hNPsD/3TqN7du3M2jQILRaLQcPHuTDDz9k+PDhdO3alVWrVjFjxgzKyspYvHgxU6dOJSgoqKNLaJbO0m+uVkajUQbjQnSk8+fP88QTTwDw3HPPWc2gKlrPnXfeyRdffMHatWt59NFHWb16tU2fkieEEKL9nTt3ji+++IKKigp69erFpEmT0Gq1fPDBB8TFxVnuYV9VVcWePXtISEhAp9ORmpqKs7MzdnZ27Ny5k7y8PHJO57B74z4AMvIOcWr5MTw8PPD29iYoKAh3d3fc3NysHu016DWZTJSXl1NWVkZpaSlFRUXs2rULOzs7Tp06xf79+6murgZqBure3t74+flZHv7+/nh4eLTr7+G+w2ouDcg8kIWLvSseHh7s37+fsWPHcujQIXr06EFycjIuLi6UlJSwf/9+brnlFubPn8+CBQuYPn06fn5+7dZeIa5EBuNCtILnnnuOvLw8evfuzQsvvNDRzbFZKpWKTz75hD59+rB27Vq+++47pk6d2tHNEkIIYQPMZjM7d+5kxYoVmM1mRo0axTXXXAPU3HK3urqakSNHWtbfvXu35f7CP/74I2lpaZjNZr7++mvUajW+vr5QpsVkMOPs4cRTLzypmCPMUDPAdnV1xdXVla5du6LX68nLy7Pcy9psNlNRUUFBQYFlYrr8/HwyMzOpqKgAwNHRkW7dutG9e3e6d+9Ot27dcHFxabM2e/q5ExwVyIlDJ9m3+RB9+vRh586dTJgwgaioKHJycnjooYdISkqipKSEpKQkQkNDmTx5MnPnzmX+/Pk88MADuLu7t1kbhWgKZbwbtCI5SqY8KpUKX19fm83GbDbj7OyMRqPhs88+w8Gh9SZmaWudMZuePXvyt7/9jTfffFMxf9C0lc6Yz9VCslEuyUbZlJrPuXPnWLZsGVlZWQDccsstxMbGAlBQUMC2bdtITEwkJyeHbdu2cfLkSU6ePAnApk2b8PDwwGw2c+211xIeHo6vry9arZbv3vgvAANG9sXf379jimukS7NRqVQ4OTkRFBRkdXq32WymrKys5sh/Tg6nT59m586dJCcnA+Du7k5QUBA9evSgR48eeHt7t2reccP7cOLQSVKTDzDluUls2rSJjIwMhgwZwhdffMH58+eZPn06ycnJbNiwgQ8//JDbb7+de+65hy+//JL58+czffr0TjURrFL7jaghs6kjs6mLjpeVlWW5xZ9oW3q9nqKiIjnVTAghRIuYTCa2bt3K+vXr0Wq1VFZWWmbfrqqqIisri6SkJMrKyizzErm61pweffLkSW666SZiYmJYvHgx58+f55FHHrHa/uyJL3P6WC4z5zzM4Ov7d0SJ7cJsNlNcXMzp06c5deoU2dnZ5ObmWg5Y9OjRg6CgIEJDQ/Hx8WnR4OXI7mP87e43cXZ34tPNb/PV11/h4uLCXXfdxVdffYVareb+++8HYMWKFezYsQOz2UxMTAwDBw5kwYIF+Pr6cu+991rdC16IjmBzh5VkAjflMRqNZGRkEBYWhkaj6ejmtCqz2Wz5hdIZB+KdNRudTmc1EL84B1vSWfO5Gkg2yiXZKJuS8jlz5gzLli0jJyeHnj17cuzYMWJjYykpKeGrr77i9OnTltnKQ0JCiI2NJSQkBHd3d7755hsCAwOJj4/nwoULHDlyhLFjx1ptP/9UAaeP5aLWqIlJVO79xWu1JBuVSoWHhwceHh6Wmearqqo4efIkWVlZZGVlsXr1akwmE25ubvTs2ZOePXsSGhra5CPUvfqG4OzuRHnxBY7uO058fDxJSUmUlpYyePBgfvzxR3JycujWrRujRo0iLS0NFxcXMjIyOHbsGIMGDWLTpk389NNP3HHHHZ1iQjQl9RtRl9FobHYuyv/payK5tZnymEwm0tPTbS6bTZs2MWzYMDIyMjq6Kc1mC9n8+uuvxMXFcfbs2Y5uSquzhXxslWSjXJKNsikhH4PBwPr16/nss8+oqqqiT58+HDt2DLVazb59+9i5cyeurq6MHTsWV1dXwsLCmDZtGnFxcbi7u5OXl8eJEycYNGgQAPv318yWXnvf8Vq1s6iHxYXi7ObUvkU2Q2tnY29vT69evRg1ahQPPPAAzz33HPfccw+RkZGcPHmS//73v7z11lt88cUXrF+/nlOnTtGYE3bVGjWxQ2sG/LvX76dPnz5oNBr27t1LREQEHh4ebN26FQA7OzvGjRtHbm4u48ePJzg4mA0bNuDv7096ejpJSUmN2mdHU0K/EQ1rSS42d2RciPZQVlbGfffdR2ZmJm+99RafffZZRzfpqmQwGHjhhRfYv38/jzzyCIsXL7bJI+RCCCFax8mTJ1myZAnnz5/HxcXFcpsxe3t7EhISCA8Pp1u3bqhUKtavX8+FCxfq3Kpsy5YtuLu7ExlZc7R77969hIWF4ezsbLXejrWpAAy4rm+71KZ0Op2OXr160atXL6DmPu7Hjh3j2LFjbN++neTkZJydnenduzcRERGEhIQ0eBp5/+tiSUnawc51qdz1p1uIjIxkz549DB06lISEBNauXcvo0aNxc3MjKiqK0NBQ1q9fz2OPPUafPn1ISkpCo9Gwa9cuXF1dufbaa9vzWyGEhQzGhWiG2bNnk5mZSVBQEG+99VZHN+eqpdVq+fbbb0lISGDJkiXMmzePadOmdXSzhBBCKIjZbCYzM5O1a9eSl5cH1AwMfX19qaysxN/fn2nTplkN/M6dO8fmzZtJTEzEy8vLsryoqIgDBw5w/fXXo1aryc/PJzc3l2HDhlnts6y4nEPbjwAwYHRc2xfZCbm7u9OvXz/69euHyWTi1KlTHD58mPT0dPbs2YNWq6Vnz56Eh4fTu3dvqw874q+NQaPTkJOZx+nMPOLj49m/fz8nT56kX79+bNy4ka1btzJ27FhUKhUTJkzg448/ZtOmTYwcOZLg4GBWr15NamoqGzZsQK1W18lQiPZgc4PxznDdx9VGrVYTFBRkM9n88ssvfPrppwDMnTu3U98ewxayiYuL45VXXuGFF17giSeeYOjQofTs2bOjm9UqbCEfWyXZKJdko2ztlY/RaCQrK4u0tDQOHDhAZWUlAAEBAQwfPhwPDw++/fZbvL29ueeee6wG4mazmZUrV+Ls7Gy5tVmtzZs34+joSHx8PAA7d+60HM292J6N+zEZTQSGdaNLj84x2WhH9p3afQcFBTF27FgKCgpIT08nPT2dZcuWoVKpCA4OJioqisjISJxdnYkeFM6+TYfYtS6VGx8ch4eHB3v27GHSpEkMHDiQbdu2MWzYMBwdHfH29mbo0KFs3ryZ2NhYvL29mTRpEn369OGHH37gt99+49y5c9x4442KfO+Q9zVla0kuMpu6EE2QlZVFfHw8RUVFPP3007zzzjsd3SRBzenqI0aMYPPmzfTv35/Nmzdjb2/f0c0SQgjRjkwmE5mZmRw8eJD09HQqKirQ6XTo9Xq6dOnC7bffjre3N6WlpZZ7gU+fPr3O6eUHDhxg8eLFTJkyhfDwcMvysrIy5syZw7XXXsuwYcOorq7m3XffZeDAgYwaNcpqG+898Qnb1+zhlv+bwJ0zJ7VL/baqvLycw4cPc+jQIY4fPw7UTJpryFOx9uvN9OobwquLnmPjxo1s3ryZP/3pTxgMBubMmcPQoUMZMWIEUHMnlo8++sjyAUztZW2VlZV89tlnFBUV4evryx133FFzj3gh2oHNfbwis6krj9FoZM+ePZ0+m+rqau68806KiopISEjg9ddf7+gmtZitZKPValmwYAHe3t7s2rWL2bNnd3STWoWt5GOLJBvlkmyUrbXzMZvNZGdns3z5ct555x3mz5/PyZMnCQkJwcnJCZVKxcSJE3n44Yfx9vamoqKC7777DoPBwL333ltnIF5RUcHKlSuJioqyGogDpKSkoNFoGDhwIFAzcVtVVRX9+1vfsqy6spq9mw4CMGBUXKvU2R6U2necnZ3p378/U6dOZfbs2UycOBG1Ws2JomMAHN17nA1rNhIeHo5er+fAgQOW12zbto2qqiqg5tKE8ePHc+zYMQ4dOmTZvoODA4899hj+/v4UFBTwySefkJycrKjvg1KzETVakovNDcZllkHlMZlMZGdnd/psioqKMJvNeHh4sGjRIuzs7Dq6SS1mK9kABAYG8u233wI1t1OxhZpsKR9bI9kol2SjbK2Rj9lsJi8vjzVr1jBnzhy+/vprjhw5Qt++fbn33nvx9/fn0KFDdO/enccff5wBAwagUqmorq7m+++/p7S0lKlTp+Lh4VFn26tXr8ZgMNSZtK2yspKdO3cycOBAHBwcMJvN7Ny5k969e9fZzr4taVRVVOPd1ZOQ6KBm19neOkPfcXJysgzMn3/pWbr2qrkE4Odvl/P555/j7OxMcnIy1dXVJCYmUl1dzc6dOy2vr50cbsWKFVy4cMGyXKvVcv/99+Pj44NWq2X9+vV8/vnn5ObmtnuN9ekM2VzNZDZ1IdqBv78/mzZtIi0tjeDg4I5ujqjHhAkT2LNnD3FxcR3dFCGEEK2ssLCQAwcOcODAAQoKCnB0dCQqKoqYmBgCAwPZv3+/5a4at956K3369LGciqzX6/n+++/Jz89n2rRp9Z6GfPz4cVJTU7nhhhtwdXW1em779u0YjUYGDx4M1MzKnpeXx6BBgzhz5gyVlZVUVVVRWVlJ0vzVAHSN8GXt2rXo9XoMBoPVQ6/XYzKZMJvNVg/A6v9qtRq1Wo1Go7H8v/ZrjUaDTqezPOzs7Or86+DggKOjo+VfOzs7m7nriJOTEyNuvoYFb/+Eh9mXceOuZceOHRQUFPDWW29ZZlFPSUkhISHBMi/AhAkT+Oijj1i1ahW33HKLZXsODg5MnTqVr7/+Gnt7e0wmE59//jmJiYmMGDECrVaGTaL1yU+VEFdQUlJimYfAzs6Ovn3lFiVKdvFA3GAwUFFRUeePKiGEEJ1DRUUFBw4cYO/evZw+fRo7OzsiIiIYO3YsoaGhaDQaiouLWbBgAUePHqVPnz5cf/31VqefGwwGFi5cSE5ODvfeey/du3evsx+9Xs+yZcvo2rUrnp6eHDhwgAsXLlBeXk5xcTH79+/HycmJuXPncuHCBSoqKgD4+eefrbZjMpg4uu0MAEa3Kg4fPoxOp0Or1VoeOp0Oe3t71Go1KpWqwYfZbMZkMlk9jEaj5d/q6mpKSkqorq6muroavV5v+behKaFUKhWOjo6WAbqzszPOzs64uLjg6OjI+fPnOXnyJB4eHri4uCj+LMAh1/dnwds/kb7rGE/0eIiBAwfywQcfoNFoOH36NIWFhQDMnz+fMWPG0K1bN1xdXRk3bhw///wz0dHRVpPvubq6ct999/H1119jNBoZMmQIW7du5fDhw9x0000EBXWeMx1E52BzE7idO3cOT0/Pjm6OuIjRaCQjI4OwsDA0Gk1HN6dJ8vPzGThwIPfccw//+Mc/bG4Wy86czZWcO3eOO++8E7VaTVJSUqf8RNuW8+nsJBvlkmyUrTH5GI1Gjh07xt69e0lPT8dkMtGrVy/69u1L7969LUc4TSYTO3bsYN26dTg4ODBx4sQ613kbDAYWLVrEiRMnmDBhAi4uLpSUlFgepaWllJSUcO7cuTrXfWo0GpydnTEajVy4cIGIiAjc3d2xs7Nj06ZNREdH07dvXxwcHLC3t8fBwYG9Gw/xnz99iU83b/697p8dchTabDZbBuuVlZVUVFRY/r30/xcuXKCsrMzyuPR0W51Oh6urK+7u7ri7u+Pm5mb5f+2joXuBt5eXprxBRmom0164k/HTRrFjxw5WrFjBk08+yYULF1iyZAmFhYWYzWa8vb2JjY0lNjaW5cuXc+bMGR577DEcHBystllUVMTcuXOxs7PjhhtuYO3atZw6dYqBAwcyevTodv+QQt7XlM1oNDY7F5sbjMts6qK1VFVVMWrUKDZv3kzv3r3ZsWOH/Gx1Ivv372fIkCGUl5czc+ZM5syZ09FNEkIIcRl5eXns3buX/fv3U15ejp+fH3FxccTExODi4mK1bn5+Pr/88gunTp1iwIABXHfddej1eoqKiiyPc+fOcfToUcttzWqpVCpcXFxwc3PDzc0NlUrFoUOHiIqKIiEhAVdXV5ydnbGzs6OiooL333+fAQMGMGbMGAA2bdrEhg0bmDVrVp0J4N578lO2r97NTQ+N464/3dq237BWZjabqaystBqcl5eXWz64KC4upri4mNLSUqvXOTk54e7ujpeXF56ennh5eeHl5YW3tzfOzs5t/oHEynm/8c0/FxHWN4S/L3rOMsv9gAEDGD16NEVFRfznP/+hb9++GI1G0tLS0Ov1BAYGkpubS58+fZg0qe6M94WFhcydOxcnJyemTp3KgQMH+O2333BycuLGG2+0mduoio5lc4PxwsJCvLy8Oro54iIGg4Ht27eTkJDQaY5Omkwm7r33XhYsWIC7uzvbtm2r82m7LeiM2TTFkiVLuPXWmj+G3n//fZ588skOblHT2Ho+nZlko1ySjbJdmk95eTn79+8nNTWVM2fO4OTkRExMDH379qVLly51BnJlZWX89ttvpKam4uDggK+vLxcuXOD8+fMYDAbLes7OzphMJiorK4mOjiYsLAwPDw/L6de1Z7rp9Xo+/fRTHBwceOCBB+qcAbdmzRp27tzJzJkzcXJywmg08v7779OrVy9uuukmq3UvlFXwaOJs9NUGXl/yF4IjA9vou9g2Gtt3jEaj1eC8uLiY8+fPWz4AKSkpsayr0+ksg/Pah4+PDz4+Pjg5ObVKu8+fLeaxa5/FbDLz/tp/4hfgw8qVK9m3bx+zZs1Cp9OxfPlyDhw4wMyZM1GpVKSlpZGamkpWVhYAPXv2ZNiwYQQFBVn9zJ09e5a5c+fi7u7O1KlTqays5JdffuH48ePExcUxduxYHB0dW6WOy5H3NWUzGAzNzsXm0rSRzxZsitls5uzZs50mG7PZzNNPP82CBQvQarX8+OOPNjkQh86XTVPdcsstvPbaa7zwwgs89dRT+Pn5MWXKlI5uVqPZej6dmWSjXJKNspnNZs6cOUNaWhoHDhzg6NGjAISHhzNy5Eh69eqFWq2mpKSEY8eOUVBQQEFBAYWFhZw5c8ZyrTZgmbgsJCTE6qisi4sLy5Yt48iRI9x5551EREQ02J61a9dSXFzMlClT6gzEy8rK2L59O0OGDLEMHA8cOEBpaSlDhgyps61d6/airzbQLbQLPSICWuPb1a4a23c0Gg2enp4NXhZae4bCuXPnLI+ioiIOHDhAcXGxZT1nZ2d8fX3x8fHB19fX8n8XF5cmHU338HUnelA4B1IOs2X5Dm5+ZDwJCQls27aN/fv3069fP4YPH05qaiopKSmMHDmSuLg44uLiOHfuHPPmzSMzM5Njx47h6elJ37596du3Lx4eHvj6+jJt2jS+/fZbvvnmG6ZOncrUqVPZs2cPq1ev5ujRo0ycOPGyP2OtQd7XlK0ludjcYFyIlnrzzTd5//33AZg7d67ltDTROT333HPk5OTwn//8h2nTpuHj48Po0aM7ullCCHFVMZvN5OTksGfPHg4ePMi+ffvo1q0bo0aNwsfHh5KSEo4ePcqWLVvIz8+3nFqu0Wjw8vLCaDRSUVGBp6cno0aNIiwsrN7rdg0GAz/88AOZmZlMnjzZanKuS2VmZrJ9+3auv/56fHx86jyfnJyMVqslMTHRUkNKSgphYWH1zsaevDQFgMSJA21mxvLm0Ol0+Pn54efnV+c5g8HAuXPnOHv2LGfPnqWgoIDs7Gx2795tuV699owHX19f/P39LY9Lr+u+WOLEgRxIOczvP29l0sPX4+XlRUREBFu2bCE+Ph5XV1cSEhLYunUrCQkJlssLvLy8uO+++/jkk0/o2rUr7u7ubN68mQ0bNhAcHExcXByRkZHcd999lgH5tGnT6NevH7169WL58uUsWrSI6Ohorr/++jqXUwhxJTIYF+IiGRkZvPDCCwC888473HPPPR3cItFSKpWKOXPmcObMGX788Ucefvhh0tPTO3zCGSGEuBqUlJSwb98+9u7dS0FBAfb29tjb2+Pt7U1JSQlr1qwBat6rfXx88PPzo1evXpbBXH5+vuWe0OPGjSMhIaHByVSrq6tZuHAhJ0+e5K677rrsNb0VFRX8/PPPBAcHk5CQUOf5s2fPsnPnTkaNGmUZBGZmZnLmzBnGjRtXd/1TBRxIOYxKpWL4zXWPmosaWq223oG6yWTi3LlzFBQUWAbqp0+fJjU11TJId3d3p0uXLvj5+dGlSxf8/f3x8vJCpVIx6Pr+fPPPReRk5pGRmknv+J4MHTqUL7/8ksOHDxMZGcnQoUPZtWsXycnJjB8/3rJvDw8PbrjhBhYvXswtt9zChAkTLKexL126lKSkJKKiohgzZgzr1q1j7ty5TJs2DTc3N6ZMmcKBAwdYuXIl//nPfxg9ejT9+/e/qj+MEU3TrMH4Rx99xFtvvUVubi7R0dHMmTOHYcOG1bvuTz/9xMcff0xqaipVVVVER0fz8ssvW72RzZ07l+nTp9d5bUVFxWU/BauPzDCoPBqNhri4uE6RTVhYmOWN++mnn+7o5rS5zpRNS2g0GubNm4dKpeLFF1/sNAPxqyWfzkiyUS7JpuNVVlayY8cO9u/fz9mzZwEsg5OqqiqcnJyws7MjOjoaPz8//P398fHxsbrmsri4mJUrV3L48GF69uzJxIkTL3u3nKqqKr7//nvy8vK499576dGjR4Prms1mfv75Z/R6Pbfccku9A6fVq1fj4eHBoEGDLMt+//13unbtSnBwcJ31Ny6pOSoePSQC3+7el/8GKVRH9h21Wm25lvziU76NRiMFBQWcOXPG8tizZw9lZWXAH0fh/f396TWgBwd/z2D9fzfRO74nAQEBBAcHs3nzZiIiInBycuKaa65h/fr1DBgwwOrshj59+pCRkcHy5csJDAy0nKp+/vx59u7dy969e0lNTcXV1ZWysjK+/PJLpk+fjoeHBzExMfTs2ZO1a9eyfPlyUlNTmThxIl27dm2174+8rylbS3Jp8gRuixYtYurUqXz00UcMHTqUTz/9lC+++IJDhw7Ve++9p556im7dujFy5Eg8PDz4+uuvefvtt9m2bRvx8fFAzWB85syZpKenW722S5cujW6XzKYuWqK6ulrx99IUra+ysrLJH/gJIYT4g16vJz8/n5ycHI4ePcqpU6e4cOGC5XkXFxeCgoIICAigS5cudOnS5bITXhmNRrZu3crGjRuxt7fn+uuvJyoq6rJHGisqKpg/fz4FBQXce++9BARc/nrtbdu2sXLlSqZMmVLvnDBHjx5l/vz53HnnnURGRgKQlZXF3LlzmTx5cp3rg01GEzPHvEhBzjn+39szGHpD3SPtonWVl5eTl5dnGaDn5eWRdfAUWWsLUGtVJD4SS2CPALRaLTt27OCee+6hV69eGAwGPvroI7y8vLjnnnusfq6qqqr45JNPcHFx4f7777caYJnNZrKzs0lNTeXgwYPo9Xo0Gg3Dhw9n0KBB2NvbA5Cdnc3y5cs5e/YsCQkJjBw50vKcEPVp8mB80KBB9OvXj48//tiyLDIykptvvpnXX3+9UduIjo5m8uTJvPTSS0DNYPypp/5/e+cdH1WV9vHf1NRJ7z2QkARCQsdQBQSpIihgWRXbqyu6lnVXXbe4a2Gbq+66FtRFUVRQwAYKEYEgAWlJiKmk90ySmWQmk8m0e98/hnO4M5lJA5IJOd/P5+beOXPvzTnz3HKec57yGNrb2wdSFRtYNHXXxWw2IysrC/PmzXPJCJDvvfceXn31VXz//fcIDQ0d7uoMKa4umytJVlYWbrvtNnz++ee45pprhrs6DhnN8nF1mGxcFyabK4fJZEJTUxMaGhrQ2NiIxsZGKJVKm31kMhliYmIwadIkjBs3rsdAd2/y6a8iw/M8zGYzjEYj2tra8MUXX6CrqwtLliyBn58fzGazzWKxWOh2e3s7cnJyEBoaitjYWHAcB57nbdbFxcWQSqV0dp3jONTU1MBisWDMmDHUTF4kEkEkEqG5rA1Z/zsNmYcUN/7+OkhlEohEIkgkEkilUkgkEpvFvox8lkqlkMlkkMvlkMvldFsmkw2J2fNIv3eMRiN+vexPaK1XYca6NHhGSdHc3EyDa4WEhCAiIgISiQRnzpyxGWwh1NXVYevWrZgxY4ZDdwTAeh+cPXsWBw8epEo5yTkfFxcHnufpgJK7uzuWLl2KlJSUS5LhSJfN1c6QRVM3Go04c+YMnn76aZvyJUuWIDs7u1/n4DgOWq22h8Lc2dmJ2NhYWCwWTJo0Cc8//zydOXeEwWCAwWCgn0kaBaPRCJPJBMBq8iKRSGCxWKi/ibDcbDbbRL+TSCQQi8VOy8l5CeRHF6bS6K1cJpOB4zhYLBZaJhKJIJVKnZY7q/tIapPJZIJWq4XJZHKpNkkkEvz1r3/Fs88+CwDYsmULnn766VElJ57nodVqYTQa6f8d6W3qr5w2b96M+vp6LFq0CDt27LB56bpKm8i9QzqIV9O1N9LvJ4vFYnPvXA1tulrk1J/7ZqS1iTCUcuI4Dm1tbVT5rqurQ0tLCziOg1gsho+PDz1eIpEgKSkJ6enpGDt2LP0NiDxIm3ieR2dnJ9ra2lBfX08VZLVajfz8fDQ0NMDLywsJCQno7OzE559/DoPBQPt2RqORLo74+uuvHZYLFd/u7m4qi4qKCqpQi8ViiMVi6HQ6GAwG+Pv7Q6fTQSwWw2g0oqurCyEhITSwHJEvz/MoPlYBAAhNCoBK3QYAVB72CxkcEMqqP9gr6lKpFHK5nPrge3h4wN3d3abMzc0NXl5e8PDwgEQigZubG53tdXSNkXuH47ge1+pIeEbIZDIsWj8PO175Am0lGjz8xydhNpuRk5ODzMxMKBQKmwGknTt3Ijg4GBEREQgPD0dERARCQkKwePFi7N+/HxEREdQKQtgmAJgyZQpSUlLw6aeform5GRUVFTh37hwUCgVNzTdhwgR8++23+OyzzzB27FgsWbIEQUFBg3pGGI1Gm3eOKzwjBisnV37uDbZNQxZNvbW1FRaLpcfsYWhoKJqamvp1jpdffhk6nQ7r16+nZcnJyXj//fcxceJEaDQavPbaa5g9ezby8vKQmJjo8DybN2/Gn//85x7lhw4douknYmJiMHnyZJw7dw41NTV0n6SkJCQnJ+PkyZPUlwkAJk2ahNjYWGRlZUGr1dLyjIwMhISE4MCBAzYXyIIFC+Dh4YF9+/bZ1GH58uXQ6/U4dOgQLZNKpVixYgVaW1tx/PhxWq5QKLBw4ULU1tYiNzeXlgcHB2PWrFk4f/68jfn+SG7T2bNnMXv2bJdoE8dxOHToEP7zn/8AAG666SZMnDgRBw4cGFVyIrEeSACdq6FN/ZXTxo0bodFokJ2djTVr1mDTpk1YsGCBy7WJcLVdeyP9fmprs3a4yb1zNbTpapIT4Wpq05WWU3R0NDw8PHDy5Em0traiq6uLdlqDg4MhEokQGBhIlYL29nZERkbC29sb7u7u6OrqwtGjR1FdXQ2TyYSysjKbWWmJRAKDwUDPWVhY2ENepPPb1NSE4OBgmM1mdHZ2QiwWQyqVws/PDwkJCVCpVKitrUVLSwukUikmTZqEKVOmoLS0FE1NTVTBTkpKwoQJE3D8+HEolUpUVlbCYDBg9erVmDhxIn744QcbOaWlpWHHjh0ICgqycZVsa2tDUFAQwsPDbWY3ly9fjvqqRhz4+wkAQOq14xAUHoAVK1ZAqVQ6lFN1dTVyc3Np5z0gIADTp09HcXExzp8/D47jwHEcgoKCEBsbi9LSUjoIwnEc/Pz84OPjg7q6OnR0dFDlggwS6nQ6G2XD0W8skUjg5+cHLy8vtLe328zMS6VSFBYWoqioiH6Wy+Uj5hkx+4YZ+OzfX6L8XBU+fPtjBEb7YdmyZcjPz0djYyMSExMRGhoKrVaLqqoqiEQilJaWIi8vj/4+ERERUCgU2LNnD8aMGQNvb2+EhIQ4bNOUKVNQXl6OkpIShISEUDeL7OxsREVFQaFQICYmBjU1NXjrrbeQmpqKlStX4scffxzQM4K8a8h6tDzLR0qb5s6dO2jL7AGZqTc0NCAyMhLZ2dk2+RVffPFFfPjhhyguLu71+E8++QT33Xcfvvzyy15TC3EcR3MC/vvf/3a4j6OZ8ejoaDQ2NiIw0Bo4g43suEabTCYTMjMzsWTJEnh4eAx7m7q6unDvvfdi165dAIB//vOf+NWvfjWgNtmXD3eb+irvbWZ83759WLx4MQ1qNtLbNBA5WSwW3H///fjoo48AAH/5y1/w1FNP0RmU4W4TuXeWL19Of/u+2nQ1yskV22QwGPDdd9/Re+dqaNPVIqf+3DcjrU2EyyUns9mMuro61NXVUZNz0uH08vJCYGAgfHx86Du7qakJLS0tMJlMdAbZvm7kf3p6esLT0xMeHh50cXd3h6enJ9zd3SEWi1FYWIi4uDjk5ORArVZj0qRJmD9/Pk031VebCgsLsWfPHoSHh2PdunXw9PTsU05ZWVk4evQo1q9fj6SkJIdy+uKLL1BVVYX/+7//oxM7lZWV+OSTT7Bhw4Ye0dmlUik+feULfLXlO6RMT8Qz/3v0ssppsNeexWJBd3c3DAYDtSggM/7d3d3Q6/UwGAzo6uqCTqeDXq+HXq9HV1eXwxl7qVQKT09PeHl52Sze3t7w8fGBl5cX3N3d4e3tDblcTgdPhut+ev037yH7m1OYvWoGHnjpTkilUpSVleHjjz/G+vXrkZCQAMCaYz4nJwf3338/vLy80NTUhMbGRjQ1NaG2tpbmRHd3d0dkZCSio6MRHh6O0NBQen2IxWKIRCLs3bsXZ8+exfz58zFjxgyUl5cjPz8fZWVlEIvFSEhIgFQqRUlJCdzd3bFw4UJMmDCBDu701Sa9Xo/MzEz6zhkNz/KR1Cae5wcde2pAyrjRaISnpyc+++wzrFmzhpY/+uijyM3NxZEjR5weu2PHDtx999347LPPsGLFij7/1/3334+6ujp8++23/aob8RlXq9Xw8/Pr1zGMoYHjOLS2tlLTnOGkoaEBq1atwtmzZyGVSvH++++P6vRlriSb4YLjOPzmN7/Bv/71LwDArbfeivfff98lAvox+bguTDauC5ONLWazGfX19aisrERdXR2USiVVvMViMeRyOUQiESwWi1Pzb8AaiI3kfibKGFHQhAp4X36xLS0t2Lt3L6qrqxEbG4ulS5f2O2Avz/M4duwYDh48iPHjx2PNmjW0A98bJSUl+PTTT7FgwQLMmzfP4T6lpaX45JNPsHbtWkycOJH+v3feeQdSqRR33313j7YZu43YdO3T6GzX4Yn/PIjpi527V44EyMCLh4cHuru70dXVRZV0nU4HnU6Hzs5OdHZ20m175V0qlcLb25sq62Tx8fGBQqGg6/5cK4Ol7Fwl/rD+r5DKpPjPoc3wC/IBz/N4//330d3djQceeABisRgGgwFvvvkmAgMD8Ytf/KJHfaqqqrB9+3b4+vrCz88PDQ0N0Ov1AAB/f39ERUUhIiICkZGRCAsLQ3Z2Ng4fPozJkydjxYoVkEgk0Gq1yM/PR15eHpRKJb1P2traEBkZiaVLl/YZcBBgzzVXh7jvDIYBmanL5XJMnToVmZmZNsp4ZmYmVq9e7fS4Tz75BPfccw8++eSTfiniPM8jNzeXPgwHArtAXQ+xWNwjn+Rw4e7uDo1Gg6CgIOzatcvpS3m04EqyGS7EYjFefvllJCYm4pFHHgHP8y6T+ozJx3VhsnFdRpNseJ5Hd3c3NBoNNBoNOjo6oFQqoVQq0d7eDp1O53AGWyKR2ChIRHnq6upCU1MT6urqYLFYEBcXh7S0NIwfP/6SI0J3d3cjKysLP/30ExQKBdatWzegoFYWiwXffPMNcnNzMXfuXCxYsKBfxyqVSuzZswdJSUlO0/AajUbs3bsXY8eORWpqKi0vKChAY2MjNm7c6PB/Ze89hc52HYIiAjFlQVq/2uHKSCQSREZG9nt/nudhMBiogm6vqHd2dqK+vh5arRY6nc7mWKlUCoVCYaOgC7fJmvi4D4SEtHgkpMejLK8SB3dk4aZNKyESibBkyRK8++67yM3NxZQpU+Dm5oaVK1di+/btyMvLw6RJk2zOExcXh5tvvhmffvopxo8fj9tvvx1qtRr19fV0KSwshMVigVgsRmhoKOLi4pCXl4fm5mbcdtttUCgUmDVrFjIyMtDU1ITc3Fzk5+cDAJqbm/Hee+9hwoQJuP7666FQKJy2aTQ910Yil6J/Djq12VtvvYWMjAxs2bIF77zzDgoKChAbG4tnnnkG9fX12LZtGwCrIn7nnXfitddew9q1a+l5PDw84OvrCwD485//jGuuuQaJiYnQaDT497//jQ8//BDHjh3DjBn9Sw9BZsZbW1upmTrDNTCZTDhw4ACWLFkyLEqOMDAZYB0hd3Nzc5gndLQx3LJxNY4fP4709HRqfmaxWAbVEbhcMPm4Lkw2rsvVIhue56HX66mi7Wjp6OhwqGwDVtNOHx8f6uscFRUFPz8/KBQKavnD8zzq6upw7tw5FBQUQK/XIyQkBGlpaUhNTaX9NGEwPOFiX0Y+8zxPF7PZjPPnz+Pnn3+GyWRCcnIyurq6MHnyZEgkEpt9yQJcfGeLRCIYjUacOHECKpUK06dPR2xsLPULJwvZVxiQzWAw4JtvvoFcLsfatWvh7u5O3ZCIeapYLKbmyg899BDNZW4ymfDf//4XoaGhuPXWW3v8vhzH4bc3/AX1ZY249cm1uOE+x5G3RxJX8t6xWCzo7OyERqOBVqula+G2RqPpcT17eXnBx8cHvr6+dC1cvL29HQ6UHPv6JF7/zXvwCVDg3wdfgpuH9ZrfvXs3Kisr8cgjj9D7YM+ePSgtLcWmTZvg7e3d81zHjuH777+3sZoQtqu5udlGQW9tbaXfR0ZGIj4+HpGRkYiMjIRCoYDFYkFpaSlyc3Nx/vx5GgQ0MTERK1ascJii+Wp5rl2tmEymQctlwDHYN2zYgLa2NvzlL39BY2MjUlNTsW/fPpr+obGx0cYx/u2334bZbMamTZuwadMmWn7XXXfh/fffBwC0t7fj//7v/9DU1ARfX19MnjwZWVlZ/VbEGa6Ps87ClUalUmHjxo1YtGgRHn3U6svlKKfoaGa4ZOOKCGNh8DyPtWvXIiUlBc8///ywvfyYfFwXJhvXZSTIhud5dHV1ob293Wbp6Oig2/Z+iZ6enpDL5eA4DgaDgbZTJpMhICAAgYGB8Pf3pzOKJpMJZrMZJpMJ5eXl9HNnZyfUajVV5sViMdzc3KBQKKDX65GdnY2jR49SBbu3gGADpaCgAIDVF3sw/PTTT/jpp58GfNy7777b5z6vv/46VdJJRHGRSIS33nqL+umSqOZNxW2oL2uEzF0KSYgFhw4dot8JA585W4YqXdlAuVL3jkQioQq0M4SWHvZKukajQWVlZY/7QiKROFTSg8f5ISDcH6pGNQ599iOW3rkQALBw4UK8/vrr+PHHH7FwobXs+uuvR1lZGfbu3Yv169f3kMusWbPQ0tKCL7/8Ej4+PlTnIf8/IiICERERmD59OgCrFUhpaSkyMzPR2NiItrY2/PjjjwAAHx8fqphnZGRg+fLlKCoqwvHjx1FaWorS0lKEhYVh9uzZSE5OtnHDGAnPNcbAGfDMuKvCZsZdF5PJhH379mH58uVDqtDs378f9957L+rr66FQKFBVVcVy0NsxXLIZCRw8eJAGmpwxYwY++OADmuJkqGDycV2YbFwXV5ENmdlub2+HWq2GSqWCSqWiSrBWq7XxtyWpp0gALABUeTYajQNOhUUCGRHlUCwW08BeJBCbj48PAgIC4OvrS9Nlkcjawtlj4eKozL68rq4O2dnZUCqViI+Px5w5c2hEdrPZjIMHD2Lx4sXUX93RLDfP88jJycH+/fsRHByMNWvWwMfHx2b23NGMOpmN/+abb1BfX48bbrgB/v7+NrP4ZFuv1+PAgQPw9vbGNddcQ8t1Oh2ys7MRERGB6OhoKgeymEwmZL1zFtpmHWJmhiFqWpDNd/ZBn5whTFXmTGEXpilzd3d3uHZzc7ssbpqucu/0BlHYOzo6nC4kJoK6TIfmUx2QeUow6740+Af6w8fHB62traiqqsLNN9+MuLg4eHh4oLCwEJ999hlWrlyJqVOn9vi/ZrMZH3/8MRoaGnD33Xf3yCzlCKPRiD179qC4uBjTpk1DbGwsGhsbUV9fj4aGBjrYExwcjMjISHh6etLo+YDVlH/ChAlIS0tDZGQkvvvuO5eWzWhmSGfGGQxXR6vV4te//jXeeecdAEBiYiJ27tzJFHHGgFi0aBE+++wz3HfffTh58iQmT56MF198EY8++uiwmq4zGIyhh+d5mEwmGnW6u7sb3d3d0Gq1UKvV0Gg01FeWRLHur0IGXMyDzXEcurq6YDKZ6LFubm7w8/ODn58fAgICEBAQQJUwR8obUcDFYjE0Gg0KCwtRWFiIuro6SCQSJCYmYuLEiRg3bly/gp8NhPr6ehw4cADl5eWIiorCxo0bbWYRAWunVSqVwsPDw2nn1Ww2Y9++fcjJycGUKVOwbNmyfteV53l8+eWXqK6uxq233uo0RS7P89i5cydEIhHuuOMOahrM8zx27NgBLy8v3HHHHQ6DeeYcycfe5qNw83TD71/9DRT+3j3OzXGcTV50+zzpzhZh9HOj0UgjoAvTwjlCmF/ckbLu7u4Od3f3HlHuyXqkxFwSiUS0/s4C/1ksFmg0GrS1tOFvd72BTrUO+gYzFL5mVFRUoL29HRzHYefOnQAu3mO+vr7Yu3cv1Go1oqKi4O/vDz8/P7i5uUEqlWLDhg344IMP8NFHH+Gee+6hLg3OkMvlWL9+PbKzs3Hw4EG0tLTg5ptvhre3NziOQ0tLi415u1KppDnESSpAEvyNBEqsq6tDXFycS1pVMAbHVTcz3t7e3qsJDGPo4XkeWq0WCoXiij48yAv40Ucfpa4Sv/rVr7B582bqA8ywZahkM5Kpq6vDfffdh/379wMArrnmGrz55ps9Ar1cCZh8XBcmG9elN9nYK9UkWrSjpauriy59KUIEYprs5uZGI40TE1oSQVosFkOr1UKlUqGlpQWNjY3QaDQArP6xJEIzWbu7uw+o/VqtFoWFhSgoKEBtbS0kEgkSEhIwYcIEjBs37pIDsTmirq4OR44cQVlZGYKCgrBo0SIkJSU5vDf6undaWlqwe/dutLa2Yvny5Zg8uf8Rykm6ztOnTzv07xVy5swZfPPNN1i/fj1SUlJoeUFBAT7//HOsW7cO48eP73Ecx3H4w/q/ouLnaqy8ZzFu/+3N/a7fpUBm/IXKubM1WezL9Xq90+uYKOtyuRze3t42yrojBZ4srj5Lu++D7/Hh5s8QFBGIl7/9M+RuMvA8j4KCAuzatQszZsyAj48P1Go11Go1qqure1ihuLu7w8/PD/7+/vDw8EBxcTHEYjHWrl2LyMjIfmVfqa6uxueffw4AuOmmmxzGLTIajXTmvL6+HtXV1TaB78RiMTiOg6enJ1JSUjBx4kRER0ePmIGUqxni9z8YmDLOuOKQF4hUKr2indaGhgaMGTMGBoMBcXFx+N///ocFCxZcsf93NTBUshnp8DyPd999F0888QQ6OzuRnp6OnJycK/6bMfm4Lkw2ww9RrEnKJfv0S0T5sF+cmXsTE23AqnCZzWYbxUUYgZzMUgcGBlJFm+S2tVgsNibLKpUKSqUSra2taG1tpbmLJRIJnY0jwalkMplNsDTiry3c5jiuh5m20WikJrpdXV0ArIq9QqGAl5cXxGJxn+bdgG3QNGdlwrVer4dKpUJXVxfkcjnNT05yL5Oc5PafSfuJeTspa2pqwvnz5+Hh4YH09HT4+fnZnIMcIzSlJ5/FYjFOnDiBnJwcLFmyhAaII/9DSENDA7Zu3Yq0tDSsWrWKlnd1deGNN95ATEwM1q9f7/A6yd57Cv/59btw93TDKwdegF9Qz2BbroojCw/7+8M+HznZz2AwODynRCKhqe3IAJSHh4dNmf1yuS0yesOgN+KJpX+Aqrkdt//2Zqy8ZzH97vPPP0dFRQUeeughGrhNqVTinXfewfjx4zF9+nSo1eoecR3IzDrB09OTWq84WsiARWdnJ3bt2oWqqirMnj0bCxYs6NPSTqvVIicnB2fOnIFGo4FIJLK5Z6VSKSIiIpCamopJkya5/ODI1QpTxsF8xl2ZK+mD1N3dbTNr8NJLL0Gn0+F3v/sdvLy8Luv/uhoZCf5hrkR9fT1+/etf46GHHqJp8QwGA+0YXm6YfFwXJpvLT2/KNdm2L3MU0MjNzQ0cx1E/aHJuomALzYWFHWqRSET9poVKIun8EqXYbDbT9eUKbCZUNvvaJnUym81UWSL5wT08PGheZ4lE4tAf25F/tlAGwrWjMo7joNfr0dbWBp1OR5VwLy8v+luRxX4AgfyGGo0Gnp6eNmVkoIS0VXj8pSISiWyUd71eD7FYTBUlYt7f2toKnU6H5ORkeHh42Pjcy2QygBdh29O70NGixbW3ZmDR7XNsvrdfExmMFHp7rpF4A46UeKEVif3iSH5yubyHgu5IgSem2cSiZLAc3nUMbz+7DV4+nnj1wAvw9rP2D3U6Hd58801ERUVhw4YNVFbEasKZ/zjP86ipqcHOnTshFouRlpZG40OQIIzCdnt5eVHF3NfXF21tbSgtLUVgYCDWrl2L8PDwPttgNBqxc+dOGAwG1NXVwc/PD15eXlCr1XQADrAODERHRyM5ORlRUVEIDAwcUdfgSIX5jDNGHQaDAW+99RZefPFF7Nq1i+YO/d3vfjfMNWNczURGRuLTTz+1Kfvb3/6GHTt24MUXX8Tq1avZS4/BuIAj5ZpERyY5iEmHnfhgO+q42yujBKLokFljApnBE6YXcnROojARs3JhADPhzKvws7Cc+IyTFGPt7e101t3X1xfBwcEICQlBWFgYQkNDaVotewXbkVLs7Pesr69HcXExiouL0dbWBplMhoSEBCQlJWHcuHHw8PDot3wGCsdxKCkpQXZ2Nurq6hAcHIylS5di/PjxA1KU7BW+oqIi7Nu3D3K5HDfccIPDQJk8z1M5k8EQsphMJmRlZaGoqAjTp0/HuHHjeqRgI8eYTCacPn0aZrMZaWlpkEql1IqhpaUFWq0WwcHBNNq8MCCbyWRCU74aHS1aSD3EaDBXYvv26l7bKhKJbIK0Odq2X/dnmyj6Q4lEIoGXl9eAJjqI1UZvyjrJKNDQ0EDdRhzNE7q7u9P/T5R0R9tkZl54Tc67MQP7PjiI2tJ67HlrH+54eh0Aq5K8cuVK7NixA7m5udQlYurUqWhqasK+ffsQGBjYw6RcJBIhNjYW9957L95//32UlJTgzjvvpHEHOI6DVqulyjkJ2tje3o76+np0dHSA53m0trZiy5YtkMvlCA4Opj7qwsXX15daYCkUCqxfv566hdTU1CAwMBDXXnstRCIRjQ1RUlKCkpISKrfw8HBERkYiPDwcERERCAwMZKbtLgRTxhkjCpPJhO3bt+NPf/oT9Qt/++23qTLOYAwlZrMZW7duRVVVFdasWYOZM2fixRdfxMKFC5lSzrgqICm0iLIszBNMlGnyncFgoAGoiPIzGMhMJVFAiNIhnJk2Go09/odQWdDpdBg3bhwCAgLg5+dHzaeFkcoH8hu0tbWhoaEBDQ0NaGxsRFNTE02vFBAQgIiICEyZMgUREREIDw+/bH7ZJpMJVVVVKC0tRUlJCbRaLTw8PJCUlITFixdjzJgxV9wyw2QyIS8vD8ePH4dKpUJsbCwNjHYpzzmdTofvv/8eBQUFGDdunNP8ysDFyPAAbPxzOY7Dl19+iaKiItxwww29+pfzPI+9e/dCq9Vi48aNiI6Opt+p1Wq89dZbmDhxItauXevweFVzO55c8ScAwB2/3YA5q2fYKOr2a7LYB28TlpNgffbB2/pjDSCRSGyirQsDuAkHmHr7nmxfqTRrIpGI/p++gp0ROI5Dd3e3Q6VdaB2jVqvpZ0fPGnsFPXF+NGpL6/HdhwcRnR6K+PGx8PLyQmxsLNLT0/Htt98iOjoaQUFBAIClS5eira0NO3fuxH333ecwCHBAQAA2btyIbdu24b333sMvfvELBAcHQywWU/cT+wCGpI0ajQatra3Izs5GZWUl1Go1eJ5HbW0tNBqNzYCEQqGAr68vdDodjhw5gsDAQMyfPx96vR75+fnYt28fvL29MXPmTKxbtw5qtRoFBQUoLi6GSqVCXV0dmpqabFIhhoWFISwsjD6zSL0ZQw8zU2dccS6HOadOp8N7772Hf/3rX6iuto5ER0RE4LnnnsPGjRuZmeggYaa2l057ezv++c9/4pVXXqGmYtOnT8dTTz2FG2+88ZJmL5h8XJeRIBtnwZ6EyrWwg0vKhAr1QKOB2882k2BPpGPs7e1Ng0TJZDKq7JN6CGeTOjo6bPxUpVIpfH194e/vD19fXzpzRGaTPD09IRKJLkk2JpMJSqUSzc3NaGpqoou94k1mmMLCwgYcYK0v1Go1zp8/j7KyMlRWVsJsNsPPzw9JSUlISUkZsoBNWq0WZ86cwalTp6DX65GSkoJZs2YhMjLyks5rNBqxfft2tLS0QCQSYenSpUhNTR2wMmgwGKjP75o1a5Camtrr/kePHsUPP/yAVatWYcqUKbScDKp2dXXhwQcfdDqQ8sqv3sbJA2eRkB6PP3/8W4glV04GFovFYRR2R9tkEEy4drTd271MXDTkcjmMRiMCAwN7ROy3V+TtI7WT9VArdDzPw2Aw9HBjceTucu6L82iv1sEjSIaY64J6xEEQi8WIjo6GQqGAp6cn5HI5zp49C7FYjFWrVlF3DJKSj6DRaLB9+3ZoNBrceuutiImJGVAbqqur8fXXX6O9vR1z585FRkYGdDqdjY+6SqVCZWUlJBIJDfhIIC4iOp2OtiE9PR2xsbHgeR7l5eUoLS1FZWUlLBYLDdBnMBjouaRSKUJDQxEeHk6XkJAQlj2mn1yKmfpVp4yzAG6ux+UIdDRnzhwcO3YMABASEoInn3wSDz/88BU1yRsNsCBUl4+mpia89NJLeOedd9Dd3Q0AePbZZ/HCCy8M+pxMPq7LlZQNMe/uTZF2FDlZr9fTdX9n1xwhkUhsOtzEf5MELyMBzEhwMGedNZ7nodPpqHmmUMkm20TJBayzNULTTHtTTdLh7M/v1x/Z6HQ6qmwT5bu1tZUG4gkMDOwxe3S5FW/AqgzW1NTg/PnzOH/+PNra2iAWixETE4PExEQkJiYiKChoSJ4BPM+juroap06dQnFxMSQSCdLT05GRkXFZ0oM2Nzdj3759qKmpwYQJE7Bs2bJBxXfp6OjAxx9/jI6ODqxfvx5jxozpdf+zZ8/i66+/xvz583HttdfScp7n8dVXX+Hnn3/G3XffjYiICIfHn/khD/986A2IJWK8tPtZxCZFDbjOw4kwEntvSjt5rgiVfkf79fZsIbP1jpT13tKuCddXKsBbW6MKT654Dt1dBmx48gZMXpxKFXWlUomzZ8/SWWjhjLw9Qiscssjlcpw/fx4ajQYZGRlISUmhg5D9ibZuNpuRlZWFY8eOITAwEEuXLrW5roXPNY7jHD5XSZBIEkOCoFAoEBAQAIVCAYvFAq1Wi5aWFhgMBri5uSEsLAze3t6wWCxoa2ujz0GxWIzg4GCEhobSJSwsjMVkcgAL4AamjLsyA00BxPM8Dh8+jGnTpkGhUAAA3nvvPWzevBlPPvkk7rrrLqaEXyZYeqbLj1KpxH/+8x+89dZbOH78OBISEgAARUVF4HneYaocZzD5uC7OZENmafpKO9RXCqLeXs1kFpr4GguDkzk6Ti6Xw8PDg3YaSYAvR36Xnp6e/Z4J4TgOnZ2dtFMo7BySbWGANblc3mNWW7jt4eFxWa5ze9lwHAeVStVD8e7s7KT1EnY0w8LCEBIScsUsHnieR0tLCyorK1FRUYHKykqYTCYoFAokJCQgMTERY8aMuSIpyJxhMBiQl5eH06dPo6WlBUFBQZg2bRrS09MvywCEXq/HoUOHcPr0aQQGBmLevHmDmg0HgNraWuzYsQMymQy33XYbgoODe90/Ly8PX375JaZOnYrly5fb/M+ffvoJ3333HW688Uakp6c7PF6r7sTTNz4PVXM7Vt23BLc9edOA6zxS6O87x2w29zpA2Nfzrbu722HwRYJEIum3Au/oe3d3d6eDcd9t+wEfvLQDHl7u+OuXf0BIVBD9Lj8/H7t378aiRYswZ84cANbnXFVVFXbu3AkfHx/MmjXL4Ux8V1cXzeJgj1Qq7eHfbr8m2zqdDgcOHEBtbS2Sk5OxZMkS+Pv7D6g/0N3djZycHJw9exatra00jzqZVSfPPuHvTUz9PT09ERISgoCAAHh4eECn06GlpQVKpZIOnnp7e9so6KGhoQgKChrVs+hMGQczU3dl+msyWFNTg48//pgGw3jjjTfwy1/+kp6DBLxhXD5GgqntSMVoNNqMhq9btw6ff/455s6dizvvvBM333wz/Pz8ej0Hk8/QYrFY+q1I6/V61NXVwcfHh84YERPv3rAPFCZUqEmQKuIT7Wx2WxjIqK+O3UCUayE8z0Ov19PgZB0dHXSbrLVarU39SC5eZzPb7u7uV3RQied5dHR0oKGhAT/++CMCAgLQ1taGlpYW2tH08fGhAdWI4u3v73/FB7va29tRWVlJl87OTkgkEkRFRSEhIQEJCQkIDQ0d0kE3nudRVVWFvLw8FBYWwmw2Izk5GdOnT0dcXNxlqYvZbMbp06dx9OhRmM1mzJ8/H1OmTMH+/fsH/FzjeR4//fQTMjMzERkZiQ0bNvQ5Q0dmxCdPnoxVq1bZtKm4uBg7d+7EzJkzcf311zv9n/96+E2cPpiHiPhQvLT793Dz6HuWc6QylO+c/jxv+1Lqe3veisVih0q7TCbHkXdPoaVShdAxQbjtTzfC08uD7nfu3DmcPn0aa9aswcSJE+k1U19fj23btiEiIgK33nqr09lus9mMo0eP4ujRowgPD8eUKVNoIEuh77vQPcgeYp1ErBBCQkIQHR2NhoYGTJkyBT4+PjbP+t58/hsbG3Hy5En8/PPPsFgsSExMxMSJExESEkLdgtRqNVQqFZqbm6HRaHoMlEilUmoxIJFIqF8/iR1Cfm/hLHpISAiCg4Ph4+MzKiYTWDR1xohFpVLhs88+w/bt23H06FFa7u3tbeMTwxQRxkjDPsgQybFLXtKbNm3CypUrcfvtt2P58uVXxPx1tCA0wezLnLu37b5maoQdO+IzSDonpLMhVKiJqafQP5v8PyH2aXycKdiXolzbYzQae1W0NRqNjQm5WCyGj48PDUoUExNj89nPz2/IZnJ5nkdnZyeUSiVdWlpa0NLSQjvnZPA2IiICkyZNQkhICEJDQ+Hp6TkkddRqtaipqaEz32q1GoA11kl6ejri4+MRExMzLO82lUqFvLw85OXloaOjAwEBAZg7dy7S09OdBlAbKBzH4dy5czh8+DA0Gg3S09OxcOFCKBQKm+uqv3R3d+Orr75CUVERMjIysGjRoj7vg1OnTmHfvn2YNm1ajxnxmpoa7Nq1CykpKVi8eLHTc2R+cgSnD+ZBKpPi4Zfvu6oV8aFGmJ98sHAcR5+x/VXqOzraEX9tGFT17WiuaMW2v32K4PSe1/2ePXuwZ88euLu70+e+v78/qqur8eqrryIpKYnGwLCflZ8wYQKCgoKwb98+ZGVlYc2aNT0ishMsFovT1I1arRa1tbX0+QYAe/fu7XEOqVTqMLo8WScnJyM1NZVGWt+1axc8PDwwceJETJo0qUf6NpPJhOrqapSUlKC2thatra003zoxk7cPmieXy2mAz59//pkO1MpkMgQFBSE0NJRmmBhNSnp/YMo4Y9jo7u5GbGwsHVUTiUSYP38+br/9dmzYsIGaqDMYIx2xWIydO3eirq4OH374IbZv346CggLs3r0bu3fvxuLFi3HgwIHhruawIDTr7m0mpC8FeyA+jKRz5evrSz+TgGIE+8jdwly6JLCOXq9HYWGhzf8Si8U2CjUJLOZMub7U/LmOfs/u7m4a8Zyk3iLbRNHW6/U2xxF/cF9fXyQkJFAlm5R5e3sPeceJ4zio1WrqB0l8GVtaWuhsklQqpR288ePHIzg4GAEBATh69ChWrFgxJMouz/Noa2tDTU0NXYjyHRQUhISEBMTHxyMuLm7YXKy0Wi0KCwtRUFCA2tpauLm5YcKECZg0aRKioqIum2w5jkNBQQGOHj2KlpYWpKSk4Pbbb+/TlLw3Kioq8OWXX8JgMGDDhg0OU58J4XkeP/zwA3788Udcc801WLJkiU376urq8PHHHyMqKgpr1qxxev8Vnz6PbZt3AgBufXIN4scPLCgX48ojFovp83ygnEg/g9ce24K2wk7ccu96TJyTTN8nXV1dOHz4MFpbWzFhwgR4eHjQd46bmxvq6uqQn58PT09PGt/DmaFxd3c3PvjgA7i5ucHHxwceHh59mtkL301ubm4wGo3IysrC2bNn4enpiRkzZmDMmDE9UkeStVqtRn19PXQ6ndOZd4vFgtOnT+PkyZNwd3dHUFAQzUtO3mHTp0/HvHnz4O7uDqVSiaqqKlRVVaG6uhoWi4Uq2r6+vpDJZNSUvr29nUaEN5lMNAOF8DcSi8VQKBTw9/dHSEgIIiIiEB0dDT8/v1EX1Z0p44wrjsViQUVFBZ5//nlUV1fjgw8+AGA1aZw1axaUSiVuu+023HrrrYiKGllBUa4GrlSgFEZPoqKi8Mwzz+Dpp59Gfn4+tm/fjo8//hiLFi2i+6hUKtxyyy1YtWoVFi5c6LKuGfYBgXpToon5tqN9+mPWbd9p8fb2RmBgoFNfQWKVQGap9Xq9jTJNPqvVavrZWXocMntDlGdi5uzu7o6qqipMmzbNxmTwSppim81mhwq2/Wf7GX5Sb4VCgaioKBsl29fXFwqFYlifA3q9niraQqVbpVL1mF0hym1ISAhCQkIcdtwuxVywP1gsFjQ2NlLFu7a2Fl1dXRCJRAgLC0NiYiJiYmIQExMzrIPKRAEvLCxETU0NxGIxxo4di7Vr1yI5Ofmy/kZmsxm5ubnIzs6GWq3G2LFjsXr1aqeR1/tzvRmNRmRmZuL06dOIj4/HDTfc0Kdrj9FoxJ49e1BcXIzrrrsOs2bN6jEjvn37doSFheGWW25xWo/mmha8/PCbsJgsmHn9FCy7c5HD/a5GRkuf4JqlU1Fyx0J89+EPeOfZD/HHD5/E2Ilx9PuxY8fi448/Rn5+PjZs2GATTE2tVuOjjz6C0WjErbfeiujoaBt3JeFar9ejpKQE5eXl6OzshK+vL0QiETQaTY/9e/McJoHtOI7D4cOHcfToUQQGBiI8PBze3t40EJsj33mS7727u7uHj3tLSwtNf1ZXV+f0/8vlcqqkR0dHQyQSUd95pVIJi8UCsViMoKAgjBkzBrGxsQgODqbB5siiUqnQ3t5OZ9E7OjpQVVVl87+kUil91/v5+SE4OBhhYWF0Rn0o42kMBVedz3hHR8dlM7NiDA6e51FcXIyDBw/i4MGDOHz4MNrb2+n3dXV19AXd2dkJb2/vYaopgzH8kNROZMZs+/bt+MUvfkG/j46OxsKFC7Fo0SIsWrTIabTfgUDS5vRnFtp+P+F2b7PR5GVKlGh7pdpZmfCzXC6nyjSZrSCKs71ybb929GqTyWR0Npp0Kjw8PHp8Fm7bp7C5EphMJnR2dvZY7JVt+9lsmUwGhUIBhUJBlW37z97e3i7Rudbr9VCpVNQ3kazb2tqg0+nofr6+vggKCkJgYKDNeriCGHIch9bWVjQ0NKC+vh4NDQ1obm6GxWKBVCpFVFQUVbyjoqKGtZPI8zxaW1tRWlqK0tJSGwV8/PjxSEpKuuwz852dnTh79ixOnTqFzs5OjB8/HrNnz76k5xTP8ygsLMSBAweg1+uxePFiTJs2rU/5q1Qq7Ny5E2q1GmvXrkVSUpLN98Q8NzIyslef3/aWDvzljpfRWNWMMamx+OOHTzLz9KsUi9mCfz70BnKzfoZfsA/++OGTCI8Lpd8bjUbs3LkTVVVVWLt2rU0AVp1ORy3errvuOlxzzTW9XqMNDQ34+uuv0dzcjBkzZmDBggU2z4vesmjYrzUaDZqamtDR0QHgolucyWTq9b1sn5ZO+L6VyWTQ6XTU+shisUChUFBrIw8PDzoAb/8O7urqcjqQQNLmeXl50dSU3t7eNpZhFosFGo2GBgDVaDTQ6XRO+xnEIsLb29sm3SUJTkrefVc6Tsnl4qpTxtVqdZ8jp4zLi16vh1wupzN4TzzxBF555RWbfXx8fLBw4UKsXr0aN910EzNBdxFIRzMoKGjUmQW5KrW1tfjss8/w9ddfIzs7m84cE1Pr9957D3PmzKE+vySKtn0Kmt7yzzqaASaIRCKnEWudKc2O9hHO6FssFhr0jKwdbTv6zpkfN0m51R/Fmqwvp1La171jsVig0+kcKtn25fbWAWKxGF5eXg6Va+G2m5uby3Q0iB+3SqWiyjZRuFUqlY2ppIeHBwICAuDv70+V7aCgIAQEBPQrBVBfDPa5xvM81Go1GhsbqeLd2NhI5RMUFISIiAhEREQgKioKYWFhw265YjabUVVVhfPnz6O0tJT6dI4ZMwYpKSlXRAHneR719fU4efIkCgoKIBaLMXHiRMyaNQtBQUF9Ht+bfJRKJb799ltUVVUhKSkJ119/Pfz9/fusT05ODr777jt4e3tjw4YNCA0Ntfn+5MmT2L9/P5KSkrB27VqnVgEdbRo8f9e/UF/WiMBwfzy/4xn4h4yeDD2jsU/Q1anHc7f9A7Wl9QgI9cMfP3wSoTEX3SosFgv27NmDgoICLFq0CLNnz6bPXYvFgh9++AHZ2dlISUnBqlWrer3fOI7DiRMncPjwYbi5uWHBggWYNGlSv35rR7LR6XQ4efIkTp48ie7ubowdOxaTJ09GVFQUTCbToKLc99VHAGCTd54o8jKZzCbLh9lspqbzJO0mUTmFMVYc4ebmBg+Pi0H1RCIRDfjX3xR75P8I36VEWXe0XKqlEMdxg75nrjplnEVTv7IYDAYUFhYiNzcXOTk5OHHiBHJycnDs2DHMmDEDAPD+++/jwQcfxJw5c7Bo0SLMnz8fTU1NWLVqFQvE5mKwaN1XBmLCTZRi4SJUmO3X9t+TaKUAelW6yGOczPwSkzUyuyx8cZJtuVzuUKF2lA6GmLjZv7SFS29KtTNTdIlEQl+4RJEm/n/Cz2QtVLiHuqNIfLGJiZ9Go8HJkycxZswYOmtPvuvs7HSYm5bIRphazFGn4HKl+LqckDRmwrRlxOyQzGYIB05IXlt/f3+bdUBAwBUPVtif55rBYIBSqaRpzsg2uVZ9fX0RGRlJle+IiAiXMI3keR5NTU2orKyk/psmkwm+vr5ITEzEuHHjEBcXd0We51qtFvn5+cjLy4NSqYS/vz+mTZuGyZMnD0jhdyQflUqFI0eO4Ny5cwgICMDSpUuRmJjY57l0Oh327t2LoqIiTJo0CUuXLrWRk9FoxDfffIP8/Hxcc801WLx4sdNnR0t9G/72wH9QX9boUCkbDYzWPkF7qwYv3PUv1Jc3IiDMH09veQTR4y66WPA8j0OHDuHo0aNITU3FypUrba6z4uJifPHFF5BKpVi+fHmfKUzb29vxww8/ID8/H8HBwVi0aBHGjRvX63O/N9mYTCYUFBTg9OnTqK+vh4+PD9LT05GWltavATJ7SOo6Ei+lrq4OVVVVqK+vh1arhUgkogPDHh4ekEqlNGCpfR+nL8WeIBKJIJVKIZPJaLYRsVhMU3eSGC6kb9XbeYCeir5YLLYZRLHHzc3NYcpPR58d9UEuxT2KKeMMh5DZNjJLceDAAfzmN7+h6U/sEaYhI75z5OU8Wh/uI4HRLhtiFkYUYeG2o8/2Zb3t3xckOrczRZkEFKuoqEBaWhoNFOPv70/3f+SRR/D555/3iE7s5eWFtLQ07N69G6GhoTCZTOjo6LDxGxMq1P0xi3P2qiDmYo6Uakdr4bazPLBDAcdx0Ov1PQLf2KefEa4d/Qb2QdmcjbpfrijoVwqDwUDN4R0p28QKg0DSmAmjqRNl28/Pb1ifJ8Lnmkgksgn8plQq0dTURAOsER9HYb7c8PDwPlNmDRUcx6GlpQXV1dVU+dbr9ZBKpYiJiUF8fDzGjRuH4ODgK3IvdXd3o7S0FPn5+SgvL4dYLEZycjLS09ORkJAwqP8plE9HRweys7ORl5cHT09PzJs3D1OmTOnzXuE4DqdOncKhQ4cgFouxatUqpKSk2OxTX1+PL774Ah0dHVi1ahUmTpzo9Hzl+VX4xy//i45WDfxD/PCHbU/YmCuPFkZzn6C9pQPP3/UvNFQ0wcPbHY+9+gDS5tgq1QUFBfjqq6/g5eWFNWvWIDo6mn6n1Wqxb98+FBcXIzk5GUuXLoWvb+9WFQ0NDThw4ACqq6sRFhaGOXPmICUlxeGAUX9l09DQgDNnzqCgoAAGgwERERFIS0vD+PHjL4tVqkqlQllZGR0Q7O7uhkwmo8+j2NhYhIWFUUs0RzPaZJvEbmlra6PpMvV6PX3XiEQi+izgOK7PmXDA+kwnv58wbSjHcb365BPIM42snf1PMntP3v1Tp07t4RrTX5gyPorheR4qlQqVlZUoLy9HcXExiouLUVRUhNLSUrz++uu45557AABZWVmYP38+AMDf3x+TJ09Geno6ZsyYgYyMDMTExDh9KY/mh7ur4+qy4TgOJpOpz4WkkXK09KZc9yfFDvF3kslkVEkWKst9fSf8LFS6e+tskpl1nU6HAwcOYNasWT1eaESpJia1bW1t0Gq19MXo7u6OiIiIPk25yL6enp5OA6IJI5Dbl/WW33QoEAZoG8hin14M6BkJXbhtHwVdLpfj0KFDQxaxe7CQGW0S4M1Z4Df7ASQSOEeobAvXrjBLLKSrq4sGgFMqlSgqKoJYLLaJV+Lp6WmjdIeFhSEoKMglfOoJer0e9fX1qK2tpQGVjEYjxGIxoqKiEB8fj/j4eERGRl6xemu1WpSUlKCoqAhVVVXgOA7R0dFIS0ujkaUvBeKHKxKJUFZWBm9vb2RkZGD69Ol93ks8z6OsrAzff/89lEolpk6dioULF9qkxzIajThy5AiOHz+OsLAwrF271unsIM/z+P7TLHz4189gMpgQPS4Sv33rYQRFBFxSG0cqrt4nuNJo1Z341yNvofj0eYglYty0aSVufGAZxJKLyrFKpcKePXtQX1/fw/eb53kUFRVh37596O7uxowZMzB37txe7xme51FdXY2jR4+ioqICAQEBmD59OiZNmmRjRTRQ2ZjNZpSWluLcuXM4f/48OI5DZGQkkpOTkZycPKgZc3s4jqOWOpWVlaiurobZbKZpJaOiohAdHY3o6Oh+x4fieR4dHR1obm6maSvJmvRlPD094efnBx8fHzrYTfojFoul176fMB6O0WiE2Wzul5LeF4mJibjtttsGdexVp4y3tbUhIGB0PkTtIcFcGhoaUF1djcrKSqo8A8CRI0dw7bXXOj3+N7/5Df7+978DsJqC/fDDD4NKg2I2m5GVlYV58+a5VKeHMTjZEEWRKMHCpS/FuD/Ks3Dpr3kTAKrwChepVNqngtyXEi3MIe3s9xA+/J2tnZmkO1v39Wi2V+7JWiaTUbPplJQUqkQ//fTT+Pnnn21mw41GI3ieh0wmg16vpwMEn376KTo6OhAfH4/o6GhERERcsZygwtl6ZzP2jszhyeLMCkEul9OZ+N4WoYI9kGAvw/lc43keer2e+p8T03j7tVarhU6ns7mWJBJJn0HffHx8XO5ZTSLyEl90+4X4pItEIvj5+UEkEtEZ4+DgYJqux5XQ6XRoamqiS2NjI9ra2gBY/eqjo6NpZzYiIuKy+NM7wmw2o66uDhUVFaioqEB9fT1EIhHi4uKo3/nlCIzb0dGB3Nxc5Obmor29HcHBwZg1axZSU1P7vN54nsf58+dx5MgRNDQ0IDo6GkuXLrUJFMdxHHJycnD48GHo9XrMnz8fs2bNcjrw2dakxv/+/DHOHjoHAJg8fyIefvleeHoPT/o5V4D11wCT0YR3/7gdWV8cBwAkT0vAfX+5A5Fjwug+HMfh+PHjOHLkCNzd3TF//nxMmjSJXmsGgwHHjx9HdnY2JBIJZs6cienTp/dpcVNXV4cTJ06gqKgIEokEqampSEtLQ2xsLCwWy6Blo9frUVpaipKSEpSVlcFkMiEwMBCJiYl0NvtyDLBaLBY0NTXRwcTa2lpoNBoA1thR4eHhCAsLo+uB9CssFgtUKhWUSiWUSiUNANrW1kb7AeT5HxgYiMDAQPj7+8PPz48uztpIgtoKlfbu7m5qHabVaum7Veh+J7QUTktLw5o1awb1u111yvhoiKZuMpmoyV1LSwuioqKoaURJSQnuuece1NfX2wSeITzzzDN46aWXAADV1dWIi4tDeHg44uPj6WhZcnIyUlJSEBcXN2ofxq4OUQAdKcLOFOT+7Nuf4/pjJiSE+AD1tfR3P0dLXwoz+a1Im3pTmh2t+9qnP49RiURiozRf6nqgirFer0dZWRmKioqoFUxJSQkqKysREBCAsrIyum9GRgZOnDhhc7yXlxciIiIwZswYfPfdd/Qa/PHHH2EwGOgItZubW4/AMY5M34XbziAm8GQW3t7MXRigzX5xZZNwe4jJvKPI8ETZtvdHt78PpVIpNZMnM/lC5Zoo3J6eni7njw6ABiQkZvJku729HWq1Gu3t7TbBf0gEXeFCov662nvLaDTamMo3NzejsbGRxoOQy+V0pj4yMhJRUVEICAi4YnIymUxoaGhAbW0tzRlsNpvh6emJ+Ph4JCQkXLbAbxqNBiUlJSguLkZFRQVkMhnNb96bRR2hu7sbeXl5OH36NFpbWxETE4P58+cjPj7exv+zoKAAP/74I1paWjBx4kQsXLjQaTBfg96I/R/9gN1v7oOhywCpTIpbn1yDpXcsHDVByxh9c/SrE/jfcx+ju8sAiUyCFRuvww33L4WXz8VBvY6ODmRmZqKgoAD+/v6YN28eJk6cSN8/nZ2dOHr0KHJycsBxHNLT0zF9+nSEhYU5+7cArNYpZ86cQW5uLtVrUlNTkZqairCwsEt6NphMJlRWVqK4uBjl5eXQaDQQiUSIjIxEXFwc4uLiEBkZedlie3R0dKC2tpbmGm9sbKSZQTw9PREWFobQ0FAEBwfTYJ4DefaQ4KFEMW9ra6OKulqttpnU8fDwoIo5eYeQbYVCMeB4LcTyTKvV0mjug+GqU8ZHUjR1juOg1WptAuCQHHzjx4/H9OnTAQBlZWW488470dLSgtbWVhuzO8BWwa6qqkJ8fLzN9yEhIYiOjkZ8fDzWrFlDzSjsUypdyXbW1tYiOjp6RL7ohIEjHK0v5bvBnm+gCjFgVQZJUAyi+EokEnAcBy8vL5vAGRKJxCaIhvCY/pSRGVpnwcDITLL9bPpA1s6+c1TW38ccUZjtzc+Fa7L0tZ/9ejAKwuW6d4S/uSPzfaIUkFn5H374Aa2trTSaOTHVJz5SQUFBMBqNfV6HZrPZJh94e3s7OI6Dh4cHjXBKTJ9JIBihSfxw+pT3hb1s7NPSCHPMChVsR9uOTOYBq08aUayFSrbwM1kPRQq2wUDcLewjyRNlmyjf9oMyZDDBz8+vh9Lt6+vb62DLcLxzzGYzfZerVCqbfOkk/RBgNf8PCwuzWa6k4m2xWNDW1oampiZq9t7c3AyO46if55gxYzBmzBiEhoZecj1IDvaKigqUlJSgoaGBzrJPnDgR48ePh0wm61U+ZrMZ5eXlNEc6x3FITk7G9OnTERsbS+uo1WqRk5ND06qNHTsWCxcudJpWTafpwqHPj+Gb/x1AR6t1ti5pyljc86fbEZPkOB/6aGOk99cuN8q6Vrz//KfIOZIPAPBUeGDpHQux+LZr4Rd0ceKvubkZhw8fRnFxMby9vTF58mRMnTqV+ozr9XqcPn0aJ0+eRGdnJ0JDQ5Geno4JEyb0OoHI8zxqa2uRn5+PwsJCdHV1QaFQYOzYsUhISMDYsWMvSWkmGSSImXllZSUNPhoUFITIyEi6BAcHXxbXBZ7naWo2oqArlUoaxwOwPv+FyrlwhnsgFkJEUSc6lv3S0dFho6xLpVIbazFHVmPkfWsPi6aOofMZJ/lgSYciICAAISEhAIDW1lZ8/fXXPToc5AW9YcMGqgifO3cOkyZNcqokPPXUU/jrX/8KwLGCLRaLaUqYO++8E08//TQAa2f666+/plFgw8LCrphpW38xGo3Yu3cvrr/+ekgkElgsFqrg2m/39Xkw313q8Zdyi4jFYqrcEmVYuHZURiJI9nffvpRkEszC0aBCd3c3jh49ihkzZtBUFMIZd/v9+xo46O8xA/0NhTPnva37s09vCvNQdz7sI4TaWyHo9XqcOHEC6enp1H9eOFs/kKU/1zH5LeyjsEskEqpoAkBycjL97tVXX8X58+fR0tIClUplY2ofFxeHiooKev4ZM2bg1KlTDv93UFAQWlpa6OdnnnkG5eXldMZdmJbE19fXJhc78WsdTDRy4nYhdB+wj34vVK7JZ71eD6VSCTc3tz7TrJCo8fbp1pxtC/OvuhLktxIOJpDFfvaevAcdWT54enrCx8eHzkb4+vrCx8fHpuxSZrcvt98rz/M0wJ39ADqZuScDWoB15j4gIICaxwcFBdHtKxVJ3mKx0FRyxG++ubmZ5gsGgMDAQERFRdElJCTkkq8xvV5PFf3q6mrU1NTAZDJBLpfTGfbExESbQX9H8lGr1aioqEBlZSXKyspgMBgQHByM1NRUTJ48mQae6uzsRElJCX7++WdUVVVBKpUiLS0NM2fOpH0xIWaTBYUnS3B832kc33cKBr3VWjA4MhA3P7IKc1f3nht6tDHafcYdwfM8zhw6h0//tQf1ZY0AAIlMghmLJ2PWyhlIm5UCubu1n61UKnHq1CmcO3cOJpMJsbGxSElJQXJyMnx8fGCxWFBeXo68vDyUlJTAYrEgPDwc48aNQ0JCAsLDw50ONHZ3d+Pzzz9HUFAQKioq0NLSApFIhNDQUOqXHR0dDV9f30Ff0zzPo62tDXV1daivr0d9fT0dvBOJRPD390dISAhdLqc1kslksrEgIjFA2trabPqNXl5edJCWKOhCZXmg7//Ozk5qjUXiqNjHWbGPLySTyXoMkqempvbQ1frLoJTxN954A//4xz/Q2NiICRMm4NVXX8XcuXOd7n/kyBE88cQTKCgoQEREBH7729/iwQcftNln165d+MMf/oDy8nKMHTsWL7744oBs74kyfvDgQbi5ufXwM5w2bRo15a6oqMDWrVt77EOWu+++G2vXrgUAnD59GmvXrqWdC3uBvPDCC3j22WcBALm5uZg8eTLNsUfC6JP1o48+it///vc0WMPcuXMhEong5uZm0xnx8fHBokWLsHr1ajp7feLECfodMTUELkYXJAtRIi/3cinnvdzjPWKxmCqnwm37z/35ztk+9rIjZtD25fbbwiiMJHrj5Rh4GOw5Bqr4OvqtnQ0AOCvrbdBAWNZfBftyKiT2v42zxdl+fZn292cf4TJQC4fBmvE78qe3V7wvx2CEvQuNxWLBkiVL6PePPPIIiouLbZSY9vZ2WCwWxMXFobKyku5LFHeRSNTjmggICMDBgwfpoMUf//hHFBUV2bRNGHBt1apVVLmurq6GwWCg16rwvnUGuQ/sB3O6u7uRlJREn8fCFzTxQXeFIHcEYYq6vhZ7/3yigDt6ppDsGfaz9Y5SuHl6el7xQYb+KBQcx9HYCvYLed8LA93ZZxFRKBQ2vojCzqGPj89ld5Mwm820Yyi0KCDmmCRbAmC9DkNCQmyC1YWGhl7SQAAxtRd2lhsbG+lsllwuR0xMDDVzDQ8Pdyrnzs5OfPnll4iPj4dSqUR1dTXa29shEokQERGBhIQETJgwAcHBwdBqtTTuTXl5OZRKJUQiEeLj45GamkrjYhB4nkdjZTNKcspRdLIUZw+dg05zMc1gVGIElt+1CHNXZ0AqGzmuLEMFU8adw3Ecftp/Fnu3ZqL8XBUtd/N0w6S5EzAhIxlJkxMQlRgOk8maaqywsBCVlZXgOA6hoaGIjY1FXFwcoqKiIJVKUVZWhtLSUpw/fx4GgwFSqRSRkZFUsQ4JCaHKtb1sOjo6UF5ejtraWtTW1tJ4E+7u7lRZDg0NRUhICI2ZMdjsBySYGgmk1tzcDJ1OR/fx8fHpYcFE9BVvb+9LmhTkeR5ardamzyDcFloeAdbBbzKATxb7AKxk6c+gN3lvEuVcGKtFGKdl9uzZmDBhwqDaOGBlfMeOHbjjjjvwxhtvYPbs2Xj77bfx7rvvorCwEDExMT32r6ysRGpqKu6//3488MADOHbsGB566CF88sknuOmmmwAAx48fx9y5c/H8889jzZo12LNnD/74xz/ixx9/xMyZM/tVL6KML1++nCaIFypOaWlpiIuLo6M+Z86c6aFYkXVERASCgoLoy7qurs7hfiQ1kVQqHVDY/CsJqZ+jRfidI0XTkWLprIz8r97KyGeO41BXV2dj9mSfOmAgo1hXYrDBfrmS9DUQMNCBhoGeR6gQA0B2djauvfZaet+QY4isrvRAT3/PTUyt+6tQO1ouB/bm/n0tA92fHENeXkeOHMGyZcuuaO7p3u4r8ps7c68Y7GdhOZmVBkDdCxz5RfcHYglAzmmxWDBnzhw6079//35UVFQ4nAU3mUw4cuQIHaC4//77sWvXLqf12LlzJ2688UbIZDLceeed+PDDD+l3YrHYxrqgqKiIWm29+OKL2L17t1OLjVdffRW+vr4wm8346quvcPLkSepWYm9Bk5GRAZlMBpPJhPr6ejpbIpRtf99PwmeF0IeOzAyTgIHEXYH47sfFxdHARCQHu6NBT7FYTC0tBoswgKRQ1sJBhK6uLvz888+IjY3t8T1RuInPoj3u7u42HTn7Th1ZLmUmiHTuSA5f4UCA8LNer6eDAcKOL2B1YfDx8UFgYCBNKRcQEIDAwEAoFIoBvVNJPYTuA8JI+2q1mgZhAqyWDcHBwQgPD6dLYGAgfWfwPE8HNUiAPZVKRU34VSoVAOuzNCQkhFryeXl5QavV0tmwpqYmanFATHPHjBlj9RfnxFA1qdHWrEZLXRvqyhpQX96ImpJ6aNWdNm30CVRg+qJJmHPDTCRNHVw6ttECU8b7R2VBDbK+OI5T3+egrVFt852nwgMxSZGISohA5Nhw+IX6QNutgbqzDQ3Keqo8Cv2lAwMDwXEcNBoNWltbUVtbS+95uVxOLWuUSiWmTZuGwMBAOolHnqddXV2oq6ujZt9KpRKtra30uS+Xy20GD4mPs3CwdCAWWTqdDi0tLQ4DahJTd4Kbm5vNM5Xk6na2kP5of7BYLD0GT+0XEnjN0TuQZJMh7zThe06YPUa4bR/sVyKRwGw2D/qeGbAyPnPmTEyZMgVvvvkmLUtJScGNN96IzZs399j/qaeewldffYWioiJa9uCDDyIvLw/Hj1sjFW7YsAEajQbffvst3Wfp0qXw9/fHJ5980q96EWX8gdseglzuZqscwvpjyy9E0eM4jubCJt9f2NlmFoZgNpkgEiSLtx7j+CLh0fPn7PETX/hMioXH8DwP8BdKehzGX9yX5y8eLywfKEM4diC8scR2v6f92l6hty+3L6PbEAEiYbkIIpHtd7bHkWMF/8duP57nnQ4aCK+DHg8OkeNyq4wBjuesMuR5cDwH8BeUInoN8OA5nnam6WdYZc9z3IXrxPZ4m+MExzj7n8AQXgY8+a3FEIutv79YLOrx+WI5WcQ2x9kPLNmXiyCCyMF5RSIxxNYTQSzqe5+L14VVhmIxyVl58V7lLsgHDmTktMxezlQmtmUcx6Gjvd3qUyYSUZlzHAee46z/+4KSZbMt+M5aPw48d3Efzu6Yy4qDAUGJzSCgBGIxGeS78L1EbKOwWfcRHC+R9DynRExlR8qsfttmGIwGGA1GOgAQFhYK7sI9UV1VhQ5Nh8BCgQwMmMFxPGbMmEHvuXPn8tDS0nLxt7tw/5Dv09PTqb9fUWEhWtvabNpzcZEgOiqKyk6n08FsMkFsN/A1EMxmMyxmM7wVCsjlckglErS3t0OpbIGFs8BitrZJaIlx/dKlCA4KgkQqxZkzZ5CVlWUzUGIVn/VZ+ctf/hIxMTHgwePYj8ewf/9+eg9abwnrfiIRsH79BsTExoDjOOTn5+No1lGIBHIRC66JGTNnIjQkFBzHoaGhgaYgsxk0lEoglUhpvAGLxQyj0Zpdob+KvMHQDYvFsdUQmZ2SyWTQdHTg+4MHre0nAeIEMpk6dQomTEgFx3FQq9X44YcfHA5ii0VihEdYIwTzHIcuvR7V1dW214FIZD23SAyHXYgLz2hSVzIA6e/nj4TEBKv1ikiMzO8zLzxLL7zHRKDnDAsLQ9K4JJgtZphNZuTk5Ni+O8UiiEViiMQip78lz9u+ezieg5enF8YmjLX2jXggNy8XHMdduAYunpt8dnROaxN5ev3KpDKIIAJn4cFbLvxPDuAsHHgLwFt4cGYenImDxcRDzInh6eYFjboThi7HcRYAQCQRwTvEAwtWzsPkaydi3KQx2PbhNtTX19vud6Ge7u7ueOKJJ2j5J598gurqaof7SiQSPPnkk7R89+7dNsEvhfsCwK9//Wt6LX3zzTcoKSlxuu+mTZtotOfMzEz8/PPPTve9//776QDYkSNHkJeX5/T3uPPOO2kspePHj+PMmTNOz3vTTTehoqICM2bMQF5eHk6ePOn0vGvWrKGByM6dO0f79I7Ou2LFCkRGWv3yi4qKcOzYMafnXbx4MWJjYwFYYyZlZWU5Pe/8+fMxZswYANagxIcOHXK676xZs5CYmAjAmov74MGDTvedNm0akpOTAVhNzzMzM3vsA1iv61BFBNprOlFythznc8th6HKc4QMApDIJfAIVkLlLAQnQbdJDJAMg5iGWXLiHJCJIZBLI3WTw8HSHSCKGxWK25t/mOXrPWyt0If2qTAZPby+4uVmVQ2JpAsCmT2jhLJBKrPpNd3e37STFhXOS54NYLIab3A3+/n6QXVA+mxqbbO5zEel3iUTw9lYgITGBuqOeyzsHk9lkM9lCtsHzEInFVtc3J2qo6MIzUyKRwMvbm+pmWq0WPM9deBcJ+4mAXO6G2NhY2meorq6GyWS82Djw4DhrH0ksFiMsLIwO6ra2tMB4IUgxz3N0v/5MCohEwPyF87H85uv73Nfh8QNRxo1GIzw9PfHZZ5/ZmJA/+uijyM3NxZEjR3ocM2/ePEyePBmvvfYaLduzZw/Wr1+Prq4uGkDk8ccfx+OPP073eeWVV/Dqq6/2eCASyCg3QaPRIDo6GjclbIRMMrw+0gwGg8FgMBiMK4cZJmi7Negyd0JjbIfGoEaHUY327jaER4ajsrKSDvJcc801+Omnnxyex9/fH83NzQCsyvaSJUt6KGoEuVyOzs6LM+833XQTvv76a6d17OrqopM7d9xxB3bs2OF0X5VKRaMxP/DAA9i6davTfRsaGhAcHAye5/H444/jv//9r9N9S0pKMG7cOJhMJjz77LP4xz/+4XTf06dPY8qUKQCAv/zlL3juueec7puVlYWMjAxIpVK8/PLLNoMU9nz77bdYtGgRxGIx3n333R6uqkJ27dqFVatWAbAOjNx1111O9922bRtuueUWAFbdYsOGDU733bJlCzZu3AgAOHDgAFauXOl039deew0PP/wwLBYLjh07hgULFjjdl7irWiwWnDp5CssWrISv3B8+bv7wlfvBU6aAp9QL7tLRmzJvNOAb5YW3vv/XoI4dkJ1Va2srLBYLQkNDbcpDQ0PR1NTk8JimpiaH+5vNZrS2tiI8PNzpPs7OCQCbN2/Gn//85x7lcl8p5C6W1mTQiESOR4wuR/kwmWmJ4HgWdiDlF8e3hq0ZI49h+p1EEPWw2iDWBIMtF1qyWL9wcm33p1x4AZFyUna5L9ZBVpGUO7KWIZYbPaoicvS7W/9ebnkI/6e17pe3XGQtvPgYs5Of83IRtVigxRcsIBz9T5EI4DjbcqvFhONyah3RRzmxrnBezttVnbTJcTl3wbLGvo48bysRYt1hP6pPTYkdlPP25RdmfXk7CwoyG0ysBWzLJXYm8bx1puVCXay1FNnV0XKxTXSWV3Rh5pWchad1t1hs6y6RWK1hSFuJZYl1f97uNxBBJpNSy5OL35O22przi8USyGRSWCwcOO6iuwux7OC4i0E+SbBJsVgCk8k2By2xDLFYLOABai1CzBq7hGbzPE9nrS0WC41ZAIggl8vQ1dVl/Q1E1n1FIhEkUgn4C9epTC6j/9NoNIKzcLBwFvJDUgUVIkAmvWhW2d3dDbFEArPZdMGax3qAVCKlsQ/Ib9ne0U5dJbq7DdSiBzyg8FEgMiISKrUKIpEIzU3NsFgs8PD0AM/x0NOgfjzkbnJMmzYVZt4MtVqFyqpKGMwGiCViKHy8YTQb0anXghdz8PLxwLpb1iFj7jVo7VBi8982A6pOeAII8Y5CYGA6Wltb0dnZCS8vL+zbtw9JSUlITk7GpEmTbKJWBwQEwNvbGw0NDRCLxdi3bx8Aa1rHpUuXggQ9JZDAWg0NDXRfwDozq1Ao0NjYaCPrqKgodHd349tvv6WyDA0NxR133AGdTgeVSkWvG3d3dwQHB6O6uppOQHl4eGDZsmXw9/eHWq3u4Z/r7u6OkydPoqWlBWKxGHPnzoWfnx+8vLzQ1NRkc+0RF6ADBw7AbDZj1qxZAKyZdiQSiU3dz549C3d3d0RGRqKrqwsZGRkArPd2WFgYDAYD9U3Oz8+HyWTCwoUL4evri2uuuYaeRy6XIyAggPrYFhcXw2AwICYmBrGxsZg/fz51FeF5nppJq1QqlJeX099YJpNh5cqVaG1ttWkTyRldX19P962oqMDixYshkUigVCohhJh5k33LysowZcoUBAcHw2g02kTylkqliImJQW1tLXJzc1FdXY309HTIZDL4+flBp9PZmGGT59K5c+eQk5tDI/N7ekrg5sVB1VGDJqMRIl6E2TPm4tZ1t8HHyxeZ332Pr/Z8BTGkkPASuLu5QyKWwqA3ABAhJDAEISGh8Pfzh9lkxulTpy/8RxGdDQdE4CwcPD09ofC2Bjl0c3NHY0OD1frRurf12SwRAzwgkUrhcSHOgkgkgk6no+8POuEusPqREX9vnke3wXDhsSN821jfTRKJFFKphH5nNJl69DtEFw4TS6zPPvIst8pWWAPh2WGtwwUrHYvZIugPiuhzinyWXHgu8bzVvY527cgfwbtGIpbQ/oLFbLZtl+ji24e8H8hvY+/2JRIBQZGDDx4+oJnxhoYGREZGIjs7m96ggNX/7cMPP0RxcXGPY8aNG4e7774bzzzzDC07duwY5syZg8bGRhrt+4MPPsCtt95K99m+fTvuvfdepzlonc2MNzY2Ur888sKhZhEXIOVmux+emDU4K7cP3kZGO+0DuzgrJy8yoVmISCSCVCp1Wu6s7iOpTSaTCZmZmViyZAk8PDyuijbZl4/UNvE8j3379mHx4sW0UzjS23Q1yYncO8uXL6f1HOltEtZ9JMvJYDDgu+++o/fO1dCmq0VO/blvRlqbCFeDnPR6PTIzM+m9czW06WqRE7l3li1b1sMke6S2qbfykdQm+/vmamjT1SQnnucHHahuQFPIQUFBkEgkPWaslUplj5ltQlhYmMP9pVIpVZqd7ePsnACoM7095MEuhPih2eMsAIuzcmeO+QMpd+Yf6KzcWd1HYptIHa6mNhFGapvIw2Qg942rt6m3ctYm1iZndRxoOam3/b0z0tt0NcqJtck12yS8d66WNgkZyW0i8ZPsGcltclY+0tpk/865GtrUn3JXb5O9cj4QBhQ1Ri6XY+rUqTSQASEzM5OavdiTkZHRY/8DBw5g2rRptCHO9nF2TgaDwWAwGAwGg8FgMEYyA3aufuKJJ3DHHXdg2rRpyMjIwJYtW1BTU0ODMTzzzDOor6/Htm3bAFgjp7/++ut44okncP/99+P48eN47733bKKkP/roo5g3bx7+9re/YfXq1fjyyy/x/fff48cffxxwgxyNjjCGF7FYjJiYGCYbF4TJxrVh8nFdmGxcFyYb14bJx3VhsnFdmGxcm0uRy4BTmwHAG2+8gb///e9obGxEamoqXnnlFcybNw8AsHHjRlRVVeHw4cN0/yNHjuDxxx9HQUEBIiIi8NRTT/WIpPj555/j97//PSoqKjB27Fi8+OKLWLt2bb/rRFKbdXR02ATpYDAYDAaDwWAwGAwGw9UYlDLuihBlXKVSwd/ff7irwxBgsVhw7tw5pKWlOfTTYAwfTDauDZOP68Jk47ow2bg2TD6uC5ON68Jk49pYLJZBy+Wqs3XoT3J2xtDCcRxqamqYbFwQJhvXhsnHdWGycV2YbFwbJh/XhcnGdWGycW0uRS5XnTLOYDAYDAaDwWAwGAyGqzPgAG6uCrG212q1TkPlM4YHk8mErq4uaDQaJhsXg8nGtWHycV2YbFwXJhvXhsnHdWGycV2YbFwbk8kEmUwGhUIBkUg0oGOvGmW8ra0NABAfHz/MNWEwGAwGg8FgMBgMxmhCqVQiODh4QMdcNcp4QEAAAKCmpga+vr7DXBuGEI1Gg+joaNTW1rJI9y4Gk41rw+TjujDZuC5MNq4Nk4/rwmTjujDZuDZEPnK5fMDHXjXKOMnv5uvryy5SF8XHx4fJxkVhsnFtmHxcFyYb14XJxrVh8nFdmGxcFyYb12agJuoAC+DGYDAYDAaDwWAwGAzGkMOUcQaDwWAwGAwGg8FgMIaYq0YZd3Nzw5/+9Ce4ubkNd1UYdjDZuC5MNq4Nk4/rwmTjujDZuDZMPq4Lk43rwmTj2lyKfEQ8yQnGYDAYDAaDwWAwGAwGY0i4ambGGQwGg8FgMBgMBoPBGCkwZZzBYDAYDAaDwWAwGIwhhinjDAaDwWAwGAwGg8FgDDFMGWcwGAwGg8FgMBgMBmOIYco4g8FgMBgMBoPBYDAYQ8xVp4xXVVXh3nvvRXx8PDw8PDB27Fj86U9/gtFoHO6qMQC8+OKLmDVrFjw9PeHn5zfc1Rn1vPHGG4iPj4e7uzumTp2Ko0ePDneVGACysrKwatUqREREQCQS4YsvvhjuKjEusHnzZkyfPh0KhQIhISG48cYbUVJSMtzVYgB48803kZaWBh8fH/j4+CAjIwPffvvtcFeL4YDNmzdDJBLhscceG+6qMAA899xzEIlENktYWNhwV4txgfr6evziF79AYGAgPD09MWnSJJw5c2a4qzXqiYuL63HfiEQibNq0aUDnueqU8eLiYnAch7fffhsFBQV45ZVX8NZbb+F3v/vdcFeNAcBoNGLdunX45S9/OdxVGfXs2LEDjz32GJ599lnk5ORg7ty5WLZsGWpqaoa7aqMenU6H9PR0vP7668NdFYYdR44cwaZNm3DixAlkZmbCbDZjyZIl0Ol0w121UU9UVBT++te/4vTp0zh9+jQWLlyI1atXo6CgYLirxhBw6tQpbNmyBWlpacNdFYaACRMmoLGxkS75+fnDXSUGALVajdmzZ0Mmk+Hbb79FYWEhXn75ZTah5QKcOnXK5p7JzMwEAKxbt25A5xkVecb/8Y9/4M0330RFRcVwV4Vxgffffx+PPfYY2tvbh7sqo5aZM2diypQpePPNN2lZSkoKbrzxRmzevHkYa8YQIhKJsGfPHtx4443DXRWGA1paWhASEoIjR45g3rx5w10dhh0BAQH4xz/+gXvvvXe4q8IA0NnZiSlTpuCNN97ACy+8gEmTJuHVV18d7mqNep577jl88cUXyM3NHe6qMOx4+umncezYMWa5OAJ47LHH8M033+D8+fMQiUT9Pu6qmxl3REdHBwICAoa7GgyGy2A0GnHmzBksWbLEpnzJkiXIzs4eploxGCOPjo4OAGDvGBfDYrHg008/hU6nQ0ZGxnBXh3GBTZs2YcWKFbjuuuuGuyoMO86fP4+IiAjEx8fjlltuYRNYLsJXX32FadOmYd26dQgJCcHkyZPxzjvvDHe1GHYYjUZ89NFHuOeeewakiAOjQBkvLy/Hf/7zHzz44IPDXRUGw2VobW2FxWJBaGioTXloaCiampqGqVYMxsiC53k88cQTmDNnDlJTU4e7OgwA+fn58Pb2hpubGx588EHs2bMH48ePH+5qMQB8+umnOHv2LLO8ckFmzpyJbdu2Yf/+/XjnnXfQ1NSEWbNmoa2tbbirNuqpqKjAm2++icTEROzfvx8PPvggfvWrX2Hbtm3DXTWGgC+++ALt7e3YuHHjgI8dMcq4o+AS9svp06dtjmloaMDSpUuxbt063HfffcNU86ufwciG4RrYj97xPD/gET0GY7Ty8MMP49y5c/jkk0+GuyqMCyQlJSE3NxcnTpzAL3/5S9x1110oLCwc7mqNempra/Hoo4/io48+gru7+3BXh2HHsmXLcNNNN2HixIm47rrrsHfvXgDABx98MMw1Y3AchylTpuCll17C5MmT8cADD+D++++3cTFkDD/vvfceli1bhoiIiAEfK70C9bkiPPzww7jlllt63ScuLo5uNzQ0YMGCBcjIyMCWLVuucO1GNwOVDWP4CQoKgkQi6TELrlQqe8yWMxiMnjzyyCP46quvkJWVhaioqOGuDuMCcrkcCQkJAIBp06bh1KlTeO211/D2228Pc81GN2fOnIFSqcTUqVNpmcViQVZWFl5//XUYDAZIJJJhrCFDiJeXFyZOnIjz588Pd1VGPeHh4T2se1JSUrBr165hqhHDnurqanz//ffYvXv3oI4fMcp4UFAQgoKC+rVvfX09FixYgKlTp2Lr1q0Qi0eMAcCIZCCyYbgGcrkcU6dORWZmJtasWUPLMzMzsXr16mGsGYPh2vA8j0ceeQR79uzB4cOHER8fP9xVYvQCz/MwGAzDXY1Rz6JFi3pE57777ruRnJyMp556iiniLobBYEBRURHmzp073FUZ9cyePbtH+szS0lLExsYOU40Y9mzduhUhISFYsWLFoI4fMcp4f2loaMC1116LmJgY/POf/0RLSwv9juVMHH5qamqgUqlQU1MDi8VCI3cmJCTA29t7eCs3ynjiiSdwxx13YNq0adSCpKamhsVXcAE6OztRVlZGP1dWViI3NxcBAQGIiYkZxpoxNm3ahI8//hhffvklFAoFtS7x9fWFh4fHMNdudPO73/0Oy5YtQ3R0NLRaLT799FMcPnwY33333XBXbdSjUCh6xFXw8vJCYGAgi7fgAjz55JNYtWoVYmJioFQq8cILL0Cj0eCuu+4a7qqNeh5//HHMmjULL730EtavX4+TJ09iy5YtzOrXReA4Dlu3bsVdd90FqXSQajV/lbF161YegMOFMfzcddddDmVz6NCh4a7aqOS///0vHxsby8vlcn7KlCn8kSNHhrtKDJ7nDx065PA+ueuuu4a7aqMeZ++XrVu3DnfVRj333HMPfZ4FBwfzixYt4g8cODDc1WI4Yf78+fyjjz463NVg8Dy/YcMGPjw8nJfJZHxERAS/du1avqCgYLirxbjA119/zaempvJubm58cnIyv2XLluGuEuMC+/fv5wHwJSUlgz7HqMgzzmAwGAwGg8FgMBgMhivBnKkZDAaDwWAwGAwGg8EYYpgyzmAwGAwGg8FgMBgMxhDDlHEGg8FgMBgMBoPBYDCGGKaMMxgMBoPBYDAYDAaDMcQwZZzBYDAYDAaDwWAwGIwhhinjDAaDwWAwGAwGg8FgDDFMGWcwGAwGg8FgMBgMBmOIYco4g8FgMBgMBoPBYDAYQwxTxhkMBoPBYDAYDAaDwRhimDLOYDAYDAaDwWAwGAzGEMOUcQaDwWAwGAwGg8FgMIaY/weYratcof9FyAAAAABJRU5ErkJggg==", @@ -168,7 +169,7 @@ " linestyle=\"--\",\n", ")\n", "for i, x_i in enumerate(xs):\n", - " xi = xi + (1 / (1 + nu)) * (Normal.sufficient_statistics(x_i) - xi)\n", + " xi = xi + (1 / (1 + nu)) * (Normal().sufficient_statistics(x=x_i) - xi)\n", " nu += 1\n", "\n", " if i in [2, 4, 8, 16, 32, 64, 128, 256, 512, 999]:\n", @@ -299,33 +300,33 @@ "\n", "\n", - "\n", - "\n", + "\n", + "\n", "hgf-nodes\n", - "\n", + "\n", "\n", "\n", "x_0\n", - "\n", + "\n", "\n", "\n", "\n", "x_1\n", - "\n", - "1\n", + "\n", + "EF-1\n", "\n", "\n", "\n", "x_1->x_0\n", - "\n", - "\n", + "\n", + "\n", "\n", "\n", "\n" ], "text/plain": [ - "" + "" ] }, "execution_count": 6, @@ -379,9 +380,7 @@ " # set the learning rate\n", " generalised_filter.attributes[1][\"nus\"] = nu\n", "\n", - " means.append(\n", - " generalised_filter.input_data(input_data=xs).to_pandas().x_1_xis_0\n", - " )" + " means.append(generalised_filter.input_data(input_data=xs).to_pandas().x_1_xis_0)" ] }, { @@ -593,7 +592,7 @@ "source": [ "# get the sufficient statistics from the first observation to parametrize the model\n", "sufficient_statistics = jnp.apply_along_axis(\n", - " MultivariateNormal.sufficient_statistics, 1, input_data\n", + " MultivariateNormal().sufficient_statistics, 1, input_data\n", ")" ] }, @@ -791,7 +790,7 @@ "\n" ], "text/plain": [ - "" + "" ] }, "execution_count": 15, @@ -1022,7 +1021,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Last updated: Mon Jun 10 2024\n", + "Last updated: Tue Jun 11 2024\n", "\n", "Python implementation: CPython\n", "Python version : 3.12.3\n", @@ -1032,10 +1031,10 @@ "jax : 0.4.27\n", "jaxlib: 0.4.27\n", "\n", - "jax : 0.4.27\n", - "numpy : 1.26.0\n", - "matplotlib: 3.8.4\n", "seaborn : 0.13.2\n", + "matplotlib: 3.8.4\n", + "numpy : 1.26.0\n", + "jax : 0.4.27\n", "\n", "Watermark: 2.4.3\n", "\n" diff --git a/src/pyhgf/math.py b/src/pyhgf/math.py index 676ee80da..42eab0590 100644 --- a/src/pyhgf/math.py +++ b/src/pyhgf/math.py @@ -1,6 +1,6 @@ # Author: Nicolas Legrand -from typing import Union +from typing import Tuple, Union import jax.numpy as jnp from jax import Array @@ -17,11 +17,13 @@ class MultivariateNormal: """ - def sufficient_statistics(x): + @staticmethod + def sufficient_statistics(x: ArrayLike) -> Array: """Compute the sufficient statistics for the multivariate normal.""" return jnp.hstack([x, jnp.outer(x, x)[jnp.tril_indices(x.shape[0])]]) - def base_measure(k): + @staticmethod + def base_measure(k: int) -> float: """Compute the base measures for the multivariate normal.""" return (2 * jnp.pi) ** (-k / 2) @@ -35,16 +37,30 @@ class Normal: """ - def sufficient_statistics(x): - """Compute the sufficient statistics for the univariate normal.""" + @staticmethod + def sufficient_statistics(x: float) -> Array: + """Sufficient statistics for the univariate normal.""" return jnp.array([x, x**2]) - def base_measure(k): - """Compute the base measure for the univariate normal.""" + @staticmethod + def expected_sufficient_statistics(mu: float, sigma) -> Array: + """Compute expected sufficient statistics from the mean and std.""" + return jnp.array([mu, mu**2 + sigma**2]) + + @staticmethod + def base_measure() -> float: + """Compute the base measure of the univariate normal.""" return 1 / (jnp.sqrt(2 * jnp.pi)) + @staticmethod + def parameters(xis: ArrayLike) -> Tuple[float, float]: + """Get parameters from the expected sufficient statistics.""" + mean = xis[0] + variance = xis[1] - (mean**2) + return mean, variance + -def gaussian_predictive_distribution(x, xi, nu): +def gaussian_predictive_distribution(x: float, xi: ArrayLike, nu: float) -> float: r"""Density of the Gaussian-predictive distribution. This distribution is parametrized by hyperparameters from the exponential family as: @@ -178,7 +194,7 @@ def gaussian_surprise( Examples -------- - >>> from pyhgf.continuous import gaussian_surprise + >>> from pyhgf.math import gaussian_surprise >>> gaussian_surprise(x=2.0, expected_mean=0.0, expected_precision=1.0) `Array(2.9189386, dtype=float32, weak_type=True)` @@ -237,7 +253,7 @@ def binary_surprise_finite_precision( expected_mean: Union[ArrayLike, float], expected_precision: Union[ArrayLike, float], eta0: Union[ArrayLike, float] = 0.0, - eta1: Union[ArrayLike, float] = 0.0, + eta1: Union[ArrayLike, float] = 1.0, ) -> Array: r"""Compute the binary surprise with finite precision. @@ -264,3 +280,8 @@ def binary_surprise_finite_precision( expected_mean * gaussian_density(value, eta1, expected_precision) + (1 - expected_mean) * gaussian_density(value, eta0, expected_precision) ) + + +def sigmoid_inverse_temperature(x, temperature): + """Compute the sigmoid response function with inverse temperature parameter.""" + return (x**temperature) / (x**temperature + (1 - x) ** temperature) diff --git a/src/pyhgf/model/network.py b/src/pyhgf/model/network.py index 4d026e29b..35285fbbd 100644 --- a/src/pyhgf/model/network.py +++ b/src/pyhgf/model/network.py @@ -1,5 +1,6 @@ # Author: Nicolas Legrand +from copy import deepcopy from typing import Callable, Dict, List, Optional, Tuple, Union import jax.numpy as jnp @@ -20,6 +21,7 @@ input_types, ) from pyhgf.utils import ( + add_edges, beliefs_propagation, fill_categorical_state_node, get_update_sequence, @@ -28,7 +30,7 @@ class Network: - """A generalised HGF neural network for predictive coding applications. + """A predictive coding neural network. This is the core class to define and manipulate neural networks, that consists in 1. attributes, 2. structure and 3. update sequences. @@ -189,12 +191,13 @@ def input_data( # this is where the model loops over the whole input time series # at each time point, the node structure is traversed and beliefs are updated # using precision-weighted prediction errors - _, node_trajectories = scan( + last_attributes, node_trajectories = scan( self.scan_fn, self.attributes, (input_data, time_steps, observed) ) # trajectories of the network attributes a each time point self.node_trajectories = node_trajectories + self.last_attributes = last_attributes return self @@ -399,13 +402,31 @@ def add_nodes( attributes. """ - # extract the node coupling indexes and coupling strengths + if kind not in [ + "continuous-input", + "binary-input", + "categorical-input", + "DP-state", + "ef-normal", + "generic-input", + "continuous-state", + "binary-state", + ]: + raise ValueError( + ( + "Invalid node type. Should be one of the following: " + "'continuous-input', 'binary-input', 'categorical-input', " + "'DP-state', 'continuous-state', 'binary-state', 'ef-normal'." + ) + ) + + # transform coupling parameter into tuple of indexes and strenghts couplings = [] for indexes in [ - value_children, value_parents, - volatility_children, volatility_parents, + value_children, + volatility_children, ]: if indexes is not None: if isinstance(indexes, int): @@ -420,6 +441,9 @@ def add_nodes( else: coupling_idxs, coupling_strengths = None, None couplings.append((coupling_idxs, coupling_strengths)) + value_parents, volatility_parents, value_children, volatility_children = ( + couplings + ) # create the default parameters set according to the node type if kind == "continuous-state": @@ -428,10 +452,10 @@ def add_nodes( "expected_mean": 0.0, "precision": 1.0, "expected_precision": 1.0, - "volatility_coupling_children": couplings[2][1], - "volatility_coupling_parents": couplings[3][1], - "value_coupling_children": couplings[0][1], - "value_coupling_parents": couplings[1][1], + "volatility_coupling_children": volatility_children[1], + "volatility_coupling_parents": volatility_parents[1], + "value_coupling_children": value_children[1], + "value_coupling_parents": value_parents[1], "tonic_volatility": -4.0, "tonic_drift": 0.0, "autoconnection_strength": 1.0, @@ -449,10 +473,10 @@ def add_nodes( "expected_mean": 0.0, "precision": 1.0, "expected_precision": 1.0, - "volatility_coupling_children": couplings[2][1], - "volatility_coupling_parents": couplings[3][1], - "value_coupling_children": couplings[0][1], - "value_coupling_parents": couplings[1][1], + "volatility_coupling_children": volatility_children[1], + "volatility_coupling_parents": volatility_parents[1], + "value_coupling_children": value_children[1], + "value_coupling_parents": value_parents[1], "tonic_volatility": 0.0, "tonic_drift": 0.0, "autoconnection_strength": 1.0, @@ -536,11 +560,35 @@ def add_nodes( } elif "ef-normal" in kind: default_parameters = { - "nus": 0.0, - "xis": jnp.array([0.0, 0.0]), + "nus": 3.0, + "xis": jnp.array([0.0, 1.0]), "values": 0.0, + "observed": 1.0, } + elif kind == "DP-state": + if "batch_size" in additional_parameters.keys(): + batch_size = additional_parameters["batch_size"] + elif "batch_size" in node_parameters.keys(): + batch_size = node_parameters["batch_size"] + else: + batch_size = 10 + + default_parameters = { + "batch_size": batch_size, # number of branches available in the network + "n": jnp.zeros(batch_size), # number of observation in each cluster + "n_total": 0, # the total number of observations in the node + "alpha": 1.0, # concentration parameter for the implied Dirichlet dist. + "expected_means": jnp.zeros(batch_size), + "expected_sigmas": jnp.ones(batch_size), + "sensory_precision": 1.0, + "activated": jnp.zeros(batch_size), + "value_coupling_children": (1.0,), + "values": 0.0, + "n_active_cluster": 0, + } + + # Update the default node parameters using keywords args and dictonary if bool(additional_parameters): # ensure that all passed values are valid keys invalid_keys = [ @@ -581,32 +629,37 @@ def add_nodes( node_type = 2 elif "ef-normal" in kind: node_type = 3 - - # convert the structure to a list to modify it - edges_as_list: List[AdjacencyLists] = list(self.edges) + elif "DP-state" in kind: + node_type = 4 for _ in range(n_nodes): + # convert the structure to a list to modify it + edges_as_list: List = list(self.edges) + node_idx = len(self.attributes) # the index of the new node # add a new edge edges_as_list.append( AdjacencyLists( node_type, - couplings[1][0], - couplings[3][0], - couplings[0][0], - couplings[2][0], + None, + None, + None, + None, ) ) + # convert the list back to a tuple + self.edges = tuple(edges_as_list) + if node_idx == 0: # this is the first node, create the node structure - self.attributes = {node_idx: node_parameters} + self.attributes = {node_idx: deepcopy(node_parameters)} if input_type is not None: self.inputs = Inputs((node_idx,), (input_type,)) else: # update the node structure - self.attributes[node_idx] = node_parameters + self.attributes[node_idx] = deepcopy(node_parameters) if input_type is not None: # add information about the new input node in the indexes @@ -616,91 +669,40 @@ def add_nodes( new_kind += (input_type,) self.inputs = Inputs(new_idx, new_kind) - # update the existing edge structure so it links to the new node as well - for coupling, edge_type in zip( - couplings, - [ - "value_children", - "value_parents", - "volatility_children", - "volatility_parents", - ], - ): - if coupling[0] is not None: - coupling_idxs, coupling_strengths = coupling - for idx, coupling_strength in zip( - coupling_idxs, coupling_strengths # type: ignore - ): - # unpack this node's edges - ( - this_node_type, - value_parents, - volatility_parents, - value_children, - volatility_children, - ) = edges_as_list[idx] - - # update the parents/children's edges depending on the coupling - if edge_type == "value_parents": - if value_children is None: - value_children = (node_idx,) - self.attributes[idx]["value_coupling_children"] = ( - coupling_strength, - ) - else: - value_children = value_children + (node_idx,) - self.attributes[idx]["value_coupling_children"] += ( - coupling_strength, - ) - elif edge_type == "volatility_parents": - if volatility_children is None: - volatility_children = (node_idx,) - self.attributes[idx]["volatility_coupling_children"] = ( - coupling_strength, - ) - else: - volatility_children = volatility_children + (node_idx,) - self.attributes[idx][ - "volatility_coupling_children" - ] += (coupling_strength,) - elif edge_type == "value_children": - if value_parents is None: - value_parents = (node_idx,) - self.attributes[idx]["value_coupling_parents"] = ( - coupling_strength, - ) - else: - value_parents = value_parents + (node_idx,) - self.attributes[idx]["value_coupling_parents"] += ( - coupling_strength, - ) - elif edge_type == "volatility_children": - if volatility_parents is None: - volatility_parents = (node_idx,) - self.attributes[idx]["volatility_coupling_parents"] = ( - coupling_strength, - ) - else: - volatility_parents = volatility_parents + (node_idx,) - self.attributes[idx]["volatility_coupling_parents"] += ( - coupling_strength, - ) - - # save the updated edges back - edges_as_list[idx] = AdjacencyLists( - this_node_type, - value_parents, - volatility_parents, - value_children, - volatility_children, - ) - - # convert the list back to a tuple - self.edges = tuple(edges_as_list) - - # if we are creating a categorical state or state-transition node - # we have to generate the implied binary network(s) here + # Update the edges of the parents and children accordingly + # -------------------------------------------------------- + if value_parents[0] is not None: + self.add_edges( + kind="value", + parent_idxs=value_parents[0], + children_idxs=node_idx, + coupling_strengths=value_parents[1], # type: ignore + ) + if value_children[0] is not None: + self.add_edges( + kind="value", + parent_idxs=node_idx, + children_idxs=value_children[0], + coupling_strengths=value_children[1], # type: ignore + ) + if volatility_children[0] is not None: + self.add_edges( + kind="volatility", + parent_idxs=node_idx, + children_idxs=volatility_children[0], + coupling_strengths=volatility_children[1], # type: ignore + ) + if volatility_parents[0] is not None: + self.add_edges( + kind="volatility", + parent_idxs=volatility_parents[0], + children_idxs=node_idx, + coupling_strengths=volatility_parents[1], # type: ignore + ) + if kind == "categorical-input": + # if we are creating a categorical state or state-transition node + # we have to generate the implied binary network(s) here self = fill_categorical_state_node( self, node_idx=node_idx, @@ -778,3 +780,39 @@ def surprise( response_function_inputs=response_function_inputs, response_function_parameters=response_function_parameters, ) + return self + + def add_edges( + self, + kind="value", + parent_idxs=Union[int, List[int]], + children_idxs=Union[int, List[int]], + coupling_strengths: Union[float, List[float], Tuple[float]] = 1.0, + ) -> "Network": + """Add a value or volatility coupling link between a set of nodes. + + Parameters + ---------- + kind : + The kind of coupling, can be `"value"` or `"volatility"`. + parent_idxs : + The index(es) of the parent node(s). + children_idxs : + The index(es) of the children node(s). + coupling_strengths : + The coupling strength betwen the parents and children. + + """ + attributes, edges = add_edges( + attributes=self.attributes, + edges=self.edges, + kind=kind, + parent_idxs=parent_idxs, + children_idxs=children_idxs, + coupling_strengths=coupling_strengths, + ) + + self.attributes = attributes + self.edges = edges + + return self diff --git a/src/pyhgf/plots.py b/src/pyhgf/plots.py index 59fc6c925..c5fbe6c8d 100644 --- a/src/pyhgf/plots.py +++ b/src/pyhgf/plots.py @@ -274,10 +274,31 @@ def plot_network(network: "Network") -> "Source": ) # create the rest of nodes - for i in range(len(network.edges)): - # only if node is not an input node - if i not in network.inputs.idx: - graphviz_structure.node(f"x_{i}", label=str(i), shape="circle") + for idx in range(len(network.edges)): + + if network.edges[idx].node_type == 2: + # Continuous state nore + graphviz_structure.node(f"x_{idx}", label=str(idx), shape="circle") + + elif network.edges[idx].node_type == 3: + # Exponential family state nore + graphviz_structure.node( + f"x_{idx}", + label=f"EF-{idx}", + style="filled", + shape="circle", + fillcolor="#ced6e4", + ) + + elif network.edges[idx].node_type == 4: + # Dirichlet PRocess state node + graphviz_structure.node( + f"x_{idx}", + label=f"DP-{idx}", + style="filled", + shape="doublecircle", + fillcolor="#e2d8c1", + ) # connect value parents for i, index in enumerate(network.edges): diff --git a/src/pyhgf/typing.py b/src/pyhgf/typing.py index 751c81132..d7a17b702 100644 --- a/src/pyhgf/typing.py +++ b/src/pyhgf/typing.py @@ -12,6 +12,7 @@ class AdjacencyLists(NamedTuple): * 2: continuous state node. * 3: exponential family state node - univariate Gaussian distribution with unknown mean and unknown variance. + * 4: Dirichlet Process state node. """ diff --git a/src/pyhgf/updates/posterior/exponential.py b/src/pyhgf/updates/posterior/exponential.py index 57f5f28c3..ca0ee4579 100644 --- a/src/pyhgf/updates/posterior/exponential.py +++ b/src/pyhgf/updates/posterior/exponential.py @@ -3,6 +3,7 @@ from functools import partial from typing import Callable, Dict +import jax.numpy as jnp from jax import jit from pyhgf.typing import Attributes, Edges @@ -49,11 +50,14 @@ def posterior_update_exponential_family( """ # update the hyperparameter vectors - attributes[node_idx]["xis"] = attributes[node_idx]["xis"] + ( - 1 / (1 + attributes[node_idx]["nus"]) - ) * ( - sufficient_stats_fn(attributes[node_idx]["values"]) + xis = attributes[node_idx]["xis"] + (1 / (1 + attributes[node_idx]["nus"])) * ( + sufficient_stats_fn(x=attributes[node_idx]["values"]) - attributes[node_idx]["xis"] ) + # blank update in the case of unobserved value + attributes[node_idx]["xis"] = jnp.where( + attributes[node_idx]["observed"], xis, attributes[node_idx]["xis"] + ) + return attributes diff --git a/src/pyhgf/updates/prediction/dirichlet.py b/src/pyhgf/updates/prediction/dirichlet.py new file mode 100644 index 000000000..03fb86e25 --- /dev/null +++ b/src/pyhgf/updates/prediction/dirichlet.py @@ -0,0 +1,55 @@ +# Author: Nicolas Legrand + +from typing import Dict + +import jax.numpy as jnp + +from pyhgf.math import Normal +from pyhgf.typing import Attributes, Edges + + +def dirichlet_node_prediction( + edges: Edges, + attributes: Dict, + node_idx: int, + **args, +) -> Attributes: + """Prediction of a Dirichlet process node. + + Parameters + ---------- + edges : + The edges of the neural network as a tuple of + :py:class:`pyhgf.typing.Indexes`. The tuple has the same length as node number. + For each node, the index lists the value/volatility parents/children. + attributes : + The attributes of the probabilistic nodes. + node_idx : + Pointer to the Dirichlet process input node. + + Returns + ------- + attributes : + The attributes of the probabilistic nodes. + edges : + The edges of the neural network. + input_nodes_idx : + Static input nodes' parameters for the neural network. + dirichlet_node : + Static parameters of the Dirichlet process node. + + """ + # get the parameter (mean and variance) from the EF-normal parent nodes + value_parent_idxs = edges[node_idx].value_parents + if value_parent_idxs is not None: + parameters = jnp.array( + [ + Normal().parameters(xis=attributes[parent_idx]["xis"]) + for parent_idx in value_parent_idxs + ] + ) + + attributes[node_idx]["expected_means"] = parameters[:, 0] + attributes[node_idx]["expected_sigmas"] = jnp.sqrt(parameters[:, 1]) + + return attributes diff --git a/src/pyhgf/updates/prediction_error/nodes/dirichlet.py b/src/pyhgf/updates/prediction_error/nodes/dirichlet.py new file mode 100644 index 000000000..cc2a50d3b --- /dev/null +++ b/src/pyhgf/updates/prediction_error/nodes/dirichlet.py @@ -0,0 +1,402 @@ +# Author: Nicolas Legrand + +from functools import partial +from typing import Dict, Tuple + +import jax.numpy as jnp +from jax import Array, jit, random +from jax._src.typing import Array as KeyArray +from jax.lax import cond +from jax.scipy.stats.norm import pdf +from jax.tree_util import Partial +from jax.typing import ArrayLike + +from pyhgf.math import Normal +from pyhgf.typing import Attributes, Edges + + +@partial(jit, static_argnames=("edges", "node_idx")) +def dirichlet_node_prediction_error( + edges: Edges, + attributes: Dict, + node_idx: int, + **args, +) -> Attributes: + """Prediction error and update the child networks of a Dirichlet process node. + + When receiving a new input, this node chose to either: + 1. Allocate the value to a pre-existing cluster. + 2. Create a new cluster. + + The network always contains a temporary branch as the new cluster candidate. This + branch is parametrized under the new observation to assess its likelihood and the + previous clusters' likelihood. + + Parameters + ---------- + edges : + The edges of the neural network as a tuple of + :py:class:`pyhgf.typing.Indexes`. The tuple has the same length as node number. + For each node, the index lists the value/volatility parents/children. + attributes : + The attributes of the probabilistic nodes. + node_idx : + Pointer to the Dirichlet process input node. + + Returns + ------- + attributes : + The attributes of the probabilistic nodes. + + """ + values = attributes[node_idx]["values"] # the input value + alpha = attributes[node_idx]["alpha"] # the concentration parameter + n_total = attributes[node_idx]["n_total"] # total number of observations + n = attributes[node_idx]["n"] # number of observations per cluster + sensory_precision = attributes[node_idx][ + "sensory_precision" + ] # number of observations per cluster + + # likelihood of the current observation under existing clusters + # ------------------------------------------------------------- + cluster_ll = clusters_likelihood( + value=values, + expected_mean=attributes[node_idx]["expected_means"], + expected_sigma=attributes[node_idx]["expected_sigmas"], + ) + + # set the likelihood to 0 for inactive clusters + cluster_ll *= attributes[node_idx]["activated"] + + # likelihood of the current observation under the best candidate cluster + # ---------------------------------------------------------------------- + + # find the best cluster candidate given the new observation + candidate_mean, candidate_sigma = get_candidate( + value=values, + sensory_precision=sensory_precision, + expected_mean=attributes[node_idx]["expected_means"], + expected_sigma=attributes[node_idx]["expected_sigmas"], + ) + + # get the likelihood under this candidate + candidate_ll = clusters_likelihood( + value=values, + expected_mean=candidate_mean, + expected_sigma=candidate_sigma, + ) + + # DP step: compare the likelihood of existing cluster with a new cluster + # ---------------------------------------------------------------------- + + # probability of being assigned to a pre-existing cluster + cluster_ll *= n / (alpha + n_total) + + # probability to draw a new cluster + candidate_ll *= alpha / (alpha + n_total) + + best_val = jnp.max(cluster_ll) + + # set all cluster to non-observed by default + for parent_idx in edges[node_idx].value_parents: # type:ignore + attributes[parent_idx]["observed"] = 0 + + # get the index of the cluster (!= the node index) + # depending on whether a new cluster is created or updated + cluster_idx = jnp.where( + best_val >= candidate_ll, + jnp.argmax(cluster_ll), + attributes[node_idx]["n_active_cluster"], + ) + + update_fn = Partial( + update_cluster, + edges=edges, + node_idx=node_idx, + ) + + create_fn = Partial( + create_cluster, + edges=edges, + node_idx=node_idx, + ) + + # apply either cluster update or cluster creation + operands = attributes, cluster_idx, values, (candidate_mean, candidate_sigma) + + attributes = cond(best_val >= candidate_ll, update_fn, create_fn, operands) + + attributes[node_idx]["n_total"] += 1 + + return attributes + + +@partial(jit, static_argnames=("edges", "node_idx")) +def update_cluster(operands: Tuple, edges: Edges, node_idx: int) -> Attributes: + """Update an existing cluster. + + Parameters + ---------- + operands : + Non-static parameters. + edges : + The edges of the neural network as a tuple of + :py:class:`pyhgf.typing.Indexes`. The tuple has the same length as node number. + For each node, the index lists the value/volatility parents/children. + node_idx : + Pointer to the Dirichlet process input node. + + Returns + ------- + attributes : + The attributes of the probabilistic nodes. + + """ + attributes, cluster_idx, value, _ = operands + + # activate the corresponding branch and pass the value + for i, value_parent_idx in enumerate(edges[node_idx].value_parents): # type: ignore + + attributes[value_parent_idx]["observed"] = jnp.where(cluster_idx == i, 1.0, 0.0) + attributes[value_parent_idx]["values"] = value + + attributes[node_idx]["n"] = ( + attributes[node_idx]["n"] + .at[cluster_idx] + .set(attributes[node_idx]["n"][cluster_idx] + 1.0) + ) + + return attributes + + +@partial(jit, static_argnames=("edges", "node_idx")) +def create_cluster(operands: Tuple, edges: Edges, node_idx: int) -> Attributes: + """Create a new cluster. + + Parameters + ---------- + operands : + Non-static parameters. + edges : + The edges of the neural network as a tuple of + :py:class:`pyhgf.typing.Indexes`. The tuple has the same length as node number. + For each node, the index lists the value/volatility parents/children. + node_idx : + Pointer to the Dirichlet process input node. + + Returns + ------- + attributes : + The attributes of the probabilistic nodes. + + """ + attributes, cluster_idx, value, (candidate_mean, candidate_sigma) = operands + + # creating a new cluster + attributes[node_idx]["activated"] = ( + attributes[node_idx]["activated"].at[cluster_idx].set(1) + ) + + for i, value_parent_idx in enumerate(edges[node_idx].value_parents): # type: ignore + + attributes[value_parent_idx]["observed"] = 0.0 + attributes[value_parent_idx]["values"] = value + + # initialize the new cluster using candidate values + attributes[value_parent_idx]["xis"] = jnp.where( + cluster_idx == i, + Normal().expected_sufficient_statistics( + mu=candidate_mean, sigma=candidate_sigma + ), + attributes[value_parent_idx]["xis"], + ) + + attributes[node_idx]["n"] = attributes[node_idx]["n"].at[cluster_idx].set(1.0) + attributes[node_idx]["n_active_cluster"] += 1 + + return attributes + + +@jit +def get_candidate( + value: float, + sensory_precision: float, + expected_mean: ArrayLike, + expected_sigma: ArrayLike, + n_samples: int = 20_000, +) -> Tuple[float, float]: + """Find the best cluster candidate given previous clusters and an input value. + + Parameters + ---------- + value : + The new observation. + sensory_precision : + The expected precision of the new observation. + expected_mean : + The mean of the existing clusters. + expected_sigma : + The standard deviation of the existing clusters. + n_samples : + The number of samples that should be simulated. + + Returns + ------- + mean : + The mean of the new candidate cluster. + sigma : + The standard deviation of the new candidate cluster. + + """ + # sample n likely clusters given the base distribution priors + mus, sigmas, weights = likely_cluster_proposal( + mean_mu_G0=0.0, + sigma_mu_G0=10.0, + sigma_pi_G0=3.0, + expected_mean=expected_mean, + expected_sigma=expected_sigma, + key=random.key(42), + n_samples=n_samples, + ) + + # 1 - Likelihood of the new observation under each sampled cluster + # ---------------------------------------------------------------- + ll_value = pdf(value, mus, sigmas) + ll_value /= ll_value.sum() # normalize the weights + + # 2- re-scale the weights using expected precision + # ------------------------------------------------ + weights *= ll_value**sensory_precision + + # only use the 1000 best candidates for inference + idxs = jnp.argsort(weights) + mus, sigmas, weights = ( + mus[idxs][-1000:], + sigmas[idxs][-1000:], + weights[idxs][-1000:], + ) + + # 3 - estimate new mean and standard deviation using the weigthed mean + # -------------------------------------------------------------------- + mean = jnp.average(mus, weights=weights) + sigma = jnp.average(sigmas, weights=weights) + + return mean, sigma + + +@partial(jit, static_argnames=("n_samples")) +def likely_cluster_proposal( + mean_mu_G0: float, + sigma_mu_G0: float, + sigma_pi_G0: float, + expected_mean=ArrayLike, + expected_sigma=ArrayLike, + key: KeyArray = random.key(42), + n_samples: int = 20_000, +) -> Tuple[Array, Array, Array]: + """Sample likely new belief distributions given pre-existing clusters. + + Parameters + ---------- + mean_mu_G0 : + The mean of the mean of the base distribution. + sigma_mu_G0 : + The standard deviation of mean of the base distribution. + sigma_pi_G0 : + The standard deviation of the standard deviation of the base distribution. + expected_mean : + Pre-existing clusters means. + expected_sigma : + Pre-existing clusters standard deviation. + key : + Random state. + n_samples : + The number of samples used during the simulations. + + Returns + ------- + new_mu : + A vector of means candidates. + new_sigma : + A vector of standard deviation candidates. + weights : + Weigths for each cluster candidate under pre-existing cluster (irrespective of + new observations). + + """ + # sample new candidate for cluster means + key, use_key = random.split(key) + new_mu = sigma_mu_G0 * random.normal(use_key, (n_samples,)) + mean_mu_G0 + + # sample new candidate for cluster standard deviation + key, use_key = random.split(key) + new_sigma = jnp.abs(random.normal(use_key, (n_samples,)) * sigma_pi_G0) + + # 1 - Cluster specificity + # ----------------------- + # this cluster should explain new dimensions, not explained by other clusters + + # evidence for pre-existing clusters + pre_existing_likelihood = jnp.zeros(n_samples) + for mu_i, sigma_i in zip(expected_mean, expected_sigma): + pre_existing_likelihood += pdf(new_mu, mu_i, sigma_i) + + # evidence for the new cluster proposal + new_likelihood = pdf(new_mu, new_mu, new_sigma) + + # standardize the measure of cluster specificity (ratio) + ratio = new_likelihood / (new_likelihood + pre_existing_likelihood) + ratio -= ratio.min() + ratio /= ratio.max() + weights = ratio + + # 2 - Cluster isolation + # --------------------- + # this cluster should not try to explain what was already explained + + # (pre-existing cluster) / (pre-existing cluster + new cluster) + cluster_isolation = jnp.ones(n_samples) + for mu_i, sigma_i in zip(expected_mean, expected_sigma): + ratio = pdf(mu_i, mu_i, sigma_i) / ( + pdf(mu_i, mu_i, sigma_i) + pdf(mu_i, new_mu, new_sigma) + ) + cluster_isolation *= ratio + cluster_isolation -= cluster_isolation.min() + cluster_isolation /= cluster_isolation.max() + + weights *= cluster_isolation + + # 3 - Spread of the cluster + # ------------------------- + # large clusters should be favored over small clusters + cluster_spread = pdf(1 / (new_sigma**2), 0.0, 5.0) + cluster_spread -= cluster_spread.min() + cluster_spread /= cluster_spread.max() + weights *= cluster_spread + + return new_mu, new_sigma, weights + + +def clusters_likelihood( + value: float, + expected_mean: ArrayLike, + expected_sigma: ArrayLike, +) -> ArrayLike: + """Likelihood of a parametrized candidate under the new observation. + + Parameters + ---------- + value : + The new observation. + expected_mean : + Pre-existing clusters means. + expected_sigma : + Pre-existing clusters standard deviation. + + Returns + ------- + likelihood : + The probability of observing the value under each cluster. + + """ + return pdf(value, expected_mean, expected_sigma) diff --git a/src/pyhgf/utils.py b/src/pyhgf/utils.py index bbc36d656..8602a78e8 100644 --- a/src/pyhgf/utils.py +++ b/src/pyhgf/utils.py @@ -1,7 +1,7 @@ # Author: Nicolas Legrand from functools import partial -from typing import TYPE_CHECKING, Dict, List, Tuple +from typing import TYPE_CHECKING, Dict, List, Tuple, Union import jax.numpy as jnp import numpy as np @@ -11,7 +11,7 @@ from jax.typing import ArrayLike from pyhgf.math import Normal, binary_surprise, gaussian_surprise -from pyhgf.typing import AdjacencyLists, Attributes, Structure, UpdateSequence +from pyhgf.typing import AdjacencyLists, Attributes, Edges, Structure, UpdateSequence from pyhgf.updates.posterior.binary import binary_node_update_infinite from pyhgf.updates.posterior.categorical import categorical_input_update from pyhgf.updates.posterior.continuous import ( @@ -21,6 +21,7 @@ from pyhgf.updates.posterior.exponential import posterior_update_exponential_family from pyhgf.updates.prediction.binary import binary_state_node_prediction from pyhgf.updates.prediction.continuous import continuous_node_prediction +from pyhgf.updates.prediction.dirichlet import dirichlet_node_prediction from pyhgf.updates.prediction_error.inputs.binary import ( binary_input_prediction_error_infinite_precision, ) @@ -34,6 +35,9 @@ from pyhgf.updates.prediction_error.nodes.continuous import ( continuous_node_prediction_error, ) +from pyhgf.updates.prediction_error.nodes.dirichlet import ( + dirichlet_node_prediction_error, +) if TYPE_CHECKING: from pyhgf.model import Network @@ -110,39 +114,6 @@ def beliefs_propagation( ) # ("carryover", "accumulated") -def trim_sequence( - exclude_node_idxs: List, update_sequence: UpdateSequence, edges: Tuple -) -> UpdateSequence: - """Remove steps from an update sequence that depends on a set of nodes. - - Parameters - ---------- - exclude_node_idxs : - A list of node indexes. The nodes can be input nodes or any other node in the - network. - update_sequence : - The sequence of updates that will be applied to the node structure. - edges : - The nodes structure. - - Returns - ------- - trimmed_update_sequence : - The update sequence without the update steps for nodes depending on the root - list. - - """ - # list the nodes that depend on the root indexes - branch_list = list_branches(node_idxs=exclude_node_idxs, edges=edges) - - # remove the update steps that are targetting the excluded nodes - trimmed_update_sequence = tuple( - [seq for seq in update_sequence if seq[0] not in branch_list] - ) - - return trimmed_update_sequence - - def list_branches(node_idxs: List, edges: Tuple, branch_list: List = []) -> List: """Return the branch of a network from a given set of root nodes. @@ -322,12 +293,14 @@ def get_update_sequence(network: "Network", update_type: str) -> List: node_without_update = [i for i in range(n_nodes)] # start by injecting the observations in all input nodes + # ------------------------------------------------------ for input_idx, kind in zip(network.inputs.idx, network.inputs.kind): if kind == 0: update_fn = continuous_input_prediction_error update_sequence.append((input_idx, update_fn)) elif kind == 1: + # add the update steps for the binary state node as well binary_state_idx = network.edges[input_idx].value_parents[0] # type: ignore @@ -356,13 +329,19 @@ def get_update_sequence(network: "Network", update_type: str) -> List: update_fn = generic_input_prediction_error update_sequence.append((input_idx, update_fn)) + elif kind == 4: + update_fn = dirichlet_node_prediction_error + update_sequence.append((input_idx, update_fn)) + # add the PE step to the sequence node_without_pe.remove(input_idx) # input node does not need to update the posterior node_without_update.remove(input_idx) + # prediction errors and posterior updates # will fail if the structure of the network does not allow a consistent update order + # ---------------------------------------------------------------------------------- while True: no_update = True @@ -400,10 +379,16 @@ def get_update_sequence(network: "Network", update_type: str) -> List: # for the exponential family node ef_update = Partial( posterior_update_exponential_family, - sufficient_stats_fn=Normal.sufficient_statistics, + sufficient_stats_fn=Normal().sufficient_statistics, ) update_fn = ef_update + elif network.edges[idx].node_type == 4: + + update_fn = None + # the prediction sequence is the update sequence in reverse order + prediction_sequence.insert(0, (idx, dirichlet_node_prediction)) + update_sequence.append((idx, update_fn)) node_without_update.remove(idx) @@ -425,15 +410,20 @@ def get_update_sequence(network: "Network", update_type: str) -> List: else: # if this node has been updated if idx not in node_without_update: + + if network.edges[idx].node_type == 2: + update_fn = continuous_node_prediction_error + elif network.edges[idx].node_type == 4: + update_fn = dirichlet_node_prediction_error + no_update = False - update_sequence.append((idx, continuous_node_prediction_error)) + update_sequence.append((idx, update_fn)) node_without_pe.remove(idx) if (not node_without_pe) and (not node_without_update): break if no_update: - break raise Warning( "The structure of the network cannot be updated consistently." ) @@ -447,7 +437,10 @@ def get_update_sequence(network: "Network", update_type: str) -> List: # create a new sequence step and add it to the list prediction_sequence.append((idx, categorical_input_update)) - return prediction_sequence + # remove None steps and return the update sequence + sequence = [update for update in prediction_sequence if update[1] is not None] + + return sequence def to_pandas(network: "Network") -> pd.DataFrame: @@ -599,3 +592,204 @@ def to_pandas(network: "Network") -> pd.DataFrame: ].sum(axis=1, min_count=1) return trajectories_df + + +def concatenate_networks(attributes_1, attributes_2, edges_1, edges_2): + """Concatenate two networks. + + Parameters + ---------- + attributes_1 : + The attributes of the first network. + attributes_2 : + The attributes of the second network. + edges_1 : + The edges of the first network. + edges_2 : + The edges of the second network. + + Returns + ------- + attributes : + The attribute of the concatenated networks. + edges : + The edges of the concatenated networks. + + """ + n_nodes = len(attributes_2) + edges_1 = list(edges_1) + attributes = {} + for i in range(len(attributes_1)): + # update the attributes + attributes[i + n_nodes] = attributes_1[i] + + # update the edges + edges_1[i] = AdjacencyLists( + value_parents=( + tuple([e + n_nodes for e in list(edges_1[i].value_parents)]) + if edges_1[i].value_parents is not None + else None + ), + volatility_parents=( + tuple([e + n_nodes for e in list(edges_1[i].volatility_parents)]) + if edges_1[i].volatility_parents is not None + else None + ), + value_children=( + tuple([e + n_nodes for e in list(edges_1[i].value_children)]) + if edges_1[i].value_children is not None + else None + ), + volatility_children=( + tuple([e + n_nodes for e in list(edges_1[i].volatility_children)]) + if edges_1[i].volatility_children is not None + else None + ), + ) + + edges_1 = tuple(edges_1) + + attributes = {**attributes_2, **attributes} + edges = edges_2 + edges_1 + + return attributes, edges + + +def add_edges( + attributes: Dict, + edges: Edges, + kind="value", + parent_idxs=Union[int, List[int]], + children_idxs=Union[int, List[int]], + coupling_strengths: Union[float, List[float], Tuple[float]] = 1.0, +) -> Tuple: + """Add a value or volatility coupling link between a set of nodes. + + Parameters + ---------- + attributes : + Attributes of the neural network. + edges : + Edges of the neural network. + kind : + The kind of coupling can be `"value"` or `"volatility"`. + parent_idxs : + The index(es) of the parent node(s). + children_idxs : + The index(es) of the children node(s). + coupling_strengths : + The coupling strength between the parents and children. + + """ + if kind not in ["value", "volatility"]: + raise ValueError( + f"The kind of coupling should be value or volatility, got {kind}" + ) + if isinstance(children_idxs, int): + children_idxs = [children_idxs] + assert isinstance(children_idxs, (list, tuple)) + + if isinstance(parent_idxs, int): + parent_idxs = [parent_idxs] + assert isinstance(parent_idxs, (list, tuple)) + + if isinstance(coupling_strengths, int): + coupling_strengths = [float(coupling_strengths)] + if isinstance(coupling_strengths, float): + coupling_strengths = [coupling_strengths] + + assert isinstance(coupling_strengths, (list, tuple)) + + edges_as_list = list(edges) + # update the parent nodes + # ----------------------- + for parent_idx in parent_idxs: + # unpack the parent's edges + ( + node_type, + value_parents, + volatility_parents, + value_children, + volatility_children, + ) = edges_as_list[parent_idx] + + if kind == "value": + if value_children is None: + value_children = tuple(children_idxs) + attributes[parent_idx]["value_coupling_children"] = tuple( + coupling_strengths + ) + else: + value_children = value_children + tuple(children_idxs) + attributes[parent_idx]["value_coupling_children"] += tuple( + coupling_strengths + ) + elif kind == "volatility": + if volatility_children is None: + volatility_children = tuple(children_idxs) + attributes[parent_idx]["volatility_coupling_children"] = tuple( + coupling_strengths + ) + else: + volatility_children = volatility_children + tuple(children_idxs) + attributes[parent_idx]["volatility_coupling_children"] += tuple( + coupling_strengths + ) + + # save the updated edges back + edges_as_list[parent_idx] = AdjacencyLists( + node_type, + value_parents, + volatility_parents, + value_children, + volatility_children, + ) + + # update the children nodes + # ------------------------- + for children_idx in children_idxs: + # unpack this node's edges + ( + node_type, + value_parents, + volatility_parents, + value_children, + volatility_children, + ) = edges_as_list[children_idx] + + if kind == "value": + if value_parents is None: + value_parents = tuple(parent_idxs) + attributes[children_idx]["value_coupling_parents"] = tuple( + coupling_strengths + ) + else: + value_parents = value_parents + tuple(parent_idxs) + attributes[children_idx]["value_coupling_parents"] += tuple( + coupling_strengths + ) + elif kind == "volatility": + if volatility_parents is None: + volatility_parents = tuple(parent_idxs) + attributes[children_idx]["volatility_coupling_parents"] = tuple( + coupling_strengths + ) + else: + volatility_parents = volatility_parents + tuple(parent_idxs) + attributes[children_idx]["volatility_coupling_parents"] += tuple( + coupling_strengths + ) + + # save the updated edges back + edges_as_list[children_idx] = AdjacencyLists( + node_type, + value_parents, + volatility_parents, + value_children, + volatility_children, + ) + + # convert the list back to a tuple + edges = tuple(edges_as_list) + + return attributes, edges diff --git a/tests/test_math.py b/tests/test_math.py new file mode 100644 index 000000000..146498b0c --- /dev/null +++ b/tests/test_math.py @@ -0,0 +1,56 @@ +# Author: Nicolas Legrand + +import unittest +from unittest import TestCase + +import jax.numpy as jnp + +from pyhgf.math import ( + MultivariateNormal, + Normal, + binary_surprise_finite_precision, + gaussian_predictive_distribution, +) + + +class TestMath(TestCase): + def test_multivariate_normal(self): + + ss = MultivariateNormal.sufficient_statistics(jnp.array([1.0, 2.0])) + assert jnp.isclose( + ss, jnp.array([1.0, 2.0, 1.0, 2.0, 4.0], dtype="float32") + ).all() + + bm = MultivariateNormal.base_measure(2) + assert bm == 0.15915494309189535 + + def test_normal(self): + + ss = Normal.sufficient_statistics(jnp.array(1.0)) + assert jnp.isclose(ss, jnp.array([1.0, 1.0], dtype="float32")).all() + + bm = Normal.base_measure() + assert bm == 0.3989423 + + ess = Normal.expected_sufficient_statistics(mu=0.0, sigma=1.0) + assert jnp.isclose(ess, jnp.array([0.0, 1.0], dtype="float32")).all() + + def test_gaussian_predictive_distribution(self): + + pdf = gaussian_predictive_distribution(x=1.5, xi=[0.0, 1 / 8], nu=5.0) + assert jnp.isclose(pdf, jnp.array(0.00845728, dtype="float32")) + + def test_binary_surprise_finite_precision(self): + + surprise = binary_surprise_finite_precision( + value=1.0, + expected_mean=0.0, + expected_precision=1.0, + eta0=0.0, + eta1=1.0, + ) + assert surprise == 1.4189385 + + +if __name__ == "__main__": + unittest.main(argv=["first-arg-is-ignored"], exit=False) diff --git a/tests/test_updates/prediction_errors/inputs/test_prediction_errors.py b/tests/test_updates/prediction_errors/inputs/test_prediction_errors.py new file mode 100644 index 000000000..f8d4f63cb --- /dev/null +++ b/tests/test_updates/prediction_errors/inputs/test_prediction_errors.py @@ -0,0 +1,32 @@ +# Author: Nicolas Legrand + +import unittest +from unittest import TestCase + +from pyhgf.model import Network +from pyhgf.updates.prediction_error.inputs.generic import generic_input_prediction_error + + +class TestPredictionErrors(TestCase): + def test_generic_input(self): + """Test the generic input nodes""" + + ############################################### + # one value parent with one volatility parent # + ############################################### + network = Network().add_nodes(kind="generic-input").add_nodes(value_children=0) + + attributes, (_, edges), _ = network.get_network() + + attributes = generic_input_prediction_error( + attributes=attributes, + time_step=1.0, + edges=edges, + node_idx=0, + value=10.0, + observed=True, + ) + + +if __name__ == "__main__": + unittest.main(argv=["first-arg-is-ignored"], exit=False) diff --git a/tests/test_updates/prediction_errors/nodes/test_dirichlet.py b/tests/test_updates/prediction_errors/nodes/test_dirichlet.py new file mode 100644 index 000000000..4f7b1d655 --- /dev/null +++ b/tests/test_updates/prediction_errors/nodes/test_dirichlet.py @@ -0,0 +1,51 @@ +# Author: Nicolas Legrand + +import unittest +from unittest import TestCase + +import jax.numpy as jnp + +from pyhgf.model import Network +from pyhgf.updates.prediction_error.nodes.dirichlet import ( + dirichlet_node_prediction_error, + get_candidate, +) + + +class TestDirichletNode(TestCase): + def test_get_candidate(self): + mean, precision = get_candidate( + value=5.0, + sensory_precision=1.0, + expected_mean=jnp.array([0.0, -5.0]), + expected_sigma=jnp.array([1.0, 3.0]), + ) + + assert jnp.isclose(mean, 5.026636) + assert jnp.isclose(precision, 1.2752448) + + def test_dirichlet_node_prediction_error(self): + + network = ( + Network() + .add_nodes(kind="generic-input") + .add_nodes(kind="DP-state", value_children=0) + .add_nodes( + kind="ef-normal", + n_nodes=2, + value_children=1, + xis=jnp.array([0.0, 1 / 8]), + nus=15.0, + ) + ) + + attributes, (_, edges), _ = network.get_network() + dirichlet_node_prediction_error( + edges=edges, + attributes=attributes, + node_idx=1, + ) + + +if __name__ == "__main__": + unittest.main(argv=["first-arg-is-ignored"], exit=False) diff --git a/tests/test_networks.py b/tests/test_utils.py similarity index 78% rename from tests/test_networks.py rename to tests/test_utils.py index f37801dd2..6d99e7434 100644 --- a/tests/test_networks.py +++ b/tests/test_utils.py @@ -5,6 +5,7 @@ import jax.numpy as jnp +from pyhgf.model import Network from pyhgf.typing import AdjacencyLists, Inputs from pyhgf.updates.posterior.continuous import ( continuous_node_update, @@ -16,7 +17,7 @@ from pyhgf.utils import beliefs_propagation, list_branches -class TestNetworks(TestCase): +class TestUtils(TestCase): def test_beliefs_propagation(self): """Test the loop_inputs function""" @@ -109,7 +110,7 @@ def test_beliefs_propagation(self): assert new_attributes[2]["precision"] == 1.5 def test_find_branch(self): - """Test the find_branch function""" + """Test the find_branch function.""" edges = ( AdjacencyLists(0, (1,), None, None, None), AdjacencyLists(2, None, (2,), (0,), None), @@ -120,31 +121,37 @@ def test_find_branch(self): branch_list = list_branches([0], edges, branch_list=[]) assert branch_list == [0, 1, 2] - def test_trim_sequence(self): - """Test the trim_sequence function""" - # TODO: need to rewrite the trim sequence method - # edges = ( - # Indexes((1,), None, None, None), - # Indexes(None, (2,), (0,), None), - # Indexes(None, None, None, (1,)), - # Indexes((4,), None, None, None), - # Indexes(None, None, (3,), None), - # ) - # update_sequence = ( - # (0, continuous_input_prediction_error), - # (1, continuous_node_prediction_error), - # (2, continuous_node_prediction_error), - # (3, continuous_node_prediction_error), - # (4, continuous_node_prediction_error), - # ) - # new_sequence = trim_sequence( - # exclude_node_idxs=[0], - # update_sequence=update_sequence, - # edges=edges, - # ) - # assert len(new_sequence) == 2 - # assert new_sequence[0][0] == 3 - # assert new_sequence[1][0] == 4 + def test_set_update_sequence(self): + """Test the set_update_sequence function.""" + + # a standard binary HGF + network1 = ( + Network() + .add_nodes(kind="binary-input") + .add_nodes(kind="binary-state", value_children=0) + .add_nodes(value_children=1) + .set_update_sequence() + ) + assert len(network1.update_sequence) == 6 + + # a standard continuous HGF + network2 = ( + Network() + .add_nodes(kind="continuous-input") + .add_nodes(value_children=0) + .add_nodes(volatility_children=1) + .set_update_sequence(update_type="standard") + ) + assert len(network2.update_sequence) == 6 + + # a generic input with a normal-EF node + network3 = ( + Network() + .add_nodes(kind="generic-input") + .add_nodes(kind="ef-normal") + .set_update_sequence() + ) + assert len(network3.update_sequence) == 2 if __name__ == "__main__":