{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Quadrature integration\n", "\n", "Quadrature methods, or numerical integration, is broad class of\n", "algorithm for performing integration of any function $g$ that are\n", "defined without requiring an analytical definition. In the scope of\n", "`chaospy` we limit this scope to focus on methods that can be reduced to\n", "the following approximation:\n", "\n", "$$\\int p(q) g(q) dx \\approx \\sum_{n=1}^N W_n g(Q_n)$$\n", "\n", "Here $p(q)$ is an weight function, which is assumed to be an probability\n", "density function, and $W_n$ and $Q_n$ are respectively quadrature\n", "weights and abscissas used to define the approximation.\n", "\n", "The simplest application of such an approximation is Monte Carlo\n", "integration. In Monte Carlo you only need to select $W_n=1/N$ for all\n", "$n$ and $Q_n$ to be independent identical distributed samples drawn from\n", "the distribution of $p(q)$. For example:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:15.986768Z", "iopub.status.busy": "2021-05-18T10:56:15.986472Z", "iopub.status.idle": "2021-05-18T10:56:17.939618Z", "shell.execute_reply": "2021-05-18T10:56:17.939847Z" } }, "outputs": [ { "data": { "text/plain": [ "J(Normal(mu=1.5, sigma=0.2), Uniform(lower=0.1, upper=0.2))" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy\n", "import chaospy\n", "\n", "from problem_formulation import joint\n", "\n", "joint" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:17.941828Z", "iopub.status.busy": "2021-05-18T10:56:17.941520Z", "iopub.status.idle": "2021-05-18T10:56:17.949339Z", "shell.execute_reply": "2021-05-18T10:56:17.949075Z" } }, "outputs": [], "source": [ "nodes = joint.sample(500, seed=1234)\n", "weights = numpy.repeat(1/500, 500)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:17.951366Z", "iopub.status.busy": "2021-05-18T10:56:17.951054Z", "iopub.status.idle": "2021-05-18T10:56:18.021780Z", "shell.execute_reply": "2021-05-18T10:56:18.022011Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD0CAYAAABgk2Y8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABd4UlEQVR4nO29fZwU5ZkufFVVd9d0Tc84MgILUUQzWL9dZBgZP45xNxsMkRDPy+vCEpkNEPdENK67QFyMa/TkTHYxBiEskNcliPuGA3og6+K6nI0uhqjJRn01onzI7j44EYIGDuAgMj3dXd1d3e8f3VVUV9dTX139NfNcv19+kZ76eJ76uOt+rvu675vL5/NgYGBgYBhZ4Os9AAYGBgaG4MGMOwMDA8MIBDPuDAwMDCMQzLgzMDAwjEAw487AwMAwAsGMOwMDA8MIRKjeAwCAM2eG6qrHvOiiKD75JFnPIVQMNofGAJtD/dHs4wfcz2Hs2DaO9jfmuQPgOOr1aRqwOTQG2Bzqj2YfPxDMHJhxZ2BgYBiBYMadgYGBYQSCGXcGBgaGEQhm3BkYGBhGIJhxZ2BgYBiBcCWFlGV5FoB5AE4DyBNCvmP6+wMAfgfASQDXAvg2IeQ/i39bBOAaACqAXxNCNgc3/MZFJCKAC4cQi4YRT2aQz2SRTqv1HlZNUe1rMBKu8UiYA0NjwtFzl2VZAvBDAN8ghPQD6JZl+fOmzWIA7iOEPAZgF4A1xX0vBbASwEpCyDcB3CnL8pQAx9+QiEQEKByPu596G1c99ALufuptKByPSESo99ACQSQiQGwV0XlJDGKraDkvq2uQD4cgxlps9/Myhma/xtWYg5t7wzA64IaWuRHAbwghSvHfrwK41bgBIeS/E0K0RCQeQLz437MB7DP87XUAcyobcuODC4ewfOd+vP7+ILK5PF5/fxDLd+4HF26InLGK4NYgma/B2DYRcSWLu7fvs93PrXEK+hrXwyhWYw5294YZ/tEFN0/ROABDhn+fL/5WBlmWIwC+CuBeL/tedFG0rokHgsCjo0MK7nghHr86drbkt18dO4tYNIxo2N0Lpag5SGIICSULUXD+BvM8h2hbi6d9/EBRc1i+bR9ef38QAHSDtGVJLzokUd/OfA3undmF+585aLsfz3NI8zx2vHYMs6+egK5xMQxzXMlxaccHvF9jI86lMlj21Nv41bGzuG7yGGzs60FHR/l5neDlWQp6Dk73xu0cg34fao1mHz8QzBzcGPfTANoM/24v/laComHfBOAhQsivDft2mfYdMO9b71Thjg4J584lAjue2Criuslj9JcMAK6bPAbxZAbKsGKz5wXva/nO/fpLuGFhD8R8jsrFRiIC0jyPZTvc7+MXnZfELA2SJIYw+FFc/818DbrGOe8XbWvBjjeO47ZrLsUDuw5emEtfD8Rc6VwqucZmiK0ilu3YX2IUl+3Yj82LZng+lpdnKcg5APb3Jp7MuJ5j0O9DrdHs4wfcz2Hs2Dbq39y4d68DuFyWZe0TfxOAn8iyPEaW5XZA5+U3A1hHCNkny/L84rZ7APTKsqy55TcCeMHFOWuOIJer+UwWGxb24MYrOxHiOdx4ZSc2LOxBPpN13NfPUp0Lh/QXt9o0UEJRcd3kMSW/XTd5DBJK6UfEfA0+OJuw3C+ezOj/lsQQZl89AQ/sOlg6lx3lc6nkGpsRi4apHjQNQVAcQc4BAOLJDPUa+5kjQ3PD8e0nhCRkWb4HwEZZls8AOEgI+Zksy48BOAvgewCeAnA1gCtkWQaAVgC7CCEfyrK8FsDfyrKsAniSEPJetSbjB5GIUFiumr3eCHx7vem0CjECbF40w7MKwu4lpHlzfvbxCw55rFnQjfufueBZr1nQDQ6ltd/M1yCpFAyZeUViNGQJJUv18M1zqeQam6EZRSsP2oxIREBIDCORUUsoDu2Z8YIg5wBc+FhYXeN4cU5u5sgwMsA1QoPselaFFFtF3P3U2yUP/Y1XdvpaktdrPLWcQyQiIB8O4Vwig8vGSPjgbAIdUhicC6PkJPvr6JCQzOZw9/Z9gc3FjdTQigrbtGgGeI5DrCWk7wcACscjkVbx4LOHLMcYDQt1pQRo8/VC9zU7rdHs4wc80TLUYOWoN+6dl8Rw1UMvIJu7MIQQz+HII3NKOORaodE5dwCISmHkeQGtLSEMp7LgciqSCe8eoNkQRcMClHQGiRywPIC5eLmWxrEklSyG02rZfq0RAXdu24en7rwB8sPWz4yazTWsYXGrqW9249js4weYcQ8Ejea5A/4SWzo6JCQzquflvddz+fn4uD3Oxr4eRHI5AAgkscfvvaXtt2XJtfjtuSRiooCVBuWP8bhuPfdaJS/5fZaa2Tg2+/iBYIz7qC8/kM9ksbEvuKBWEEinVSjDCgY/ikMZVly/9F738ZNEE5Q22+o4y4qBU7/zN8NvEJG2XzQioH/3YQAcNlYQCK1VAtZISPRi8I9Rb9zTaRUdLWFsXjQDRx6Zg82LZlSNzmg0+DHUQakuaqHesFOP+Nlv4HQcr78/iJXPHEAqm8Oj86YVnpnFvRDzhRWHouZqnoBFg5fzGNU/ippjH4ARgFFv3DUE4Sk2G/wYWJrhG05lPckC/RpeL/ArNbTab/X8bjz+ciFF41fHzuJTF0cxRgpj6HwSSjwFoBBsXbrNPgMXqM2Hzct5zB7+0m37mIc/AsCM+wiG0RuT2lrK6rr4MbA0g/mjV496WvpbHWdjn7XhjUQESG0tiLZFC+OPtbgyPOm0CjGfc1yVmTXrAEr2e3TeNKx9kWD3gRMl18joCHjxkqnXPZUNNN/C7f0dyeUyRjNGfUAVGJkBGGPAcny7iJWz5RJt+oaFPZAEIKHCc3DUGKQbTmXxo1ePYt3eC+kLbgPSVmoZ833QpJdxJVs2/iDoM6cAsdsAshfVleUx+3ogCjy+btbOe5ij+Xry+Zyr+9toirFKMRLfZxqYWsYBI/FhMCo+9qz4LPp3H7ZUd+Qz2YpUG0EaBqv7ILaKOJvIUHXl2gfEyrDlON5xXmKsxVFX70Zx4lWZU3LMVBZqPo9PirkDA6fjePzlAZwZUlyrtmgfIUmA43VoibXgLotr8MTiXqSKlFMzYSS+zzQwtcwohJFvtcv6rFSZUm3uPBYN47Ixki13bOaMt752DAkVjiqRSERArCXkyEu7uUZe+X3tmPFkBltfPYqsmsODzx6C/PAL6N99GCtvkTG+XXTNw9OolRzHO45dEgWsnt9dFmOQRMa5NzOYcR9BMPK18VRWN7oDp+NVM8BB10cxI57MONalMRu22VdPcMUhc+EQjg9Sjp3KeuK8NX5/y5JeT6qrWDSM2VdPKKsN9MCug1gx66qye0SraVNJkHY4lcVz73yI/rlTQVbNQf/cqXjunQ8xnKqfHJihcjDjPgKg1ccxeqq5fF43upteGcCaBd1VMcBug5Z+kc9k0SGFbcdvNmx2KxUjYtEw1u89gvUWypitrx71HFBMp1WIAo+h84Uqp23tUcfAaDyZoY53UqdUco/sdOuVrKDymSwWXj8J/bsP6yuHhddPqmuuB0PlaNpwOGtPdgFcOIRlBr739fcHcc9Tb+PJJb0lhbs2L+4tqZUS1PVKp1WgSDMEjXRaRQTAxdEwtiy5FpIoIJ7KIp8uGC2xVdQrVWrz11YqTkWy4skMTp1X0NkaQf/cqegaF8PA6TjWvkjw/KGT+PPPT/E1J4XjsdyiqBiNzhnmOOvxpkrvERcOYbnpPi/fuV+PnTgVZqNBL2BWfD4SigoOeWSVnOe5jwSMFNvSlMZdDx65fIEaEX7S/mnbUzMqi7XSjQZKaQDxg9e5p9Oq/vfk0IVjaM/A+HaxpFLlnndPujJ0mkH88OOkZcDZD22lqDmdEgJKDTAs5phOq4iGeGzo6ymrp6N9wDTYUS+Dw0rFFSZTGRV3b9/XtO9UEBgJtkVDUxp3Ow/G6gVqNHh9gJy291Kytt6wm4sXmJ+BXB54dN40TOqUdLWMk6HTPNbWWAQb+3pKC6/19YDPe/dcJTGE8e0i9qz4rL4S2PTKQEnJYvPHTc1kIcLZMDvd50pWUM3+TgWFkXQdmpJzb/bGA16TRpy2D6o+Ti16bAaVMGN+BnYfOIFZ634OoJBtnEwUDKTWqIILhyznk06rSAylEOULtIQWUNz5xnEkVHi+Bqm0ipWz5RL+euVsGUmlcC9ovLk2bq+KnE2LZgAcV/E9a/Z3KiiMpOvQlJ57M3mqVvDaXMNp+3RaRUeHWNGS3MtqwpzEJIkChlNZV+e0m4uade8pOz0DXldHOY4v07u//v5ZTx5bJCJAzefLesXe/8xBPLGkV79ufj1DqwYoVqWJ7VaANDqs2d+poDCSrkNTeu7Vlt9VG16VDW63r0Sv7tajNnued23fh99+nMLW145B4XhEpbCt9x+ULt7pGfC6QrD76LjxhrXrIkWsdfNSJASlmExUiWdo1NyrubynwmB2FSKb/Z0KCiPpOjSlca+2/K7a8PoA1eKBc2t0rIzmA7sOYvbVE7DzzeNI5OyTh5zm4pYaMj4DZNUX8cTiXlzSJur0i1cjSvvoHB9MuKqVo10XWk7BwOk4lu/cj2FD/oHx7348Q6s5jm8XwfN8yfWLRATwkTA6YyL6507Fl6ZNKPsQeH2nakHh1QPNbluMaEpaBqiu/C4o0JbBXntnBt1r0wpul6M0o9k1LgZcPQHLd9grRezmIkneZISaTHIwmyujJkJK1nI+CUWF2CqWXT8rKeHq+d1Y+yLRywDY0SbadXn85QGsnt+NB3YdLDvOr46dhSQKviWLZpjv2dzpE7Fytoyl294qaReYzuaw3KCCWT2/GwDw/KGTJVSg9k5FwwVD3dYetXzWRpKixArNYFvcoGmNe6PD6QXw+gBV+4Fzq5OmfQQGTsc9Nbe2motXGSFgo25Y3Fsyn2U3d+GrN12BaETAB2cVdEhhRHDBGGkfnS1LrkU0Iuh6990HTiDEc47NxrXrolWO3LRoBtpawiXHufHKTgwrKqSIgE2LZqA9Gsb5ZAYhnoPqQ1Nuvmf3feGqMr7/nKkuj7bS6p87FWeGFMsMWKeG8SNJUTKSwYx7ldBsL4Db1QHNw33unQ8x4aYrqMk4kYjg6NVJIr3OC82wUumXlhDy6Qw2L5qB1pYQBuNpfN3gva5Z0I2Li7V1zAHiRU++YbuCsVqRGa/L84dOomtsq571+atjZ3X6KcxzGEykSypcGsfiBeZ7ps3dCFpdninjYpYfb6uEOPNz61UQwFAfMONeJTTjC+BmdWA2KJpa5o7PTAafU8uScbRU/q/9/hWOiUsJCpVix0cbVxJzp0/EvTO70DUuptdFUYYV8DxftiK4/5mD2LLk2rIV1rKbu2xXMNQVWT6n15aRxBBVa6/mYamm2bLk2orvmdgqll0/rS6P+ZoOp7OWXLKb59avomSkZH42C5oyoNoMqEWnoXrBqNhIxVM4Ozisa8u1VH5NL772RYL3PxpGIqNi62vH8N6pOFpbwuAiYUSl0uCmKPCeA8eax3zfrClYecsFffld2y90E5JEgdITlS8LNK7b+x52vnkcTyy2LgBmp8LRastoiqVkIlOmYKKNJYgKjFbB6oulMNYumF5e8ZGyknLz3PoJ8LN+rrUHq+eO6tR/dtvkISg0Sg1rq7rmr6z8HHbv/y1uu+bSkkDjhr4eiLkL16OjQ0IioXj27qJSGHlBwF3brOuyg+PKNOz3zZqChTdMKltlaHVlaPXonerXO90Hav34xb16u75KYPaOBZ7D3//yKGZfPUHPmN3z7knc8ZnJ1DrzaZ4v59xNz61XL9xrvftK0CjvQiUIop47o2WqhFooXBoRVpz8pE4Js6+egAd2ldIRy3eUxyD8BI5zHI9Wir48Fi30OTWP6Y6brigxsk6BRg1eKAlLbj6dsaZ90sGs6MzXLxIRsPD6Sa7VOW4T4rzeJzu6p9KGMQzWYMa9ihgpkiqgaKgiYVPVwAy1XovROAynsq6VNH4Qi4bx3il6JUirMdGadHQZAo1OgVM7Y2nF5d9x0xWItYQQqmKFTjP8OhnKsBLoc0v7KCaV7IiWVdYTjHNncITOl27fh6seegFLt72Fj5MZ5G3qtWhccyGImEdcoSTvpLJ6VqsQ4n0lxMSTGex592RZN6ENhobb5m5KNG45oRQCjQCoNWDcJLkYufkvTZuA2665VL9+d27bh1RGxdD5pK9sYjOcEooq7bYVBGg8vZqH71pDIzWRKigwzh3Nx9FZeZSSJFZtDjS+9NF50zBGsve8tX1nTx2POdMmYIVJQnn4xDn0Th5TVu7WS2wiKoWRUIGdbx7XueVhJYsQl8fwkPXYnGIifjli7VkycvN2PWwr9Y6rEdup1vsQiQgIiWHkwel1+VtFAfLD/+q5B6/dvKv5LtQKrIfqKIRdVcFqgcaXXjZGcqyJou3b/7//HY/85D90Jc2WJddi7YsEV45t07NajZ4bL4YLHlmsxdEjy3F8iWEfOB0v9CW1yQsqSzNf3IvWiKB3T2p10VvVDsaVgdvOUH4QVJXNWmE4rWLptrcKz+72fRgcTmPZzV0l27hRlTXbvOsBZtybCHY1QhS1el1zaBTGB2cTji+hcd/dB05g9vpfYNGTb+CTZBr3zuzClPEx9M+dirnTJ+r7aIW2NAPgJJmLRcPY+NIAZq//BT79recxe/0v8P5Hw2U1VszQ6Iqh80mkMiru3LZP/2BSjY7LvqJGGqKaPWybqUStpUHesR933HSF57pJzTTveoEZ9yaB5rHftX3fhTrht8iYO31iUSddPY/Fii9ds6AbHVK4rMenmQO12vcHf9IDgEP/7sO46qHSuQAXyhm49cjMHx9jjRU3mmq3RmfNgm7k8nlX3K5xZdA1rhUbAqi3r8F4nYdTWV+eb7Vgx4PbZRN7LdQ1kvNIgoIriyDL8iwA8wCcBpAnhHzHYpvbAXwXwHJCyL8Yfr8fwGQAHwGYAuBrhJBk5UMfXbAqZ2CU7iWUYCpEWvH5ANASFvD00huQUFTwHJDPl6plnDI3dbVGKguOA/7if9FliJreXIOTqsZNjRVjCr15jnZG59F503DZGAkDp+N47F/dFRHTUJI9GhECkcXSrjMAbHxpoKJCZJWiko5hXtU5tqqlMAusAi6MuyzLEoAfAphKCFFkWd4ly/LnCSE/M2xzBQqG/wPTvr8D4EEAlxBCcrIs/zMKH4mng5zEaIBdNcYNC3sgCjwq/WLSXs5IiC+py6IZbaAQMNXKECw36cY1g6oMKyWS0M5LrDnoKeNj2LRoBra9dkwvwAU4e2RuaqwYNdXmOW5e3EutIDlr3c+RzeVLShskix8IL8Y5KFksrWbRE4t78eefn1JXnbhTPaVKmnib4VRdlMEdLXMjgN8QQrSn8lUAtxo3IIQcJYS8bLFvAkAaQHvx3zEAh32OdVTDjXSvUtCCVOcSmfLAVSRcEti1K/rlZS4AMK/3Us8UhlHuZ7dkN89xbJuIZEa1bFPIIY/rJo8p0DyG0gZLt71VRvNodIRfOadb0D7yrS0hS6ljLeWCTjx40LXSG0Hi2chwQ8uMAzBk+Pf54m+OIIScL9IyP5Zl+SSADwEMmLe76KIoOI6q6Kk6BIFHR0fjf+7NTZw39vWgJcQDITGQOQghnqqKMf8WawmVZHhqQcNy7zdrOS6ruUTDAnK5PCICjy1LroUkCoWPl1DwQbx4ZLTjG+eoGe0VO/djfLuoN9g2nnNjXw+GFbU8u3bnfmxZ0osOSQSAQplcw2pgY18POjpE9wN2CVpxNdp19jMuv8+Sl7Gp2VyhbnxYCNzTbpb32Q5BzMGNcT8NoM3w7/bib46QZbkHwP0AZhBCsrIsfx/AtwF807jdJ5/Ul4JvFp17xIK31cYdxBysqgpqqhgjNMrC+CGwalKxYWEPVCWDc5R+nlqWppbxmsvlS+aQLLoUXp4OjU8f0xrBE4t70WrIAj13LlEyx3tndpUY7ef2n8B9s6bgT2+6AkKIRzyZQZQHOjuty+ZKYsFbFltFLDM1KVm2Y39V6qZEItbNPqyus99x+X2WvIytmmiW99kOHnTu1L+5Me6vA7hclmWxSM3cBODvZFkeAyBLCDlvs++nAJwlhGhr6pMAJrk4J4MFqlXOQDOIrS2hspK9Gud+45WdJd4fzwHLbu7Cur3vASjIHLvGtpYZVLulciqj4m4Dl7+xr8czl22eBy2oqx3TyPua9edzp0/EbddcirtM8QWu6JGObRN13v2DswkkizRSLcs7eyknYB6XFjfQaJKgufnRWk+pUeEqQ1WW5S8A+GMAZwBkCCHfkWX5MRQM9/dkWeYAPATgawB+CeApQsgeWZYFABsBpACcA3A1gBWEkJPG47MMVWc4VeGrxNsyZvoZa6AY1TJaXZnjgwms33sEp84r2LCwBzvfPF6i0nDLoVajSqDbY2rXkud5LN32lr49LZP0iSW9yOXyGFKyJU02tPly4VDNKh56gfF6aBSUeWVldb+a4X2wQ7OPHwgmQ5WVH0DjPwxuUsz9zsGtQaRt58VTN8KpdK4f2B1z6HyybGzm60pWzYH8cPn+ZNUcnDiXxDf/8aDlddIVODtLOX4pLFgWV6sVjPPrnzvVdQmERn8fnNDs4wdYyd9Rg2q27HNLKTipNLzCbzcfI8yrmXjKOqB3fDABKSKUVRq06ipltf+Jc0l86uIo9ToNDiuF4xRjCMcHE3jkJ/+hr27qVeHQPL9m6wzGUBlYhmoToJqp1nayQXMmZJAZgVaZqxv73GuerWrs5PL5smOunt+NdT89Qs1yNcrpcsVa68b91y6YjpYQj+ODCdv5p9MqkM/jK1vewOfWvoLn9p9oiHon2vxYRufoA6Nl0PjLODfUSVCcu0b5SAKQUFHCxVs1fahEp2z2vKNhwfUcaNfkySW9yINDNCJg4HQcj788gN0HTrimfMzNshMZFSt27sfYNtGRs64G1RQUvFSPbPT3wQnNPn6A0TKjBkFm9plBUzjkuBCW77xgPDVVjF+OnXZuo/on6kHXS1vNRIvNqRc9+ZYvysc4pkhEwNj2Ah2jGez+uVP1LFVVSZfMPwiqqVpgSpbRB0bLNAGCzuyzOr4508/KeG58aYCaCVlr2NEMfho4WyGdVnUeHyitapnLlV//fCZrmelajzovVmAZnaMLzHNvEtS6ZV8jeaFe293RvFTgQi0ct54rn1Mtz8NblHxw23/U61yZEWbwA8a5Y3RxdIA7A1KNDj9OYxLEMCSxlPKxGweAC/NIZSFw0GkZ85z8zkdsFbH1tWMljUD2vHsSd3xmsuWHtpJnqdbX3Hhev7GPRsRoep+Zzt0Bo+lh8GJAauVF0sbUGhEQFQvSwnU/PaJXijQHk93MyW/SlNcgaSXPUjUSu5xgde029vUgYkE7uT1evVceo+l9Zm32mgC1qt5n157MPAYANeFoaWM6E0/jqodewIPPHipp5mGWgbppueZWTmq+BklaY+8q0FP16C5kde2W7fAn36S1gGSNq+sDZtwbALV8KeySkbhIGK0tYbx3Ko6trx2r2Ytp16NVMzgP7DqIe2cWOg6Zjasbo+hG5211H4bTKjYtmlGVLkrmj3g9tOhBflBYX9PGAjPuDYBavhRWBmTZzV0YHE7jbkMLv9uuuRQ73zxekxeTZtQGTl+gPbTGJFbG1Y1RdKOgod0HnuM8KZVoxtvpI87nc4G243ODID8orK9pY4EZ9wZALV8KPp8rM3J33HQFlu8oNWoP7DqI2VdPqMmLSevR+vjLF0r/Xzd5DJJp1dK40gw3n8/phpYLhyAJsDXSdu323NBTkYiAc6kM1Xg7UWIJFdj5xnH0z50KsmoONi/uhSRUt3SB5bXrs1YDOcHth6KWDURGM9h6qQFQS9lhjuOx881jejLOwOk4YpQuSl3jYjWRPmrSxS1LegtqmVQWuXweZ4YUhHjOUBc8bWnorKSPfD5XzLA1l//NYpASnKz0PnDhEJbZ1ACKRcMY3y5iz4rP6td+0ysDhTEDev0gLWFMC6ZWE4W2dDw2L+5FqxjCwOk4dr5xHAuvnwTRY/llN8l2Tn1WGYIDU8ug/tH1ICRwbudgpf6glbrdvLgX+XRtqhqapZB8Poccx/tWXfhRnkSlMHK8oBu5Pe+eLBg5l/fBSVkjtbXg42SmpGzwmgXduDgaRlQMuVblBK1ICVKl4zS2WiiC6v0+BwGmlhkhqHYGqhFWS+c9754s53r7esDn1JoZdoXjsXTbPp3OSKgFT3DofKEPU1t71NMS3ivVpdEixrjDwhsmeaJFnGiJXB64/5mDJbTM/c8cRC7vjdIIOvgeJC3olAXLePnagRn3BkGtUsOtONaF10+CxJv46FwOyURtslGpXLSpCbcXQ+Y1UGg5hh37kePcvyJO5QdaW6zpr9aWkOuSCdUIvtdSpcOqU9YOjHMfZaCl5icThY9JPWp72wUyjU24vdSx91pszcmjdEN9OJUfSCgqpYG06rqwVzXqslezMF09zzXawYz7KETQdWoq5YBpgUxzE27A2ZAZxxJSsnoDDadx0cbw3qk4+ncf9hT0y2eyiAN6oDSCwn4c8lizoLuMc+dQ4Nnd3JdqBN/NH5aEkoVapQ5SrDpl7cBoGYaKEAQHTKMkOOQ9LeHNY7lz2z6kMiqGzicdqS6rMayeX5BjGmkiN6Bdj6ySQUwM4dF500BWzcGj86YhJoaQVdwb5qAqXpphpAVFga+qsWXVKWsDppbB6IquB42geqtaFQ4DgHw4hHOJDC4bI+GDswl0SGFwlONVqsQwev3vnbrQ6AMo78WqbdvaEkJCUfW5hgQeX/ufb1HHEITSpdr1W5r9fWj28QOsWQeDBSp98b3uT+OApaK0z62OOZ1W0SGJJbK/SERAOpvDg88eKpWIehyLWz5ao0UAlElDtV6sY6QwIih45ztfO4bbrrm0pDvTxr4ejG8vHaFxDEFQYrUu/8zQnGC0TBPBKbOvUorEz/52pQMqVXN4VYYEpcSwUr2snt+N9XuPIBYN6+OaffUEPLCrVNq4bMd+rJh1VcVjYGCoFMy4NwncGN5KZXJ+9rfjqjX41TF71URXykdrH8+29ihaIyGs+/J0kFVz0D93Kta+SHDqfKHRtDaurnExy/FN6pQathsTw+gBo2WaBFw4pKenA9aywEppCT/7m9UPw6ksfvTqUZ2rBvx7rl6VIZUoMWhp8Y+/9B42vjRQItmL48LqxGp8v/04qZd3SChZ5GqU5VsJjHRcQskWKLEGHzODPZjn3iSopKztcCpbUeKP0/5G9UMuncHC6ycF4rn68cT9KjG4SNhy1fKnN11RljWsjWvPuyexen532aplzR6C2et/AfnhF9DaEmp4I2leFS7dto/VYR8BYGoZNEd03UkJ0tEhIZFQymrUrJ7fjefe+dBVjRSrGjde9jcex09Q1+o+1KKzTyQioK096qnjkpVaxqljVDUQxPWpRweoaqIZ3mcnMLXMCIXXhtAaNFriicWF6ooDp+NY+yLB7gMn8Pr7Zx0zO93sHwEcjUmQao5aKEO4cMEw0xKpxFaxbJ7mceXSAqSIUFbJsppce1AVFquR9cpQfzDj3mCgvrD5nF5czMmwdlp4oW5fVrv9W1tCGMzmRly51lg0jG/vPozV87tLZI1rFnTj4ecO4dR5xXGeZWWLA1pl2HnmbuIwblDLktMMtQPj3D2imo0GIhEBfCSMzpiI/rlT8aVpE0oUKxrfqyk2tCYPZlQqCaTtn1DUEdlGLZ7M4NR5BWtfJOifOxVHHilkjz72rwTP7T/hep7ptApR4APLvHRSSAVVYbFaWa8M9QUz7h5QzV6n2rHvKpac3fPuSfzN/301fv3dL6F/7lS0toRcn7/Sl5W2P62qoR+ZYyQiQIy1oPOSGKJtUUhtLZ6P4elcNh9kbb5nhhTcuvHfkM8D6/cewb0zu/Dr734Je1Z8FuPbxZqXpXWSpgal6zeXnN6ypLdqJacZagcWUIX74IWfwJPbgJfx2HOnT8TKW+QSimBDXw9awwLu3LbP8vzRsFAyh2pkqnKRcEmVRv38i3uhxFOu524VuF2zoBttYggwbON1/Fb7AHBshBKJCOAiYcSKwdEQD3w0nLZsqpEYKp+nEUEG85yafwTR5MUKzR6QbPbxA6xZR83hpwGEW0/feOx7Z3aVZT4u37EfeXCuz19pcSar/QUOWLOgVPq3ZkE3BIvHy27uVh7p/c8cxMeJjE41eV0hWe7D8+ApEkc+EkYkIlzYb/u+ogzwLZxXsnh234dl41Nr7II4eea1bPLC0HxwZdxlWZ4ly/LfybLcL8vy/6Bsc7ssy7+WZfm/mn6Xi/s9IMvy87IsXx/EwOuBQBpAULhb47FpmY+SKNS10UFUDGHtHqI3cO6fOxVr9xBExfL52M2d9pG8bIxUkt7vhdunNdugUUmSGILC8QiJFsZ/R6G0gHmfWEttYwt29Joxm5bjODTCCpyhseBo3GVZlgD8EMA3CCH9ALplWf68aZsrAJwG8IHpdwHAOgB/TQhZDeBrAI4GM/TawyuX7cXTNx5by3w04rrJYxBP1TfwpQUeZ6//BT79recxe/0v9JR8M+zmTvtIfnA2UZLeb7WvBjOPTjPiWoMMI5bd3IWhVAadMRE5wLLQV9e4WNn43H5Egwq6p9MqJKEgTT3yyJyCRLV4KOMq5a7t+/Dbj1PY+toxlnzEoMON534jgN8QQjRS+VUAtxo3IIQcJYS8bLHvdQA4AH8hy/KDAP4vAB9VMN66wusy2Iunbzx217jW8p6mC3uQT2fqugz38nGzm7vxOLf1TMQrKz+Hp5fegDGtEfD5nK1aR2wVEZXK2+8NDqex7Oausn0+SaZLioDdN2sKFl4/Cfc89Tbkh1/AXdv2YeVsGXOnTyzZb1jJlsxz06IZAMe5MtgKH0zQXevreleRMrpr+z4kVOts2gd2HcTsqyeMCPUSQzBwDKjKstwH4HZCyG3Ff98J4HOEkEUW274CYC0h5F+K/74dBa9/MiHkE1mWnwKwlxCy1bhfOp3Ncxw1LlB1CAIPVc1V5djnUhks23Eh4LWxrwcdLS6bPqg5SGIICSULUbD/DldzDn7H5TT3TC6P4XTWchvzvmsWdGPtnkLxrs2Le6mB3bu37yvLrr3zD66Ems9DioQwlMrgHoug+LovT8cfPPZyyRi0eabSKnWcZiSzOcuxbVnS63gPra7zUosA+tNLb7AMtJJVcyA//AKOPDIHatb/s8DzHJIZ1fWz12io5btQLbidQzhsFfEqwM0n/jSANsO/24u/ucF5AP9JCPmk+O9fAvgcgK3GjT75JOnycNVBNaPrkYhQlnjk5VzJoeL/O2xXa4WAm3EZ555UslDzgBDiL6hvwiEs27G/JAln2Y5CEk4kdyFh6/hgAo/9K9HT+ltFiiSzJYTNi3vRWsyu1comZFJpAMBgWkVnTLTcd/xFLfjPv/kiUpkcVCWtX8vkUEHJRBunUSVVoIfo9e2tyhjYXbu29qgt1WROOtLovHgy4zuzNBIRkOb5kg+ZGwVOLcpEuMUoU8tQ/+bmk/w6gMtlWdaIyZsA/ESW5TGyLLc77PsGgM4i9w4AlwM44uKcIwajuaWYNveh80kMp1VdkaJRFXa6eW1fAJi17uclVSapMYlkBvl0BsOpDKaMj+GOz0zWjZJGeyWUrOW+xwcT+PDjwqfKnBzmNnYSi4Ztx+YWmoJHK4lgPhaHvGWZ5T3vnqw4BmP84LoNZtupm6qZ9MdgD0fjTghJALgHwEZZllcBOEgI+RmAvwLwZwAgyzIny/LDKBjv22VZnl3c9yyABwCsl2X52wDGAvjbqsyEoWFBU79YBTvNhtCKf9/z7knrmETRW6R9TNNpFbl0xrIRx8+PnMYlMRHRiICziQzyBgMfT1l/EOKpUiOaVLLokMJlx9/Q583gatdr3U+PlFWd3LCwB1mlNPbyxOJefOrilpKPmV/4yXqlqpsi5bERFvCtHVgSEypbxjXKcrSRl6J2yTgfDSmOCUZWiTqSAOQ43td1l9pacCaexmVjJAycjuP1X3+EWb83npq0JLW14ONkxjapyTjO8e0iVsy6CpM6JcRTWfA5FcmEe8/deL3mTp+Ie2d2oWtcDMm0ClVJU+dZiwqRVuewq6j5lS1v1LzaZCO/C27BqkLWGUFV5RvpsCtMZSyGllCyUJVMmbdt1YAjmShs48dIZJUMpIiARU++gV8dO4u99/0h7n/mYAmnfv8zB7FlybUACvr+h557V2/AMXA6jrV7CNbd3oNEMfZgLuL13P4TepDXi2E3X6/dB05g94ETulG0M+xBPItai0Ez565p663OESpSXVYVNVm1yfqhucLgDYZK29qNFthJKI00iijw4MKhMn7Wjmrxw+maJa2TOiVq0lgkImA4lcX3v9wDAPjGj/db6vupdIaPxCc/tYGCehbTaRUdLWFLuS3tHGoeluPlkK9r0t1oB7NCFYDVwbaHcQkfUrLYvLgXsRbrcriRiFAmf3TyPCvxVo312MVYi6XnmUyrheOb5JVdY1ux8PpJJcbWb9lcGpXitV1g0M+iMqyU7Wf3AcunM2XjzSo5xx4EDNUDM+4VoFnrYNequ5GV4c1T+oly4RCWeaxNHlQ983w6Y2mE1Fxe91S14z+w6yCeWNxb1hfVTTMVt9dI/zh5aFKiPYtj20Sdo//gbAJJxfr8fp4Bu+edNl63H6lGiV2NJLCAKvwHYKpVlc8P3M6hVmP2WkHTqQKiFfzsQ4PXQCGt9Z4ghl036wiyvV0kIiAfDiGuZEsCv1b31ukZoD1L1Xp2gj4uC6gWwDj3CtCMVflqFSfwKqmjlRxIKlkqpx5UPXPAmtf3enytWYdTMxUNQTXb0M7Nc5weGLa7t36fgWo97yx2VR0w414hmi1JKUiDYoQ5sJmkJAvRDKOm0jDXcxlOq1SddLU7CHk9vhY3cJvME+THCQBiLpupVPIMVON5r9YzOdrBPo2jDG7iBF75Txp3vGnRDNxj5twphrFAB4gl/Cw4rozzNnLqfgKPXuD1+LS4wZNLejGcVsuuj5S3CDj2FTpeAfA8F7cxoEaLFTXaeEYKGOeO0cXROfGb2t93vnkcs6+egK5xMQwr9ok4VO54cS+Qz7s2vOY5BMmpVwo3HzzaePd/+xYs3faWJbeu1diJRcOIp7LY+upRbHxpwBfv7Ja79su5VwuMcy8HS2Ji8Awnb5QLh7DztWO47ZpLy9r8iRHB8mWzk8hpy3c/8OLRVVNt4VZySRuvJAr02vYoXL/homFft/c9AP6UP/q9LUpOE4oKDnlklZzldk8u6UUeHCRRQDyVRT5dn0qK1V6BjVYwzn0Uwo43jUXDmH31BMs2f7QAV9DcsQa3nHc1G5cD7gN+VnGDDQt7bGvTGBtu3HbNpSV15f3yzqmMiq9seQM9f/0i7ty2j3othtMqlm57q3DNttO3qwWaLXbVDGDGnaEE8WSG2uaPZmiqFdjU1RnFTkRbllyLVgvjU221hduAnzm788klvWgJC2gVhfJCZ3092Prq0bKGG/fOvNBwxM8H0u21YAqVkQ92JxlKkM9kMcxxngJclSyrjXRKQinULzHvl8qoJU04jJRIJCKA53k8decNGDgdx+MvD2D3gROBZgp7DfgpwwrymWwJj7zs5q6SDN3WlhA2vjRQsp/W3i/Ec76zOd1mqtptV0mjD4bGAfPcRxCCqJ2dTqvgcyq1pK7dfnbLaquxmemUpRYUgp2Hqe2/dNtbkB9+Af27D2PlLYWWeUGqLYKo9bJu73u4e/s+vZHGMIWqSSjZijTkbimyalFpDI0DppbByImuf5zMBKY4CDJASVNDtEYE3GnRRs6YoWmnmIknM5YqnUfnTYMUEQJNKHN7PbRnyUnpU+9sT7vtJEls6vdhpLzPTC3DAKDQb9NOE+4VXmub2IFWA2bLkmup1IBmRJNKFstu7tJlmQOn49jz7kk9C9Rq/0mdEobOJwMNyvmt9UKjcqqlEHF7XLvtJKmiITA0CBgtM0Ig0fqKNkCWH80IS6JAbXmnqV6GMyr+2+9fATFUeFTFEI8ln5kMPp+zpRbc0EJ+4PY4bqicailE3B6XKVRGNphxHyGg9QZtBA6VaoRT5QZwzYJurPvpEZ1f3/nGcSTSKh589hDkh1/Ag88eQlzJghME31LJra8dAxcJezb0XiSXjVJ3iPUwHb1gnDu8cXSNWpo0aM49SNjxuwAMahkVDz93CM/tv9AMe8+Kz6J/9+EyXn3LkmuRHEq6uh/GDNq50ydi5S1yaYKWy+vkpopjI/G9fnn9RpqDHzT7+AHGudccjd5Wz9iyrpE+PFb8Lp/PIcddMMpD55MQxDBOnS/ltGmae0kUkBxy5sLNUsmYKGClqaWe29hEUA0xauUgBFXvnqE5wWgZD2j0xI9G5lCNY8tnskioKKM3RIEvo1mGaXRTyln/bSWVjIR4jG8XS7ZzG5sIQj5YrWxaK/ol6GqLjOJpLjDj7gGsNGkwoH0kFTVXxlPzOdWaV087G1Sr8yzbsR8rZl1Vsp1bAx1EJm41HATaB8Nr2WU/52AGvnHRGC5nk2A0liatBoVAV8+EMDiUKqNZxIjgq1WbnVTyxis7Pff19CJfpF23avTdpdIvi3sD62HKKJ7mAzPuHuCnT2YzgxZjkCQeOY73bfBpH8kEpd+nG4251Vg3L+61/hinsiWVEwUOyOTy6GyPOs7H71i02IxXB8HNx9Vr42o/H+dKPkqNKkIY6WC0jAfUUt7WCPwmjULI8UJFy3MavSEKvOWc3VwLq7FuffWoZRkFPqfqlRMffu4QBhNp3L1934X58Dyikn+qjYuEqdSLF2rHLRXipPcPIg7jN97A6Jz6gUkh0XjSKT8StmrMgZZCT1bNwae/9bz+m5+mzmZvjs/nkMwBy3aUzlkSgIQKx2tBH+sXMZzKlniNXDikSxppUsvNi3uRT5cnQzmho0OCEOIdSw+48WTdNtCuRbMLv+cIsgl4JeNvNrAG2SMUjaLKoXlrA6dLuyD5CSqbPcocx2PZjvI553nB1bWgjXU4lS3zXGPRMMa3i9iz4rOYMj6G/rlTy+qot4ohX9dbUXM4PpigFAVTIbYWlDpuvGkvpYarvaL0ew4mQqgfmHFvQDTKC2FJIfT1YM+7J0u2CyKoTJtzq8umz1Zj3bRoBsBxZXROUsli5WwZ/bsP46qHSqtJavMZOB3Xz0Gjhax+l8QQ1u89gtXzu0vGsrGvBw8/d8gTLeGFCnFLv1RC9/mheFj1yfqBBVQbELSg23DKut55tUBLPlp4/SS8/v7ZQIPK2pzHtom4d2YXusbF8MHZBFJp1VUA0jzWpJItNKU20wgRQM0D95sSmR7YdRD9c6fizJCC1fO78dw7H2JC+2TboHJCRdnv4bSKU+cVrH2RoH/uVH0e6WxOz7x1qzLxE8C3o3zqkYQ32kQIjQTGuaN+HB3tRbTiNzWDs/D6SZbL4VrOoRrqh0hEAMIhDClZ3P9Mae/WiMDjHrNBqoDrjUXDVE78vVOFqpPadTby8wAwd/pE3PeFq3BJTLRser1lSS+yuTzOJTK4bIyED84m0BoR8Dc/+Q/sPnCi5Fxumnx7udZOvLhb/jvoZ6nWahnGuRfAPPc6wcmLEiPAE4t7IYkhDJyOY+2LBLsPnMDr75/1pC2uxosVZDlg4zGjYrjMo16+o6DX9irns6O2qCsjJYsp42OYeNEV4HIqkgkVne1R/TjGujRP3XmD5fGjEQGD8TQefPaQfl+//+XpJdt5oSW8XGs7vTscNPZaoNmuI5ZfVON5YXCGK+Muy/IsAPMAnAaQJ4R8x2Kb2wF8F8ByQsi/mP4WBfAGgBcJISsrHnUDoFKj6ZQUkk4XDIvZw/SS8FLNZbiV2qUS7TtgU7a4JaTzvG7HpjWlttK4g+Pw9NIbcHwwgfV7j+DUeQUbFvbgR788io0vDRiuk1DyIbh3ZpfeOPz/fJLE3vv+EJeNkfT2fmeGFCSKVJDxvv7lPxzAo/Om4flDJ6tKS9jr3QXqRy2pZBu6ZhKDPzgGVGVZlgD8EMA3CCH9ALplWf68aZsrUDD8H1AOswrAO5UNtXEQhHbXTdC00mBUtVQ3VvNPqMDW145VpGW2K1vsJRDIhUPY+urRsqDmhr4e5PJ5XdP+4LOH8NCtv4u//+q12Pnmcazb+56tLl0rYlYIvHJ6GeL+3YfxzS/K+MGf9ECKWH+gJnVKVc+NoD0vxwcTthp7NY+GUGcxBAs3apkbAfyGEKK5Ta8CuNW4ASHkKCHkZaudZVleXNznaCUDbSQEYTTdGO5Ka5lUS3VDm//sqydUZBysCodtWNgDPp/z9DGNRcPY+NKAHtQkq+agf+5UxCIh3FNcLRlrzeTysGxWHYuGSySAyWJw996ZXVj5zIGS49z/zEEk0zkMnI5T72u1C7rlM1lsNCVtrZ7fjfV7j5TNxfihiblUJBlR6+YnDN7h5u0bB2DI8O/zxd8cIcvy7wH4XULIt2RZ7vYxvoZEEPVB3KgIKm3FVq1aOLT5d42LlfzbT70Uq7LFOS6E5Tvd1zXR5r37wAk9iHnjlZ14eqk1T651hLJriYdioHvDwh50xkTL40zsiOIbP96P1fO7y+rF10Idkk6rkNrCeHTeNJ0uWvsiwZkhpWwuxvvipySCRuOMbxexYtZVmNQpIZ7Kekr8avQS2s0ON8b9NIA2w7/bi7+5wR8BSMmy/FcAfh9ARJblFYSQ9caNLrooCo6jBn2rDkHg0dHhvnGkRh9Y1UbxchwJwJYlhaBpQslCFAoLKaselmo2h2hYAMKC5d9pc9jY11OS9bmxrwfRsICoxbaKmisbixVo8x84Hcfc6RN1KeOwx+shCDwkSSybrxDiqR/TaNja07Oat919c3udJACJtPVx/s8nSdw7swufujiKTYtmoK0lhGRatb2v1UCrKGDRk2+4uucavD4ny7ftw9g2Efd9obTxyca+HnR0iBZnKId2nPLeur3okNwdwwpe3+dGRBBzcJRCFjn3gwCmEkIUWZZ3Afg7FDj0LCHkvGHbVwCsNQdUi3/rBxCzCqg2mxSyWp3rKwFtDm4Cv17nQ9t+32/OYurEDl9djuzm4CeF3WreAFx1hPIjOdy0aAYSaRV/+Q8HbOcehHrJ6Rh+zmHcJ6FkoSrWfWi1bd47FYcUEfDNfzzo6b4YQSsZ4VYmSgOTQhbgSucuy/IXAPwxgDMAMoSQ78iy/BiAs4SQ78myzAF4CMDXAPwSwFOEkD2G/ecDuBdABMDjhJAdxuM3m3EHGq/SXSUPtBhrwd3b91VkPPl8DnlewF0ej2OeQzKjutL9+/2YBnXfpLYWnImndfqjIxrGih/vL5v7E4t70doS0q+Rmzo5duOshWPhtrbM00tvqMg4V6vuDDPuBbAkJoyuh8GMSERAm4Xk0o8HVYknFokISPN8WeEwzWg12sfUPNdff/dLkB+2LrImP/yC7t0n0yrGtbeUyCe9FAIzJ1UBwRfisnqWrAzxKys/hwefPeR7LNX6UI2m95kVDhvB0NQGQoj3pTbgwiFqoSs35VyNSodKOv9w4ZBl4TBNcZNOq8hnsnoTDi4cqkhZUalKw6x2oqlkBk7Hkc3lMbZNRFzJ4r5/OKDLJ1feImN8u1imSrFTY9Wr7pDVedfvPVKmzvESPK5lCe3RCCZkbWIEoTaIRcP49u7DZQqPjX3ONUyszr1p0Qxse+0YZl89QQ+q8jlnysFJgRSkssLNsZxWCma10553T5apn1bP78baFwkA4N6ZXZb1bB6dN62sZpCf7NpqF+KyOu+p8wqksLsuWTSw7NXqgXnuTYyg9PbGQldk1Rw8Om8apLB9+jnt3GGBx8LrJ6F/92HID7+Au7fvQ0JFSSVFP/0+K52r0VPnbZpp2I3R6N2bvc47PjMZklBQPx15ZA6eWNyL5975UJdiaglQRmjJTT969WjJ8e1yIILo4+oHtPNmlUzDNmUf7WDGvYkRxBJde2nPDCm4deO/YdGTb0CKCMgq9p4g7dx5h2xHmpFW87Bd4lcyV7OxppY5KB7L7YfEXAI3mchAFHgMfhRHLp3Bwusn6fP54Kw19fXbj5NYt/e9kuPbGfAgqAw/lBSjUJoPjJZpYlS6RNeoh0ui4WKRMgHDqayrpTXt3JIo2NIrdvVPolmeusSvZK7GOj5zp0/EUCoDsmqOHtTcfeBEybFaW0J6uV5tm+cPnbRMyjLTNxqsShBb0TZr9pCya0RLXgMKgU3tt6HzyapQUrT9GimgzeAMZtybGJXUyqa+5EXNtxhr0RtIc8gja9I9W567rwenPknZGmE7Ix0NC1CGFUv+1Xy+ZTd34Y6brkCs5YK3aycn1GrCrLxFLikfvHp+N7rGtmLh9ZOQzxS478HhNPp3Hy7bxvwhsbqGG/t6dP7czCeLkQv89PHBhF7p03yNgHIuOqiYg1PBOiuwTNLmBJNCormlU349KprG+MklvWVNLtYs6EZMDIFzSJYJ8UCe45FIZ20ljTT5mySJtvdBO19rSwiDw2ksp5xD21Y7T//cqejffVj/fysteq6YNk/VXlv0VPWr0/YjAQxKE+5Grmp+H+rRB7USNPP7rIHVc2fQPTzN63ULKmcOrqxk7f3PFFQdY6RwiXdn9C4jEQFxjsfynYV6I4/Om4ZJnRKGU1ndcGr70OrlOKXna+cDx2H5jv1l3ueWJddCDOf02uSah/r4ywNYPb8bn7o4ajnn1pYQBuMp2+sSM2zjdA2dauro12Bxr2l1lKPu4+ZctKxc429JSgkGO3oriFpKDLUHC6iOUtAUGTTO/LIxkm3w0hiEfG7/CXxu7Sv4ypY3kM/ny7xRP704NUQiArWKYTQi6MoWY+/V3QdOYO2LBEMp50qcXsosV1qSOZVR8ZUtb6Dnr1/Endv2lShmzEFPrT497VxWCp98OFT223BaxaZFMzypbVgf1OYEM+6jFDRFBs2IfHA24du7s4PRiClqzlG5YZd0NXA6rnvxCUUt2Wb3gRPY9toxbHBIuvEiNbTa1ik/wDgPK0UOHwkjKoXLjHIun7cdl9XxziUy1ufgOE+ql3rJLxkqA+PcMbo4OiPcFteice5G+C3u5ZV77rwkhvt+vL+sGqGWMLT7wAmdQ/5oSCk7tiSgpGOUVQcpwFsRMeO20bDg6j7QuG+yag6GlaxlrZ8nl/RCzeUtx2V1PFpJBKdyELTaMs2ilhlN7zPj3BkskU6riACIo+B5x1Hw0sR8zoIPtq/T7Ue540e5YU66mjI+RlWeWNWGTyYKx9XiBIlinMD8caGpdqyMnHFbc4lcmlGkqYYGTsepCU9Rkd5u0Op4mra+1tmsDI0BRss0IYLsgmOViQkASjyFwY/iSA4lkRhKOXppfpJc/FA55qSr/+dn70GKCDgzpFgm/Nhx+16zXr22V7Tb3orqWD2/G4+/PGDbzcnpuhiP1yGF9d9u65mIV1Z+Dk8vvQHgOE/PTBBtJRlqD0bLoLmWcTQq4+Jo2PMcrKiU+2ZNwZ/edIVepraay+9KpIRBNOf2WsXSTiaJfL6sFrrT/CIRAXwkDEkMlSRU3TdrChbeMMlW6unmumirppAYRiKjUuWpZjApZP3BqkKOQtC8TUWly+hoMHvOc6dPxG3XXIq7ig2kq+2h+Q3UWaX9+1HfeFWB2Mkkt752DO+diiMaCYErBkWdVibptIpcOoPBuIL+3Yfx/KGTuPHKTiy8fhIkHtRVEG3lZrVSSadVqLm8bcVNJ9SrEiVDZWCce5OB9qJJYgjJIcpOFJh52ntnduGBXQc9ceCVwKx5p3UAqha8xgloPPmZIQW3XXNpaQeqvh5wLjTlNN2/MTZghJ9s0Up16vWqRMlQGZjn3mSgeZsJxbsszew50wJ51fTQjN6mKPA1VWB4jRPQpI/pbE7/KOqe8Y5CMTQ3KxMvun8/1TG9rFAiEQGKmtNXBVEpDIHnKqrbzlAfMM+9yUDzNkWBR9Ljscxe43DKe/Zis8NLPXErL1vgOYxpFal0TT6dqajeuRl+vHC3KxR9VbBtX8l2O18/ivc/GtazjuOpbFkpBobGAzPuTQa7ioF+j2csIWBVDKzVpjjXaIJVwDILYDgP6kcx6GYUfigSu5IPRtCkqf1zp2L2+l/guf0n9EDqaH4OmgXMuDchrAyGU10Wt8ctMQKpLLa+ehQbXxrQe4CKkTBiNVDSNBrsqmjyORUb+nrK1C1WH91Kk4H8VgJ185GhrQq6xsVK/s1qyjQHmHFvItQiS1AvzgWUyN+0HqD3P3PQdSBvJMEu4So5rJSU86XdmyBK51p54QUpaAid7dGKngu7xCrjv0cyTTeSwAKqTYJaJ5KYvThjD1C/Lf3sYA7k+ZlXkMld5uO4kTUqwwrUbI4aFA2iLaLxXIMfxZHPZJFQEchzQZOm7nn3JAukNiGY594k8JOqXwnMXpydkqbSJTotkOfFo7XzigFvtWKsjsMFEGyuRuncIJ8LbVWwZUkvJDGkrwru+Mxk/Pnnp4w6Kq7ZwTz3gBGU92hGrRNJzF4crQdoEEv0IDxa2jFCYnmFRTvPltr4m+ccK0o6IajSuSXNvnke49vFkr9X8lyk06reB7aSBDGG+oN57gGimu3Iap1I4qYHaFBLdOOHa+70ibh3Zhe6xsWQLHZucnPtvDQfMXu2xliGtp8R49tFxNNZ7HzjuN5bdVjJgs+perKRG1TSFlGD1TO2ZkE3cnlQW/YxjE4w4x4gqkmdBGEYvMKuByitabPXZXskUmjKTVbNwYlzSbSEeCwza/gNH0daUJmm0bdr2A1AVwW9OnAGV45tw5TxMey97w+x7qdHdGO5YtZVWLZjP8a2iZhdPMbZ4TTGxiIA3BtRt5JEO1g9Y1qnrOcPnazJc8HQHGDGPUBUsx2ZV8NQDWVN0E2b9f237yvxQse2iSWUiPZxtDufJApYPb+7rMZ7glYCIJXF3cXjLLu5Cwuvn1RWw57ngFPnFUzqlDC+XSyrIW9shu33GnoF7Rmb1CnhyCNzGC9eRDPVn68WGOceIKrdjsxtmnpQyhqn+EGlXLnV/vc/cxD3zuzStzF62XbnG05l8dw7H6J/7lSQVXPQP3cqnnvnQ3CwKAHQ14Otrx7VjzP76gmW41h12zRsWXItfvtxEitmXVVWYmDZjuDUQm5h94xZPRd+Y0DViBvVEqxEMTPugaJR2pEFEaB084GoNMjrJmnG+HG0O18+k8XC6yehf/dhyA+/gP7dh7Hw+knIKpmy+jGdrRFsfGlAPwZNCSSJAnK5HP7p7Q8xqVNqiMqIXp4xPx/5SETAuVSm7oaxEmGCouYCkZw2O0bXbKuMIDjVIBAEPeQmflBpkJe2/wdnEwjxXBl/bHc+x2tvokKMx9GaY4xtE/Vg7gdnE0gqWf2jceqTVEPU3fHyjPmJAXHhEJbVUHJrhUrpPkm0bqA+2jJrmeceMLxU+KsWgqCH3Hjlla5UaPuPbRNx5JE52LLkWrQaPDan87m99ubj7Hn3JDYtmoFvflHWPf8Hnz2E4eL+Yj6H9pZQw1RGdDtPPyurRqjdXunKU4uzGDEaFUTMcx+BCEJZ48Yrr3SlYlXPHaqKYSWLO//nW5ZeWxArI6vjcBynZ+ACpR6rZkAjLkoMNBL8rKwaoXY77QOjFbBzuv6iwNdcWdaIcNVmT5blWQDmATgNIE8I+Y7FNrcD+C6A5YSQfyn+9mkAqwC8DeBSAIOEkL8278va7FUO8xwqVQvQ2vk5tXqrBB0dEpIZtS4t3by03LO7tnbPUq0VHH7uYSQiIM3zrlvyVQO09o9uWw92dEhIJJSmVssE0WbP0XOXZVkC8EMAUwkhiizLu2RZ/jwh5GeGba5AwfB/YNp9DICdhJB/Lm7377Is/4QQss9x1AwVwU5y58bI1Ct+UE05qR1oHutwKovOS2Ilun4/fLAbHjlo4+/nHqbTKjo6xLquUKxWnnfcdAXu3r7PdSwg6FLLzQg3nPuNAH5DCNGu0qsAbjVuQAg5Sgh52bwjIeRXmmE3nG/Y72BHM9yoB9xu41ZBUY/4QbXlpDTQ+PwfvXq05DpxkbAvPtiJR65WYTi/97CecSOrDlmxFnqQlMEabjj3cQCM3TnPF3/zBFmW/wjAHkLIf5r/dtFFUXAcdXVRdQgCj46OAAqiVxHnUhksM3h9G/t60NFxoaYIz3OF5bTNNkBRJrat3APasqQXHVLptrWGIPCIAtjY11NCC2zs60E0LCBa5Xsk4ULRrGElix/98ijW7X0PwIXr9PTSG6hGJhoWqM+SEOJt96v1fVHUHCQxhISShSiU+niN9D6o2RyiYYGajJZQsmVjbaTx+0UQc3Bj3E8DaDP8u734m2vIsjwTwEwAK6z+/sknXhvEBYtG59zFVlHvXg9AT6Ax8tDRthbHbYACt0xrsG3kluuR4afdB6vAZS3vT3KocJ2MWnigcJ0SikoNOCrDCvVZEltF2/3c3pdK4YaHb8T3wbJL2MIeqEoG5yw490Ybv1d44Nypf3NDy7wO4HJZljX34SYAP5FleYwsy+1OO8uyfCuA2QCWA/gdWZZvdHFOBgPcyNPstL1G0GiPpJLVKR2praWMIkjzPKS2Fp0mqFb1S6Bx5KQ/6OvB/m9/Ae8/+iXs//YX8IO+HnDIWzbJFnjO9ho4yThrRUcFVVO+1vDazJzBhXEnhCQA3ANgoyzLqwAcLAZT/wrAnwGALMucLMsPA7gcwO2yLM8u/t4L4McA/guAlwH8MwC5GhMZyXDz4rvV9loZmU2LZmA4rerG/Ew8XWYAlu3YjzPxNBSOR1S6UEb3vh/vx9lEBm3tUYixlppnMrqBnw9RiAd6Lx+De4rX5J6n3kbv5WPA8xw6YxE8saQXRx6Zg0fnTcMjP/kP3LltH9I8j0wuT41f2BknPp+ruKSwG9RDxx6UI9AIH/1mgispZLXBpJD2cLuU/jiZcSV7M1Mu4LgSJcKvv/slyA+XywLJqjlY9OQbeGJxL+7avg9j20SsvKW0mFYlsrlq3Ae/ks6WWAvuMlwTAPqH8NpVe/XiYo/9K9GrR954ZSfWfXk6WsWQpz6z2hh3vnkcs6+eYCopHKznbiUzNEtNg7wP9ZLUNvL77AZBSCFZhmod4NWTcbskdbtsNXtAZiWClo5vhNZLU0sm+dWxs7h3ZldZMa1GW+L7pSFaKeqM9mjYtsjZ+ItacPf2fbjqoRew9bVj4CJhx/usjXHd3vcwe/0v8OlvPY+7t+9Djiu8nkFSYLWuf9SsNNBIALvCNYbfuhludLt+tb1mjffjLw9gzYLukmbYq+d3Y+2LRNd+Xzd5TEWt96wCtrTf/Xp4kYgAnufx1J03YOB0HI+/PIDdB06UjdHqnLT68OcNNNevjp3FlPExzJ0+EbsPnMB1k8fg+GACr78/iLnTJ+K2ay7F3dudWwc6FUQLsgFMrfMX6pW3wMA895ojJIaRSKt46s4b8JNlf4CxbWLdPRmzN3dmSEFMDOGJxRd45XU/JTgzpBT6ieZUbFjY47v1Hk3TDQRXqlU7x9Jtb+lVIlfeImPu9IklY6SNReDKA6cbFvbguXd+WzLX44MJrLxFxn2zpmBjXw/W7z0CAK5XNcZmJXtWfBZzp08suY7V8HxryV3XK2+BgXHuAGrH0Vmldq+e3411PyVYd3tPRZK3SudA85jtfg+JYSQyqudUdRrv+8SSXty1rZznfmJxL1o9cNh253h03jRIEUEfox0HzedzyPMCWlsKuncOwFJDE29tNXNmSCkEpZUsVhbr09DiFsZyBlZ89Or53XjunQ+x8PpJEPM5tLVHXZdFCApBvg9RKYxEDq7KBgQFxrkXwGiZGsKqnOoDuwot0oyeTD005jRKx+53P8W07KiSVoqcUxJDuOqhFzxREnYdi4bOJ/X97WiDgvHMIBW/MPYtS65FNCJg4HQca18sBFNDPIf2aBhtLSFs6OvB8h379biFXQEuq5K8D+w6iCeW9CKnFMoYN0IhL7+IRAQkVGDnm5X1nmXwB2bcawgngwNUt8l2NeCF59fntu2tEk8VAM4MKRimZCEOnI5btt2zg1PtdzfbWc1VDOew6Mm3qMlIF3dI2LxoBloNhp5WmZBa/VAMYUgpnL8evXODgvHjpWX6aqsihupjVHPumgpBCPE1aSdG5R9TF7zdWqkLqpmERIPV3B7YdRD3feEqbFo0A7lcHmsWdJfw3GsXTMfjL1/IFNW8aqfxu1WFeFWPuNleGVZwdnAYYs5evRRPWecmHB9MgAuH9BXcJW0inljcC7Lqi7YqKKdrUut73gi14UczRq3nXg8PmeqFpS94ibVQF9RrdWC3comnsri7qJ3XlvAfnE0gJHC6jhy4kE3rNH63qhCv6hEv2zutagSuvI6OFoP5/penYzCbK5+jjWG3uyb1uOfNTCmNBIzagKqbZI5qwIlP9zsuL0Gkes3d7ryxaJgaOPzKljdKDFJLWChJuqr2+J3umfHvCSULVcm4Mpidl8Rw8lwSGTWPy8ZIegzizJCiJ4q5naPTPfVyz4MKSNYjgQlgAVUNo9Zzr5f+1smbqwXHWq+5280tDlh7ealsmZccaxE9j99soPl8DjmOd/S+q+kRx5MZ7Nr3IW675lIsetLwAevroSZR0ebodE/rcc8bpadwo6LawolRy7k3qv62FgWS6jn3SIjHo/Omgawq6OcjocIjmM9krXuUpjMlmmwAeoKR2/GbtexbXzuGhApXenqnGEglMRKt+fZz73yI/rlTQVbNwebFvZB47/fIaft63XNWD8Ya1arfb8SopWXqtWQ0nj/Ir7aXpWi95u5EDWht9uzoD60Gy23XXOq6po35vHtWfBb9uw+7oiic2u95ac9nBbs8Ai/3yGl7L8drdlqjGcbv5l1gtIxP1HPJWG+5Yz1b6I1vF7FnxWfRNS6GgdNxbHploIQaUIaVkrIAYqtYUuBseZGHHjgzrAdeE0oWuTSd5zafN5lWMb69tAEGjaJwCgpWGjS0yyMIMtDLKJLGQi1oslHruRtR6y99NQKazeCtSG0t+DiZKalZs2ZBNy6OhpEYSpXMwcrTfHrpDb68ZKvzbljYg+cPnUT///53APTr79YjrkVFx1qhGZ4lOzTD+GvhuY9azr2eGKn6XycdtZoH7n+mtN7K/c8chGrxabfiso8P+qtlY3Xe5Tv3449mXIrbeiY6attpcQKjDv2Om67AnndPQn74Bdy9fR8SKhqytr0Z9ch3YKhNdc5RS8vUEyNR/+uGaqI2OW4JQTE53lYfwPV7j5Tpwt28ENTziiGsum0acrkclaLgwiF83cLDenJJL4bTasl8V8/vxsCZYew+cMJ1Jm09Sk0Yz211z5q7+2hzoBY0GfPcA4IXD6jWNbVrATeqES+KDattT51XIIUFz0oi2nkHTschiYX71NYetbxvtFVWHpxltq1W393NSqwWigk70O6ZouZqcv7RjmoriZhxDwBeX9KR2A/SDdXk5aNG2zarZDy/EPlMtqyF3er53djz7slCZqzpvkWlsP6hHk5lsezmrpLjXTd5DCRRsJxv17iYvo3TSqzejSxo90wS2YJ+JIDdxQBgVd3PaVnut7FGo8IN1eQ1dZ+2rVcqI51WIUk8Ni/uRasYwsDpeKGs7g2TsPXVo+X3bXFvWZMNANj40oD+7zilmcfA6bjrlZgXxUQ16BvaPUsozbmC1K6RVitqtKuBmOceAJo5QBpUQM2tV+5lKWq1rV8qI5nIIJ/OYDiVwZTxMdzxmcnobI1g40sDJdtpVRnN3vSf3nRFySorn86Uz7evB13jWgsfkYiAtvYopLYWiLEWy+vrlqaqFn1Du2ei0Hxmod4UVyOCee4BoJkCpCUeYCqLXD6Prwegt6+VjtrPKsk4RvNqieZ9GzG+XQTHlSrOaPPNqXmkMiqW79yP8e0iVs6WyySY2vV1W2qikjk7XQ+rOXhFPYPCGqp1jZoZzfeJbkA0S4C0zLvZvg9xJYuxbWIgnK8Xr9zviiHIVZJVyYP1C3uw592T+jZzp0/Eytkylm57q8wjtJqvouZ0Hv2ez3VZSjC16+s29lLNlWGlQb1G8ZibefVcLTDjHgBqHSD1axitAnj3P3NB4QHU5oWw66HqhCBrpKTTKqSwoGvY++dOxQuHTmJe76W6wb/vC1fZGmgzJEM3KbsG4sYxWBlX4z2mBXUbYWVY76CwhkatFVVPMFomINQqQFqJNpnm3WgKD6A2LwRtCb1lSa++jZk+EjggKoYQT2WxadEM3PPU2xjfLmLFrKswqVPCcCoLqa2lsI0HaiCrZCBFhJKqjLdd8yk8saQQfAXgKU08YegmZddqz47KoN1joDSo2wgrw3pVGDXDSHEZn4t4KquvskYbmOfeZKhEm0zzbj44m6gpneQkwbOijz5OZnDfj/fj7u37kM7m8P9+9Vo8dOvvYvf+3+K9U3FIYgjZPLD99WOeqIGyVdfiwgfmrm37cNVDL3jOihUFXqfoNr0yUNZZasPCHvD5XPnKhS9IMAH6PTYHdc1F1eqRadooHrN2H59c0ouHbv1dPPjsIf3ZGa2BVWbcmwx2htHpAabFBsbGIjXV29MMgibBo9FH93yuSzd02VweO94oVIfs330Y8sMv4Ovb9+GLV0/A2DaRSg1YGUEjNYJ8HvcUVxXZXB7rfnrE0kDbfQC1j8W623twcTSMzYt7ceSROXhicS86YxHkeQE73zxearx37EeOFxCJCPTeqi0hS268nrx3I8Wb0mkVai6PZTvqTxM1AkbfjJscNGXO8cEExkhhW2UATR2RGMogMVSL0RdAU4mIAo8knOkjzdDNvnoCHth1sITeWbFzP9Yu6EZcUXVuu6yMrl2JBNO5dx84AZ4Dtiy5FpIouNbVGym6SETAR0UFjVWpAn1OYgjD+bxn9ZWdUiRS/Hu1lCyNVm2yUWiiRgDz3JsMVgqP1fO7sX7vEVeB0CBSniulAGgBaA125QK0/x5WspYBy/HtIiIhHv27D5d5sX5LJJw6ryCXy/m+ZrTG4PfO7MLc6ROxZ8VnQVbNwVAqg9aWEPh8ztIb5vPW1Judp18Lj76RGnI0Ck3UCGDGvclgpfBY+yLBqfNKbTopBUQB2BkEq6X+mgXd2PTKgP4xi4YFDCvlHZlWzLqKuiwPukSC03XSPoA8z1vWj//02FasvEXWaaV7nnobg8Np5Iq0jdadqX/uVOx88zhyFoqiSETAcCoLsmoO9qz4LOZOn4i50ydi731/CI7jkEirgUld3c7X7Qe/GnECakevBgg+1xqMlmlCWCk8NvbV5gH2kyzip1yAcal/6pMUcvk8vv/lHr10wB2fmVzwcPt6sNxQJXJSp0Q14EGXSKDBiv5Zs6AbuTx0Gua6yWOQzKhltNLyHfvx9NIbsPGlAazb+55+zBDP4c8/P6WEWtDPYyiVsHbBdIhhDn/xv0opIKBw7mpQFH6az1SrYU06raKjQ2wYmqieYMbdBo2QeWcFKwMUDQs1aVDgldP0+xJrvHU+k4XA87hv54GyLM5kWoUYEUquA63mi9YU21hfZs+7J7Hw+kmWJRL8yFojEQGKmkNbexQfDyZKPOb7nzmIR+dNw/OHTupzaBWtSxEnFNUV5271oV35zAE8Om9ayW8P7DqI/rlTsfvAiapQFH4++NXOKDV29BqtYMadgnq3wnOC2QBFO4Krwm33UQsy2OfW05eK/DytfZw5eGkVrOXzOSRUYPnThoJgfT2QeCCZqPx+6s/Ltn1Uj3lSp4Qjj8zR50C7lhzyrsoS0D60l42Ryn7rGheriKKweyb8BDFZ4LP6cGXcZVmeBWAegNMA8oSQ71hsczuA7wJYTgj5Fy/7NiJGa60Kp4+a23ooGgLz9PNZDLp46ak1X7gQlu98u4wC2bxoho+rVA6r56XMY05lkTf0eo0Altcyq2Qgwpkaon0cPjhbuoK7bvIYJNMqNi+a4Wv16fRM+Kmt1Ez1mJoVjgFVWZYlAD8E8A1CSD+AblmWP2/a5goUjPcHXvdtVIzWWhVOihKvpRa8qheCSGe3CtZW+37ayTe1IPDWV4+WzMPuWrpRoNCCvx1SuOw3VUn7VrI43RM/QehG0sePVLh5Y24E8BtCiOY2vQrgVgA/0zYghBwFcFSW5f/hdd9GxWj1LNx42l446Wp7+m5R7ftJO34yreqKpucPnSwLilZStoJa1TFDp7H8wOme+AlCN5o+fiTCjXEfB8CY4nK++JsbuNr3oouiZSVVawlB4NFhwVmb+3Vu7OtBNCwEym8HBdocvMJYG0WDlj3q9/gSgC1LeiGJISSUrF4vXDIdThD4qpxfQ7Xvp/n4axZ041v/dEhXyNx4ZWcg87CCms0hGhaAsGD5m/la28H8LHm5J37O6XecbsffjAhiDm6M+2kAbYZ/txd/cwNX+37ySdLl4aqDjg7JUmkSMSkx8plsTRQpfkCbg1fQApKqksG5Cr2qZPEzT7vbHR0SVCVTtfNX+35GIoL+EdNq5Z8ZUhDiuUDnUW2Yn6VqPhPVQFDvQj3hdg5jx7ZR/+bGuL8O4HJZlsUivXITgL+TZXkMgCwh5LzXfV2csyEw0lrhuUG9l8vVPH+172c6raJDEgs1amD9MWlG2qHezwSDPzgad0JIQpblewBslGX5DICDhJCfybL8GICzAL4nyzIH4CEAlwO4XZblDCFkD23fKs6HIQDU+6NW7/MHhZEyD2BkzWW0gMvn8/UeA86cGarrIEbTMq6RwebQGGj2OTT7+AFPtAw1WMlqyzAwMDCMQDDjzsDAwDACwYw7AwMDwwgEM+4MDAwMIxANEVBlYGBgYAgWzHNnYGBgGIFgxp2BgYFhBIIZdwYGBoYRiFHVrEOW5d8BsArAdELIdRZ/51GoST8EYDKAvyeE/H81HaQDXMzhDgD/BcCvAcwA8ANCyGs1HaQNnMZv2O4PUage2kMIebdW43MDN3OQZfmrADqL/5tOCPmvNRyiI1w8RxcD2AzgAICrAPySELKltqOkQ5blT6Mw/rcBXApgkBDy16ZtWgCsBfBbAFMAfI8QcqTWY6XB5RweAPA7AE4CuBbAtwkh/+nm+KPNc/99AP8MgJbV9WUA7YSQRwA8AGCbLMvBtoqvHE5z+BSAFYSQNQDWo/CCNhKcxg9ZlscBuB3Ah7UalEfYzkGW5T8AcDkhZB0h5CEA36rl4FzC6T7cBeBk8V1YAeAHReenUTAGwE5CyBpCyHIAC2VZ7jVtswLAcULIowD+FsDf13iMTnAzhxiA+wghjwHYBWCN24M30s2qOggh/4jSEsRm3IpCsTMQQs4CSAGYWoOhuYbTHAghjxBCUsV/8gDiNRmYSziN37B6eqhmg/IIF8/RVwDwsiwvl2X5uwAazUFwM4dTAMYW/3ssgP2EkFzVB+YShJBfEUL+2fATD2DYtJnxfT4EYLosy+01GqIj3MyBEPLfCSF5w99dv8+jyri7QCW16xsKxWJuywHcV++xeMRfAdhCCPm43gOpAJcDmEQI2YACLfBPRZqjmfAUgIgsy38H4Ak0cDVXWZb/CMAeC7qiad5nmzlof48A+CqAh90ekxn3UlRSu75hUDTsawBsJYS8Xu/xuEWRI70awExZlv8KwEUAvtYsrRkNOA/gDUBfAf4fANPrOiLveAzAPkLInwGYA+ARWZYbahULALIszwQwE8A3LP7cFO+zwxw0w74JwEOEkF+7Pe6oCqhaQZblVgASIeQMgJ8A+CyA7cV69S0ADtdzfG5gnEMxRvC3AHYRQn4uy/J8QsiuOg/RFqZ78CeG37+OQlC7oQKqVjDN4WcAPl38nUchIPZ+HYfnCqY5XAbgIAAQQpKyLH8CQKzn+MyQZflWAH+Awgp1gizLlwMguNBn4icotPr8N1mWpwE44NB/ouZwmkOxD/XjANYSQg57eZ9HVYZqUYGxBMAXUfgSfh/AfwMwjRDy9eKL+CiABIBJKNADjaaWcZrDOgB/jAvG5NOEkMvqMlgLOI2/uE0YhYD2XwLYDuCHhJB/r8+Iy+HiHkQArAYwiEKA+w1CyNY6DdcSLubwuygoOd5Bgco4W2xy3xAoBh5/DuCt4k+tKBjB30NhrN+TZTmKAi12EkAXgO82mFrGzRyeRWE1e0Lbxk5lZsSoMu4MDAwMowWMc2dgYGAYgWDGnYGBgWEEghl3BgYGhhEIZtwZGBgYRiCYcWdgYGAYgWDGnYGBgWEEghl3BgYGhhEIZtwZGBgYRiD+f2Y0L0HGZ39NAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from matplotlib import pyplot\n", "\n", "pyplot.scatter(*nodes)\n", "pyplot.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Having the nodes and weights, we can now apply the quadrature integration.\n", "for a simple example, this might look something like:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:18.024240Z", "iopub.status.busy": "2021-05-18T10:56:18.023969Z", "iopub.status.idle": "2021-05-18T10:56:18.043137Z", "shell.execute_reply": "2021-05-18T10:56:18.043442Z" } }, "outputs": [], "source": [ "from problem_formulation import model_solver\n", "\n", "evaluations = numpy.array([model_solver(node) for node in nodes.T])\n", "estimate = numpy.sum(weights*evaluations.T, axis=-1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How to apply quadrature rules to build model approximations is discussed in\n", "more details in\n", "[pseudo-spectral projection](../main_usage/pseudo_spectral_projection.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gaussian quadrature\n", "\n", "Most integration problems when dealing with polynomial chaos expansion\n", "comes with a weight function $p(x)$ which happens to be the probability\n", "density function. Gaussian quadrature creates weights and abscissas that\n", "are tailored to be optimal with the inclusion of a weight function. It\n", "is therefore not one method, but a collection of methods, each tailored\n", "to different probability density functions.\n", "\n", "In `chaospy` Gaussian quadrature is a functionality attached to each\n", "probability distribution. This means that instead of explicitly supporting a\n", "list of Gaussian quadrature rules, all feasible rules are supported through\n", "the capability of the distribution implementation. For common distribution,\n", "this means that the quadrature rules are calculated analytically, while\n", "others are estimated using numerically stable methods.\n", "\n", "Generating optimal Gaussian quadrature in `chaospy` by passing the flag\n", "`rule=\"gaussian\"` to the\n", "[chaospy.generate_quadrature()](../../api/chaospy.generate_quadrature.rst)\n", "function:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:18.045643Z", "iopub.status.busy": "2021-05-18T10:56:18.045353Z", "iopub.status.idle": "2021-05-18T10:56:18.163031Z", "shell.execute_reply": "2021-05-18T10:56:18.163318Z" } }, "outputs": [], "source": [ "gauss_nodes, gauss_weights = chaospy.generate_quadrature(6, joint, rule=\"gaussian\")" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:18.165417Z", "iopub.status.busy": "2021-05-18T10:56:18.165082Z", "iopub.status.idle": "2021-05-18T10:56:18.228118Z", "shell.execute_reply": "2021-05-18T10:56:18.227679Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD0CAYAAABgk2Y8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABMUElEQVR4nO2deZgcVbn/P7V396xZJvu+VfaQsCUkBBAQFFGBK8giisgmeMF9vfeCP3cEvYhsgnIVEBBQBEQ2WUMIkEASSFLZ922yzdbdtf/+6Jk4hEmmZ6aru2ZyPs/DQ6a76tRb1ae+deo973lfKQxDBAKBQNCzkEttgEAgEAgKjxB3gUAg6IEIcRcIBIIeiBB3gUAg6IEIcRcIBIIeiBB3gUAg6IGopTYAoLa2oSjxmFVVSerqMsU4VKeIu30QfxsLZZ8sS6iqjCRJhGGI74f4flAACw+faxglwsYcNTUV0sG+y0vcTdM8BTgb2AmElmXdcMD33wYGANuAo4D/tixrRfN3FwHTAR9YY1nWnZ05iUIgSQe9DrEg7vZB/G3sqn26rhAoMrqmsmJbPfVZl5SuMqZfOZoEsh9g215JbYyauNsHwsZ8aFfcTdNMAXcAkyzLsk3TfNQ0zZMty3qh1WblwNcsywpN0zwPuBE40zTNIcA3gOnN371lmua/LMtaFcXJCASdRZYl9KTO9gabO19ZxROLt2J7/x6pyxKcMK4fV8wdxdQhVdhpG88rzEheIIiCfEbus4ANlmXZzX/PA84A9ou7ZVn/1Wp7GWhs/vdpwELLslrcLvOBjwFC3Hs4siyh6iqKppDQFPwwJJP1kPwAx+nayLfQyLJEojzBbS+u5rcvrWlzmyCEF62dvGjt5JQJ/bjls9ORMg6u6xfZ2kOjaQqhIpNMaGiKhOPl3jRCzxcPo8OMfMS9H9DQ6u/65s8+hGmaOvB54OqO7FtVlSzKK4yiyFRXpyI/TmeJu32Qn41BGOIDjyzczJ/mb2DjnjSGJnPC2BquPHE0I/uUocnR/N6duYZOEHLLC6u485W1eW3//PKdXP6nhfzu4iMp60S/jep3doKQuozLXS+u4akl26jPuvQp0zl7xhC+OGcEZSklr+veU/phqSm1jfmI+06gotXflc2ffYBmYb8d+L5lWWta7TvmgH1XH7hvsSZGqqtT7NuXLsqxOkPc7YP2bdQMlXo34Nw732B7fXb/57YX8MSSbTyxZBuXzB7BN04dR7bJJggKO5fe0WtoGCordqfzFvYWXlu9i3tfX8/5Rw7BzbqR2pgPRsrguRU7+dajS/BbXdOtdVlufXE1d7+2lrsvPoopAytxMk7R7Ss0wsYcNTUVB/0un1DI+cBw0zSN5r9nA0+ZptnbNM1K2O+XvxO42bKshaZpntO87TPAkaZptgwXZgFPd+IcBN0AWZbQDe1Dwn4gf5i3ngfe3ISilz5YK1QUbn+pY8Lewv+9voGEoRXYoo5jJFRW7GzgG48s/oCwtybrBlz6f2+zN+uh60qRLRSUgnbF3bKsNHAVcItpmj8CljRPpn4H+HLzZveRE/3fmqb5UvN3WJa1Gfgl8CvTNG8C7haTqfkhyxJ6Sqe6dxlGmYGqxn9JgqarPPbOlkMKewu3v7yGREKjlAEFqirjBAEvr/zQi2hebK/PsmDdbgyjxA8pVeHm51bSXoJX2wu45V+rCJT4i7ssS2jN/T9RnkDT4m9z3MirV1qW9Rzw3AGffavVv88+xL73kRN/QQcwUga/fmEVf35zIyeN78cvzpmK35Bp9wYuJbKm8Kf5G/Ladk+Tw7zVuzhmSCXZbGkmWFVV4bkVtXTFM/SP97YzZeDBX42jRlFksm7AG2v35LX9U0u28ZOzptCYdWLdl7Skzh/f2MAdL6/lqBG9uO2CGfh+tuBuvJ5M/IeDh0DTFMoqk8gRTc6VClmWcIKAu19bR5Pj8+SSbazb1YSqxnv0YmgKG/fk72Ncs7MRSSpdF5RliV2NdvsbHoK6jEtI6fqfLEts2Zf/nJXtBexLO8hyfG99SQJDV/jlsytptD1esmpZuGFv7Pt/Z9DLDIxENK69+P7CeeB5Ab7r9bineRiGlOkqNeW5aY6EJjOgKhH78/SDkISWf5dKGQpQunMKw5BUF/3PRgzcZYkO2qCrMnEu0hOGICMxpFcSAE2RGNG3LNY2dxbJ9/EiCqct/YxWFwjDkGymY5EK3YEwhEzG4an/nMOzy7Yze0xflDDELdDy96jI2C4nmv346ztb2t1WkuDk8f1x3dKdUxCEjOvfNZfK6JpyVAkOHX8SHZ7nM6JvOb1SGnvT7d8LY/uVk1AVmoJSWZwf6Sabv18zh6eXbuOoEb2o0BXcdqJ8uiN2hC7J0g87BG3i2h56EPDpyQPorSkdDrcrBbIfcOUJo/La9rjRfUjpCp5XukVAjuMxY1gvBlQmOrW/LMFFxw4jKOE5hCHYtsv5xwzLa/tL54zEi9kisrZwHA/J8fjU5P4MSGo9UtijRoh7jHFdn0zGjd2KzoPhOD6DqhJcMffQAt+nTOeXn5kGbmnPKwwha7tcNDM/YTyQj4zvhwwlX/kZuD5XnziGiQMrD7nd8WP78slpg3C7SX/yvO7V/+OGEHdBQXHSDl/5yBiuP3Pi/jmDFiQpJzBP/uccklLuYVBqAtfn88eNYHB1skP7JTWF735sAnJQeleZ7wd4tsPDV8zk7BmD0ZUP3tYpXeGS40Zw50VHYqed2M/dCAqDFIdJimKl/I37qra42wf52ShJEqqhkjA03ly3m9W1jaR0lRPNGpKaAq4XmbB35hrqhkqDF3L27a+zs6H96BlDlfn9F45m8oCKdld7FsrGfFBVBTQFVZF5fvkO9jQ5DKhKcNL4fji2R+D6eaUt7in9sNQUaYVq11L+CgQdIQxD3KyLZ7scMbCCGYNz7gLX9XGa4veK7dgeFYbGP649nu8+tpQXlu84aOz7jGG9+NGnJzG40uiUsEeJ5/ng+YSKxClj++7PRZ9pEPHhhyNC3AWRkZvsi5+Yt4Vju+i6wo3nTMEPp3Dv6+uZv2Y3DVmPlK4wYWAllx0/kr7lBoHr4cQ4SitXWCS+9gmKgxB3gaAZx/GRvYBEQuOyOSP54uyRyFIu3a8E6BLYWafkE6gCQT4IcRcIaM6Dripoqsxrq3exYN0eVu5oJOP66IrMiL4pjhreixPG9UPXJSTf7zZvJYLDEyHugsOalsnfjB/yy38s58kl2z5QgamF+Wt38+c3NyFLMHdcDV8/dRzDe6fwsm7BaqsKBIVEiLsgMiQpl9hKliXCMLciNE5CqCgyRpnBw29v4mdPr2hT1A8kCOElq5aXV9byuZnD+fbp4/Gy8YrFlmUJRZGRJAgCCIJATKgehghxFxQcVZWRVIVEQmPznjT1WRddVehf2ZwkyfdxbK+kWQkVRSZRZvDNRxbz1NLtHd4/DOGP8zewcMNe/nzZTHQoucAbhkogy8iKzJa9adKOT5mhMrR3Csf1kX0/FmsLBMVBiLugYMiyhJrQcYKAe15bx0NvbfpQvpNjR/bmsuNHMmdsDXbGwSmR31pP6Xz3saWdEvbWvL+1nvN/9wZ/uXIWvu/j+8V/Ymmagp7UWbmjgTtfWctzy3bgtRqpG6rMmdMGccXcUQyqSuBkxKTw4YAQd0FBaCkyfecra/jNv1YfdFS+YN0eFqzbw9DeSR740kwqElrR8+boCY0XrVoeX7y1IO29v7WeXz5jce1HxuKnu5ZCuKPouoKka3zh3rd4c13bOd1tL+CRhZt5ZOFmTp88gJs/My2Wxb0FhUWkHxB0GUmSMMoMfvbPFdzywsGFvTWb9mT4xG9eo97x0YtYyUhRci6j7/9taUHb/cPr69lSlylqVSZNU5ANjbNvf/2gwn4g/3xvO1+49y30pI6iiNu/JyN+XUGX0QyVJxZvzbsKUwt1GZfP3vUGekIrWsEVWVO4740N1GcK6w4KQ/jfF1YRFFEw1YTGlfctZPXOxg7t9+a6PVz/xPsopS4PKIgUIe4xRtNyk5Jxrx+pGyq3vrimU/turcvy5JJtaEUqlp1IaPzpjY49hPLl+eU7CZGKMiLWdYVtdVnmrd7dqf0fW7SFoEi2dhZVlUkkNFHQu5PE95c9jGlxczSF8NyqXey0fZLlCRQlfuUEDUNl4Ya9HSr1diC/n7cOpQjirqoyuxpsttW1X8C7M/hByBtrdxflYRwoMne+srbT+3tByB/nr0eO4cBBknIT3q6i8PyqXWxudElVJmP9IIoj4r0shqgJjYcXbuZHTy3f/9nFs4bzzY+Ow+9izc9C48syf3h9fZfaWL6tga37MtQYSqSTfKqqsDhP33RneXP9HmaN6BXpMSQJkobGE12cEL5/wUauOnE0dTErBKMaGs+vqOUbjyzeP3/ziakD+fk5U0nXd34QcbjRrR+FsiyRSGpI8RvQdhpZltA0hRufsT7w+R/nb6A+66HGoGZnaxRZYv2upi63s25XU+R+d0mSWLG9IdJjrN+VJuoYFFmW2dvk5LXo6lC0pDeO0/0jSTnX2Q+fXPaBifknl2xjw+6mHueiMSJ0u8ZLKTqIYWhUlCdi75PuCLIsUdtgt3njrq1til3VekWWyBRgYUyj7SFFrDIh4EUch16MOreSBNkCveHk+ll81F2SJNKOR10bWTdX7miMXf/vCrIsUVmRgIgGbN36SmUyDvv2pXvUqjvfDxhQmaCm4oNVjAxV5ohh1SWtOdoWnh9Snui6d686qUVe3V6WoMyIdiBQZii5p0iEhGFIqkCRLklNify6d4QgCEloCqP6ln3gc1mC2WP6xK7/d4UgCNm7twkvoiLZ3VrcgR63ECMMIZt1ufvioxjSK1f6rabc4LcXzMB3/djlCAmCgKOGd83HLEswdUh15KsmfT/giKHVkR5j/IAKtIgHwr4fUpXU6HfAAKCjTBpUWbA3gELiZFzu+tyRjK7JCXx1SuOmz0xDl6Uet7LW84LIHq5iQjWGuLbL8KoEz311LhnXJ6WrZDNu0Vdy5oPkB1x2/CgeeHNTp9s40eyHKkfv0nBdn0mDqiI9xjEj+lCMpDnZrMuFxw7jV8+v6nQbX5ozkiCG4u46Hn2TKk98ZQ62G5DUFbJZFzdmla/iTrcfufdUXNulsS5DmHVp2JfGteMn7JATzL7lBjOGVXe6jSvmjkIugq86CEI0WWLakGgEvjKhMmN4dVHchIHrc/GsEaidnISuTKqcPnlAyXL7tIdrezTVZQiyDo11adysW9JEc90RIe4xJ25umLYIXZ/rz5yE3ok45Llj+zJ1SFXRCl8Ens9lx4+KpO1zjhyCbXtF8WH7foBMyBVzO3cu3//4hKLZ2hWCIBSi3kmEuAu6jG27DO+V5LcXTkfrwEKrGcOque3CI8k2Fe9127FdPjKhH2b/ioK2W5lQ+cpHxhIWccLPy7p8+aQxfObIIR3a79qTx/DxyQPwYvo2KCgMQtwFBcHJOBwzrBcPXj6T0TXlh9xWV2Q+e/RQ7rv0WHzbKWoERBjmRPHWC6Z32qXRFv/vU5ORg6CoE35BEGI32Vx/5kS+87HxVCYPPYVWU2Fw02em8aU5I3HSjhgR93DEhKqgYDgZhzG9kjxxzWys5tzir6/ZRUPWQ1dkBlUnOf+YoZx/9DA8P8BO2yWJfrBtj35lOj87ewrfeGRJl9u78NhhnDyhH5mGaNIaHArfD8g22Xx2xmA+P2s4/3xvO394fT2rmuu/lukKkwdX5XLoj6kha7vYTbYQ9sMAKQ4+t9rahqIYUV2dYt++dDEO1Snibh/kb6NhqPiyTMpQ0VUZPwixPR/P8fFdL7KiFvnal8tfYvDs8h1897GluJ2054uzR/D1j5pkG7N5z49E9TtLkpRLn6zKJHUVVZZw/ZCM4yF1oPpVT+qHpaQYNtbUVBz09VOM3AWR0DJBWh/T8LUwBCdtc6rZjyOvm8vVD7zDsm31ee/fr8LgV+dNY+rgqg4Je5SEYYjdHC5rE68cRILiI8RdcNgShjlXUm9D5dGrZvHa6l387tV1hyx8MbqmnEvnjODTRwzGdTyyMUvkJhC0kJe4m6Z5CnA2sBMILcu6oY1tzgN+AlxrWdaTrT7/JjAC2AWMBS61LEukdhPEBsf2cB2PY4dWc+zFRyJLEtaOBpZsrqMh65HQZCYOrGTioEoSqoLreGRiMloXCA5Gu+JummYKuAOYZFmWbZrmo6ZpnmxZ1gutthlJTvg3HbDvAOC7QF/LsgLTNB8n95C4v5AnIRB0lTBkv0tDliXGVCcw+5QhSbnvgiDAzbo0BfF0MwkEB5JPKOQsYINlWS3vn/OAM1pvYFnWOsuyXmxj3zTgAJXNf5cD73fSVoGgKARBiOP4ZDIO6bRDJuNg254YqQu6Ffm4ZfoBrZNg1zd/1i6WZdU3u2UeMk1zG7AZWH3gdlVVycjTvUKuOHJ1dSry43SWuNsH8bcx7vZB/G2Mu30gbMyHfMR9J9B6OV9l82ftYprmEcA3gRmWZXmmad4E/Dfwrdbb1dUVxwUf9/CpuNsH8bcx7vZB/G2Mu30gbGyhpubgK63zccvMB4abptmSX3Q28JRpmr1N06w8xH4Ag4E9lmW1JA7ZBiTyOKZAIBAIukC7I3fLstKmaV4F3GKaZi2wxLKsF0zT/AWwB/iZaZoS8H1gOHCeaZquZVnPAP8EPt48Yt8HTAaui+ZUBAKBQNCCWKEaI+JuH8TfxrjbB/G3Me72gbCxBbFCVSDIA0mScmkTgFCSSOgKiiwRBGB7PoEfoEmIyBlBt0CIu+CwR1VlJFXBMDTmra5l3prdLN1cx/rdTdhugKJIDK5OMnlwFceM6M1HJ/XH94JcvpYeVL9X0LMQ4i44rNETGoEs89uXVvPw25upy7Sd43xf2uX9rfU89NYmjL/KfHzKQK47ZSzVSR3PdsVIXhA7hLgLDktUVUZL6vxrxU5+8Lf3qO9ABXrbC/jrO1t4YvFWvvKRMVw2dxRuxhGjeEGsEOIuKDiSBIah4QGSLKMqEmEIruejSuC7fkmr2KuqjJEy+MYji/nH0u2dbscLQn71/CqeXbaD+790LDoSjlO6mqSKIqHpKm4IqiIjyxJ+EOL7AbqYKzjsEOIuKBiSBKqhkTA0Xlu9i+eX7+D9rfXsTTsossSovuUcMaya844aQrJMI3Q8XLe4o11FyQn7Vfcv4uWVtQVp8/2t9Zx12+v89cvHoYVhSc5J1lUUReaxd7awYN1urO2NZF2fMkNhwsBK5ozpyxlTB+I4Hr4Q+cMCEQoZI+JuHxzcRk1T0BIaT7+/nZ8/bVF7iFS4kgQnmf34+TlT0CRwD+LnLqR9LRhlBj95egUPvrXpoNt0lmlDqvjzZTNJN2QPWXi6kL+zbmiohspP/7Gch9/ejOMf/I0opStcOnskV504GjfrHvQtozv3wzhR6lBIUUNV0GV0XUFNaFx1/yK+8ZclhxR2yGVZ/NeKnZx440u8tWEfekovjp2GyrLtDZEIO8DizXX88Y31qAktkvYPREto7Ei7nHLzy9y3YOMhhR0g7fj85sXVfPK388gAmiFe3HsyQtwFXUJRZLSEzufueZNXVu3q0L5Njs8V9y1k/ro9aBELoiRJ6Amdrz70bqTHufnZVWS8AE1TIj2ObqjsSrucc/vrbKvrWO3W1Tsb+eSt83CQ0HUh8D0VIe6CLqEmNH769HLe2bSvU/sHIXz1ocXYAZEKjW6oPPP+tg4LYUdx/IA7Xl5DqER3aymKhGpofPH/3qbB7twE7s4GmyvvW4SW0IqSkVVQfIS4CzqNYahs3pfhvgUbu9ROxvX56sPvIkco7rKmcPdr6yJrvzWPLtpCwlCR5WhEU9ZUbv3Xajbu6Zo/d9HGvfzt3S3CPdNDEeIeU4yESqI8QU1NBWWVSYyERtwGWIEs85t/raYQc/Kvr9nNniYnEneGosikHZ/3tuRfALsrNNoe89fujuRcWlIk3PfGhoK0d+cra9FjKO6GoWKUGdTUVFBelUSPYf+PO0LcY4ie1Flem+bCexYw8rtPceat83h13R70lNH+zkVCkiBhqDy3bEfB2rxvwQaIwJ2hqjJLt9QVvN1D8ebaPYQRqJFhqLy6alen3TEHsmF3mvW7miKfI+gIWkJjQ73NpX98m5HffYrTfv0qz6zYSaI8IQS+Awhxjxm6rrKzyeHCuxewZHMdYQhrahu5+oF3WLRxH0aRIjHaQ1UVVu9owCtgvPS7m/ZFIoihJPHmuj0Fb/dQLNlSRxTR7j4wf+3ugra5YN1eFDUeUqBpChk/5D/umM9b6/cShrBxT5pvPbqU55btQDfi0f+7A/H4RTtJbgm5FplvsxT4ssQdL69pUzRve3kNgRyPn0xRZJZvb2h/ww6wamcjyQhcBH4I2+ujnUg9kB31NkoEbyE+uWiXQrJ8ez1u6RYMfxBF5p7X1mG3sYL59pfXoMToDaMQaAktsremeChFJ1FUheryRKxeKQvB1n1tC9HWfRlUJR4PMkmCbIFXYrpegBLFyB3wi7wi0wsC5AjORQLcduLZO4pTwlQQB+KHsHlf22U3t+3LosXkDaMQKIpMdUUCVCHuH8LOuuzbl8YukP8xDsjAcWP6tPndzFF9ir60/WCEIfQpK+wcQFVSa3chTmeQgESRBwAJVSGIYvV3mLtOhaQ6pRGXl19Vgrlj+7b53bGjepPuQfe67wfs25fGtwu3Qrs13VrcgdiIXaEIXJ8vHDeC0TVlH/i8ptzgW6eZSBGIX2fwPJ/Jg9srodsxJg2uwo4g8ZYmw7j+5QVv91CM6VeOH8GIWJNh8qCqgrZ55LBexOXd13M8zpo+mCmDP3iOlUmVH5wxETkm/b9QuK4fWZ6f+MVAHeb4foBsezzxlTn87Z0tLNywF7N/BecfOwy/BIm2DobnBfSvTlFTbrSbbiBf5o7tiwoFn4gM/ICjR/QucKuH5oih1WhyNOfykfH9uOm5lQVpT5Ja3gjjMSIOghA77fDwFTN5+r3tzFu9i5F9y7ho5nDwA+xsNKPcnogQ9xjiOh6+53PGxP6cPqEfiiThph38mI1aMlmX848dyi0vrO5yWwlN5uwZQ8g2Fn7i03UDxg+oQFfkSNw+bTFnTB+CCI7lOD4j+5Yxtl85qwowsXrC2Bo0WSradckH1/Xx6rOcMrYvJ47pgyJJ+Fk3dv0/7nR7t0xPJQhCnKyLb3u5/8ewYweuz5fmjKJXqus+4MvnjorsFTUMQ2zH5/TJAwredluMrilneJ+yyIp3uI7H9z4+ocvtyBJ852Pjkfx4vA22JgxD7Jj3/7gjxF3QaXw/IPR8fnHO1C61Y/av4Mq5oyObWAKQ/IArTxgVWfut+eLsEbgRFu1wsi5HDe/FmVMHdqmdq04YzcBKo0cFJAj+jRB3QZdwsi4zR/bm+50cSQ7tneT+y47FyTqRFpBwHI9hvVN8ZHy/yI4BMLxPirOmD8aLuCKTl3X4+TlTmTWq7ciq9vjE1IFcfdIYPOHD7rEIcRd0GTvt8Nmjh3DbhTOoTOY/jXOS2Y8nrpmDFgS4Rag/6mVdbvrMtA7Z2BEkCX5z/nTcIhTM9rwAJ+Nwz+eP4pLZI/Jelq8pEt/86Dh+cc5Usk22qMjUgxHiLugyYRiSbbQ5bkQvXvnmSXxpzshDCuiMYb24++Kj+M35R4Dj4RTJLeC6PlIQ8Ktzj4gkrvvak8cyvFeqqOeTbcry1ZPH8uQ1czh1Yn+Ug5yYocp8+ojBPP+1E7jomGFkGrPCj93DEWX2YkTc7YP2bVRVGUlVSCQ01tU2smjjPnY0ZNFkmUmDK5k6uBpDlQk9P5KwtnyuoZHSed6q5Rt/WUyhBq6XzB7BN04dl9doOIrf2TBUAkVGUWTe31LPu5v20eR4VCc1pg/rxfgBFdiuj+z77U709oR+GAdKXWZPhEIKCornBeAFeLbLwJTGpyb3ByQkKTcB67ke2WxpR4xOxuEUsx+//8LRfPWhd9mb7vxDRldkvn26yXlHDy2pm6NlUtSXJSb0TTFlQPn+VMy+77db11XQ8xDiLoiEMMy5DeKy6Ko1YQh22mb6oEpe+uZJfOfRJTz93vYOtzNtSBW3XjCDSl0h22jHQjyDIMRxPByn1JYISo0Qd8Fhi5N1UVWfX5wzhW+fbnLXK2v527tbSR/CbaHIEqdM6MeVc0czfmAFTsbByQglFcQPIe6CwxrPC/Aabao1hW9+1OT6T05my940izbuZcX2Rhw/QJUlhvdOceTwXozpV07G8ZD9gMa6trMXCgRxQIi7QEBzAjrXpy7jUK3KfHRcDaeOqyEgl1VSkXIPgoa6dEHKCgoEUSPEXSA4AM8LchPDAkE3RsS5CwQCQQ9EiLtAIBD0QPJyy5imeQpwNrATCC3LuqGNbc4DfgJca1nWk60+N4HzgQxwAnC9ZVlvFsB2gUAgEByEdkfupmmmgDuAr1qWdT0w1TTNkw/YZiQ54d90wOcKcDPwQ8uyfg5cCqwrjOkCgUAgOBj5uGVmARssy2optzMPOKP1BpZlrbMs68U29j2aXLDBV0zT/C5wJrCrC/YKBAKBIA/yccv0Axpa/V3f/Fk+DCf3cDjfsqw60zTvAxzg3tYbVVUlkSKoFH8giiJTXZ2K/DidJe72QfxtjLt9EH8b424fCBvzIR9x3wlUtPq7svmzfKgHVliWVdf892vAiRwg7nVFWgwS92RDcbcP4m9j3O2D+NsYd/tA2NhCTU3FQb/Lxy0zHxhumqbR/Pds4CnTNHubplnZzr4LgD7NvnfIjeQLU9lXIBAIBAelXXG3LCsNXAXcYprmj4AllmW9AHwH+DKAaZqSaZo/ICfe55mmeVrzvnuAbwO/Nk3zv4Ea4FeRnIlAIBAI9iPyuceIuNsH8bexUPapai43egtBEOJ5fkFSDxwu1zBKhI05RD53gaAdZFlC1VVCWabMUNm8N82GHQ1kPR9NkRlYlWB0TTmOF+B7PqHnixQFglgjxF1wWKMoMrKuomkKjyzczKOLNrNiWwNOGyXoZAlG9i3nY5MHcMnsEWg6ENOc9QKBEHfBYYtuaKiGyi+fsfjzWxvJuoceiQchrKlt5NYXV3PbS6s5deIAfnLWZHRFxrVdkS1SECuEuAsiRZLYv4ahVCXoDkSSQE8ZrNzZyLUPvsuWfR0PxQ1CeOb97cxbvYsbPjmJ0yf1x07b+H5czjFX2jAMQ/HQOUwR4i4oOLquEMgyqqqgazIZ20eRJQxNoTHrooYhruOVROwlCYyUweNLtvFfj7/XZeFrtD2+/pfFvLluCP9z5iSyTTZ+Gy6dqJEk0A0VD4lUQsMPQhwvIKHJhEDW9lDDcH+tVUHPR4i7oGBomoJqaOxosLnntTUsWLebtbua9gtohaEyeUgVnz5iEGdOG4Rje3hFdmdoSZ3Hl2zjB397r6DtPvT2Zvww5IZPTiZT5GLUekLDMDReX7OLB97cyOJNddQ22vu/H1SVYMbwXlw8azhTBlfhZF0cIfI9HiHugoKgJTR8SeI/H1jEK6vaTh/UYHvMX7Ob+Wt286Mnl/NfZ07kjCkDcdIOnhf9pKRuqKzZlea/Hi+ssLfwyMItTBxYxdlHDMTNuJEcozWyLKEldZZvb+BrDy8+qHtpa12WrUu28eSSbYzrX85vzp/OgAoDN+MIl00PRuRzF3QZPamzbEcjJ9z40kGF/UAabI9vPbKEqx9YhJHS0TSl/Z26gKLIqIbGfz74TqSC9vN/rqDJDdD1aMdNsiyRKE9w60trOO+uN/KeN1i5o5GP3/Iaf313K0aZQRFSOglKhBB3QZfQEhpWbSOf/8ObNHbiVf8lq5ZL/+9t9JSOLEenNJKucvOzFpv3RpvHyPYCrn3wXVQjWnHXUwY3PWtx1ytrO7yvH4Rc/8QyHlm0BS2pR2CdIA4IcRd0GlWVkVSFK+9bhNuFKJH5a3dz58trUQ2tgNb9G7l5MveBNze1v3EBeHvDXrbsy6Lr0byN6AmNtzfs5ffz1nepnR89tZztDTZ6xA8iQWkQ4i7oNJKu8t+Pv8eeJqfLbf32xdXU2V4k7hlVV3nsnS1kirjY6M5X1hAohT8XWZbQDZVvPrK4y235Qcg1D7yDnojmoSooLULcY4phqGgpnfKqJHqZgRGzG1BRZAIknlyyrSDteUHIb19aQ6gUvkuGssRjizYXvN1D8dTSbZRH8JupusoTi7eyq7HrD1SAVTsbWbypDiNmo3fDUFGTuf5vNPd/MT/QMYS4xxA9qbO+zubahxYz9xcvcen/vc07W+sxyoz2dy4Sqq7w0Fsb8QsYq/73d7dSliy8IJYZGu9vrS94u4ci6wZsq8ugqgW+xRSZ+xZsLGiT985fjxcj5dSSGtszHt96bClzf/ESn/v9m8xbvxejLCEEvgMIcY8ZhqGyfm+Gs29/nZdX1lLbaPP2hr18/g9v8drqXbF5hfZCeGv93oK2mXF9NuxOo6qFc2eoqsy2ugx2CZJ8Ld68r6DnIkmQ1FWWFfhBtWjDXpIxGbnrusqejMcnb32N55btoLbRZvHmOq66fxF/X7w1snmZnki3FndZlkjEROwKRaDI3PiMhdfGiPjGZ1eixeQmNHSVFdsKPxp+b0tdQUe7siyzcU9pUsOu3NFY0PKRiiKzZW+6zb7RFXY22PhBWJRSl+0RKDK/em5lm3l+fv38qh53vxuG+oHU0oWkW4u7ZqhUVCR61A+e0FXe21rX5nfrdjWhNOcMKTWqLNHoFH6VY33WK+j5SRIlGbVD7rgBhRNiSZJocqKZFM66fiz6laJILN3Sdv+vbbRpsr1IQ2aLiarKVFYmUSIasHVrcfccn/rGbI/Kl+G4PiP6lLX5XU1Fzuceh1WFfhCSKKDLoYWkphT0/MIQNKU0YqApckFvsDAMSRTah9+MFtHosaP4fnjQ/l9uqJQbalFTO0SJ5wXUNWQJI4riiscv2kl8P8DOuD3mxwbAD7j6pNFtfnX58aPIZKNf1p4PWcdjbP/ygrc7aVBlQRNvBUHIwKpkwdrrCMN7F/a4vh8wpFeq4CPsqqRGQldikbVTCQKu+ciYNs/xc7OGk872rNTKTtaNrB5Atxb3nohju8wc2ZtfnTuNIb1y4tC3XOd7Hx/PBccMxY/AFdIZFAmmD60uaJuqLDG6X3lB88x4ns/wPmWoJXiVnzGsF247OeI7QhiCGwSMrinsQ3XakCqaMvHoV7btMa6mnLs+dySj+uZG8NUpjWtPHstXThpDEJP+3x2Ix+ycYD9hCHaTzUfG9uX0yQMIAUWWyGRcsk12LEZXAKEXcP4xw7j1xTUFa/OUif3J2l7BR2aZ5reM5dsaCtvwIVBliWF9yqjb21TQdgPX5zNHDuGnT68oWJvnHzMMlZC41JOy0zbHDK3mya/MASl3LdMZl0xjNjb9vzsgRu4xJAxzr2uNdRnS9Rnq96Zxs26sOrbr+lQYKnPG9C1Ym1efOBo5KPzkZ+gHfHzKwIK3eyhONGtIR+BC812PC44ZRrJAK3n7VRicZPbDsePh7mvBybo01ef6f10M+393QIh7zImzfzF0fW78j6kktK53o88cOYQRvVORTI4Hrs/nZg4vqmvmirnRPKh8P8T3fH7wiQkFae+mz0zDjnGJwLja1R0Q4i7oNI7jkVQlfvzpKV1qZ3RNOf9z5kTciCaLfT+AICza6H10TTmTB1dGFsXl2S6fPmIwHxnfr0vtfG7WcI4YmiveIeh5CHEXdAkv6/LRCf348acn05mB8dh+5Txy5Sx8x4u2PJ3nc8MnJ1GZjHaaSZLg1+dNi1QwwxCctM2t50/vtMCff/RQvnv6eNxMYXLUCOKHEHdBlwjD3ATYJ6YM4O/XzGFEn1Re+0kSXDpnJI9fPRvZ8yIv++a6PlIQ8NOzuvaW0R5fnD2SYb1SkZ+P5wXYzQL/o09PztsHX53SuOOiGXz/jAlkm7KxKegtKDwiWkbQZXIjSYehFQb/vHYuL1o7+cPr63ln494P5XmvqTD4xNSBfGnOSMp1tagC42ZdThhbwxeOG8G9r68vePszR/Xm66eOI9uULXjbbeF5AX5Dhk9M7s+Z0wZy/xsb+cvCzazb9cEIHVmCcf0ruODYYZwzYwie45FpKI6NgtIhxF1QMBzbxXVcZo/oxXGj+pAyVLbsS1OXdlEUmSG9kiRUhaztgh/gpO32Gy0wdtrmm6eZeEHIfW9sKFi7s0b14e7PH4WTcYo6Gg5DcDO563vhUUP4/HHDkZDYsCeN7fqkdJVhfVK4XkDg+WRFOOFhgxB3QUEJQ7Cb/c1O2qZKlamuzKVN8DMODSUWliAIyTZm+e7pJpMHVXLDE8u6VMRDluBLc0Zx3aljcdJOZKsN28P3A3w/wMm6SJLEwKQKSZUwhMa6tIg6OQwR4i6IFK9ESbsORRCEZBqznD6xHyeN78fXHn6Xeat3d7idsf3K+dV5RzC0OkGmIT4j4jAMS/aQEcQHIe6Cw5IWd0ZCV7jjwiPZ3WRz1ytr+fvibYcs9K0rMieNr+GKuaMZP7ACN+vipEXEiSB+CHEXHNY4jg+OT7Wm8M2Pmlz/ycnsbrRZsqUOa3sDWddHV2WG9U4xfWg1Q3qlaMy6KEFAU12m1OYLBAdFiLtAQC5UEtenLuNgKDKzhlYxZ0QvgjBEliSCIMTzfPY154oR6asEcUeIu0BwAC2TkwJBdyYvcTdN8xTgbGAnEFqWdUMb25wH/AS41rKsJw/4LgksAJ61LOsbXbZaIBAIBIek3RWqpmmmgDuAr1qWdT0w1TTNkw/YZiQ54d90kGZ+BLzTNVMFAoFAkC/5pB+YBWywLKtlxck84IzWG1iWtc6yrBfb2tk0zc8177OuK4YKBAKBIH/yccv0A1pXOahv/qxdTNOcCEywLOt7pmlOPdh2VVXJolReVxSZ6ur8cp+UgrjbB/G3Me72QfxtjLt9IGzMh3zEfSdQ0ervyubP8uEsIGua5neAOYBumuZ1lmX9uvVGdUUKKauuTrFvX7oox+oMcbcP4m9j3O2D+NsYd/tA2NhCTU3FQb/LR9znA8NN0zSaXTOzgdtM0+wNeJZl1R9sR8uyftzyb9M0E0D5gcIuEAgEgsLTrs/dsqw0cBVwi2maPwKWWJb1AvAd4MsApmlKpmn+ABgOnGea5mmt2zBN8xxgLjDTNM3zC3wOAoFAIDgAKYxBRqHa2oaiGBH3V7m42wfxtzHu9kH8bYy7fSBsbKGmpuKgk5WiWIdAIBD0QMQKVYHgACQJVFVBUf499gmCAM8LYpP5USBoDyHuAgEgyxKqriKrCpoqs2ZnI2tqG8m4PposM7R3ivEDK1AkCcfxCD0/lumMBYIWhLgLDmtkWUIxVDRN5ZGFm/nzmxuxdjQctLjFkF5Jzpo+mEuOG4GuQ+B4QuQFsUSIu+CwxTBUVEPjrlfXcsfLa8i67Yv05r0ZfvOv1fz2xdV8+ojB3PCpSciuj9NcfUogiAtC3AWRIcsSqirvX33s+0FsRrl6UmdXxuXKe95k5Y7GDu8fhPDYO1t4ddUubj53GtOGVOGk7ViUs1MUGUWRkaRcURLP88VcwWGIiJYRFBRJkjASGsmKBInyBGv22SzYXMeirQ00BFDdqwwtqaPrSsls1JMaq3c38YnfvNYpYW9NbaPNxX94k6fe246eMgpkYcdRVRktoVHZK4WnKry3s4kFm+pYvqsJKaFRUZVCT2gfmCQW9GzEyF1QMHRDRUto/GPpNv4wbz3vb/3w4uWUrnDG1IFcfeJoeqUMvKxT1FGlbmhsqXe46O43u1QYuzVhCN99bClJTeHEsX1wM8Vz0UgSqIZGIMv87tW1PPz2JnY1frjs3+DqJBfNHMZFM4ejCDfSYYFYxBQj4m4fHNxGPamzrcHm6gcWsaa2qd12ZAm+dPworjt5LE7GKVhB50NdQ0WRMcoMPvqrV9iyr/D5jAxV5l9fP4FyORdR0xkbO4KiSBgpgyeWbOOHTy4j7bR/DXuX6fz87CnMHNUbJ932g7U798M4IRYxCbo9ekrnrY17OfPW1/ISdsj5rO96ZS0X3fMmakJD06J306gJjZ/+Y3kkwg5gewHXPPAOWkIj6iSnsixhlCW44cllfOexpXkJO8CeJofL/rSQe+atx0gZkdspKB1C3AVdwjBUNu7NcOV9i3D9jr+ALdq4lyvvW4SW1CMVGk1TqM963LdgY3QHAd7ZtI+XVu5EN7RIj6MldG57cTUPv725U/v/+vlVPLl0G2oiWjsFpUOIu6DTyLKEamhc/cA7eF3wm7+2ehdPLN6KGqEghorMna+sjaz91vzulXVIanRvIrqhsmlfhttfXtOldm54Yhm2T0kntwXRIcQ9xhiGSiqlYxhqLF+fVV3l4bc3sXFP1/2KP316BYmEFknRFkmCZELjr+9sKXjbbfHOpn3sSTuRuZpkTeV//v4+XZ2Hzrg+P/7HckIlnuKu67n+nyiCm6snIsQ9hmiaQlllktV7s9z/9mYWb2+kvCoVuxGWbqj84fX1BWmrLuPy3LId6EbhA7hUVWFtbSON9sEnOQvNq6tqUdXC316aprA37fD2hr0Fae+f721HUWRkOT7qqaoyqcokmxsd7n97Mws21VFemUKLoG/0ZMTVihmyLKEndS7709vMW717/+cTBlbw8OWzUIMwFguBFEWmPuuxYXfhogGefX87s0f1Llh7LaiqzMICiWG+LNq4j49PHFDwdhVV5tmlWwvWnuMHvLl+D0cPrsQu4sPvYEgSGCmDr/1lMf98b/v+z0f0SfHYl2ejaWHBIqt6Ot1+5B6nEUchUHWVx97Z8gFhB1i+rYH/fWFVpL7cjqCqMu9tqStom0u31KFH4MpwAnivjZj7KLG2NyBHsGDIDeHdTfsK2uZb6/YQxsTvoRsaL6+s/YCwA6zfnebHTy0jjEn/LxRR6le3Fnc9odGnTzm63nNeQHz4UMdu4bnlO1CLEDKYD7Issa0uW9A2dzU66BHdvOkij0qbbA85ortrV6Nd0PZ2Ntp0ItApEtwQnlq6rc3vnl++k1Si59zrsizRp085WkQRS91a3EM/oD7t4Puld1MUkvKDdOCKhBqrHCFqgUcdsgxRLKoLIZKJ2kMhS1LuwFG1XUAUSYJ4DNyRJSg/iG89bv2/q4RhSF2TTRBEo1/dWtxd18dusnuUuGvAJceNaPO7i2cOJ4zJufp+yOia8oK2Obx3GZk8F+N0BEXKrcwsJr3LdPwIhEgChvdJFbTNUTVlaHER9yDkktkj2oyOufDYYbGYFygUYQhO2sGPoM9DNxf3nkg26zJxYAU/PmsyvVK517UyXeHak8dwxtRBuDHp3J7nM2FgRUHbnDakijCCUYwCHDW8V8HbPRSTB1chRfAWogJHjyjspPOxI3sTxGTQ4DgeAysS3HLedGoqconYEprMpbNH8PlZI/APkdZB8EF6jgOrB2E32ZwxsT/nzBjCvrRDVVLHdlyyjdlI3BadIQhCgiDkuNF9eH3N7vZ3yINzjx6KHMFo1/N8pg2tLni7h+LoEb1QgEKn53Icj1Mn9sdQZewCRE3VlBtMHFhFY3188rTYaZvjR/fm1W+d1Nz/NWzHJ9tk9yi3TNSIkXsMCUNwsy6NdWl0P6CpPo2bcePXsT2fq04cXZCmJgysYExN+SETbnUWzwvoW64X3J1xMHRFZu7YGly38OcSBCG+F3DmtEEFae+iWcPI2m4s8tC3xs26NOxr6f9Z3EzPm1uLGiHuMSYMcwUu4nbjtWDbHtOHVnPKhH5dakeW4FfnHoFrR5eG1rE9vnCQuYxCc/rkAfh+gB9RCEroevzXJybud9t1lpF9y7hszqhYuzpy/T+mN0DMEeIu6BJe1uWmz0xjSK9kp9v41ukmAysNnAjnEzzH49yjhlIR8SpHSYKrTxqNFOEo0/MC8HxuvWAGSicjllK6wh0XzcC1vfi9EQoKghB3QZdwXR88n79++ThG15R1aF9Jgm+dZnLRMcMjL3ARBCGe4/HDT02K9DifmzmcgZWJSNxLrXGyLpMHVHD3xUeR7ODah6qkxoOXzWRAmY4T4duSoLQIcRd0Gcf20IOQJ66Zw5fmjMwr/n1k3zIev3o2Fx0zjGyTXZRXb9d2+ejE/pwwriaS9of0SvKd08fjFanKkZNxmD64khe/cQLHjswvgua0Sf15+ZsnMqw6Iaox9XBEJaYYEXf74NA2yrKEYmggS9z3xgZeXlnL+1vr9xeSGNY7xbQhVZx/zDCmD6vGzboFj1tu7xpqmoKS0PiP2+dj7Wgo3HFTGk9cM5tKVW7XvVTo31nXVVRDZfO+DPe+vp5FG/axurYRPwjRFZlxA8o5dmRvLjluJFVJFd/2Dpmfpbv3w7hQ6kpMQtxjRNztg/xsVBQZRVMIJInyZG7ST5Ykso6H7fpohGSz0bgt8rFP1xXQNS66ewFLC5Afp6bC4JErZlFlyLh5nFdUv7OuKwSyjKLIlCU0gjBElqAx6+UWv/lBXkm3eko/LDWlFncR5y4oOLlIkdyE4r7Mh4s1lzqnn+P46Eg8dPlMfvviau54ZW2nV5N+fMoAfnb2VALXy0vYo8RxfMDHI7dWQnB4I8RdcFjiOB6e53P58SP55BGD+OnTK3hlZW3eBTCmDqnia6eO46hhvXCzTizSMAsErRHiLjhsCYIQJ+0wMKXxv+cegRME3L9gI4s27GXpljr2pv894ZjSFSYOrGTa0GrOP3ooA6oSBK5PprGwmTEFgkIhxF1w2JOb1PVQVZlLjh3GxccOoyyh4Qchrh+gyBK6qtCUdZHCECkIyDQIURfEGyHuAkEzucVBOfdKfdZFkiQkKZdjPBODwAOBoCMIcRcIDkIYhrFN/SAQtEde4m6a5inA2cBOILQs64Y2tjkP+AlwrWVZTzZ/Nhr4EbAIGALstizrhwWyXSAQCAQHod0VqqZppoA7gK9alnU9MNU0zZMP2GYkOeHfdMDuvYEHLcu60bKsa4HPmqZ5ZEEsFwgEAsFByWfkPgvYYFlWS+DsPOAM4IWWDSzLWgesM03zf1rvaFnWWwe0JQNNnTdXIBAIBPmQj7j3A1qv065v/qxDmKZ5FvCMZVkrDvyuqipZlBqXiiJTXV2cnN6dIe72QfxtjLt9EH8b424fCBvzIR9x3wm0rqdW2fxZ3pimeRJwEnBdW9/X1WU60lynifuS5bjbB/G3sVD2KYqEoshIkrQ/r36hikUcLtcwSoSNOWpqDl7qMp+skPOB4aZpGs1/zwaeMk2zt2male3tbJrmGcBpwLXAANM0Z+VxTIGg6GiagprQqKxOEeoaK/dkWLStgfdrm8jKMlW9ylCTOkbEOeEFgkLQbi+1LCttmuZVwC2madYCSyzLesE0zV8Ae4CfmaYpAd8HhgPnmabpWpb1TPPk6UPA28CLQBnwW3IPDIEgFqiqjJrQ2N3kcte/VvLP93ewp+nDOXHKdIU5Y2u4Yu5Ixg+oxM06zflcBIL4IbJCxoi42wfxt7Gj9mkJDUmV+c6jS3n6ve157zdtSBW3XjCDCl3ByzodiofvadewFAgbcxwqK6Qo1iGIBEWR0RJazo1RZqClDNSEFiuXhp7UWb07zYk3vtQhYQdYvLmOU25+mWeW78BIGRQhHiAvdF1FTWioqdx1V5M6akJD62C1JkH3Jz53mqBHoOsqoSoTAL9/fT2vrdpFXcbDUGUG90ryuZnDOXpELxzbw3NKV79TS2hYtY187p43sTuZ0dH2Ar796FIyjs/Z0weXNM2ukdBQdJUNu5v43avrWLWzgYzjU2aoTBlcxeXHj6J3eYLQ9bFFab3DAiHugoKhJTSa/JDrH13KCyt2fihH+rJt9Ty3bAcDqxJcfvwozj1qCHbaKVgUSr7ouoIDXPKHtzot7K25/ollTBxUidm3rOil6yRJQk/pLNy4lxufWcmybfUf2mbJ5jruX7CRI4f34gdnTGBU7xROG3n2BT0L4ZYRFAQ9obG1webjt7zKs8t2HLL4xba6LDc8uYzv/fU9EmUGch41VwuJmtD5zwffpamAk6HXPfguqq6iKMW7pSQJjDKdh97ezCX3vt2msLdm4Ya9fOaO+by6Zjd6Si+SlYJSIcRd0GUMQ6XRCzj/dwvYl85/5Pr44q38+B/LMcqM9jcuEImExrsb9zJ/ze6Ctru1Lsvv561DLqJvWzU0nnl/Bz/+x/K89/GCkGsfepf3tjWgGVqE1glKjRB3QZcJFYXv/+096jIdd0ncv2Ajm/dm0PXieAgDWebOV9ZG0vaf5m8gYWhFmVyVZQlNV/nB4+91eF8/CPn6w4tJJIpjq6A0CHGPKbquoKV0qnuXYZQZsYoyaY2qyjhBwCsrazvdxp2vrCUogjtDliUUReLV1bsiaX9ng83ybfVFiUxRdZVHF20m63ZuzmB7fZY31u0u2kO1o+i6gpbM9f9EeSK2/T/OCHGPIZqukg4l/vPBd5l6/bN84d632drkxvI1WlIVfv/a+rxrj7bFP5ZuQ1XlyH3vqiqzYntDpDna31i3G7kIDyrNULl33voutXHXK2sJ1fiFSOq6gqcofOuxpUy74VnOv/sNVu/NoiXi1//jTLcX92IkHCs2RlLjonsW8MqqXWRcn0Ub9/LZ372BkVDjd76yzMtdGLVDLqRw0Ya9qBELjaoqLNqwN9JjLN1ShxNx8I8sS6Rtn7W7upZg9fU1uymPoWAqhsZlf1zIs8t2kHZ83ttSzwV3v4GiKihKzPp/F4nydu7W4p5IaPTtW96jFmioqsyuBps1tR+8cfelXd7dtC9256oqUqd87QeyN+1GPnIPQj5Q9DoKGrIeRKw/sizRWKBY9azrx2rAIMsSXhCyaOMHH8JZN+DllbVoWs9xz0iSRO8+5ZG9kXRrcfc8n7pGu+hx0lESBCGVSY22dK5vuUEc0kW0JgxBK8BoylDlyM9NIoz8AaIUQSjDENQCuX4UWQLi06fCMCShKRjqh8+vX0X8+n9XCMOQukabMCL96ubiHuBknJKtcoyCIAghhAuOGfaBz08YV8Pg6iSuG69EVZ4fMLA62eV2BlUnI79xwxBG9Y02v/ag6mTkN1UYhvRK6c3C3HkqDBVVlmJVJzYMIWt7XHHCqA98Pn1oNVOGVGHbXoksiwYv6+JFdE/3nHecHoRvu3zv4xM4dWJ/Xlm1i+nDqjl5fD/sdPxWFSphyAXHDOtS3PiQXknG9S+nMeK8/p7nM31Yr0iPcdTwXqhAlBIUBCGu6/OR8f14btmOTrdzzpFDaIxhn/Jtl8uPH8Xs0X15dtkOJg6s5GNTBohVtR1EiHsM8f2ApvoM0wdWcMSgylxtwvpsLF9JHdvj1In9qU5pHVrA1JqLZw3HKcKIzPMCBlenqEpqBZknaItZo/vgedG/XclBwBVzR3VJ3C87fiTE0KUZBCGZhiwT+qYYf8IoZCDTkO1Rb+jFoFu7ZXo6tu3hZF2yWTeWwg45F0E263LxrOGd2r/cUDn/6GH4bnFetzNZl3OPGhJJ2zNH9aZMV/AKkK+mPWzbY9KgSsb1L+/U/ieOq6HCUGPn5mtN6/4vhL3jCHEXdBnf8bhy7mhOGFfTof00ReLeS47G93x8vzg3b+j5XHb8KPQIYtGvPmkM+MUTSyfj8qcvHkufso7liRnRJ8Ut508ncHqW/1rwQYS4C7pMEITYaYfbL5zBxyYPyGufckPlgS/NZGzfMtwiZlL0vABDlvj6R8cVtN3TJw9g+pBq7GzxBNNxPAwp5O/XzGZo7/wmtScNquSvX55N6HixHrULuo4Qd0FB8DwfO21z02em8fAVMzl5Qr82wzn7Vxp8/aPjmPftkxjbtzSpZ92suz+vfCEYWJXg5+dMxc2W4Fxsj3JF4tnr5vLzc6Yc1E0zfWg1t10wnUevPA7Z83HEqL3HI8rsxYi42wf52WgYKoEi4wWwYN1u9jQ5JDWFYb1TTB1STdZ2CVw/kvUJ+V7DlmLYF969gMWb6zp9vP6VBo9ddRzlipT3pHAUv3NLIjG1uWDHiu0NNGQ9qlMaUwZXUVNuEHg+ju22G/rYU/phqSl1mT0RLSMoOC2xyKoqc/yIXsjNsdRBENJYn45FXLXr+kgSPHDZTP7fk8t48K1NHW7juNF9+M3501GDoCjRPociCELsrIuddRlUpjFkbF8kKRc37vsBmcZsSe0TFB8h7oLI8LygKJEjncVxfHzf5nsfG8/ZMwbzk3+s4N1N+9rdb3ifFP/5kbGcPqk/btbFiZnv2nF8IF42CYqPEHfBYY3vB/hNNmbfMu774jHsaLD52ztbWLKljpU7Gkg7ProiM7KmjCmDq/jY5AFMGlSJY3tkGrOxeAsRCNpCiLtAAM21T1366AqXzx6OJ0kYqtLsUgqx3QApDJCCMPKVtAJBIRDiLhA0oygyKDKaofHq8h28vGoXjVmPpK4wdXAV5xw5BMf10DRFhBEKYo8Qd0GktBThCMOcDz4O0VltoesKakLn1hdX89Bbm9jT9MGwxkcWbuYnTy/nzGmD+PZp49Hl/KNjSoGmKfsnVMWD6PBEiLsgEoyEhqQqNNgem3enMTSFcf3LsW2P0PNjNdGq6ypoCufeOZ/3t9YfdLusG/CXtzfz0opaHrx8Jn0TaqwEXpYlVF1FN1TW72qiPuvRp0xnUK8knu3h2F5sH66CwiPEXVBQJAn0lMF7W+v55bMrP1B0oTKh8h9HDuGrp45Dtt3mqI7SIssSWlLjs3e9cUhhb01to81n73qDf1w7JzYuGlWVMVIGf3pjA/e+vp5tdf8OfRzVt4wrThjFJ6YMxE7bRUv1ICgtYoWqoKBoSZ1nl+3gwnsWfKiaTn3W4/fz1nP2ba8j6VosqkqpuspfF21hSQcXMtU22vzs6RWEbRSVKDayLGGkDL7y4Dv89OkVHxB2gLW7mvj2o0v56dMr0FNGpKXdBPGh9D1TcFBkWULXldxEXzdA0xSa3IBvP7b0kCGCq3Y28v2/LoUYFGfWdJXfd7LQ9JNLtqE1R9SUElVX+fObG3lh+c5Dbnffgo3MX7MbXY9f3dS26G79P26IqxZTNF0hUZ5gzT4bNamj6d3Ag6bI/O7Vtfh5pGd9+r3tIEslvXF1XWVNbSNrahs7tb/tBTyycDNqiX8b3VD5w+vr89r2jlfWImnxv+11XSHZ3P/lhIZmdIP+HzPi/yu3g94dRK8TJFIGZ902j7Nvf52Tb34ZxVBjX/ndMFSeWrItr229IOSZ97ej66UbvSuKlNeK1EOxZEsdpXRhq6rClr0ZNu/NL/Z+4Ya9BEglf9toDyNlcOHdCzj79tc54caX8GQZNQYusEKjadG9+XXrq6XrKkZSi73odRRZlmjMeqzckRtR1jbYrN/VhCzH++fSFJmGDqS83Zt2kUroAJYkqUP2tkWT7Ze0vLQkQUMHI3bStlfS694eLaa90/zgrcu4vL+lrke6Z7SEhh7RW0m3vlqO45FuzPa42f8gCEnpCieaueIXkwZVMrqmPJIsioXE9gL6VuRfOGJwdaKkFXbCMKQq1TX/c0VCRSqhvIdh2KFiHbIElUkt1pWNWpKdfXLaIABG15QxY3ivWIXPFopsk002opKP3d6n0dOEvQU77XDbBTNw/ICEppBtsmN9QwK4jse5Rw3lpmdXtrttSlc4eUJ/Mg2ly1boeQFzRvftUhtzxvRBlUqXpsvzAvpUJpk4sJJl29oP5TxpfD9cz499vLuddvjp2VP44acmkdQUMmkn9oObzhDlPZ2XuJumeQpwNrATCC3LuqGNbc4DfgJca1nWkx3ZV/BhPM/Hq88gyxINTfG+EVsIXJ/PzRzOHS+toamdGPbzjh6K7fglfWC5rk/v8gQzhvX6UNhmPlQmVU6bNICm+tKm0/Udjy+fOJpr/vzOIbeTJLjmpDHI3UAkfT8g3dz/67tJ/48b7bplTNNMAXcAX7Us63pgqmmaJx+wzUhy4r2po/sKDk3cR+ut8f0AyQ+495JjSB4ihv3EcTV886MmYQyqAYWez+VzR3Zq33OPHEo2Bqs+HdvlRLOGK08YddBtJAl++MlJjO5btj/ffnegO/X/uJGPz30WsMGyLLv573nAGa03sCxrnWVZL3ZmX0HPwsm6jKtJ8fzX5nLBMUNJtYqGmTiwkps+M43bLpzRvFKy9CNIx3aZO7aG0yblV/u1hXH9y7nulLHglX51ahiC3WRzzYljuP/SY5gzpu/+SUlFljht0gAe//JsPjV1EE7aPnRjgh5DPm6ZfkBDq7/rmz/Lh7z2rapKFmX2XlFkqqtTkR+ns8TdPsjfxrIEfO/jE/ifT06iyfbQFBlFllAlCQnQyhMlte9Afn3eNK57KOSZ93e0u+3EgZU8cNmxGIqM0YnziPJ3PnZkb6YOrUaRJTKOT5mh4noBSU0hDEOSevvH7Un9sJSU2sZ8xH0nUNHq78rmz/Ihr33ripQfO+51F+NuH3TcxqyUCzl0KM4rdmevoarK/OrcI3jJ2snvXl23PwyvNSP6pPjCcSM496ihuFmHfU2dGwUX43eWJAlJgsasQxhCR2YFemI/LAVFqqF60O/yEff5wHDTNI1m98ps4DbTNHsDnmVZh5qib3Pf/E0XdHfCkJL7pPPB8wL8hgxzRvbm+LE11DbazFu9iz1NDhUJlSOG9mLCwIr9FZji7gsOw1BUiTrMaVfcLctKm6Z5FXCLaZq1wBLLsl4wTfMXwB7gZ6ZpSsD3geHAeaZpupZlPXOwfSM8H4Gg04Qh2M0VmXprCmdNGYAk5Sox+b6owCToXkjdYVQlEAgEgo7RrVeoCgQCgaBthLgLBAJBD0SIu0AgEPRAun1umQNpL91B82raXwJvAUcAD1iW9fci2jcA+BEwzbKso9v4XiaXxqEBGAHcY1nWG8WyL08bvwDMBNYAM4DfWJb1epxsbLXdCcALwBGWZb0XJ/tM0/w80Kf5v2mWZX2iWPY1H7+937kXcCewGBgHvGZZ1u+KaN/oZvsWAUOA3ZZl/fCAbRLk7uctwFjgZ5ZltZ/cqHj2fRsYAGwDjgL+27KsFcWwr0eN3PNMd/Atcp30Z8DPgZuKayVzgMeBg63aOheotCzrx8C3gT+aplnspOft2TgYuM6yrBuBX5MTgGLTno2YptkPOA/YXCyjWnFI+0zTPB4YblnWzZZlfR/4XjGNa6a9a3g5sK25L14H/KZ58FEsegMPWpZ1o2VZ1wKfNU3zyAO2uQ7YaFnWT4FfAffEzL5y4GuWZf0CeBS4sVjG9ShxJ790BzuAmuZ/1wALi2QbAJZlPcIHV+0eyBnk1gdgWdYecutPJhXBtP20Z6NlWT+2LKtlXYwMdK6UURdoz8ZWb0DfL5pRrcjjd74QkE3TvNY0zZ8ARa9akoeNB94r71qWVbScEZZlvWVZ1uOtPpKBpgM2a32/LAWmmaZZGRf7LMv6L8uywlbfF+1e6Wnink+6g5uBY03TvBn4b+APRbItX7qS7qGoNK9vuBb4WqltaYPvAL+zLKvj6R6Lw3BgmGVZ/0vOrfDXZjdInLgP0E3TvA24ixIuQDRN8yzgmTZcGrG4Xw5hX8v3OvB54AfFsqmniXs+6Q7uBe62LOtrwFnAQ82rbeNCV9I9FI1mYb8RuNeyrPmltqc1zX7YycBJpml+B6gCLo1ZRtJ6YAHsf0PbDkwrqUUf5hfAQsuyvgx8DPixaZpFfYsEME3zJOAk4KttfF3y+6Ud+1qE/Xbg+5ZlrSmWXT1N3PenO2j+ezbwlGmavVu9qg0lN7kBsBcIKPF1ME2zzDTNltffp8i5l2h+6CSA90tlWwutbWyeA/hf4AnLsv5pmuY5pbUuR4uNlmVlLcu6wLKsnzXPrdSRm5gu6eroA37nF4BRzZ/L5Cbd1pbKthYOsHH/vWJZVobcdTQOtm9E9pwBnEbuDXGAaZqzDrifW98vU4DF7aREKap9zfOAdwI3W5a1sJj3So9boWqa5qnAfwC1gGtZ1g0tqRIsy/qZaZpzyE3CLAJGkhuZ3FFE+04ALgZOJ/c0vwn4IjDFsqwrm2/0nwJpYBg510Kxo2Xas/Fmcte4RYxGW5Y1NE42Nm+jkZuU/jrwJ+AOy7KWxcG+5tHcz4Hd5CaoF1iWdW8xbOuAjRPIRYO8Q87Vsac5UKFY9h0JvAy83fxRGfBbYCL/vp+T5Nxa24AxwE+KGC2Tj32PkXuL3NqyzaGiuwpJjxN3gUAgEPQ8t4xAIBAIEOIuEAgEPRIh7gKBQNADEeIuEAgEPRAh7gKBQNADEeIuEAgEPRAh7gKBQNADEeIuEAgEPZD/D0PVQV73H/vXAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pyplot.scatter(*gauss_nodes, s=gauss_weights*1e4)\n", "pyplot.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since `joint` is bivariate consists of a normal and a uniform distribution,\n", "the optimal quadrature here is a tensor product of Gauss-Hermite and\n", "Gauss-Legendre quadrature." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Weightless quadrature\n", "\n", "Most quadrature rules optimized to a given weight function are part of the\n", "Gaussian quadrature family. It does the embedding of the weight function\n", "automatically as that is what it is designed for. This is quite convenient in\n", "uncertainty quantification as weight functions in the form of probability\n", "density functions is almost always assumed. For most other quadrature rules,\n", "including a weight function is typically not canonical, however. So for\n", "consistency and convenience, `chaospy` does a small trick and embeds the\n", "influence of the weight into the quadrature weights:\n", "\n", "$$\\int p(q) g(q) dq \\approx\n", " \\sum_i W_i p(Q_i) g(Q_i) = \\sum_i W^{*}_i g(Q_i)$$\n", "\n", "Here we substitute $W^{*}=W_i p(Q_i)$. This ensures us that all quadrature\n", "rules in `chaospy` behaves similarly as the Gaussian quadrature rules.\n", "\n", "For our bivariate example, we can either choose a single rule for each\n", "dimension, or individual rules. Since `joint` consists of both a normal and a\n", "uniform, it makes sense to choose the latter. For example Genz-Keister is\n", "know to work well with normal distributions and Clenshaw-Curtis should work\n", "well with uniform distributions:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:18.231575Z", "iopub.status.busy": "2021-05-18T10:56:18.231020Z", "iopub.status.idle": "2021-05-18T10:56:18.241273Z", "shell.execute_reply": "2021-05-18T10:56:18.240901Z" } }, "outputs": [], "source": [ "grid_nodes, grid_weights = chaospy.generate_quadrature(\n", " 3, joint, rule=[\"genz_keister_24\", \"fejer_2\"], growth=True)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:18.245319Z", "iopub.status.busy": "2021-05-18T10:56:18.244821Z", "iopub.status.idle": "2021-05-18T10:56:18.337442Z", "shell.execute_reply": "2021-05-18T10:56:18.337156Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD0CAYAAABgk2Y8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABjcklEQVR4nO2dd5wdVfn/32f6vXdreq/AACEBAgFCbwqKoqCISFFUVOx+f9avfm1fReULKopIEVHEBiIdaaGHhIRQQhIy6b1syva99079/XF3N1tuzd2F3c15v15Kdu6ZueeZmfvMM+c85/OIKIqQSCQSydBCeac7IJFIJJK+Rzp3iUQiGYJI5y6RSCRDEOncJRKJZAginbtEIpEMQaRzl0gkkiGI9k53AGDXruZBnY9ZXR2jsTH5TnfjbUPaO3Q5kGyFwW/vyJGVItdnMnLvA4TIeX6HJNLeocuBZCsMbXulc5dIJJIhiHTuEolEMgSRzl0ikUiGIEVNqNq2fTZwIVAHRI7j/KjH598CxgDbgWOB7zuOs7L9s8uAo4EAWOs4zi19132JRCKRZKNg5G7bdhy4Gfia4zg/BGbZtn1Wj2YVwH85jnMtcC/wf+37TgC+DnzdcZxvAp+2bfvgPuy/RJIVRRFlTZYJAaalo5o6QRShKOUdq5z9JZL9oZjIfS6w0XGcdPvf84HzgHkdDRzH+Z8u7RWgpf3f5wBLHMfpSHVcALwHWF1OpyVDAyGgr0VJFUVgxAz8CExNwfUCvKRbYr8EVoXJ86t28djynRw1qYaPHjuRdFsa3w9LOpZuasRiBs0pn4Sh7tcxJJL9oRjnPgpo7vJ3U/u2Xti2bQAfB75Qyr7V1bFBnZKkqgo1NfF3uhv9iqIIwjDjifvC3gBQhCAIQrQiotoA0BSBH4Soee6VdBBy2wvruGHeGkxN4Y+fmMPsSTVoJdxfQRRx72tb+e59ywB48I1trNjWxA/edzgVWvHTVEJAXYvL2b94mt0tLmceOorffWw2Wp6uCAEZ1y8QRT75QgGihPb56Li2Xa/3UGYo/3aLce51QGWXv6vat3Wj3bH/Hviu4zhru+x7UI991/TcdzAvIgCoqYnT0NDWJ8fqzx+VYmoIBEHaK2k/zdIJEKhE+CmvbHs1S+fG59bx4OvbWPCdM6nf25q3vW7p/OGljdz/+lae/n+n0VCf+7uHDa/gxmcyt1/aD7nx6TX87pKjCFLF26zFDP71ypZu2x54fSvXfngWe3a35NirN0ZM5475G9jdknlzeHplHTubUlSr5IzeVcvgew8uoyXl85uPHoXXlv+tQ4sZfO2eN0i6Abdefgxuazpv+0LU1MRJpz1SCExVkGop7XiWpRGqatn9yEXHM7qv3vj68rf7TjByZGXOz4oJQxYAk23bNtv/Pgl4xLbtYbZtV0HnuPwtwC8dx1li2/aH2ts+Dhxj23ZHrDIX+M9+2HDA0J+1U7QIlLD0IQERhlSYGgR9M5wgoogPHjWOL5wxnaQbFNX+/UeO5ctnHkSqQHs/CJk0bF8kNmVEnJLfCaOIicO6R3MTauOkvcJ97YoCHDQq0fl3TFcZnjDyPrwVIi49fjKfOGkKUTEP+Sji8hMm84kTpxDtx7XNhu+HJAwVr0R7ATwvQAlK368UZH2h4hDFVGKybftdwIeBXYDnOM6PbNu+FtjrOM7Pbdv+N3AEsK19l4TjOHPa972MTAZNAKzKli0z2OUHBvvTv1T6wl7T0gHwXL+oNxXL0kEI0imPfPesYWqkEfz6qdXUxnW+cOZBpFtLG+c2TY0kgo/dtpC1u1oZnjC49YpjOHh4gnQJ4/dCCGKVFv9asoXl25q4/IRJTKy28Aq8RVjt5yZV5NtG5tzQfm6K7l5W5L08uMgnP1CUc+9vpHMfXAxke3VLRzdUttQn0RTB+No46aSLm/aL2r9jMvXVTfXMGFdNS9qnJm7w+uZ65kweVtKEqK6rmHGDXc1p2tyAUZUmmiJw29y8D6h3koF8bfuDwW5vPuc+IITDJAcWuq6iGBqKIvDdAK/AHICmKQhDy0y6+0FOR22YGm/VtfDxOxaTbnfA46otHvrSyeh6ccMMuqnxwOvb+O79yzBUhXE1Frua07S6ARcdM4HvvudQ8AtH70IIzLjBZbcvYsnG+s7tP3jf4Vxw1LicGTyGoYLW/rP0fdwCw1CGoaIYmUg/dH1ct7iHmGToI1eoHkDouopWQrZHB5qpMXJkJVr7cEE5CAFm3OC/7nmD82+cz56Uj2nmjzHMhMkPHlrBpbe/TEsQYRjZ2+umznfvX9bp2AG2Naa4/olVRGpxdpuWzu+fy0zIukHIhj1ttLY72Pte24qiKkXlrBumyqNv7ujm2AF+/thKNC37MRRFoFkGX/zHa3z2r0tQDB01T79VVUExdK744yI+9oeXCbX9u7490Q2N4cMrMOJGyfsqish5fSRvL9K5DzA0Xe23Y1txE80s3UFXJkwuvGk+NRVm4cYFUFWFuuY0T71Vx7rdrdz9ymbI4yxVVaE56fPA69tYtrWJe1/diprDgWmqwtpdvTNvlm9r3JdmUQBNEWypz5695YcROxqTKErhn02A4LVN9b22p/2QzXvbsjptTVNZtrWRF1bvZuG6vSzasDevs9Y0hUUb9vLa5gbe3NrIC6t3oWnl3z+RpvCFv70KQuR9uGRDN3Wqq2P9tmhLzfFglPRGOvcBhr8fGQrFkmxNFZzMy0Zza5q7PzuXhuby09uCIGRkpcn5R45jxrgqLj1+EuSZUA2CkEpL42PHTeK4qcP46JyJBH72c+QHIYeMrui1/ahJNUWnWHhByOTh2fOeDVVhbE2MsIisFBU4buqwXtstXWHS8ARBlswj3w+YOb6ac2aM5gx7FMdPHZZ3fN/3Q46fOoyTDhrO8VOHcbo9Cj/HuSkF4Qf85qNHE4VR1n7mw0t7NDa29Vs6b+CHB0T+fV8gJ1T7gME+KVMqZee5dxlDj/wAr8Bkp6pm2iNACcI8Y+466xqSXPaHlzuHUqYMj3Pf508C1y/K8RmWzpPOLr7+r6W9Prv0+El8412HFLXiVQhBosriqjuX8OKa3UDmBeWaC2Zy7uGj8JLZH7K6rhK1R99KEBQx5q5BR7RexBh9IeS9PLiQ2TL9zGC/QUplINubyZbReM7ZRVVM49gpw0gl3YIPkA6EACth8eDSbfxm3hp2NKWoNDUuOX4SXz3rYFKt6aKjWV1XMWIGa+pa2LCnlROnj0BXwEu6AzZXeyBf2/5gsNsrnXs/M9hvkFIZ6PZ2TOpFUYTr+iU7UiEERkzHMDTSfoilK6TSPkHaL3mYAjLRtaIIfD8Y8LoyA/3a9jWD3V6ZCik5YFBVBUVX8RUBkcAwddy0X3ReuRCZ6D8dRvzxubXsbEqRMDQuOW4SoypNSKYJguKfFqapEQiFSAE0FV2I/Vr5KZGUinTukiGDbmhopsbtL67nqbfqiBkqHzl2Au89YmzRwyl6zOA/y3fy3/e92W2e9w8vruey4yfxnfccSrIlXfBhIYTATBgs3drIbS+sp64pzTGTa/niGQdhmkrB3H6JpFykc5e87Rimhm7q6KpCa9LFT+dfNt+xUEdVBK7rZ834URSBbumcf+OL3dIhF63fy3POLn524UySzan8/TI0djSn+c59b2btz10vb2LWhGrOOXQUboGsI83SuGfJVn788IrObSu2N/HAG1t59EunUGloWRccGaaG2p4nHrh+wZW1uqVhWQaCjFxBoX5JDhxkKuQBhBCi2HTvbmimRrwqtl858j0xDBVPKHzwpvkc85MneXHdXtQ8x9V1FXSNT/xpMWf/8jmW7WhGz9JeMzT+uXhz1jz3h5ZuZ0djKrP6Mw+hIrjp2bV5HzS3PL8evcAinY4x/+uecHp91pT0+cVjKwmz5I8bhkajF3LBTS9x/o3zqUv6GHkWeOmmxo4WjzOue5aTr32GDQ0p9AILwopBNzTiVRZarPRFTCALkwwUpHMfYCglLhophZraOFbCKmkfVRUomsr7f/sikVL6opZeKAq/f24tq3a20JTy+eFDy4nlce6qpvDH+RtYsrGebY0pvnf/MtQsC7184FmnlxJ1J0+9VVdwgY+qKbyxuSFvm7W7WtqrPOVuo2kKK7c30ZYjLfHFNbuJZXHCgRBc98QqVte1sG53K9c8upIgzxeFQvDzx1ayvTHFruY01zz6FmERC6wKUVFh8pFbFrJxb1vB1cM90Syd4cMryr9PciDKrLB1ICGd+wCjP2/b5la3ZO2RMIwwNIUrT5pCzFDLXkCiCJjaZZHQxNo4fp5FQQKYPnKfbO7k4YmsfRAR1ORZLj+iwig4Th5FkCjgzDRFoCkib3QfRVAdy/3AqokbBFlsVgTdFlAVkisWwKTaLudyWLxP9HBTXsCVJ05h2siKkrODRBDS0JIuaqHX/iBgwIquDTTkmPsAY39S7YrF349JvCiCZEuaC2dlJiXL/WG5aZ8LZk9AUxW21Cf59ClT8xYPSaU8zjx0FL+/dDab69v42HGTsrbXiLhi7mQefGNbr89iusp7Zo4l1ZJ/zJ0g5ANHjWP5tqacTd49YzQtORYgdeB5AWOrY9ijK3F2Nvf6/JI5E/GzRPWh53P1adMZXxPDD0MunD0hf7EMP+Db7z2UcbUxgiDkypOmkm4rfxVxujXNuYeOwk97Jaduel4A/ZgNJFenFo/Mc+8DBnuubKmUa68QAsvSiBAEflAwNVCIfRrnnpc7V9yqsPjzgg38+qnV+O1OoMrSuOXyYzhsVEVB6QVVFZgJi3N+/XxWfZmYrvLoV05hmK4UfAPSDY2WMOKS2xayee++Y71/1lh+8aFZJFtSWR2Vooj2oRBBOu0VdGaqqqC3zyV4blB2cCDv5cGFXMTUzwz2G6RUBqq9QoiMcqUieGH1LhKGxtzpw0mns2fYZMMwNHxV4Vv3LuXplXWd6ZBHT6zhmgtnMrbCKPpYuqFhxQyWbmlge2OK2ZNrqYlp+KnSI+K3i4F6bfuLwW6vdO79zGC/QUploNurqkpGoyWKcN2g5KEkw8jou0QItjYkGVFhkDA0Qq9wamI2TDOjoxME4YBfwDTQr21fM9jtlStUJQcUQRCWNTzhugG4AaoqmJDQScSNshxAej8eCBJJuchsGYkkC6apga4RaSptXoBp6fu1RkAieaeQkbtkSKFpCkJTQRFEgBJG+EUW4YZ9So4rtjdxx/wNbGtIUh3XufjYiZxhj8rUYy0ynVSIjAxxIDK52VEYooRh2bK8EkkxSOcuGTIYlo4vBLc+v44X1+zG1FQ+cNQ4LjpmAm7KxSvgVDVNwYgZXPWXV5i/Zk+3z551djF1RIK7P3sCRg7pgJ7HMuMmz63axV9f3sSe1jQzxlVz9WnTGZEwcNsKa8JLJOUgJ1T7gME+KVMq5dqr6yqRqhCRiS7ShXRa2qPxMMrks2cbwzZMjZ1tPh+++SWaUt0/t0dX8q/PzcVPuXmzVPS4wfcfXJE1V76Dg0ZV8NAXT6alMbf9QkC8MsbX7n6dx5fv7PaZqgh+e8nRzJ1Si5/F7o5zAxQV5RuGRtC+3F+Lsp+bUpD38uAi34SqHHOXFETXVWJxI6PzUiaGoSIMjV88sYqv3/smq/e0oecpvJ0ZJjG5+YX1/PjRlexM+lk1bjRD42t3v97LsQM4O5u5Yd7qfRWLsqBpCn4IDy/N7dgB1tS18MqGvXmX5RumzrOr6no5doAgjPjGPW9gmlovDRbd0Ag1leufWs0vnlhFWih5dWw6cul/8NAKvvfgcva6YZ9oy6iqQixmlCw9IBlYSOc+wOhP0SXN1NBKrEyvqgp6zOD2BRvRLL1szZBIVfl/97zBPxZv5umVdXzsDy+j6mpOuyNV4fsPLuPm59dx/+tb+dDNL2FZWjd9EU1T2NPq5l1Zes+SLVTmkSfQdY0H39iWr5xrJ/96dQt+HmGAQAj++vKmnJ+3ugFPrdiZKZHXBdXQ+PSdr3DXy5v4x+LNXHb7yxh5HnxGTOeKPy7ioaXbefTNHXzstoVY+yn21ZVYhcldr2ymJaR0bRlDRe/HyWcpSlY80rkPMEQ/3bxCQG1VjESitB+/qiq8tb2JXz21mje2NKBp5d0yuqawcse+JflpP2RbQzLnQ0NRFd7avq99U9Knvs3r9iMXQtDQln9opzHpoeTzOCLTphha0j5RvkMJwd7W/GPqdc3pXgJYcVNjxfZ9D6i1u1ox8pxvXVVYt3ufCua2xsyq13KFtZJuwLWPO9z/2taS70ehqtRUWgVF2vaXCOngi0U69wFG0E8rF6MI6utbaWkqoK/SA8/zsUdXsui/z+LICTUlC4/1Pl7AR+dM7Pz7kNEVTB6WyFm8OgpCLj9hUuffc6bUUhPTu+WxB0HIlBEJdDX3j/7QMZUkvdx9j8IIe0xFUTZMHZ4gn+uKwpAZ46rzHmP2pNpeufjNbS6XHrfP1g/NHk9rlmGmDtpSPu+fNa7z73cdPpogDMvW/1EFLP7u2Vx50pSCk9A9CdIejY1t/bZYKwojqS9TJHJCtQ8Y7JMyxaCqSqczKsdeRRGYCZN1u1rZ2Zzi5ING4CbdnBOHQgjMuMHuNpfdzS4zJ1STak33ch56zOB/H32Le1/dmvU4v774KM48aHjOyVshBImqGCf8bF7BCH7+t84gIcjpwAxDZXc65N2/fp4giyOyR1dy3xdOpLWxu36NqgrMuMXGva0EYcT0kRV5K0hpmoKVMHl9cwNBGHHM5FpSrW7OB2UxdFxbVVUIw3DAFvLuKwb7b1dOqErKpq/UKsMwItmcYnKVyfETqmlrTuXNCImiiFRrmlpN4ZBhMVpyRIWR5/PjDxzBmYeO6rZdUwRfPvMgzj5sFG4e9ckoikinPf7nvMPy9v/S4yZSZWl5I1PXDRiR0LnxkqNJ9CgQYo+u5K5PH4eX5QESBBFtzUnGJ3QmVRq0NiXznnffD2ltSjJjZIJZoytobUqW5di792XoO/ahjozc+4DB/vQvlYFqr66raKbO7tY0zzi7iOkq5x4xBsKIoAiFRSHAjJvMc3ZxzaNvUde8Tz63wtS48sQpXH36dFKtqaKKZGuWjmlqzHurjrqmFLMn1WKPrcRLemUPb/UXA/Xa9heD3V4pHNbPDPYbpFQGur26rrZP6GWEw0p56xACNFPHsnRe3VjPxj2tDK8wOeXgEaRdnyBd/GpX2FdyTwhBGIYDXmdmoF/bvmaw2yuFwyQHFJ5XWCM+F1EEXsrDT3scMSrBrDEVWJZBY2Pbfk3khWFEShatlrwDSOcukWQhivapOZqmLjM0JIOOopy7bdtnAxcCdUDkOM6PsrS5GLgG+IrjOA932f4NYAqwGzgY+JTjOL3L3EgkEomkzyiYLWPbdhy4Gfia4zg/BGbZtn1WjzZTyTj+zT22jwG+A3zJcZwfAAkyDwmJZMCiqgqGpVNRHWP4iAoikZkc7Qv5BYnk7aKYyH0usNFxnI7UgfnAecC8jgaO46wH1tu2/YMe+7YBLlAFNAAVwPIy+yyR5MQ0NUJFwTK1THqj64NffAUkw9KJVIU/vbSBu1/Zws6mFJWWxnkzx/LZ06ZTFTfwkm5RaYKqqqDoKrqhoSmClBcgghA37ck0Q0m/U4xzHwV0LeHe1L6tII7jNLUPy/zTtu3twBZgTc921dWxspdMv5OoqkJNTfyd7sbbxkC11w0j1tS18Ltn1rBow14MVeHdM0bzpTMPprLKQC+wbN0PI9bubuVjty3sJkBW3+Zx18ub+PvizVx/0SzeddgYjDyrYSGzTN4NQm6fv4F7l2yhIelyyKhKrjp1GqccPAJ9gN7vA/Xa9hdD2d5inHsdUNnl76r2bQWxbfso4BvAbMdxfNu2rwe+D3yza7vGxsE9BD/Y06lKpVx7NU1FaAohxcnUqqqCqquEZIpvZMsRN2IGDy/bwXfvX9Zt+10LN3H/a9v452dOYHylgZfju1RVQbH0Xo69K0EY8V93v8F9n69garWZs9+6rhLqKh/43Xw27913b7+ysZ5X/rKEjxw7gR+8bwZtzb3v+45z0/6FBd84DEMlVDokgrOfm1KQ9/LgYuTIypyfFbNCdQEw2bZts/3vk4BHbNseZtt2VYF9xwN7HcfpuOO2A1YR3ykZQCiKwDQ11ALRajHohoowNX773Fp+8NAK1jakMPIoGRqGihE3uXPRZq5/ajW700EviWBNU/CB/3lgWdZjtKR9vvC3VzGySAV3oOgqf35pQ07H3kEYwQ3zVhPmUceMNJX/uX95N8felbtf2cJrm+uxetihmxqRrvL7F9Zz43Pr8BQlr4Svbuq0hPDzJ1ZxzWMODV7fSP4Kkbne5YrESd5ZCt4JjuO02bZ9NfAb27Z3AUsdx5ln2/a1wF7g57ZtC+C7wGTgYtu2PcdxHgceA97bHrE3AEcAX+0fUySFMNvHoUsp86aqCrEKk/lr9nDi9OGkWtOFd8qDbhlcfOsClm3NqB/+Z9kOnvn66VTqatYoVRgan7trCS+u2Q1k5HZf+vZZ3bRuFE3lzws35pXr3bCnjVU7m5labWWP/E2Ne5ZsKcqGZ506NFUhUESvFElFEWiawhMrduQ9xu0vbuCXF83qts2yDN796+fZtDcTST60dBvPfv0MfNfvNUYvhMC0NM76xTPsaVegfHz5DhZ+5ywCLygrdTNWYbJkYz1HTqzBUPyS7hdNU9B1lWSRCpuS/qOox7zjOE8CT/bY9s0u/46An7T/r2ubAPhC+d08cNAMFb8famwqiqCqKkYUReze3VL0foahcv/rW/nOv5fxg/cfzkVHjiu8U54++GHU6dghM8795IodfOSo8Vmdu6VrnY4dIOWFvLh6F2dMH97p3IMI1u1q7bVvT9bUtTC9Npb1M1NT2VmkYmYYQUPSIy4EmdH1fSiKoK4pjVdAnmDDnlZUVdDxmNE0hbrmVKdjB9jZlGbD7hbGxvVe50bTFNbUtXQ6doCmlM+bWxs5dHh8v4dnhID1e9q44o7FXHTsBL7zbhso/n5UDJ2KhEEQlD9ElPX4moKIoqLkHw505HvXACPsJ8nfMIyob0rS2Fya5K/vh7z3iLFcfdp0PnjU+LKEqcIwwlAVxlV3H5k7dsqwnBIBfhBij943rigEHDWxplt7RcCIisI69aOrrZxyuF4QUpWnMEZPEoaW9VhRFFEbL3ycERVGt+g6CCJGVppUx/btG9NVJg1PZD03HTLHZpehE00RHDK6siyRtyiC6SMTfPmsg/jcqdNQKc2Jhp5PQ3MKL4+8clmE0rEXi3TuA4z+XAnpp/2cE4q58LyAyPW5+pQpKL5ftk53OuXy508exwnThjF9ZAU//eARTB2eyBnluSmPO66cwzkzRjN7Ui2/v3Q21T1VGYOQy0+YnPd7hycM5kwelvN7kimP82aNLcqGoyfWoGUZkoGMk1YVwdzpw/Me49LjJ6N02T+KIlIpn79ddTynHDyCudOH89dPH4/nZteyCcOIwAu44xNzmDWhmhnjqrjtimPao9ryAoRUa5pPz53MSEvDLfF+8b0AL9V/qZ5ypXDxSOGwPmCwz7iXSrn2mqZGqCqoikLgB/hpP2+BCcPQCBWBUAQijPCy5IlbFRa/eGwld2UpbycE/O5jszlxSi1eDp0XXVdpjeC0/3sWv4ADueXyY5g7qSanNrxpatQlfT7wu/m0ZhliO3ZyLX/51HG0NqV62W1aGoHIxFxqFBUsHm5YOkJVQEDkh7hl6tjIe3lwIfXcJQOKdNrHa3NJtaTao7z8ztR1ffyUh9fm4uaICr1kmu+89zB+fP4MJg7bN64+e1Itf/308Zw8bXhOxw6ZN5QKXeX6i44kXzr8p06eyknThufVhk+nfUbEdR7+8sm86/DRqO0HrI3rXH3adP78yeNIt7lZ7U6nfPyki590Czp2yLzZpFvTpFvSZTt2ydBCRu59wGB/+pfKQLVXUQSaoWGYGikvRFMFRBD5QVGOEsCIG2yqT/Kbp9cw762dnRk4syfV8rnTpnHi9OGkW9NFDQ90vKGYukrSC0gYGsmUR+iVJkP8djJQr21/MdjtlZK/kgOCMIxwUx5uKlNA26X0MVq3zWVylcn/fWgmunYUDW0uCVNDFYLID0i1pIoeT+5Y5OQLgRDQ1Ja/aLZE0pdI5y4ZkpQz8dbhlANFEBOChK6WFd1FUSS1ZCRvO9K5SyQ5yDwgpFeWDE7khKpEIpEMQaRzl0gkkiGIHJaRSHrQkXWDqqCpCqkgRLd0Ij/A76cVxBJJXyOdu2RIoWkqaAoIQURmIZCfY5VnNnRLRzc07n11C3e/spndzZlsmffOHMMnTpyCroOfKq5YhxACw1TxUcjI0ER9IssrkRSDdO6SIYEQoMcMkn7IH55bx5KN9eiawjmHj+biORM7l8XnQ48ZrNndypV3LKa5x7L73z69hpueXctPPngE5x0xhnQBdUzD1DAtgydW7OCeJVtoSnpMH1XBp0+eyqTaOOm24nLlJZL9RS5i6gMG+0KIUukLe01TQwiBW2RUbRhqe/sg68pOM27y+Fs7+da9S3tJ/9bEdf5x1QmMSRh4OVaWmqbGrlTAeb99gZSXf+jlpo8dzUlTh+VcEWoYKqGm8eGbF7B2V28FzivmTuZb59gkW9JZbTGMTMzleb2lfnvSob0OmRTOcn/P8l4eXEj5AcmAQdMUKqpjvFnXyrw1ezATZt4CE6qqEK+Ksa3V442dLVRUx3q1NwyN7c0pvpnFsQM0tHlc+oeXMS0tZznHUFH4xWMrCzp2gGv+sxIjX58Nnc/85ZWsjh3gzgUbeWTZjl526LpKZXWctQ1JVu1to6I6jmHkLsqtGyqJKouFmxuYv6GeRJWFbsiXcUkGeSdI8pIZ7jCpjOm0pDzctvKKdRhxky/87TWeXpmp1Fgb13niq6ei5yjWoVk6P3poOXe/kimkMaE2xqNfPgXVDzuX8Ieqwk3Prs0b5e5pdfnPsh2cffCIXlIEqqoQCcG8lUVVj2RLfZI3tjRy+Ih4r1J7uq6ysznN4g31eY9xy3PreN/Msd2ifzNu8PE7FvHy+r0AzJpQzd2fmYvntWUt1mHFTD7wu/k4OzMljqeOSPDol08m8Msr1qFbOpUJk5Trk25zS5ZKEAK5aGsAICP3AYbaT6XNhICqmjiJquzFKnJhWjovrt3N9O8+yjynDrMEzfOeaJrCruZ0p2OHTPHp2+evhyxl61RVwQujTscOGcf6z8Wb0btEtJahMr9LQY9cPL2yDi+L01FVweqdzQQlOMRXNuxFyaIwpmlqN/tysXZXC64fdh7DMDRWbm/udOwAS7c0snjD3s5hmq6YZsbmDscOsH53K48v35G1fbEIAbvaPA7//mPc8PQahJ77zSEbesygdlhFv5XoU1SR8+1L0h3p3Aca/RbxiC7/Xxp+mFk+X36RBJF1TDjMExhmiwDD/QwL+zKazHesqMiL2LNVNrty2yoIsrbPOOhyCMKIMIrwyxI36ycHLN8IikY69wFGf6kFRlFEc2MbLU3ZizbnIp3yOOOQkaz6ybmcc/ho0gUKSOfD9wNGVVmccvCIzm1VMY1PnjwFclQbMjXBB4/aV9pvdJXJR4+bhNdFJz3lBpwwLX9xDIDT7JHoWXxOEEQcNKoir9RvT2ZPrs1RrCPg9ENGFdx/6ogEpq50HsN1fQ4fV8XsSTWdbQ4fW8XxU4dnTZ1Mp31OOXgk00cmOrdNHBbj3Bljeg0VlUIUwZhKk+U/Opevv9smKrE4i5d0qd/bWlbFrnyEYVT2pPGBgsyW6QMG+4x7MXQdRy3HXl1XMeMGC9ftYUdjmvfOGkPkhTmzWFRVwUqYrNrZTF1zmpMPHoGb8rpVlDIMja2tLuf99sWcEXVNXGfBt8+krTmV1SnrcYP/96+lzHur8JDK+JoYT/3XabQ0Zj8HsUqLK/64iFc3NeQ8xjUXHMH7Dh/dbfzfMFTMuMmrG+vxw5Djpg4n3ZbOWaBaNzTMmM6TK3YShBHnHDEGL+mVlUffcW0PlHHzwf7bzZctI517HzDYb5BSKddeIcA0dYQA1y1O29wwNITIFNXI5pzNhMmDS7fzPw8s6+WUqiyNv111AhOqzJy57qapsb3N4/wb55MusAr1houP4vSDhudMhdQNFV9V+dDvX2Ljnt7n6eJjJ/D9980g2dK7EpMQ+1IhXbdwKqSiiG7ty82dl/fy4EI6935msN8gpTIQ7RUCjJhJfcrjtufX8eqmBnRV8O4ZY7jshElQRAk6I2awfEczn77zFdqyRMtCwPffdzgfnj2edGs6r+PVTQ3T0nlo6Xb+tWQLzSmPaSMq+PQpUzlkVMV+ZaG8HQzEa9ufDHZ7pXPvZwb7DVIqA9leXVeJVAWhZOQHlPZC0sU6Ut3SUTWVfyzexD1LtrCrOU2FpfGeGWP41MnTsDSBlyxOfkBRBLqhEYhMhkcUhpm6qGWMifc3A/na9geD3V5ZiUlywOB5AZQ4Cdht/5RHqPp85OjxfHTOJHRVIYgiPNcn8gPctuKj7TDsXeB64Lp1yVBDOneJpAdBEBEEHuCRIhPdtcni05JBhkyFlEgkkiGIdO4SiUQyBJHOXSKRSIYgcsxdIsmBokgdE8ngRTp3iaQLQmQKbaCqRELQlHSJ0j7xqhi+65dU1anncQdA1rHkAEI6d8mQoSOvXDM0dE1BEdCS8hF+UFRuuaIIzITJko313Pz8Ohas3dP52aFjKvnkSVN435HjcNvcrPLEPTEMjUhVqIwbBGFIFEEy5RGWkHcvkewv0rlLhgSapmDGTR5auo3bX1zPqp0tCAEnTR/BF8+YzoyxVaTb3JyiU4oisBImP/vPSu56eVOvz1fuaOab977Jv1/byu0fn0MUpfMWyzZiBnuSHjfMW8l/3tyBG4SMrjK57ITJXHniFIK0l1MzRiLpC4paoWrb9tnAhUAdEDmO86MsbS4GrgG+4jjOw12228AlQBI4Dfih4ziLuu4rV6gOLsq1V1UVVF0lAkQRBaMVRXRWPgr8sFfUrCiCWIXFl/7xWk7hr+++91AuPnZiztqnuqVz7+vb+N9H3irY//fNGstPP3gEbp5jLd/ZwpV3LMbNEqEfNKqCf199ImHa62WLqirt1ZQivCJ0d3RdRahKRgk36H1uSkXey4OLssrs2bYdB24GvuY4zg+BWbZtn9WjzVQyjn9zj+0q8Evgx47j/AL4FLC+VAMk7zxqlmIa+4NuaGgxg78s3swvn17D9jYPI27kbq+rxCosHllRx+0LN5FEoMe6FwzRDI17X92SV9Hxp4+uZE+ri56l+IQQYJk6Nz+3rigbHnlzO24QZi1IoSiZmqaf+8uSrI4dYE1dCz98aDmR1r0vuqGhxw3++dpW/rpkK4ql5y1BqFs6rqJw84vruen59ZlzU0Yxla6oqlK2LrzknaWYX+xcYKPjOB1hynzgvK4NHMdZ7zjOM1n2nUNGtf9Ltm1/B3g/ULhkjqRf2J/sD0URxCtjuKpCoiqWtfpQsQgBVszg/Btf5LonVnHngo2c99sXWb+nDcvK7sTMuMHFty7k2/9+k18/tZrTr3uWhlTQrbaobmrcMX9Dwe+/9fl1RFkeUqapsWDdHna1FFdCMIrgjvkbEFrvB4VuaDy8dDvNBcb4H1m6HU1TOs+noggMS+ddv3yenzzyFj//z0rOuv45VF3L+mBVVQVUhbN/9Ry/f24dt76wjrN/+RxuRNlVkMyESaCrJKpiWR+G+RCi7wIBSXkUM+Y+Cmju8ndT+7ZimEzm4XCJ4ziNtm3fBbjAn7o2qq6ODeqUM1VVqKmJ98mxgjBCLcOB5iMSgogIpYRBsDCKuHPhRn766Eq+cc4hfPLEqfttrxDwyqYGNnSRwQ3CiNtfXM9PL5hJTY3Rq/3a3a28ubWxc1vaD/nj/PV86xybeNwEwA0j1u1uLfj9SzbVo2kqFTXdo1shYPm2xhx7ZWdNXQsiy3lI+SEL1+3Jsdc+0n7ImroWZo6rIooyzv0/y3awoynV2WZPq8sDb2zjkjkTe2XoCAF/WrCRpuS+h0irG/D3xZu5+tRp+10HSVUV3trcwMW3LuR9s8ZyzQUzqUmYRe+fDkIsXYUo6pfsID+M0Prw99GXv92BRjHOvQ6o7PJ3Vfu2YmgCVjqO0/HLeRE4nR7OvbGxtOpAA42+HLfLVSi6XISAiqo4QRjS1pwqvEM7pqUzuspCVQSjKy1c10dTxH7ZaxgqVpao0jJUAj+goaV7v3RdxczSPm5kzlFbcwohKLourKoIgjDq1fd43EDZjzeawA9p6DHurpp60W83ihC0tKTxvADL0olliZIThkoy6ZJMdte2icUMEmb29um0Ryq5f1o4tbUJqmI6pqYwusokCHpfl3xolo6lqzQ3p/JOOO8vmqYQBGGfPTiGwJh7zs+KeX9aAEy2bbvj8X0S8Iht28Ns264qsO/LwPD2sXfIRPKrivjOA5b+cOyQGUpobmwrybEDuGmP0w4eweqfvod3HTY6Z8Wkoo7lBthjKpkzpbZzW4Wp8cUzDkLJkjvueQEjKyzOmTGmc9vICpOrTplG1O44Omq7HjmhuuD3n37IqKwFW4Mg5Pipw0qyZfbEGtQsPlwTEeccPrrg/lWWxiGjKzodYDrtcdJBI5g5fp8dh4yu4NwjspfNS6c9PnDkeCYP3xd1jqu2uHjOxG4lCEsliiLGVhi89eNz+drZhxCVeCw/5bWX2eufVE/f7zvHPtQpNlvmXcCHgV2A5zjOj2zbvhbY6zjOz23bFsB3yUyYvgjc5TjO4+37XgCc2b7vJOBLjuN0C9Vltszgoq/K7O1sSvOeI8YQBWHOCkkdKY6r61rY3ZxmbnsFpK5l9kxL56WN9Xzurldzfq+hKsz/9hlofpDV8VRUx3jPDS90GzLKd6wl3zsbL5nuVTS84w3pjOufZXtj7gfpVadM5QunTutmt66rWAmT1zc14Ichx04Zlr/MXntBkHlv1eGHEe8+PPPwdcvQi5f38uBCFuvoZwb7DVIqfVNmT0MI0Sdl9oQAq8Lihnmrue2F3slYhqrw+8tmc+zEGtykm/34ls5LG+q5+q+5HxAdfPbUaVx96jS8HMfSTY296YCLbl7AntbebU49eAQ3X3YMqdZ0L9tlmb23l8Fur3Tu/cxgv0FKZSDaqygCM26ybncrtzy/jje2NKApgrMOHcVVp04jpio5HXsHmTqs2/ifB5bndKofnTOR77/vcJIt2Qttd2BYOoqm8teXN3Lfa9toSftMGRHnkydN5cTpw0m1uvj+wFvENBCvbX8y2O2Vzr2fGew3SKkMZHtNUyNQFHRNyYzH+wGRHxblSIUAPWbQnA647YV13Pf6VpqSPqamcPZho/nsqdOYNjKB29Z7OCYbHYu1hJpJefSDECUMcdOFI/J3ioF8bfuDwW6vdO79zGC/QUplqNvbUYe1Mm50pug2tbkoQVhwNe1gZ6hf254MdntlDVWJpAQ66rDWt0921tTE8QsM6UgkAw25lEwikUiGINK5SyQSyRBEOneJRCIZgkjnLpFIJEMQOaEqkWRB11UUTSGIIIgiDEMb8pkykqGFdO6SIYdhqKhqJs/dzyE3kG9fdI1WN+DuhZvY3eKSMFTOP3IcU0Yk8Etc3t+xglQICMOoqHJ/EklfIJ27ZMigmxq6qbF5b5IlG+vRNYUz7JEYhk7k+QVF2XRDw1UEX/3bq8xf012295bn13HomEquu+hIJlSZObVwOlBVBdFey/XplXU0Jj0OGV3JrAnVpNN+wf0lknKRzl0yJNAtnT1Jny/d8Qortjd1bhcCzj5sNNdfdCSGIKcIl2GoeIrCeb95gbrm7EU7Vu5o5kO/f4l/fuYEptRY3cTLuqKqClbC5P+ecPjby5tId3lzGFdt8dMLZjJ7YjVum8ydl/QfckJV8rZjWDoV1TFqhyXQYkbByj2mpWEmTOJVMXSrt166Yai0+SEX3DS/m2OHjCTwkyt2csltCzFiRs6iMELX+K+7X8/p2DtI+yFX/WUJlmXkLENnxAy+8+83uWP+hm6OHWBbY4pP/XkxS7c2YVi9i4YYlo5VaRGrtDAsPW+pOyEyZfWqa+NU1yaynhvJgYt07pKC6LpKLG50K223vxgxg5c3NXDOr1/gqB8/ya/nrcZKmDmdkm7pOHuSXP7HRZzz6+f5y+LNWBVWNycdqiq/eNyhKZV7PHv5tib+s2xHZ6HtrmiaQpsX8MKa4ipA7mpO8+KaXVmPpesqdS1pHnhjW879wwh+9PCK9kLY+zDiJi+s28uFN73EB343nydX7caM566CZCZM7lu6nZN+/gwn/Gwef1+yBbOEqkm50DSFWMzAzFO/VTLwkc79AMIw9awOKe8+hkqoa9wyfwOeopbl4DVNwYvg6r++yqa9bbSkff68YCN/W7QJzejdL0URaLrK5be/zBtbGtlSn+T6J1bx5Ip9TloIiJka/3lzR8Hvv2vhxqw1VDVd419LtpQk5nXPki34WYrZRarCn17aUHD/NXUtbK5v66xRqusqe1pdvvSP11hd18LaXa18/Z432FSf7JT07Yppaqyua+FHD61gV0uava0uv3jM4fXNDWU7ZTNucsfLm9jR5qObpRXc1nUVs4+KdEvKQzr3AUbPaK6vUFWFeMKguipWUlX7SFH4zbzV/PbpNfzyyVWEyv7fMpqm8PK6PQQ9pHKfW7WLbDG3pqm8ubWRlNd9aGPeW3X47YdQFEFj0sMtQhN+a0MSLYtzD6Moq+56Pva0uGQ7kWEE2xqKKxu5pT7Z+caiaSrPOHW9HjCPr9iBmq3gtSJ4csXOXpuffKuOsIx6xIoieGzZDq5/chXfue9NohKHeVRTx4rp/Rb1q5oqC3AXiXzvGmAE/aTxHQQhLW1uyYWLlSjiirmTqWtO8YkTp6CUoSIaBBFHZCmHN3N8NQrQ0z0HQcghoytRRMZpdnD0pFpUAQGZMfUKU+vVJhu1cSNrYRBFQKJEZ1RhaWQ7kUJAbcLIskdvhicMOlRZgyBk9qTaXm2OmzKMMEufRQTHTqmF57pvP2ZSTVkRWxRFnHnYKM4/chwXHD0OUeL19l0P31cIyij1l48wCBkISraDAfkIHGCUW0knH37Kwy8xzzqd9hluqvzk/TMYHdfKKuHmeQG1MZ1vnWt3Fr4+cfpwPn/6dKIsaYpBEKIJuPbDs6iO6SgCzj9yHB89biJ++4KiMIzw/ZBTDxlZ8Ps/cuwEoiyOMgxCPnDkuJJsOW/m2KyRkRZFfOy4SQX3H1NlceiYqs7sHdf1OWhUBV856yAsXcHUFK46ZRpHT6rNWUN17rThfHTORFRFoAj44FHjOfvw0WVdoygC4fn88LzDOHpcdckpm74b4Ke8fnPA0rEXj9Rz7wMGuyZ0qZRjr6IIVFPHNDTSfoAqBKHr5UxRFAI0Uycey4zjtqV9ItfvtjDJNDU2N7t88Kb5eDmKaIypsnjqv04j3Zq9glKs0uJjf3iZpVsaC9pQaWos+u5ZtDVnP1aiKsbH71jE4g31OY/xiw/N5NxDR+F2cZ4d5ybRPmbdmvII0l7OB76qKqimhqaphFFEFEYEaa+s4tTyXh5c5NNzl5G75G0lDCO8pEtrU5Ig6ZJqSeV07JCJJL2UR2N9G00NbXhtbi/nlU77TKi2uP2KY6mK9Y6np4+s4N+fP5HAze0oA9fnug/PIqYXnjC+5sKZpNO565W6SZc7PjGHkw8a0eszU1P47nsP47wjxuKlu0fFHeemob6VhvpWvKSb900uCELctvZz2Jom3Zouy7FLhhYycu8DBvvTv1QGqr26qWNaOk8s38HiDXvRVIX3zhzDzPHVuEmvoDaMbulsa05z1Z1L2JplUrTK0rjmwpmcetAI0q358+F1XUWzdHY2pbj7lS20pH0OGlnBh4+dQOCH/Tp0UQ4D9dr2F4PdXllmr58Z7DdIqQxke4UQGKZGAAgyE8Kl6LnopoZp6byyoZ5/vrKZPS0uCVPlvTPH8p4jxpQsHaC311ANI1AF+J5fVP3Vd4qBfG37g8FuryyzJzlgiKKIdBm6LV7ax0v7HDmmgsPfdzhCtOvEBCGtTamSo+2Okn1A1nRPiaS/kGPuEkkW0mmfMO0RpDziukp6gA6jSCS5kM5dIpFIhiDSuUskEskQRDp3iUQiGYJI5y6RSCRDEJktI5H0QFEEmqERKQJVUUj6IbqlE3pBVm0aiWQgIp27ZEihqgqKrhIARKAL8Nzcq0l7olsZWeT7X9vKv5ZsYU+rS8LQOGfGaK6YOwU9ivBTblHia0JkpBE8BBGZ12QRhAXL/UkkfYF07pIhgRCgWQaREPx5wQZe29SApgrePWMM7581Fjftd9NxyYYRM1i2o5nP/mUJLT0WPq3Y3sRvn17Df7/3MC46ZgLp1lReB28YKmbcZOG6PdzzyhaaUx5TRyT41MlTqa0wcdvySwtIJOUiV6j2AYN9lVup9IW9mqYghMD3g6KiYFUV7e2zD4uYCZMnVuzk2/9+E7+H0xyWMPjbp49nTIWRc3WpaWqsbUjxkVsW5BQf6+CnHzyC82aMznksw1AJNJWP3rqQVTtben1+5YlT+Po5NqmW7MJjqqogBEXrxHQU/OiLNwJ5Lw8upHCYZMCgqgqxCouWSLCp2aWiOp610lAHiiIw4ibCNGgOM4qLPatBdVQl+sa9S3s5doC9rS4X37qwvdBD9t+CYmh8/4HlBR07wM/+sxLT1HKWBjRiBpffviirYwe446UN/GPRJtQedmuaglVh4WsqaaEQq7TQshXqaEfXVRJVMXalA3YmfSqqe58byYFLUcMytm2fDVwI1AGR4zg/ytLmYuAa4CuO4zzc47MY8DLwhOM4Xy+715K3Fc3SUVWFMAhL1vfuiRk3+MFDK/jXki0ATBwW44HPn4SmqfhZCpXoMYNbX1jPTc+uIYwyhT3+ftXxqEHUObkZqSo3zFud9w2gMenx98WbuPjo8QRBdxsMQ2V7Y4o3txaW+wVoSfs89MY2zj10VC+pA9PUWL6tieXbmnLsneGW59dx6QmT8dNeZ7+tuMk3713KQ0u3A3D2YaP47SVH09qU7GWbEAIzbnDlnxazcN1eAGZPquGuTx+P72d/IygW3dQQmkoYRgRyZe6gpWDkbtt2HLgZ+JrjOD8EZtm2fVaPNlPJOP7NOQ7zE+C18rp6YJAvUisXI26gx4qrEtSBaems2tXKh25ewLIdzRhl1MfUdZVtjalOxw6weW+S3z6zBrLY3VG4+sZn1nRWWXpzayN3LtiI2j4UIQRYhlpUceuH3thOlKVMoKqqPLOyriRbnl21Cy+LzwsQ3Pfa1oL71zWnWVvXgqZl7DBNjaVbGzsdO8BTb9Uxf81uzCx1TE1T4+mVdZ2OHeDVTQ089Ma2vG9ChRACWoOIS257mQeXbkctsUKVamhYFVbON6RykSX2iqeYMzUX2Og4TofG6XzgvK4NHMdZ7zjOM9l2tm378vZ91pfT0QOF/oqRFEWQsHQq46U55wh4YfVuVte18Nyq3WX1T1FE1vqi2xp6R6aQiU7rmnpL627a29alpJ7ADcKixu2TXpC1fqwQmc9Kwc0xHh4BbUWWmOvaHyEEm/b2HvvdsKct6/CPogjW727ttX3jnjaisvyqYPm2JlZsb+LplXVEJdZj1XQVy1BR1f4ZHgqJcg6HSbpTzGN5FNDc5e+m9m0FsW37cOAwx3H+27btWbnaVVfHEGUU9X2nUVWFmpr4O92NwrTXTy21r1844yBOnD6coyfVoCsCRdl/e+dMHUZNXKehbd9wxkXHTMw8dGK9HzyHJEzG18Q69dWFgI/OmUhl3CBqbx8KGFVpUtecX2PdHl2JEKJX3xVFMH1kRUl2jK+NYRkqRo8+hxHMGFdVMHpXBEwbmSBuaMTbu/PuGWP44YPLaW1/OJiawgeOGodpar2icSHgwtkTuGHe6s55AlURXHTMBOKWTixLtF8MiiI45eAR3P3ZEzh0TBUxTYESo3cRRcTjBvF4aW+J7wSD5re7HxRz1eqAyi5/V7VvK4YLgJRt298GTgYM27a/6jjOr7s2amwsrlr8QGWwz7gXQlEEM0YmOsvKlWOvbmo8+IWTuO6JVexuSXPp8ZM4fuowGhvaskbfuqHxwBdO4rdPr2Z3i8sVcyczbXiC+vp9Uath6Vx6/CR+9dTqvN991SlTwfNpaOme5qgogrMOG0WFqfVKgczFJ+ZOwUt5tPWI+FVVcPGxE7n2MQc3z4Kn0+1RKNDtPOqWziNfPoXfPbuGIIy4+rTp6EB9ffZznYgZ3PPZudz4zBqCED5/+nSGxfWc7YuhpiZOa1MSe1gcP+nSMMQXbQ323+7IkZU5PyvGuS8AJtu2bbYPzZwE3GTb9jDAdxwn58yR4zg/7fi3bdsWUNHTsUsGPmEYFaxiVCxe2qfa0PjR+w4jInMD5ssZ91wfXVf56pkHEUWgEZFu6x6hB17Ap0+ZxuPLd7Jie/bb8eJjJ3DQqAqSzalen4VhRCrtc8lxE7nthcKjh8dPHcaICpNUS+9jBUGEEoT86PwZfOe+N7PuPzxh8NMPHoHoMbTjpTxqDY3/PsfOnJsC2vRu0mVajcW1F84EBBoh6Ta3YP8LEUX02fWWvHMUHHN3HKcNuBr4jW3bPwGWOo4zD/g28HkA27aFbdvfAyYDF9u2fU7XY9i2/SHgVOAE27Yv6WMbJIMM1/XxUxmt9IxOev72nhdk2qe9rFWVgiDET7nc/dkTuPKkKVR2GUaYUBvjh+8/nO+/fwZuHscXeQFfO/uQrHVPuzJleJybLzuGIJ3b6fopl/NmjuGPHz+WI8ZXdW7XVcH7Z43l0S+fQkzJ7kC7nZsi3iLSaZ8g5RGkXNIp6ZAl+5CLmPqAwf5qVyoD1V5NU0BTsUyNnU0pVEVhRIWBm/bxi5Ag0HUVI2Zw58KN3PnSBrY17ovMqyyNDx8zga+efQih6+MVEdmalo6iq7S6AW1pn5GVJp4fogRB3qLg7yQD9dr2F4PdXllDtZ8Z7DdIqQx0e4UQqKogiihZ6EtRRCadz9RZVdfM7pY0VZbOjHFVpNM+kR8UvXK0g44Vp2EYDXjJgYF+bfuawW6vrKEqOaCIogjf3z8nGoYRYcrDS3lMqjCYUmUSjxs0NCT3ezGPVJKUvBPIFQESSQ48L8hE6xFylaZk0CGdu0QikQxBpHOXSCSSIYh07hKJRDIEkROqEkkWTFPDR4AiaPMCTEvHTfty7F0yaJDOXTKkUBSBbmh0JMtoImp3ysXtr5sapqXzyoZ6/rZoE7tb0iRMjfOPHMd7jhhDOu2XJHtsmhqBEIQRaCKzmlZmz0jeDqRzlwwZdEtHNzTue20ri9bvRVcF7505lpOmjyCdcnELrPg0LJ36dMDlNz/P5r3d9Y6edXbxw4eW87tLZjNrXBVuMv8y/44FUavrWvjn4s00pTwOHl3B5cdPRgd8qZMu6WfkIqY+YLAvhCiVviqzpygCzyu2zF7+0nMd9U8/c+crnaqKHUweHuefnzmBmCCngzeMTPWjc379PPVtuSNzTRH8/aoTOGREPGdNVl1X0SyDq+58hQXr9nS3QxF8972HctHsCaTb0lltl2X23j4Gu72yzJ5kwKCqorPM3oamTJk9PY+k7L4yezqNQabMXocz60DXVVq8gE/+aXEvxw4ZjfOLb12IGTOy6rkDhKrKzx9bmdexA/hhxPfuX4aWrzSgofH/7nm9l2MHCMKIHz/8Fos21mP0kOVV1UyZPVdVSJZQZq8uFbCjvcxez3MjOXCRwzKSgmiWjmFoeG5p483Z0GMGP39sJXe9vAnICHE9+MWT0YIwa6SqWTp/XriBXz+1mjCCoybW8LdPH08Q7CslF6kKNz21mnSeSHfjnjZeXL2L4yZW9xLYUhSBrik80qUKUj6cnc2s393KxAqjl/iXpmW0ZJ5YsTPvMX4zbzV//fTx0OV8mgmT7923jPtez2jBnzNjNL+++ChaskhiC5EpWXjVnUt4sb0K1Zwptdz5yePws5TlKwXd1NBNHc8LCNLegJdMkGRHRu4DjP4qsycExCotzLhZ0n6WpbOyroV3/+p5lm5vxiqjzJ6mKbSkg07HDplKQ3e8tB5F6x1xKooAReGGefvK7L2+uYF/v7YVo0u0bxoaT72V35kCPLR0O36WMkWaprByR3Peh0NPnl+1K+u10nWVx5bvKOhc39jSSBjRWVVI11U27WntdOwAjy/fyfJtTVnL5hmGxqubGjodO8DiDfW8sDp7Wb5i6Siz997fvMC9r29FKbFkn2ZqVNbE+60cXseQlaQw0rkPMPorSBJCYOoqllnia7sQrNrZwtaGJG9tbyqrzJ4QgqYskf/eVo9sblUIQVvaJ+hxUna3dNdzV4TorEaUD9cPc/a/53cUPFaQ/VhCUPRDwu+SNSOEoCHZ+9zsaXXJUvYVIQT1rb0ndXe3pMt0foKdTWk27mlj2dbSr7eiKuiq0m+l8MIoGtRV295O5LDMACPspzS5MIxobsxe7Sgfbtrjw8dM4OzDRlEbN0i2pIjtZ/TueQGThsU5fGxVZ1ENQ1W44oRJKFmcaxCEVMUNTpw+nJfWZsavK0yNS+ZMIujiQFOuz5ETanjGyV8g7NjJtWgCerrEMIyYNKy0UmszxlYRZemz74ecMHVYwf3HVFnEDY3m9qwbz/M5emItU4bH2bAnM8E3rtrilINHZC0w4ro+Zxw6ipEVJrvaH3Y1cZ33zRpXMJMnH1EUcdDIBAu/cybVMZ10m0sp07Re0qXJ9UtWziy6f2HUb3WGhxrSuR9ABEVEtz0Jw4i2piQxVaG1qfxyiG7S5Z7PzeXvizaxoynFpcdPotbScLNErZCpTnT7x4/l8eU72NaY4qJjJqADXmqfy1HCkM+eOi2vc7d0hYvnTCTd2rvOqu+HJCp0jp86jJfX7y1ow7CEwckHj6S1qXeWhev6HDK6kukjK1i7qyXnMa6YO5lUl4IfUQSppMsjXz6Ff7+6BT+MuOiYCbjJ7GPeYRjhpT2e+Nqp/GnBBsIg4uMnTiH0/bLz6NOtaYz2611qMBBFxWf5SPoXmQrZBwz2dKpSKddeVRVoupaJwMKwYOEKRRGYpgZC4OWICmMVFr95eg23vrCu12e6Krj5smM4ZkJ1zglh09JYuTuTVVPoJ/Hf7z2Ujxw9PuexdEOj3g248PcvZc2+OeXgEdxy2TEkW1K9HLeqis4xdtctvOBJ05TMfIWA0Ctda74n8l4eXMhiHf3MYL9BSmUg2qsoAjNu8ua2Rm5+bh2vbqxH1xTOOmwUXzj9IGpjWt4yewBm3OSJlTv51r1v5hyDv/LEKXz93YeQak3nzSIxTJ1QVbjthXXc/9pWmlI+00Ym+ORJU3j34WNIt7n4/sCrxjQQr21/Mtjtlc69nxnsN0ipDGR7LUsnUAQxU8sMdaR9RJFl7YTIpGo2pHxueW4d97++lTY3QFUEZx82ms+eOo1DRlWQbsvv2DvQNAWhqRimhqYopDwf/HBAa9QM5GvbHwx2e6Vz72cG+w1SKkPdXl1XiTSVqriOH0ToqkJT0kUJwqKKVg9mhvq17clgt1eW2ZNISsDzAvAC9iZdhIDq6jhegSEdiWSgIfPcJZI8DIAXW4lkv5DOXSKRSIYg0rlLJBLJEESOuUskPejIconHdHRVIQgjtJiOEhTOyZdIBgrSuUuGFB2pkAlLJ4wikqniUyEhU/AjEII/vLief72yhV0taSpMjffMHMPVp01neNzAS7pFjcV3PCRMU0NXFZKuTzTAUyElQwfp3CVDAiEEZtzgrR3N3PTsWhZt2IuhKbzr8NF85cyDqYzpOSUOOtAtnXV7k1zxx0W0dEl5bEn73PPKFu5dsoWfXTiT9xw+hnRbbxmDbscyNIShcctza/n3q1tpSLocMrqST588lbMOG026NS3L7Un6FZnn3gcM9lzZUinXXl1XiVSFCNCIeumr96QjAg4BLYqy5pobcYN/v76NHz20otdnMV3lH585gcnVJl6OPHVVVcDUOfXaZ7o59p4IAfd+7kQOqrVy5rwbhoqnqJx/44tsa+wt+vXBo8bx0wtm0pZFq0fXVWg/NyIIC1ZXMgyNoF2BMde5KQV5Lw8uZCUmSVmoqoJl6X2i0a0bGqGmcv1Tq/nO/ctYubsNI27kbq+raDGD217ayM8ed9je5qHHuqtSapqCG0b878O9HTtA0gv44t9fxcyjZqnoKn96aX1exw6Z1MjfzFtNqOaWTo5Ule89sCyrYwe4//VtLNmwN6OX0wXd0PBUhRueWcsv560hichbpUo3NRr8kP99dCU/fPgt6pI+ehl6+x0oisCydLQsGvuSwYN07gcQpqllLfyQD1VVsBImz6/fi5Uwyy4mols6l92+iLte3sTjy3fy0VsX0pDyc5aHU02Nz9y5hN89s4Z7XtnC+TfOJ0B064eiqdy5YGNeLfzNe5O8tb05p/2mqfGvV7YUZcNzq3ehqiKrZrmiCFRV4akClZj+OH8DQRehdiHAihlccNNL3PHSBv6ycCPn3zgf0zKy6pcLITAtnQ/8bj73vbaVB9/Yxgdvegnd0MrWUrcqLBZsagBDwzBKc/CaphCLlf+AkZSPdO4DDL3EH1OxKIqgqipGVZVV0n6GoXL/61v5wt9e4x+LN6Pr+z9NoygCPwg7tdwhU5zkqRU7c0aJpq51q0Wa9kPmr9ndrX0QZcroFWLdrpacjs/QVHa3FLcKNYqgMelldbqKIqhrTuMX0J7ZuLcNVd23v6oq1DWn2FK/b6hmV0uaDbtbsj5QNU1hzc4WGrqoTrakfd7c0lhWxC1E5jx95i9LuObRtwhEaS5CMXQqKqySg4hiUTWl23mT5EY69wGG7/VfsY6GpiQNTdmHCnLheQHvnzWOr5x1MB+aPb4sJcMwjDA0hQm1sW7bT5g2POfkoheEHD62qvNvRcDsSbXd2isCRlUWLh84ttrKmaXiBSFVseIdUsLQsh4riiKGJXIX4u5gVKXZTXwsCEJGVlrUxvdFvQlDZfLwRNZzEwQhU0cmsPR9P2FdFRw6trKsidoogumjKvjGOTZfPPMg1BLn5CLPp6E5hef1jwZPGET7VZfgQKSou9m27bOBC4E6IHIc50dZ2lwMXAN8xXGch9u3TQd+ArwKTAD2OI7z4z7q+5CkPye4c00m5sP3Q0Ta46q5k/Bdv+AEXyHSKY+/fPI4fvDgcna3uFx50hQmDYuRasmefeKlXP505Rz+9+EV7G5x+cypU6k01e7yvUHI5XMn84cX1+f83pEVJrMn1+YsOJJKebxv1jj+9NKGgjYcM7kWVSFrdB4EEYaAk6aP6FbftCeXnzC5W/WpKIJ0yuUfnzmBax93CIKI//fuQ3BdP2exDt8LuPOTx3HtYw5BGPH1dx8CYVR2Fk6qJc3H50wgDEufoO3Q5ekvBkICyGChYLaMbdtxYCkww3GctG3b9wI3OY4zr0ubqcAU4AfAdV2c+xxgnOM4D7T/vQK43HGcJV2/Q2bLDC7Ktdc0NQKlvc5mEOKlvbx54x0ZIUII1CginaVIhpUw+eW81dwxf0OvzxQBt1x2DHMm1eQusKGrtCE47f+eKViP9faPH8uc8dWk0zkKf5gae9yA8387n+YszvGEacO44xNzaG1K9XJWHeeGCFTCgplEpqURKioCEGGY9dyUgryXBxflZsvMBTY6jtMRWs0HzuvawHGc9Y7jPNNzR8dxFnc49i7f11rEd0qGMOm0j590cVvTuKn8jh0ypeuClIefdHM6Lzfp8vV3HcIvPjST6SMTndvnThvO3Z+dy3FTanM6dshEnAlN4YaPHo2WZ0Ly86dP57gpw3Dd3MdKp31qTY1Hv3IK580ci94+RjyywuSrZx/M7R+fQ7rNzRqFdpwbP+UWdOwA6ZSP15bGbUuX7dglQ4tihmVGAc1d/m5q31YStm1fADzuOM7Knp9VV8cGdUVzVVWoqSmtwPJgZiDb+8Ejx3P+keMJwghVFQRBiNmewhkvYpLvtINH8MiXT+E381bz+PIdnUMvc6cN5+rTp3Ps5Fp0RWBVF7a/IgY/v3Amv7r4KNJ+gKEpBEGErgiMROE5gneCgXxt+4OhbG8xzr0OqOzyd1X7tqKxbfsM4Azgq9k+b2wsv/DyO8lgf7UrlcFgbyZYiIgiKPXuGhvX+OkHZnD9RUfSlPKImxphEEIQ7FeR8LTI9CdVRPWmd5rBcG37ksFu78iRlTk/K2ZYZgEw2bbtjlDjJOAR27aH2bZdlWc/AGzbPg84B/gKMMa27blFfKdEUhZRFO23Fns67eMlXVqakmi+j6kI0q3pooZJsveFosrySSR9SUHn7jhOG3A18Bvbtn8CLG2fTP028HkA27aFbdvfAyYDF9u2fU779mOAfwInAM8ADwB2fxgikfQ1USTT7iSDF6kt0wcM9le7UjlQ7NU0BSEEFRUm9fVD3144cK5tB4PdXllDVXJAoSgZeYIoAt8PSh6eMS0NNI0212dnY4pqP2JUVQzP9QumbWZD11WEyEzuSiVIyduFdO6SIYOiCFRTR9dVlm9txNRV7NGVpNJe3jTIrhgxg/X1SX7yyGss3lDfuX36yAq+dvbBnHbICNKt2dMYe6KbGoaps6U+yd5Wl4NHV2AIQeSVvxhMIimEdO6SdwRVzSxK8v3iIlmlfRFTrshXUQSxCosb5q3mzws2kGqXcRhZYfKD9x/OqQePIN1aQIPd1HhrZzOX/3FRr4VMa3e18MW/v8b/nHcYH549vvsK2SwYls6WpjRfvm0Ra3e1ZPoo4NwjxnDth2YB5HTwHfo3xU7CqqqCEBR9LiUHBlJbRvK2IkQmOlYsg6YAElWxgkJXuqUTq7CIDI1YhZVVelg1NG58Zg23PL+u07FDRnzrS/94jeXbm/JK/gJYlsHX7n4j7wrVa/6zkhCRV/5Y01TSEXzklgWdjh0yImmPvrmDz961BC1LX4QQGHEDI26ixQ2MhJlX4VFRBFbCxFNVUkLBqsjfXnJgIZ37AURH9FsqesxgxIgKjFhu3fVi0UydZ9fs5pifPMVp1z3L5//6KmYePXfD0nljWxNH/++THPOTp/jhwyt69UNRBIap8eccujBRBL96ajXkccimqfHKxr1sz6HB3kEQRty5YANKDonijJEKtzy/ltYcpf3mr9nDjsZULzldzdR44I3tHPW/T3L0j5/kry9vQjVzP5BUU+fWF9dz/M/mMffnT/PreavRrL65RsOHV+z3w0KqNg4MpHMfYCh9UBAjG0JATW2CRMmSvxrbm1LY33uMDfXJXgUmSkXTVW54ajVB+5DDs6t2saMplTN6jxSFXz21mnT7kMO/lmwh6QXdImdVVVhb15LTmQIsWr+XijyRu6oqLNlYn/Pzrry5tZF8GZJCESxctzfvMZ5btQu1R8GPirjBdU9mRMDCCG6Yt5qKPNroiZjOLc+t6/z7Ty9txGyfvC0HVVeZ9aMneGjpdowSr7cWMxg2rKJs3f9cCCW7jr6kN9K5DzD6676NImhuc0mVqPIXhiGjqyzOmTGacTVW2YtxgiBiRBd5XkVAdUzPOUEZRREjKvZFo6amEDd7y+1a+SJpwNIzpevyYRapg64XeABHUeH+VFga9OiR64eMrNh3boYljM6HYDbCMKK6i0RwhakhRPnKiboqOPeIMcwcX03JqUFBQFOb22+LthQhF4QVi3TuA4z+nBTzki5+ieJSvh+iBCE/ev/h6GHhmp4FCQKuv2gWx06uZdKwONd+aBaaIPdEaRDyswtncuaho5g5vpqbLp2N53WXwfW8gAm1MSYNy60R8v5Z42jKU9Ta8wLec8SYokx4zxFj0PM8hJUo4kOzx+f8XFcF7zliLG6PNw037fH7S4/hhGnDmDOlllsvP4ZkMvfEbTLpcvOlszlsbCWHjK7g95fOJtkH4mHpNpf/PtdmUrVFqsRVub4bkG5N95sDlovKikcuYuoDBvtCiFIp117D1EBT0VQFz/Xxi5D8DdV9EsFuFgemWxpLt7fwyT8t7qWzPqLC4NEvn4JR4OFkVlh89i9LulV+6snICpPnv3kGrU3JnBGyoghilRYfuWUBy7Y29fr8G+fYXDpnIl4Wx22YGpGqAAIlCArqqRumhqJnIvbAC7Kem1KQ9/LgIt8iJunc+4DBfoOUykC114gbbKpP8qunVjN/zW5MXeF9M8fytXcdgh5FBYuV6LpKpGt8+OaXWLurtzJ1dUznn585gbEVRkEnqusqumVw6wtr+efiLextdZkxroqrT5/OidOGk2pND8jCEwP12vYXg91e6dz7mcF+g5TKQLbXNDNRfoWlE0XQkvIQflD0cJJhauiWzr9f3cpdCzeytT5JdVznQ7PH84kTp6KE2d8csqGqCoquYlk6uqrQlvaI/BB3P1a5vl0M5GvbHwx2e6X8gOSAoWMYo77AIqNcuGkf3wt4/4zRfODIcRi6QhBGuGmfMO3hliAf0CE3UOzqWImkL5HOXSLpQRhG7dG5h6sIqqpitEoHLRlkSOcuGXIIQXsOebRf2UcdwynxmE7aD4mEQIvpCH//soXU9sngIAhlGp/kbUM6d8mQQjd1LEtnc30bpqYwvCqGn/Zwi8zvNwwVzTK48ek1/GPxJurbPExN4f1HjuPb5x6KLgSeW/yxhK6RDkLqWlwmDY9n6sGmfenkJf2OdO6SdwTDyKTveV5QlKPrkM31vCBnlokRN1i4vp7v3b+MXS2ZnPZDx1Ty+8tmM8zUCjp4RREYMYOP3fYyr21u6Nye9kP+tWQLL6zexaNfPgVdVwtG8IahgqHzxb+9yvOrdwMQN1SuPn06V544hVRLKuekaocsQcbWvF+T0epprw3rurnPjeTAQy5ikhRECNHujPtm+ayZMNnYlGbBpgZiFRZ6gdWcRtxgjxuyYncr8Uor69J2w1DZ1eLy+b+92unYAVbuaOYjtyzEsIyC/dcNjYeXbu/m2LuysynN9U847Xno+VEMnavvWtLp2AHa3IDrn1jFkyt2omfRjBEio2y5ocllXWMaq8IqKBwWq4ixdGcLr21vJl6Zv30pGIYmNWIGOdK5DzD6UzdDNzV0o/SXtVilxdqGFPFKi3L9u2XpLFi3h/N/N5/P3vUqX/rHayh5+mRZOst3NHP2r57jktte5n8eWIbI0j5UFG57cX3W5fq7mtM84+wsqIsTKgp/X7Qpb5v7X99GZR6hM8hUcGpKeby0NvtiqD+8uB4li9SBZmrc/8Y2PvC7+Vxw00v8ZeHGvOdGMTT+OH89l92+iI/fsZjfPL06b/tiMeIGW1rcnAqc+dB1FaOA+mY5SF2Z4pHOfYChKP13SWqqYugl/vBUVSHlBVx86wKaU37JP/aeBMCCLqJar26sx8gXuYuMimLHaMOiDXvRsvQhBLY1JHMeZsPutoKOQVUEe1rzp1C2uQF+EOV9yCmKYEcedckt9Un0LG8fQQQvrd0X6S9ct5d808ERmfPRwSsb6in76QtUxQ0u+cNClmysL10ATFOoriz9oVA0QgqHFYt07gMM3++/Cj319a0FC1b0JAhCNOCtH5+LpRZfXCMXShRx6fGTqLIyEeanTp6aV8wsDCI+dtxERlaYCAGfO3V61vFuFTh6Uk3O4xw/dVjBEnd+EDJ9ZEXeNmOqrHZxrtxtgiDkoNGV6DmGNY6aWEMqy6SsSsQnT5pK3FAxNYWrTpmKmueLlCjT3lAVNEXwyZOnIvpgorahOcWS772LoybW9NK/KUSQ9mloaOu3coKhzDgqGrlCtQ8Y7KvciqGrQyvXXt3SiVk6KT9EiSLcZH4VQd3USCRMwjAi5fp4SbeXc9U0BWHqnHn9c+ztEX0fM7mWuz51HC2NuSN7yKxuXbGrjY/etjBnm2+cY3PpsRMKLkzSYwY3PL2GO3pozGuK4L7Pn8iUaotUVo0cnUTcIIqgLeVl1Z/pihEzMrIJZAKDQhWiCtFxbQs9wIYKg/23K1eoSsqmL3/oXsrDT3sIIYqKwry0T0Paz+twfD9EV0Me+uJJXPPoSp5eWUfMULnw6PF89exDcAs4Scisbj1ifBWfOnkqt7+4vtfnJx00nCtPnEKyJX9BD4DQ9fjGOTbja2P8ecEG6prSzJ5UyzfPtZlUEyOVoz9eyqOhhAVTbtLFSwlEH0vhHgiOfagjI/c+YLA//UtlINvboSBZGdMJo4i2pEfkByXVajXjJs7OZm59YT3rd7cwstLksuMnc5o9ErfNLXohk6IINEPDMDUMTaEl5SOCgHSJMrpvJwP52vYHg91eGblLDhjc9rHsZi/z31LnCMIwItmS4vDRFfziQzNRhCCKInQByeYkYQmH65AxCDyftKLg+4Xz1iWSvkI6d8mQQjc0zJjOtoYUcUOlosIgSHslyQYYcYNlO1r45ZOrWLa1kQnDYlx1yjTeffjokgpRqKpANQ0iATubUkwaliCd8vDSUqdG0v9I5y5521EUgWpoBBFoUdQZbedCCIFmaoQRKHkKbhimRnMQ8cEbXmDDnsyr9qkHj+Dmy44hDNNFZXCYlsab25q4/I+L6PDhb21v5r/ufoNvnmNzybETCIscEzcTFv/3+Er+snATQRhRE9f5/aWzOWxkRU4H36FrEwGRFxTss6Yp0JEzX8Lwk2ToI1MhDyAURezXKlPdUNHjxn4tgMqGlTC5a/FmfvzoW0S6WnCFqpkweeDN7Vz31Go0S8+Zey10ja/98/VOxw7w/Ord3PrCOpQC39GBomv8/DGHbMH5Tc+uxTSLW6lrmhqvbqznTy9t7FxY1dDmcfVfX8W09Kzp6EIIrITJH17awO+fX0+swiy4QtVMmPz66TVc99RqzLjZJ6tKdT1zvbX9LIYuV7YODKRzH2Co/VQ1HiBeYWEUWF2ZjUSFxbf+vYyqSqvsPmiawt5Wl+ueWMUDr2/jlufWoeSxuWMR1Q8eXMHfFm3izy9tyPmQiRkqizfU99r+rLOLqMiHWsxQeXNrY9bPWtI+W+uTRTmvQAgeX7Gz1/aGNo+Ne1rRsqxQ1XWF5duauPGZtdz6wjoWrNuT98Gn6yovr9vDnxds5K8vb+LplXXoeh88gHWV659cjaprJS9GMiyd2tpEvy00EvsZoByISOc+wIj6cYFGOuUSeqVnarQmPf7vw7NoLiKdsBBhGDGi0mT6yASGqnC6PRKRx+QwDKmM6RwxvoqEoXK6PYowx1BFGEaMrjJ7bZ80LF60oFYYwrBE9gegEFCbMIqaFFWB6SMSvbYrAkZXWVnH7YMg4uDRFYyrthhZaTJzfHXegtBBEHHE+GqGJQyqYzpHT6rpk8VDURDyjXNs1DyFy3MR+gHNLf1XIJsIKY5WJDIVsg8Y7OlUxaCqSucPvexFTIZGPJFZcdqadPGS+cewDUPFiptoqqC51c25gEi3dJ5evZv/uvuNzm0JQ+XhL5/CMF0pOLbfcYw/LtzIb+at6fXZ6YeM5IaPHoVbxCpfVRWYCYtzb3iezXv3LZ66Yu5k/t/Zh+C2ZT+GbmgkKkwE0NKaLlz3tX2BF0Brq1v2ZG3HtVVVkffBMlQY7L9dWUO1nxnsN0ipDFR7hQAjbrKlIcndr2yhKqZx+QmT0aHoUneqKrASFt9/cDn/fnVL59j7CdOGccvlx4LrFb0k3zA0NEvnn4s3s25XC+86fDTHThlGujU1YB3nQL22/cVgt1c6935msN8gpdIX9uq6iqKIzpqnfdXeMFSMmMG6Xa3omsKk2hitzamShglMUwNdo9X1Wba1icnD44yttlDCkLYCwmI9iVWY7G3zaEx6TB4eJ/SCgtF4Rs9dFPWmkWmf0cYv9lzmQ97Lg4uyFzHZtn02cCFQB0SO4/woS5uLgWuArziO83Ap+0reHjRNIYqikqNGXVeJNBXhB/tVZq7X8QyNNIL1u1s5bHRFQT0Uw1DxVZWN9Ukm18by6q3olsFlty9iycbMxOqXzzqIK0+YXHT6ImQmQ697bCVLNtYzcVicva0ufhDy108fX/QxICNXvGhDPZ/88ysAjKwwee6bpxPkKVCimxqNXogXhAyzNPwihmXqkj5+EDG20uiTYtyqqqAYGlEQ4hf5gOlACIGmKX1yn0jKo+CEqm3bceBm4GuO4/wQmGXb9lk92kwl47w3l7qvpDv9mkamaful9x2pCl4YFVWkohhCIbhh3mo++afFVMYKZ++EisKNz6zhyj8tJlFAstgyVF7btC9jZtH6vUVnynQgFMGbWxtZuaOZJ1fsZMnGet7c2kiixNTAEJjfRdN9V0uauqZU3kySUAi+/8ByvnnvUkQR5ztSFL53/7Ki2xeDYai4YYS2H/eKZWkly0qXgqL0iarxAUExd8NcYKPjOB0zQPOB87o2cBxnveM4z+zPvpLu9OdYbJD28Pcjsou8AFNk/tsXqFHId887jKf/32k0FTE5KdqzN5762qm05JiI7KA15XPezLGdf1949HiUUocew4jT7ZHdNp168EiaCkz89kQBPnjUOMz2VM9Dx1QyuiqWNwNFCSNu/NjR/OkTc6CIBUkiDLn1imO585PHEfWRXHQq5WOKzP1SKsmkR6pEWelSCEMpalYsxTyaRwHNXf5uat9WDEXtW10dG9S5q6qqUFMTf6e70f+0R2R9Ze+IjrTCIiM9vSNyLhDtX/vhWXzqlKlUWTqjqywMRRT9HR1cfdpBRBE8vbKOw8dV8T/nHU5cV4jrpdk9Tdd44ZtnsLm+jcPGVqEJqK4u8hiaASWtS1ALnpuCR1AVqqtjZR1jMDGUf7vFOPc6oLLL31Xt24qhqH0bC+hsD3QG+6RMqQx0e4WAqVVmRhO9Kcn+9FRVBZ+aO5lPnTQFVVEI0h71LfsXGeuawkE1FsnmFK0DvNDEQL+2fc1gt3fkyMqcnxUzLLMAmGzbdsfqkJOAR2zbHmbbdtX+7FvEd0ok+00UgeuWN/kbBBFeysNrc7HKnCD0/RDXzT2JKpH0BwWdu+M4bcDVwG9s2/4JsNRxnHnAt4HPA9i2LWzb/h4wGbjYtu1zCuwrkUgkkn5E5rn3AYP91a5UpL1DlwPJVhj89ubLc5faMhKJRDIEkc5dIpFIhiDSuUskEskQZECMuUskEomkb5GRu0QikQxBpHOXSCSSIYh07hKJRDIE6ZuKxwcIheSLbdv+BPA5INW+6XbHcf7ytnayj7BtewzwE+BIx3HmZPlcISPx3AxMIWPrwre1k31IEfaeDvwaaGjf9IjjOP/3dvWvL7FtezoZW18FJgB7HMf5cY82FnAdsBU4GPi54zir3u6+9gVF2vsJhshvtwPp3Iuki3zxDMdx0rZt32vb9llZVtx+1HGcDW9/D/uck4EHgKNyfP4RoMpxnG/btj0MWGjb9mGO4wxWIe9C9gJ81XGcZ9+W3vQvw4B/OI7zAIBt2yts237EcZwlXdp8FdjkOM61tm3PBG4HTnn7u9onFGMvDJ3fLiCdeynkki/u6dy/aNv2DiAO3Og4zt63sY99huM4/2qPVnNxHvBEe9u9tm2ngBnA0rehe31OEfYCXG7b9rFkBPBucxxnc4H2AxLHcRb32KQArT22nQf8d3v7N23bPtK27SrHcZrejj72JUXaC0Pkt9uBHHMvnmLki58DfuE4znXAK8A9b1Pf3gnKkYIejKwA/rf92v4TeLJ9aGpQY9v2BcDjjuOs7PHRkLy+eewdcr9dGbkXT0H5Ysdx1nf582ngQdu21UE8VJGPcqSgBx2O49R1+fdy27ZrgInAxnesU2Vi2/YZwBlkhmB6MuSubz57h+Jvd9BHHm8jBaWPbdv+mW3bHQ/Mg4ENg/nm6Ilt2wnb7ixR9AiZoSrax9wtYPk71bf+oKu9tm13zC102GsAO9/J/pWDbdvnAecAXwHG2LY9t4eMd9frOxN4YzAOyXRQyN6h+NuVkXuROI7TZtt2h3zxLtrli23bvhbYC/wc2AH83rbt9cBM4LJ3rsflYdv2acDlwNh2OefrgU+QsetzwN3A0bZt/wCYBFwxmH8MRdi7HrjBtu0VwOFk7E3lONyAxrbtY8gMLb0CPAMkgN8BF7DvXr4BuK79XBwEfOqd6W35FGnvkPntdiDlByQSiWQIIodlJBKJZAginbtEIpEMQaRzl0gkkiGIdO4SiUQyBJHOXSKRSIYg0rlLJBLJEEQ6d4lEIhmCSOcukUgkQ5D/Dwe9v9mCto7EAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pyplot.scatter(*grid_nodes, s=grid_weights*6e3)\n", "pyplot.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a sidenote: Even though embedding the density function into the weights is\n", "what is the most convenient think to do in uncertainty quantification, it\n", "might not be what is needed always. So if one do not want the embedding, it\n", "is possible to retrieve the classical unembedded scheme by replacing the\n", "probability density function as the second argument of\n", "[chaospy.generate_quadrature()](../../api/chaospy.generate_quadrature.rst)\n", "with an interval of interest. For example:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:18.339659Z", "iopub.status.busy": "2021-05-18T10:56:18.339341Z", "iopub.status.idle": "2021-05-18T10:56:18.348733Z", "shell.execute_reply": "2021-05-18T10:56:18.348426Z" } }, "outputs": [ { "data": { "text/plain": [ "(array([[-1. , -0.70710678, 0. , 0.70710678, 1. ]]),\n", " array([0.06666667, 0.53333333, 0.8 , 0.53333333, 0.06666667]))" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "interval = (-1, 1)\n", "chaospy.generate_quadrature(4, interval, rule=\"clenshaw_curtis\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Also note that for quadrature rules that do require a weighting function,\n", "passing an interval instead of an distribution will cause an error:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:18.350984Z", "iopub.status.busy": "2021-05-18T10:56:18.350666Z", "iopub.status.idle": "2021-05-18T10:56:18.421194Z", "shell.execute_reply": "2021-05-18T10:56:18.421430Z" } }, "outputs": [], "source": [ "from pytest import raises\n", "\n", "with raises(AssertionError):\n", " chaospy.generate_quadrature(4, interval, rule=\"gaussian\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Smolyak sparse-grid\n", "\n", "As the number of dimensions increases, the number of nodes and weights\n", "quickly grows out of hands, making it unfeasible to use quadrature\n", "integration. This is known as the curse of dimensionality, and except for\n", "Monte Carlo integration, there is really no way to completely guard against\n", "this problem. However there are a few ways to partially mitigate the problem,\n", "like Smolyak sparse-grid. Smolyak sparse-grid uses a rule over a combination\n", "of different quadrature orders and tailor it together into a new scheme. If\n", "the quadrature nodes are more or less nested between the different quadrature\n", "orders, as in the same nodes get reused a lot, then the Smolyak method can\n", "drastically reduce the quadrature nodes needed.\n", "\n", "To use Smolyak sparse-grid in `chaospy`, just pass the flag\n", "`sparse=True` to the\n", "[chaospy.generate_quadrature()](../../api/chaospy.generate_quadrature.rst)\n", "function. For example:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:18.423609Z", "iopub.status.busy": "2021-05-18T10:56:18.423345Z", "iopub.status.idle": "2021-05-18T10:56:18.434162Z", "shell.execute_reply": "2021-05-18T10:56:18.433902Z" }, "scrolled": true }, "outputs": [], "source": [ "sparse_nodes, sparse_weights = chaospy.generate_quadrature(\n", " 3, joint, rule=[\"genz_keister_24\", \"clenshaw_curtis\"], sparse=True)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:56:18.436567Z", "iopub.status.busy": "2021-05-18T10:56:18.436300Z", "iopub.status.idle": "2021-05-18T10:56:18.492627Z", "shell.execute_reply": "2021-05-18T10:56:18.492846Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD0CAYAAABgk2Y8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxBUlEQVR4nO3daYBcVZ338e/dq6rXdHYSshDCCZAEErbEICggiOgjBBWZweVRZ9xHmWdm1BkfR+ZxHMdxZhxnVFxwRllcEBcEBdkVkrBnJweyJ2TpbL13Vd3tedGd0Ft1V6e6u5b+f1511V36nKpbvzp17rn3GHEcI4QQorKYxS6AEEKIkSfhLoQQFUjCXQghKpCEuxBCVCAJdyGEqEAS7kIIUYHsYhcA4NCh1rIej1lXl6S5ubPYxRgzUt/KNZ7qCuVf38mTa4xcy6TlPgIMI+frW5GkvpVrPNUVKru+Eu5CCFGBSqJbRohSEoYh27dvZcuWjXR2dpJKJZg1ax5KnYXnecUunhB5kXAXFcd1LWzbIo7B9wOCIMp72507t/Hwww8QxxG+7594ft++faxe/QfOP385S5dekPfPedM08DwbwzAIw4hMJhh2fYQ4GRLuomI4joWTcNjb1Mljej9Jx+KaRdNxXYMgnSWKBj9vv23bKzzyyO8Igv4BfPy5559fg+9nWLbs9YPuyzDA9hwc1+Z3G/fT2JJh6ewJnDOzjkzax5eQF6NMwl1UBMexsBMOH7r9eVZvO3Li+Vt+s4n3v24O/+dKRbotnTPgM5lMzmDvKQgC1q9/kblz5zN16rSc67kpj99vPsjnf72RtP/aL4e5k6q444MXUu3ZZCXgxSjKK9yVUlcAK4FGINZa39Jn+WeAacB+4HzgC1rrLd3LbgKWACGwTWv9nZErvihHQeCzadMG0ulOzjjjTCZMaCh4n4Zrc/PP1vUKdoAohh88tZPZE1Ncu2g62bQ/4PZab8r7f4VhyNq1z3HVVW8dcLnn2ew40sFf37Oevjdd3XG4nRu+u4aHb74UPxv0Wz5cnZ2dbNq0Hog566xFpFJVhe1QVIwhR8sopVLArcDNWusvAouVUpf3Wa0a+Eut9VeBe4B/6d52JvBXwF9prf8G+JBSav4Ill+MINe18BLOqP6POI659957WLPmjzz//NP8/Od30tzcVNA+bdskHUQ8/NLBnOt87487cL3cbZmNG9cN2Wo/Lo5jduzYShiGAy4PTZNvPrY1Z3DvPdbJqm2H8bzCXmvf97n77jt47rnVPPfcGn72szvIZjMF7XMormfjDfI6itKRz1DI5cAurfXxo+Yp4JqeK2it/6/W+vihbAJt3X9fBTzfY9lq4OrCilzZLLs4o1MNwyBZ5bGrOU0yOXoBn06naWw8cCIYwzBi587tBe3Tskw27G0etBW891gnfhhjmgOfCO3sHN6FLIZhkM1mB1zmOhZr9zQNuv2aHUeJCxxiffhwI5lMhiiKiKKuE8AHDuwvbKeDsG2LwDRpi7oaAsVgWEbO91D0ls9X8BSgtcfjlu7n+lFKucD7gI8PZ9u6umRZX0xgWSb19aliF6NgfhizYHotRhwP2qospL7V1S6maRJF0Yl9TZkysaDXzzCgPjX4F5JtGniOSbI2OfBy2yIzjEZvHMc0NFQPODQyE8bUJh0aW3PvsD7lkEw4JAtovcfxROI46vE4Ytq0SQW9lkO9tyFQ7XW14FLlf8hXzGd3IPmEeyNQ0+NxbfdzvXQH+7eBv9Nab+ux7el9tt3ad9tyvvwXoL4+RVNTR7GLMWYKre+VV76Vhx66nyAImDdPMW3arIL2Zxhw9im1TK72ONQ2cKBecdZUOtMBfufAre1p02awffsr5DszWVVVNR0dAZ2d/btm3ITDO5bO4CsP6AG3NQ1453kzaWtND2uYZl+GkWDZstezatUTAJx//jJct7qg11KO5fIyeXJNzmX5hPtqYLZSyuvumlkBfEsp1QAEWuuW7n75bwJf01pvUkpdr7W+B3gQ+KRSyujumlkO/GehFRLlbc6c0/jQhz5BHMeYZuHdUHEM6UzAP61cyIfveIGwz4iYhiqXL7z1LIwwd5Cee+557Nq1Pa9+d9t2WLLk/Jy/NkM/4D3L5/Drdft4aX9rv+WfuOx0XNMgW0CwH7d48RIWLToXqOxL6cXwGfm0VJRSbwLeARwCfK31LUqprwJHtdZfUUr9AlgI7OvepEprfUH3tjfRNYImBF4eaLRMud84rNy//YerVOvrplx2HevkPx55hae2HibhWLxl4XRuftN8nBj8zMAjZaCrm+U3993D3r17Ic4dujFQXV3Nn9z4fhzHzV0W18LyHL73x+387Lm9HG3PcvYptXz40nmsOG0imY7MkOPui6FU39vRUu71HezGYXmF+2iTcC8vpVxfz7OJLJOqhEMcx7SnA4wgxPcHHtlynONYHE1n+fdvf486I41j9A/4IDbIYnHhldfzpkVzyOQYVnmcZZmYjoXn2TiWSUcmgDAim/ELHgI5Wkr5vR0N5V7fwcJdxjSJinL88v6mjoH71nOJLZMfrN7D/ZkFnGEdYpF9EM8IiDEwiIkw2BxM4aVgKlteOMxVS+YBg4d7GEaEYYQ/xJeAEKNBwl0IwHNtHtnSSITJlnAqW8Ip1BoZEkaAH5s0xQni7pHD6/Y2Extd940pxa4VIUDCXQigK6g7e3XdGLTECVpyZLc/AidDhRhNcj93IYCsHzJnYn6X7le5FlWenfewSSGKQcJdCIAw4r3LZ+e16rVLZtBZwidFhQAJdyEAyGYCrl44ndMmDd56r/FsPnnZ6SDdMqLESbgLQdc492xnlp99ZDlnTK0ecJ2GKpef/PkykqYx5NBKIYpNTqgK0S2bDXA9m3s/cTHP7jjK7Wt2caAlTW3CYeXSGVy9cBqZjJ/ztsFClBIJdyF6yGYCspmAc6ZVc/Z1CzGMrrsQGmFEe0taTqKKsiHdMkIMIJMJCNM+QWeWlGORSfsS7KKsSLgLIUQFknAXQogKJOEuhBAVSMJdCCEqkIS7EEJUIAl3IYSoQBLuQghRgSTchRCiAskVqmLM2baJ63YdekEQkc0OPSn1WPM8m9AwMC2TziDCSzglOT2e41g4jkUcxwRBJPe8ESdIuIsxY5oGdsIhiOGnL+6jLRNwyfxJnDG1Bj+dJZstfjCZpoFX5fHS/la+/+R2dhxuZ1K1x00XzeKNC6aQ6ciWRIDatomdcGnq9PnNs3sghqsXTmNyTQK/M0sgd60c92SC7BFQ7pPsDtfJ1NcwIFmd4FuPb+PbT2yj5+x058ys40cfuBD8oKgBbxiQqE7wzw9s4fY1u/stX3ZaAz943wWk2zOEYfHC07JMElUen7lnPb9Zv7/XsivPmsrXbzj3pMsox3J5GWyCbOlzF2PC9Rx+/9JBvvl472CHrjlJP3rnC5iuU5zCdXM9m2d2HB0w2AHWbD/Ktx7fhulYY1yy3kzH4huPvNIv2AF+v/kgX/ndFowil1EUn4S7GBOWa/GdJ7bnXL5q2xGaOnycIoZSZJp878kdg65z1zO7SSaK9yVkGJBIONz1zMBfQAB3P78Xz7MxzZyNOjEOSLiLMZF0bLYcaB10nY37mosaSJ5joYco49H2LB3ZsGjlNE2Tw20ZWtK5T0J3+iH7mjol3Mc5CXcxJqI4pjYx+Pn7iVVuUUejhFFMXXLwVrllGiQdq2jljOOYmoSNMURu1yackhvZI8aWhLsYE20dWa4/b2bO5VNqPBafWo/vF29YZOCHXL90xqDrXLZgCmk/KNq93aMohgguPn1SznXOnz2BpGMW9aSvKD4JdzEm4iDk5ivOGHACascy+No7zyHdWdxx5JEf8r7XzeHUhuSAy1OuxeeuXoBZ7NAMAr583UImpPr/yqhN2Pzz9YuISmC4piguGQo5Asp9ONVwnWx9HdfG9mx+uHoXv3zxVTqzIRfMaeCTl81jUsol25kdhdIOj+vahLbJF369iQc3HcAPuw7N5adN5B/efjZTUk5JzKHqeDahafLtx7fx0EsHieKYyxZM4eNvOB3XAP8kyyjHcnkZbCikhPsIKPcDZLgKqa9lmViOhWlbGAb4QYgZxiV1larjWGBbWLZJY0ua+pSLYxoQhGQypVNO2zYxbAu7e4SR74cYYWFXqcqxXF4GC3e5QlWMqTCMuvuCi9/6zcX3Q/BDQtOgwTGp9uySDIAgiCCITrqVLiqbhLsQOURR3HUCU4gylFe4K6WuAFYCjUCstb5lgHVuAL4MfEprfV+P5/8amAMcBuYDH9RadxZedCGEELkMOVpGKZUCbgVu1lp/EVislLq8zzpz6Qr+PX2enwZ8Dvik1vrvgSq6viSEEEKMonyGQi4HdmmtM92PnwKu6bmC1nqH1vqxAbbtALJAbffjamDTSZZVCCFEnvLplpkC9Lwmu6X7uSFprVu6u2V+qpTaD+wFtvZdr64uiTHUJXclzLJM6utTxS7GmJH6Vq7xVFeo7PrmE+6NQE2Px7Xdzw1JKXUu8NfAUq11oJT6V+ALwN/0XK+5uby74Mt9ONVwSX0r13iqK5R/fSdPrsm5LJ9umdXAbKWU1/14BXC/UqpBKVU7yHYAM4CjWuvjg4P3A4k8/qcQQogCDNly11p3KKU+CnxDKXUIWK+1fkQp9VXgKPAVpZQB/B0wG7hBKeVrrR8EHgDe0t1ibwIWAp8enaoIIYQ4Tq5QHQHl/tNuuKS+lWs81RXKv74yE5MQQowzEu5CCFGBJNyFEKICSbgLIUQFknAXQogKJOEuhBAVSG75K8QADANM08Q0u0aamaYht/8VZUXCXYgeHMcitkySnsOhtgztnQHJMGZSdeLETEelNGuUELlIuAtB15R1dsKhuTPge09s5Z4XXqWtx5R6jmVw5VnT+PAlp3H6lGr8dJZsViahFqVLwl2Me65rYSdc/vJn63hw04EB1/HDmPs37Of+DftZOKOW/3n/hTiegV9Cc6oK0ZOcUBXjmm1bWJ7Ljd9bkzPY+9r4agtv+cYf6YzB8aR9JEqThLsY17yUy8fveoH1e5uHtV1ja4Z3f3cNjuecOOkqRCmRcBfjlufZvHywlSdePnRS2+860sHPntuD7UrrXZQeCXcxbkWmyXf+sL2gffz3UztxpWtGlCAJdzEuWZYJpsFDmw8WtJ/dRzvY8GozngS8KDES7mJcsiyT9XubCEbgwqTH9SEM6XcXJUbCXYxLpmnQ1OGPyL5a0gGhXLwqSoyEuxiX4jjGs0fm8PdsE2m3i1Ij4S7GpTiOOaU+OSL7mlGfxJJ0FyVGwl2MS9lsyOlTqpk5obCAt02DlUtnyK0IRMmRcBfjVjYT8N7lswvax5VnT4UoJgyjESqVECNDwl2MW6EfcOMFs0g4J/8x+Oil8zAjCXZReiTcxbgVhjFhEPLtPz2PkxnJ+Okr5jOnIUVGbh4mSpCEuxjX/LTP0lPr+M8bl+AM46zoJy47nT+7eC6Zjuwolk6IkyfhLsa9bEeWi0+byH2fvJg3L5yGNUgz/sK5Dfzwf1/Ah18/l3R7hjiWAe6iNBmlcHAeOtRa/EIUoL4+RVNTR7GLMWYqtb6eZxOZJrFhcMeaXWx4tZnWdEDStZg9McX/ft0cJqRc4iAkkx6ZC6BKTaW+t7mUe30nT67J2RKRG2II0e1437llmbz3wlMJ4lMxjK7HcRhBGJFuSxe5lELkR8JdiD7CMOo1tLGmPkVTe6aIJRJi+KTPXQghKpCEuxBCVCAJdyGEqEAS7kIIUYHyOqGqlLoCWAk0ArHW+pYB1rkB+DLwKa31fT2eV8CNQCdwKfBFrfUzI1B2IYQQOQzZcldKpYBbgZu11l8EFiulLu+zzly6gn9Pn+ct4N+Af9Ba/zPwQWDHyBRdCCFELvl0yywHdmmtj48Fewq4pucKWusdWuvHBtj2AsAAPqmU+hzwNuBwAeUVQgiRh3y6ZaYArT0et3Q/l4/ZdH053Ki1blZK3QFkgf/puVJdXRLDKN/ZDizLpL4+VexijBmpb+UaT3WFyq5vPuHeCNT0eFzb/Vw+WoAtWuvm7sdPAm+gT7g3N3fmubvSVO6XMA+X1Ldyjae6QvnXd/LkmpzL8umWWQ3MVkp53Y9XAPcrpRqUUrVDbPs0MLG77x26WvIv5/E/hRBCFGDIcNdadwAfBb6hlPoSsF5r/QjwWeBjAEopQyn1ebrC+wal1FXd2x4FPgN8XSn1BWAy8O+jUhMhhBAnyF0hR0C5/7QbLqlv5RpPdYXyr+9gd4WUi5iEEKICSbgLIUQFknAXQogKJOEuhBAVSMJdCCEqkMzEJEQPURSxe/cO1q59nmPHjhIEPo7jkExWsXjxEk4/XeE4TrGLKcSQJNxFRbFtE9uxCOOumxqZBmQzAVE0+GjbKIpYu/Z51q59ljAM8f3XJsD2fZ+Ojg6efPIx/vjHx1iw4CyWLXs9rusOuk/DAM9ziIAIsA2DMAjx/bDwigoxBAl3URE8zyayTMIY7nhuD/ua0tiWwZnTanjrOaeQzYbEfkAQRP22DQKf3/3uXvbvf5UgCHL+j+OBv3nzBvbs2cW1176LqqrqfuuZpoHl2iQ8hye3HmL19qNk/JD6Kpfrl8xgcnWCKAjJpP1+2woxUuQiphFQ7hdCDNdY19eyTBzHwjC6Jq/OZnu3fJ2EQ0s25Jb7NvPYlkb6NtJTrsV1S2bw2asXEGUCstnXAjyKIu6//5fs2/cqYZg72PsyDIOamlre+c6b8DzvxPO2beKlPH7w1A5+uGoXh9r6T6y95NR6Pnf1AtTUavzOLD0/grZtYdtdp8LCMBrzVr4cy+VFLmISZcm2TZyUi+E5PPzKYe7b3Mir7T6p2iSO1/Wj00k47G3JcM03nuSRl/oHO0BHNuTOp3dz3bdWgWPhuq/9YF237nn27x9esAPEcUxbWyuPPvrgiecsqyvYP37Xi3zt9y8PGOwAL+5p4t3fW8Oj+hBOsqtrx3Es3CqPDuB3Ww7x2y2NtIQxXnUCx7UG3I8Qg5FuGVGSHMfCSbh87pcbuH/DfsIeqb1gWg1fv+FcZtR6BBjcdNvTtGaGDuetjW2877+f5cd/dhHZbNDdz/7coF0xgzl+8rWjo51UqgrTtfnyb1/iMT30TVOjGP7q5+u55yPLmT8pRTqCT931Ik9u7T3dwYVzG/jPG5fgeQZ+HnUU4jhpuYsT4jgminr3STuOhZ10qa5LYidd3BFoRYZhwO7dO9i27WXa2lr7LTcM8FIu7/vvZ7h33b5ewQ6w5UAr1397FZkIbntyB00d+fddr93TxKZ9LXiezZ49O0862HvauHEdlmVimAY/eXbP0Bt0C6OYHz+zm9iyuO5bq/oFO8AzO45y7TefwujRXdNTJpNhx45t7NixlXQ6XVA9jrOTLlW1SdyUh+f1bv8NdIyI0iQt93HEMCDXKZaNG9fy1FNPEEURs2bN4cor30pVVQIcm8/9YgPP7TzG0tn1/NPKxQWVIZPJ8Itf/Lg71A3iOOYtb3k7M2fOOrGO69ms2naE53Ydy7kfyzTwHJMfP7N72GX4zh+28y/XL2Lt2ud7jYo5GWEYsmHDWlZccim3r9lFMMSonL4WTK/l9tU72X00d7/v/uY03/3Ddv5sxRzocUK4ubmJX/ziJye+oEzTZOXKdzNhQsNJ1cWyTII45t8ffoUHNh3gtMlV/OO1i5iYsMmmfVat+gPr178IxJx55iIuvfTyASfZMYyu91UUl7TcxwnHtaipT1Fbn8Jxere+DxzYx6pVfyAMQ+I4Zu/e3axa9QSxbXHzz9bxu40HONSW4cFNB/mLH79IJjz5ltvq1X+gubkZ3/fx/SxB4PPAA/f2ag2GhsmPVu8adD9nTq9ly4FWjrRnh12Gx7Y0Upt0aWrK/eUxHNlsBj8KeWjzwWFv+9bFp3BXHl9QP3l2D1XJ3kMvH374d6TTnfh+Ft/Pksmk+f3v78uxh6GZjsV/PbqVH63ZRWNrhjXbj3LTbU+TSLi88soWNm1aRxxHxHHMyy9vZuPGtb23Nw2qapPU1L12TkQUj4T7MI32dIBuwmE0/kVkmnz6J2v56oOa2Or9th86dLBXSysMQ/bvf5WapMuTr/TuKnhq22GqCvjgNjYeJIp6jwAJw5DOztdarqZp8GrT4LNzVSdsWtMn16USRDF+FBEEIzMU0TRNQj97UuWZWOWyr2no7pQj7dmugfs9HDt2tF8LuaAvLNPgiVcO9Xpqf3Oaw+0ZDhzY16sLKwgC9u17tde6rmvzuG7k8n97gkRy8GsATpblDtw9NVLKeLbPfuTrdZhG++dmkA1ydp0Uwoxj/vHahdiWSZTtHWp1dfUYhgl0ha5hmEyYMJGOjM/CGbW8sLvpxLpnTa+lM3vyw/Pq6uo4cuRQv9cxkUie+DuKYhqqBg+HjmxI6iT7/w0DHMvEsmxg4BEtwxFFMZblkDyJ8rRlAhpSbs6RNcelXAuzT/BUVVWRzfbeLpWqGnYZjoujmEUz6tj4asuJ52qTNhOrPOrrJ2Db9omAtyyLhoaJvbb3/ZA3LpjC4pn1pEdpDH8cRP3OwYzo/iuoN0la7iVmqCspT5afCTD8gKAzi98nnE89dQ5KnYVlWdi2Q21tLa9//WVEfsjXbziXM6Z2Xahz+pRq/vPGJdh9U2YYXve6S/E8D9O0AAPbtrnkksuxrNeC0Ypj3nX+zEH3s7WxjbOm11J1EoF6wZwG2tI+1dW5558cDssysSyLZacNv6/7oZcOct3SGUOud+25M2jt0wV12WVXYdsOhmFiGF2v5WWXXTXsMpwQRHzu6jNZPq8rtKfWenzrT5aSyficffY5TJ8+A8uysW2bKVOmsWTJBb02D8OIjtY01Sb4oxTuo/X5qERyEdMIKPcLIY5ra2vF933q6uoxza7vfdezsb3urqIYshmfqqRbUH07Ojp45ZUt+H6WmTNnMW3aKb2WG4ZBVW2Cq//jj+w8MvD/sUyDJz/zRv7r0a3c+fTwTqp+/73nceGpdWxYv4Ennni4oJOqhmFy9tmLuOyyN5ExTZZ/5dFhtf7eed4MPnv1mVz2r0/Q3DlwOao9m4duvoQkcb+Lmpqbj7Ft2ysAzJkzr19rerjq61N0BiGuZYEB6XQWv7u7KY5jWlqagZja2vpR76IcC+X+2ZWLmEReqqtrmDCh4USwQ9d9WTpaOulo6aS9pXNExlqnUinOOWcp55+/rF+wQ1eIZNI+P//I6zj7lP5zsNcmbG7906XUuBYffcM8XCv/w/jUhiQXz59MNhMwb958+nVkD5NpGixevJQgiEg6FpcvmDKs7a88exqeafDzjyxn5oRkv+VTaz1+8ufLSFrGgFer1tVNYOnSC1m69MKCg/24TFuGtpYOWps6TgQ7dH3p1tXVU1c3oSKCvdJJy30ElPu3/3CNVX1d18ZJOLzS2Ma96/aRDSKWzqrnzQunkc4EBGkfN+myZucxPnbXC0P2xTZUufz64yuotQ2y3V9Sq1Y9wYYNawnD4Z9HMAyDqVOnsXLljUDXNQFWwuGdt65my4H+4/f7+sQb5/HhS+aRaU9juzZewuHZncd4dEsjMXDJ/EmsOH0S6c7smF3AJMdyeRms5S7hPgLK/QAZrrGur+fZxIZBTNdPzWwm6HVC1k25bD7Qyt/+ciM7DrcPuI9lpzXw9RvOJWEY+Jned3z8+c/voqnpGHE8vCGenufxrne9h5qa135dOK6F7Tn87S838tsN+wcc9z6p2uXmK87g7eecQqYj06sf2fNsou5WsUVMJjM6J9hzkWO5vEi4j7JyP0CGqxTr6yYcXM9m8/4W7lizmwMtaRzT4PQp1XxgxVzqkjZRNuh30zHoOgfwq1/9lNbWlrxa8IZh4Loub3/7u5g0aXK/5bZtYTgWmAa3r9nF87uOkfEj6lIO1507g0vVZNJpnyDjl9zojFJ8b0dTuddXwn2UlfsBMlylXF/Pswm6W74xYAHkcXdF38/y+OMPs337K4Ax4I3ETLNrVMqUKdO44oqre7XYB2JZJubxe8t3Xx3sGF0npUvgYzegUn5vR0O513ewcJdx7qKiZPr0Tefb0eI4Lm9601vo7OzkpZc2sG7dC3R2dmCaFnEcYds2CxYsZPHiJdTV1ee1zzCMCPtczVv4qHoh8iPhLkQPyWTyxOiTKIrIZrNMmlRLW9vwb3MgRDHJUEghcjBNk0QigW1LG0iUHwl3IYSoQBLuQghRgSTchRCiAkm4CyFEBZJwF0KICpTXMACl1BXASqARiLXWtwywzg3Al4FPaa3v67MsCTwN/F5r/VcFl1oIIcSghmy5K6VSwK3AzVrrLwKLlVKX91lnLl3Bn2t24C8BLxZWVCGEEPnKp1tmObBLa3384rqngGt6rqC13qG1fmygjZVS7+neZkchBRVCCJG/fLplpgA971/a0v3ckJRSZwFnaq3/Vim1ONd6dXXJsr4/tGWZ1Nenil2MMSP1rVzjqa5Q2fXNJ9wbgZ7zkdV2P5eP64C0UuqzwMWAq5T6tNb66z1Xam4efDLkUlfuNx8aLqlv5RpPdYXyr+/kybmniswn3FcDs5VSXnfXzArgW0qpBiDQWrfk2lBr/Y/H/1ZKJYDqvsEuhBBi5A3Z56617gA+CnxDKfUlYL3W+hHgs8DHAJRShlLq88Bs4AalVK9ZepVS1wOXAMuUUjeOcB2EEEL0IfdzHwHl/tNuuKS+lWs81RXKv74yQbYQQowzEu5CCFGBJNyFEKICySwEQvRhGOB5DhgGYRwTxjHJpEMmExBFZX16SIwjEu5CdLNtE8O2SHgOf9x6iBd2NdGaCUg6FvOnVvPWRdPJZAOMPCbcFqLYJNyFANyEQ2yZ/OCpHdz19G4ODzBn6hfv3cS1557CRy6dR13KJdsh86qK0iXhLsY9J+GwpyXDTbc9TVOHn3O9jmzIXc/s4e7n9/Kv7zyHN5wxmUx7Juf6QhSTnFAV45qbcDjQluUdt64aNNh78sOYv/jJWh5/+RBuyh3lEgpxciTcxbhlmgauZ3PTbU+T9qNhb/9/7l5HczrAcaxRKJ0QhZFwF+OW49rct37/gP3r+fDDmFuf2EZsS7iL0iPhLsYty7X5wVOFTTPwq7X78BwL0yzfW1aLyiThLsYlx7HY39zJS/tbh155EB3ZkPs27MfzZGyCKC0S7mJcMk2DnYdH5oZRrxxsQ65tEqVGwl2MS4Zh0J4JRmRfndkAuaRJlBoJdzEuxXFMXdIZkX3VJBzsMp4mUlQmCXcxLoVhxILpuacoG44ls+qhBOZFEKInCXcxLgVBRMqxeN28iQXtZ1K1yyXzJ5PJ5HcBlBBjRcJdjF9ByJ9fclpBu/iTi2aRzvjScBclR8JdjFvZbMCy0yayYNrJdc/Upxw+sGIucSCnU0XpkXAX41YcQ6Yjw50fuojpdYlhbZtwTO744EUYYUQQDP/WBUKMNgl3Ma5lsyFWGHHfJy/OuwU/scrlVx9bwcxaj2xa+tpFaZLL6sS452cDXNfiVx9fwZrtR/juH7azatuRfuudOb2GD6yYy1sXTyebCSTYRUmTcBeCrha873ewdHoNt/7pUtqzIfpAK02dPjUJm1kNKU6pTxJmAzpa0zLdnih5Eu5CdItjyGQCyAQkbJOlp9RgGAbJpEtbW5qOls5iF1GIvEm4CzGAIHjtRGki4cicqaLsyAlVIYSoQBLuQghRgSTchRCiAkm4CyFEBZJwF0KICiThLoQQFUjCXQghKlBe49yVUlcAK4FGINZa3zLAOjcAXwY+pbW+r/u5ecCXgBeAmcARrfU/jFDZhRBC5DBky10plQJuBW7WWn8RWKyUurzPOnPpCv49fTZvAH6itf4XrfWngHcrpc4bkZILIYTIKZ+W+3Jgl9Y60/34KeAa4JHjK2itdwA7lFJ/33NDrfWzffZlAu0nX1whhBD5yCfcpwCtPR63dD83LEqp64AHtdZb+i6rq0tilPEEw5ZlUl+fKnYxxozUt3KNp7pCZdc3n3BvBHre6Lq2+7m8KaXeCLwR+PRAy5uby/uGTPX1KZqaOopdjDEj9a1c46muUP71nTw59xwE+YyWWQ3MVkp53Y9XAPcrpRqUUrVDbayUuga4CvgUME0ptTyP/ymEEKIAQ4a71roD+CjwDaXUl4D1WutHgM8CHwNQShlKqc8Ds4EblFJXdT9/HvBTYBnwGPBrQI1GRYQQQrzGiEtg2vZDh1qLX4gClPtPu+EaT/U1DKirGz/1HU/vLZR/fSdPrsl5slLu5y7GlGGA69kEdB2TFhD6AWFYOt/vlmViORbJpEMcg2EaOEkHgqik7utuGAauZ+MDBl0fZj8byCxRApBwF2PIcSy8lMuqbUe4+/m9dGZDlp3WwJ9eNBszCPFLYE5Sx7Fwky7ff3I7d6zZTWNrhirX4rolM/jLKxWOaeBngmIXE8e18ZIO92/Yz283HCCO4aqzp/L2c08hk/ZLooyiuKRbZgSU+0+74TqZ+lqWQaIqwQd/9Byr+0w+XZ9yuPvDy5mctIsaSqZpkKxOcNNtz/DC7mP9lk+rTXD/X1yMFYRFbcE7jgWuzfXfXs22Q229ls1qSPHLj70OOwzJZodfRjmWy8tg3TJybxkxJkzH5rYnd/QLdoCmDp8P3/48bsIpQsleY7s2v3jx1QGDHeBAS5ov//YlYtsa45L1FtsW//dXG/sFO8Duox185p71RS+jKD4JdzEmUkmHO57elXP59sPtbGtsw3WLF0qmbXHX07sHXee+9ftJeTbFuubONA1cx+KBTQdyrvPolkYMw8A0y/fCQFE4CXcxJmzT5GBLZtB19hztKOqVyo5lcqAlPeg6mSCiPRsUrZyGYXCsI4s/yAnoKIZDrRkJ93FOwl2MiWwYMXNCctB15k6uppjngLJBOGQZU65FlWsXrZxRFNOQcvHs3B9d2zSYWpuQUTPjnIS7GBOZtM97l8/OuXzBtBpmNaRO6iTgiAkj3rssdxkBrj33FDoyAcX6DorjmIwf8rZzTsm5zpVnTyUMIwn3cU7CXYyJIBtw00WzedNZU/stm1Lj8b33nk9Q5KGQ2UzAWxZN55L5kwZcPmdiis+8eQEUe6y7H/L3bzuLRTPq+i1SU2v4p5WLISid8fiiOGQo5Ago9+FUw3Wy9bVtEy/lseVAKz99bg+d2ZAVp0/kf51zCtl0QDZT/HHutm3hVbnc8/xe/mfVLnYcbmNStccNF5zKn7/+NMJsgJ8t/hhy17W7rhnYepj7NuwniuAti6ZxyRmTyXRkTvoXkBzL5WWwoZAS7iOg3A+Q4Sq0vl6PK1Rto/SuqjRNA9u1sV2blGuRDSM6O30iPyQMo2IX74QTV/vGBjHgANmsX1CXkRzL5UVuPyBKSqbHhUql2HkQRTHZtE827dNBVwCUwtWzfcUxZNKvvZaDj0US4430uQshRAWScBdCiAok4S6EEBVIwl0IISqQhLsQQlQgCXchhKhAEu5CCFGBJNyFEKICyUVMQgzAcSwMyySMwY8iXNcmWwK3HRAiXxLuQvRg2xaWZ9OaCbjzqZ0caElTk3B453kzmTupCj/tl8S9ZYQYioS7EN0cx8JJuPzl3Wt5cNPBXst+uGonC2fU8oP3XYDn2WRlAmpR4qTPXQi6bsLlpVw+8MNn+wX7cRtfbeG6b63CcGzsQSbLEKIUyBEqBF13qlyz/ShP7zg66HqvNnXy3T9sx5AJqEWJk3AXAghNk+8/uT2vdX/8zG6SCWeUSyREYSTchQA8x+Kl/a15rXukPUt7NpAJqEVJk3AXgq57o1t9wtogxiXApP8EHZYhwS5Km4yWEQLI+CHL5jbwm/X7mGa2stA+wAyzhRgDg5ijcYoNwTR2hfXMnVSDY5lkS2j2KCH6knAXFcW2za6TnaZBHIMZxwR5TONnhhF/9vq5NL/0FLOsJiwiuhryXdtNMjq42NnJIttj6YVvJ8hjrHvXNHgOoWFgGBBHMWYUnfT8pkIMh4S7qBhOwiHA4PtPbufJrYfxbIu3n3sK71g6Ez/tD3qFaTYb8PKzj3Oa04wRDzxPqmNENBhp0psfI73gRgwjd6/m8cnAn3j5EHc+vZtjHVnOml7LRy6dx6Qql2xHtuD6CjGYvMJdKXUFsBJoBGKt9S0DrHMD8GXgU1rr+4azrRCFcjybV1sz3PCdNbT1uMDo+V3H+NHqXdzz0dfhxDG+P3Cr+dixo2zZvBkjHrxVbRJz7OhRtm59hfnz1cDrmAZelccn7nqRR7c0nnh+074W7nlhL//2rnN5w/yJ+J2lNy+rqBxDnlBVSqWAW4GbtdZfBBYrpS7vs85cusJ7z3C3FePP4cOH+OlPb+eHP/wu69a9MCL7dD2Hj9/5Qq9gP25rYxtf/u1LxFbuw339+heJooFb7H35vs+LLz6Tc7nt2ty3bl+vYD8uiuFvfr4ey7KwrMJPyu7evZM77/wBd9xxGzt2bC14f6Jy5DNaZjmwS2t9fHL1p4Breq6gtd6htX7sZLYVvRlFHIVhJxy8Kg9rkBAsVBRF/PrXd3PkyCHa29t4+ukn2bNnV0H7dF2L7Yfb2XmkI+c6967dR1XSIdfLu3PnNuIc3TEDOXLkML4/cMvbsExuX7M757bZMOIXL+7FdQvrFW1vb+OBB+6lubmJlpZmHnrotzQ3NxW0z8EYBrgpFzflFe04Lebno9zkc3RNAXoOAG7pfi4feW1bV5cs6zfNskzq61Mjsq8ohmINn44Ng5f2t7DwlNpBT0AWUt/29naC4LVQjOOYzs6Wgl4/wzA4uG/wMeqdfkjaj6irG/j/hOHw7hVjmhaplE1VVf/9BXHMwZb0oNvvPZbGdm3qvZO/GKql5RCWZREEXWXv+ruD+vpTTnqfg723hgFN6YCMHzC1PklchMFCYRT3G7JaiJH87JaafMK9Eajp8bi2+7l85LVtc3NnnrsrTfX1KZqacrcay4Xt2cyqT3DsWAfxIJ/cQuobx1BVVU1rawtxHGMYBvX1Uwp6/RzH4rRJVYOuMyHl4Nlmzv/jOC7p9OCB3FMUhaTTIb7ff39OymXe5GoaWzMDbNnlrOk1+JmA1vTJ97u7bk2vgI2imGSyrqDXcqj31vEcLNvk2LHyP96h/D+7kyfX5FyWz+/v1cBspZTX/XgFcL9SqkEpVXsy2+bxP0URBJmAIO0PGuyFMgyDlSvfzRlnnMmpp87hzW/+X0yZMrWgffp+yOQajyWn1udc508umk3nIEGq1JlYVv73i5kx41Qsa+C2kRl1DavMpS7p8OaF08gUeGfJRCLBypXvZu7c05kzZx7XXfcuUqnBv+QK5Wd8ggK+kMTYMfL5ICul3gS8AzgE+FrrW5RSXwWOaq2/opQygL8DPgg8CdyhtX4w17Z993/oUGtZXw1S7t/+w1WK9XVci4xh8o5bV7HnaO9fgpfMn8StN51Huj1DGA7cr97e3sYdd9xGGA49Bt22Hd785rcxa9acAZcbBiSqE3zjka1894+971dT7dnc/sELmTshiV+CIVmK7+1oKvf6Tp5ck7OPKq9wH20S7uWlVOvrejZe0uXhlw7y6JZGErbFO86bwYJptWQ6sgTB4MH94ovP8uyzq0/0YQ/Etu3uXxxvG/Q8kWkauCmXI+1Zbl+zm6PtWc6ZWcc7zz8VPxuUZLBD6b63o6Xc6yvhPsrK/QAZrlKur2EYuJ5NEHe1oK04zrv7I47jEwEP9GrFG4aBaVqcdtrpXHbZVXl34TiOBZZJFINlkNfVssVUyu/taCj3+g4W7nKFqqgocRyT6dEqHk6vtmEYLF16IWeccSabNq1H681ksxkcx2HWrLmcc85SGhomDas8vh9C94VTctMBMZYk3IXoo7q6hosuWsFFF60Ayr91J8YnueWvEEJUIAl3IYSoQCVxQlUIIcTIkpa7EEJUIAl3IYSoQBLuQghRgWQo5DAMNfGIUur9wEeA43eguk1rffuYFnKEKKWmAV8CztFaXzDAcpOuyVlagTl01XXNmBZyBOVR3zcAXweaup+6X2v9L2NVvpGklJpHV11fAGYCR7TW/9BnnQTwNeBVYD7wFa31y2Nd1pGQZ33fT4V8do+TcM9Tj4lHztZaZ5RS9yilLtdaP9Jn1XdrrXeOfQlH3MXAr4Fzcyx/F1Crtf6sUqoBWKOUOlNrXa7X6gxVX4BPa60fH5PSjK4G4Cda618DKKU2K6Xu11o/32OdTwO7tdZfVUotAm4DXj/2RR0R+dQXKuezC0i4D0euiUf6hvsnlFIHgBTwX1rro2NYxhGjtf55d2s1l2uA33eve1QplQbOBtaPQfFGXB71BXiPUup8um5d/T2t9Z4h1i9JWutn+zxlAu19nrsG+Nvu9Tcopc5RStVqrVvGoowjKc/6QoV8do+TPvf85TPxyBPAP2utvwY8B9w9RmUrhkImcSlHm4H/1/3e/hR4qLtrqqwppa4DHtRab+mzqCLf30HqW3GfXWm552/IiUe01jt6PHwUuFcpZZVxV8VgCpnEpexorRt7/L1JKVUPnAoUNkdgESml3gi8ka4umL4q7v0drL6V+Nkt+5bHGBpy0hKl1D8ppY5/Yc4HdpbzwdGXUqpKKTW5++H9dHVV0d3nngA2Fatso6FnfZVSx88tHK+vCxwsZvkKoZS6BrgK+BQwTSm1vM8EPD3f30XAunLskjluqPpW4mdXWu550lp3KKU+CnxDKXUIWK+1fuT4pCXAV4ADwLeVUjuARcBNxStxYZRSlwLvAaYrpT4P/Cvwfrrq9RHgZ8ASpdTfA7OA95bzhyGP+u4A/kMptRk4i6765j8vXwlRSp1HV9fSc8BjQBXwTeA6XjuW/wP4WvdrcTpdE/GUpTzrWzGf3ePk9gNCCFGBpFtGCCEqkIS7EEJUIAl3IYSoQBLuQghRgSTchRCiAkm4CyFEBZJwF0KICiThLoQQFej/AyahesEQQA5/AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "idx = sparse_weights > 0\n", "pyplot.scatter(*sparse_nodes[:, idx], s=sparse_weights[idx]*2e3)\n", "pyplot.scatter(*sparse_nodes[:, ~idx], s=-sparse_weights[~idx]*2e3, color=\"grey\")\n", "pyplot.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that in sparse-grid, the weights are often a mixture of positive and\n", "negative values, as illustrated here with the two different colors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note here that the flag `growth=True` was used. This is because a few\n", "quadrature rules that are nested only at particular orders. Passing the\n", "`growth=True` flag implies that increasing at the order which is nested. This\n", "affects the four quadrature rules: Clenshaw-Curtis, Fejér, Newton-Cotes and\n", "grid." ] } ], "metadata": { "jupytext": { "formats": "ipynb,py:percent" }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.5" } }, "nbformat": 4, "nbformat_minor": 4 }