{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Generalized polynomial chaos\n", "\n", "Generalized polynomial chaos is an advanced polynomial chaos method for\n", "dealing with problematic random variables. The problems it deals with include\n", "heavy tailed distributions (like Log-normal, Cauchy, etc.) which breaks\n", "premises for using chaos expansion as approximations, and stochastic\n", "dependencies, which there currently does not exist numerically stable method\n", "for creating. Let us consider an synthetic exponential model than encompasses\n", "both issues by using a multivariate log-normal distribution for its\n", "uncertainty:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:14.662194Z", "iopub.status.busy": "2021-05-18T10:57:14.661828Z", "iopub.status.idle": "2021-05-18T10:57:14.669124Z", "shell.execute_reply": "2021-05-18T10:57:14.669393Z" } }, "outputs": [], "source": [ "import numpy\n", "import chaospy\n", "\n", "coordinates = numpy.linspace(0, 10, 1000)\n", "\n", "\n", "def exponential_model(parameters):\n", " initial, rate = parameters\n", " return initial*numpy.e**(-rate*coordinates)\n", "\n", "\n", "distribution = chaospy.MvNormal(mu=[10, 1], sigma=[[1.0, 0.09], [0.09, 0.1]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are interested in the mean and standard deviation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Monte Carlo integration\n", "\n", "As a baseline we can solve this using quasi-Monte Carlo integration. It\n", "requires no modification compared to the stochastic independent case. It\n", "consists of generating samples:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:14.671634Z", "iopub.status.busy": "2021-05-18T10:57:14.671321Z", "iopub.status.idle": "2021-05-18T10:57:14.961289Z", "shell.execute_reply": "2021-05-18T10:57:14.960818Z" } }, "outputs": [], "source": [ "samples = distribution.sample(10**5, rule=\"sobol\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "evaluate model for each sample:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:14.964481Z", "iopub.status.busy": "2021-05-18T10:57:14.964040Z", "iopub.status.idle": "2021-05-18T10:57:17.161797Z", "shell.execute_reply": "2021-05-18T10:57:17.161507Z" } }, "outputs": [], "source": [ "evaluations = numpy.array([exponential_model(sample) for sample in samples.T])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and performing analysis on samples:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:17.164282Z", "iopub.status.busy": "2021-05-18T10:57:17.163970Z", "iopub.status.idle": "2021-05-18T10:57:17.453255Z", "shell.execute_reply": "2021-05-18T10:57:17.452918Z" } }, "outputs": [ { "data": { "text/plain": [ "(array([9.99998, 9.89953, 9.8002 , 9.70195, 9.60479]),\n", " array([0.9999 , 0.9815 , 0.96431, 0.94833, 0.93353]))" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mean = numpy.mean(evaluations, axis=0)\n", "std = numpy.std(evaluations, axis=0)\n", "\n", "(mean[:5].round(5), std[:5].round(5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also plot the final result:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:17.455620Z", "iopub.status.busy": "2021-05-18T10:57:17.455300Z", "iopub.status.idle": "2021-05-18T10:57:17.515914Z", "shell.execute_reply": "2021-05-18T10:57:17.515560Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD3CAYAAAAE2w/rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtcElEQVR4nO3dZ5Ak533f8W+HyTO7szlcTugDQAQGgAkkCFISZdNyCQpVpqskucx3VlGibakUSpZliipJpeiykmnRoiRSJC2RFEVRjCIJAcQhHO6AwwF3ffluc5jd2ZmdPD3tFz2zN5tnJ8/s/1N1tbs96enb2d8++/T/eR7Ftm2EEEJ0HrXVDRBCCFEdCXAhhOhQEuBCCNGhJMCFEKJDSYALIUSHkgAXQogOpe92B8MwRoGPAQ+ZpvlI8ZgX+F1gCjgF/JZpmlca2VAhhBDrVdIDfwz4EqCUHfsIcMc0zd8E/gD4RP2bJoQQYie7Brhpmn8HxDcc/gBwpnj7q8BDhmH01L95QgghtrPrEMo2hlkf6rHisdjGO775179pf9+9w3cfGPLwM+89VeXLthdNU7GsQqub0TByfp1Nzq9zuVyasvu9qg/weSBU9nVP8dgmS4ksiVQWl+Z09icjeSbnYgQ91b50+wiH/USjyVY3o2Hk/DqbnF/nGhoK7X4nqq9C+QrwdgDDMB4AXjFNc1PvG8AGFhPZdcduL3Xnf7oQQjTTrgFuGMbjwE8AY4Zh/IphGD7gfwJHDMP4FeC/Ah/a6TkWVtcH+E0JcCGEqNmu4ximaT4FPLXFTT9d6YvMxzcEeCRV6UOFEEJso+ETeTRVIZ7Jk85Za8fimdymYRUhhBB70/AAHwq6gS2GUSKJRr+0EEJ0tYYH+EiPF4D5TQEu4+BCCFGLJgS4B4CF1Qzlu//cWU6RL8huQEIIUa2GB3jY58KtKaRyBRLZu+PgWavAxLJczBRCiGo1PMAVRWEoWOqFrx9Gub4o4+BCCFGtpiwnW7qQOb+aWXdcAlwIIarX1ABfWM2uGwePJLNEk7lmNEEIIbpOUwI84NYIuDVyls1yan1gX5NeuBBCVKUpAa4oCsMhZxx846xMCXAhhKhO07ZUGykOo8zF14+D315Kks1355KQQgjRSE0L8MGgG0WBpWRuXWBbti2TeoQQogpNC3CXpjLg37oa5crCarOaIYQQXaOpu9KPhErDKJvHwQu2zMoUQoi9aHKAly5krp9Wn8pZMitTCCH2qKkB3uPV8eoq6XyBWDq/7jZzXoZRhBBiL5oa4OXlhBurUcx5KScUQoi9aGqAw/bj4PFMjumVdLObI4QQHavxy8mGvOu+Hi4ubBVJZslZ6+u/L8swihBCVKzhAX56NLTua7eu0ud3YdubVye8PBdvdHOEEKJrNDzA7xvv2XRstDgOPhtbPw4eTeWYickwihBCVKLhAT7e6yXkca07NlbcpWd2QzkhwKU5GUYRQohKNGVDh3uGAuuO9Xh1fC6VTL6waXXCS7MyjCKEEJVoShXKPcPBdV8risJocbPjjcMoK+kck1GZ1COEELtpSoAf7vPh0bV1x9aGUTYEOMDr0gsXQohdNSXANVXh5KB/3bHBgBtNVVhJ50mWbXYMzji4rI0ihBA7a9pEnlND64dRNFVhuLhG+OyGWZmJbF6WmBVCiF00LcBPDAbQFGXdsdG1YZTNpYOvzcgwihBC7KRpAe7RVY70rx9GKdWDL6xmyRfWz8q8srAqO/UIIcQOmroWyj3D68sJvS6NPr+Lgr15r8ysVZCp9UIIsYPmBvhQEIX1wyhjxV74zBbVKK9Ox5rSLiGE6ERNDfCgR18b9y4Z6707Dr6x8uTOcorohok+QgghHE1fTnbjpJ6QRyfo1shaNpHE+mEUG1t64UIIsY2mB7ixYVq9oiiM9zqzMqdXNg+jXJiObVovRQghRAsCfDDooa+4O33JWoDH0pvCeiWd49aSTK0XQoiNmh7gsLkXHvY5i1ulc5sXtwJ4eWqlWU0TQoiOodfyYMMwfh44CiwCp4APmaa5a3f5nuEgz91eXvtaURTGe7xcjySZXsnQv6GHfmV+lWTWwu/WNj6VEELsW1X3wA3DGAV+CfiwaZr/HQgAP1LJYw/0evG71v/uuDsOvnkYxbJtLkxLL1wIIcrV0gNPAlmgB4gCQeC1jXfSNJVw2L/xMA8d6eOlO3d74aN9Gh5dJZG1SBWcYZVylxZTvP8hH8qG6fittN25dQs5v84m59f9qg5w0zRjxSGUzxmGMQNMAtc23s+yCkSjmxemOhh0cSabX3dstMfD7aUUdxYT+EfWlxvOZPO8fCPCsYH2+YaFw/4tz61byPl1Njm/zjU0FNr9TtQ2hPIw8PPAB0zT/A844+C/Wunjj/X7canrX368uMnD1MrW+2Kem4xW01QhhOhKtVShHACWTNMsdaNnAG+lD9Y1leMb1ggfCrrRVYVYOk88k9/0mKsLCWJpmZkphBBQW4B/DXjNMIzfMwzjvwGPAL+5lyfYao3w8eLU+qno5l54wbY5PykXM4UQAmobA7eAn67lxU8NBVAVZd0aKAd6fdxZTjMZTXN6wzg4wPnJGO88PoCuts/FTCGEaIWWTOQp8bk0DoZ9644Nh9y4NIV4Jr/lcEkyl5c9M4UQghYHOMA9G2ZlqorCgWJN+OQWwygAZ+9EG90sIYRoe20Q4JuHSQ6G7wb4VgtZzcbT3FmW9VGEEPtbywM87HcxHFy/RvhgwL02qSea2lyNAvBC2VR8IYTYj1oe4LB5jXClbBhlu5rwqwsJlpLZLW8TQoj9oD0CfMM4OJQPo6S2HEaxsXnxdrTRTRNCiLbVFgE+2uOl1+tad6zf78LnUknlCiwlt568c2E6RipnNaOJQgjRdtoiwMGpCS9XPowyEd36gmWuUJCKFCHEvtU2AW4Mb65GOdzn1IhPRtNYha23VTs7ESVnFRraNiGEaEdtE+CH+nz4XOs3bOj1uej16uQsm7n45v0yAVI5S3bsEULsS20T4KqibBpGASfYgR3rvp+7Fd22hy6EEN2qbQIctp7Uc6hYjTIbz5DJbz1UEs/kuDgTa2jbhBCi3bRVgB8f8OPW1jfJ69IYDrqx7e1rwgGevbW8ZbmhEEJ0q7YKcF1TOb7Fjjuli5kTOwyjLCezvD632rC2CSFEu2mrAIfNszIBxno9aKrCUjLH6hYbPZR870ZEeuFCiH2j7QL85KCzRng5XVU5UNzoYaeLmYuJLJekFy6E2CfaLsC9Lo2j/dsPo9xZ3npqfcn3bixJL1wIsS+0XYDD1pN6BgNu/G6NVK7A/Or2i1gtJDIyFi6E2BfaMsDvGQ6gsH4YRVEUjhR74beXdl4L/OnrkXXbtAkhRDdqywAPuPW11QjLlQJ8OpbetiYcYCmZ5eKMbLsmhOhubRngwJYbGvvcGiMhD7a9c0khwDM3IjI7UwjR1do2wI3h4KZhFICj/U4v/NbSzhczo6mcrJEihOhqbRvgPV4X472bh1FGezx4dJV4Js/yNuuElzxzY4nsDkMtQgjRydo2wAGM4c2LW6mKslZSeGuXi5mJbJ4XZb1wIUSXausAPz0S2vJ46WLm5Ep617XAz9xaIpmVXXuEEN2nrQM87HMx1rN5GCXk1RkMuLAK9q4XM7NWgWduRBrVRCGEaJm2DnDYuhoF4Fhx0asbkZ0vZgKcm1yRHeyFEF2n7QP83uGth1HGe71rFzMjiZ0vZhZsm+9eXWxE84QQomXaPsDDfhejoc3DKKqicKxYUngjktz1eS7Pr+64EJYQQnSatg9wgHtHtx5GOdrvRwGmV9KkcrtfqPznKwuy0JUQomt0RoBvM4zic2uM9Xiw2X19FICZWJoL07L1mhCiO3REgIf9W1ejABwfdC5m3owkK1rA6rvXIjuuoyKEEJ2iIwIc4N5tqlEGA26CHo10vsDMSmbX50lk83zv5lK9myeEEE3XQQEe2nJtFEVR1vbRvLaYqOi5Xry9TCQhZYVCiM7WMQHe63NxYIu1UcCZmenSnD0zK6n3tmybb5oL9W6iEEI0lV7Lgw3DMIAPAingceDXTNN8oR4N28p9oyEmVzZfrNQ1laP9fq4uJLi+kKT/iHvX57oRSWDOrWJsMzQjhBDtruoeuGEYGvD7wEdN0/xt4EPAzXo1bCunR4KbNjwuOT7glBROraQrXvvkm1cWdl1LRQgh2lUtQyiPAArwYcMwfgn4IaCh0x2DHn1tJcKN/G6NA71ebCqb2AMQS+d4+oZc0BRCdKZahlCOAG8HPmia5ophGJ8CssAny++kaSrh8OZd5qv11pNDTL88teVtp8d6mFxJc2spyQMHe3Fpu/9+enkmzjuNYUa2KVPcSb3Prd3I+XU2Ob/uV0uAx4DLpmmWtr15BngPGwLcsgpEo5X1iCtxwK9j5SysLWq+Q26Vfr+LpWSOa3NxTgxuXk98K5997jY/+chBlG2GZ7YTDvvrem7tRs6vs8n5da6hoa0nL25UyxDK88BAcSwcnB75lRqeryJel8bJoe2DuXTbtcXKJvYATK2kOC/brwkhOkzVAW6a5hLwC8AfGobxq8AQ8Af1athO3jDas+1t4z0eAm6NZNZiKpqu+Dm/czVCPJ2vR/OEEKIpaiojNE3zi8AX69SWip0YCuDVNdL5zdUmiqJwz3CA85MxriwkOBj2VjQ0kslbfMOc50cfGm9Ek4UQou46ZiJPOV1Vtt3oAeBw2IfXpRJL55mN7z69vsScX+XyXLweTRRCiIbryAAHeMPY9sMoqqpwqngB05xP7GkJ2a9dWqhoaVohhGi1jg3wQ2EvYZ9r29uPDvhwawrLyRyLe1j3JJnL8/VL8/VoohBCNFTHBriiKNy/w8VMXVXXygivzFe2yFXJ63NxzLnVmtonhBCN1rEBDvDA+M61kscH/OiqwvxqluXkzvtmbvTVS/MVT8kXQohW6OgA7/e7Odi79dR6ALeucrS4b6Y5v7cedTKX56uX5mpqnxBCNFJHBzjAA+PbD6MAnBoKoCkwE8vsuRduzq/yqmzBJoRoUx0f4PeNhtDV7U/D69I4Vtzw4VIV49pfvzxPNLW34BdCiGbo+AD36Cqnh3de0/ue4SCaqjAXz7C0x514slaBL1+cld3shRBtp+MDHODBXYZRPLrKiRp64RPRFM/eXK6qbUII0ShdEeBH+n071oSDMxZeqkjZS114ydM3IkxGN+8GJIQQrdIVAa4oyq69cLeucnKw2Auf3XsvvGDb/MPFWdIyS1MI0Sa6IsDBGUbZatf6cieGArg0hcVEloXVytdIKYmmcnxVZmkKIdpE1wR4j9fF8YGdd+dwa+raGimvzaxWdWHy0lycc5PRapoohBB11TUBDvDGg7273ufEoB+PrrKcyjG1Uvl64eW+ZS4yG6vusUIIUS9dFeAnhwKEPDsvca5rKvcWl6J9bXYVq7D3Xni+UOCLr8p4uBCitboqwNUKLmaCU7US8ji79tyscAf7jZaTWb5wfuvNlYUQohm6KsABHj7Qu+vFTFVRuH/MWQjr8vwqWatQ1Wu9PhPjzK2lqh4rhBC16roA7/W5KtqNfjTkYTDgJmfZe15uttxT1yLcqrIXL4QQtei6AAd406HdL2YqisIbir3w64uJqpeOLdg2X3x1RtZLEUI0XVcG+IkB/64zMwH6/C4Ohr0UbLg4U/1emKmcxedfmSZX5VCMEEJUoysDXFEU3lRBSSHA/aMhNAWmVtJVTe4pmYtn+MfXZP1wIUTzdGWAAzx0oHfHZWZL/G6Ne4qrGV6YjlOoYdXBS3NxvndDLmoKIZqjawPc59K4f3TnLddKTg0FCLg1Yul81WWFJf9yPcLlueqHY4QQolJdG+AAbzkcruh+mqrwQPGC5qXZVTL56ifo2Nh8+eIcMzJTUwjRYF0d4CMhD4fC2++ZWW60x8NIyE2uYPPaTG070ucKBf72/DSxtFSmCCEap6sDHODRw30V3U9RFB4Y70FR4PZyas8792y0ms3zufPTZPJSmSKEaIyuD/B7hgMVlRQChDw6J4uTgM5PxWq6oAmwsJrh869MV7XeihBC7KbrA1xRFN5yKFzx/U+PBNcuaF5dqH6GZsmtpSRfeW1O9tQUQtRd1wc4OOujeHStovvqqsLDB5wFsS7PrbKaydf8+hdnY3znWqTm5xFCiHL7IsDdusobD+y+SmHJcMjDoT5nhubLU7G69J6fu7XE87dlY2QhRP3siwAHeORwH6qy8yqF5R4Y68GtKSysZrmzXJ/NjL99ZZEL07G6PJcQQuybAA959bXFqyrh0VUeKK4t/upMvC6bN9jY/NPrc5jztZUpCiEE7KMAB3jbkb5d1wovdyjsZTjoLDlbr6GUgm3z9xdmuBGp/QKpEGJ/21cBPhj0cHJo97XCSxRF4Y0He9FVhZlYholofWZXWrbN51+eqdvQjBBif9pXAQ7wjmOVTewp8bu1tW3aXpmKVb1u+Ea5QoH/d36KyaiEuBCiOjUHuGEYPsMwLhiG8bv1aFCjHej1cbTfv6fHHO7zMtbjIV+wOTe5Urea7qxV4LPnpphakRAXQuxdPXrgHwPO1+F5muYdx/r3dH9FcWrDS1Upta5YWM4J8WkJcSHEnim19CYNw/gJIAE8CARN0/y5jfcpFGzbasOdaj7+9A3uLO0tiCeWkjx7YwlNVfiBe4cJ+901T7cv8egaP/X2Ixze418HjaRpKu34vasXOb/O1s3n53JpFVVb6NW+gGEY9wH3mqb5y4ZhPLjd/SyrQDTafpv+vmU8xLXZvdVkjwTdHAp7mYimefZ6hO+/bwSrhqVny2WzeT7+3Wv8+MPjHGmTEA+H/W35vasXOb/O1s3nNzRUWclzLUMoTwJpwzB+EXgMeNQwjI/U8HxNdXwgwIHeypaaLffQgR4Cbo2VdJ5XJlfq2qasVeBz56e5viglhkKI3VUd4KZp/oZpmh81TfO3gGeAF0zT/MO6tawJ3nVib2PhAC5N5ZHDYRQFrs6vMr1S340b8oUCf/fyNOacTPYRQuysHlUoPwq8G3ibYRgfrL1JzXN8IMDBKnrhfX7X2nZt5yZX6lZaWGLZNl98dYZXpurbwxdCdJeqx8BLTNP8PPD5OrSlJR4/OcCnX5rc8+NODvqJJHPMrKQ5eyfKYyf697TWym4Kts1XXp8jmbN4+9G9/6UghOh++24iz0ZH+v17rgsHp7Tw0aN9eHWVSDLHazON2cj4O1cX+aa5IOuJCyE22fcBDvCek4NVPc7r0njkSBgFuLaYbNisyhfvLPP3r86Sl519hBBlJMCB8V4vxnCwqscOBtw8MF4cD5+INWwj40tzcT7z0iSpOqyKKIToDhLgRY+fHKx6DPv4gJ9DYS+WbfPcrSjZBk0umIim+MsXJlhK1rbhshCiO0iAFw0G3GuLVu2Voig8fLCXXq9OImvx0p36rZey0VIyy1++MMHtPc4iFUJ0HwnwMu8+MYBbq+6/RFcV3nokjEtTmI1nuNigi5oAqZzFZ85N8XKdJxIJITqLBHiZoEfnrUf2ttxsuYBH561lFzVv1XHRq40Kts0/XZrjG5fn67YeixCis0iAb/C2o32EPNWXxw8FPWu72r88FWNhNVOvpm3p7ESUz7w0RSKbb+jrCCHajwT4Bi5NrbqssOTogJ+Tg35s4PnbUeKZxobr7eUk//e5ibpP6xdCtDcJ8C28YSzEeI+35ucY7fGQs2zO3Fyuy6bIO4lncvz1ixOcm4w29HWEEO1DAnwLiqLwA6eH97QB8lbP8cihXsI+pzLlzK1lcg1eu9iybb52aZ5/uDhLNt+d6yQLIe6SAN/GeK+36rLCEl1TefvRPgJujWgqz/O3o1hNmE15cSbGJ1+40/DxdyFEa0mA7+A9pwbw6lpNz+F1abzjWB8eXWVhNctLE42rES+3mMjyF89PcF5KDYXoWhLgOwi4dd5zaqDm5wl6dN5xrA9dVZhaSXNhOt6UEM8XCnz10hyff2VapuAL0YUkwHfxxgO9jNV4QRMg7HPxtqNhVAVuRJJcnGlOiAOY86v8+ZnbDa1LF0I0nwT4LhRF4V/fN1KXtb6Hgh4ePeLs5nNtMcnrs6tNC/F4Js9nzk3xTXOBfJduBCvEfiMBXoGRkIdHD4fr8lxjPV4ePezM1ryykOByE7dOs7F58c4yn3jujtSMC9EFJMAr9O4TA/T5XHV5rvFeL2853AvA5flE0/e/jCSz/NWLE3z76qKsMS5EB5MAr5Cuqfyr+0Zqqg0vdzDs482HnBB/fW6Vy3PNG04BZy2V524t8Ykztxu2EYUQorEkwPfgaL+fhw/WVhte7nCfjzcddEL80txqUy9slkSSWf76xUm+fmmejEz+EaKjSIDv0ftODdHrrc9QCsCRfp8zJl68sPnyVKzpIW5j89JklI8/e4vLc41bBlcIUV8S4Hvk1lU+cH/9hlIADoS9vO1IH6oCt5ZSnJ1YackSsfFMni9cmOFz56aIJhuzNZwQon4kwKtwtN/PW+pUlVIy2uPhncf60VWFyWia525FW1budz2S4ONnbvPPl+el5FCINiYBXqUnTg4wEqp9gk+5waCbx47349YU5uIZnr6x1PBVDLeTLxT4jjnP/372tgyrCNGmJMCrpGsqP/bmg2h1mOBTrs/v4t0nB9YWwHrqWoRYunWbNaykc3zhwgyfOjvJbExqx4VoJxLgNRjr9fLEqdo2f9hKyKPz+MkB+vwukrkC/3I9wuJqa3eiv7Oc5C+en+DLF2eJpWV8XIh2IAFeo0eP9HFyMFj35/XoKo8d72esuCnEMzeXWr6WiY3NqzMx/ux7t/n21cWWDe8IIRwS4HXwQ28YIeSpX2lhSWmn+xODfmwbzk/FeHkq1vJNjPOFAs/dWuJPnrnFmVtLcqFTiBaRAK8Dn0vjyQdH67Lg1UaKovDgeA9vOtiDqsDNSJJnbiyRybe+95vOW3zn6iJ/8swtzt6JyrR8IZpMArxODoZ9vLcB4+ElR/r9vOt4P15dJZLI8Z2rEZbbpFZ7NZvnG+Y8f/rMLV6akCAXolkkwOvo0SN93DsSatjz9wfcPHHKubiZKl7cvL6YaPrMze3EMzm+ftkJ8rN3WlfHLsR+IQFeZx+4b4ShoKdhz+91abzreD/HBnwUbLgwHeeFO9GGb5i8F/FMjm+Y8/zR07d49mbratmF6HYS4HXm1lV+/KFxfK7a9tLciaYqPHygl0cO96KrCtMrGb7dRkMqJclcnu9eW+SPnr7Jt68uEm9hPbsQ3UgCvAHCfhdPPjjWkIua5Q6GfTxxaoCwTyeZtXjqegRzfrXlVSobZa1S1cpNvnxxlrl4ptVNEqIrSIA3yNF+P+8/Pdzw1wl6dN59YoDjA06p4euzqzx9fYnVTPv1di3bqSP/xHO3+fTZSa7MN3cNdCG6jV7tAw3DOAF8DDgHHAQipml+tF4N6wZvPNjLcjLLc7eXG/o6mqrw0IEeRns8nJtcYSmZ49tXIjwwHuJovw+lwX8JVOP2cpLby0l6vS7edKiXhw/0NnTYSYhuVEsPvB/4rGmav2Oa5s8C/84wjDfXqV1d44lTg5werv9Mza2MhDy879QgB8NeLNvm5akYz95cJplt34uIK+kc37m6yP/6l5v8w8VZ2R1IiD2ougdumuaLGw6pQKK25nQfRVH4t28YJXFuiokmhJNbV3nkcJixnhQvT8WYX83yLXOR+0aDHB/0N3xcvlr5QoGLMzEuzsQYCnh4+GAPbxjrkV65EDtQ6jEGaRjGk8B7ij3xdQoF27baqMStnjRNpdJzS2Ut/vx7N5lr4op+qZzF+TtRJpadXxx9fhePHO2jz++u6PGqorT0gqiuqtw7FuJNh8KcHA7WfShoL9+/TiTn17lcLq2iN3vNAW4YxhPAk8BHTNPc9L+Zy1l2NNraRZgaJRz2s5dzW83k+asXJ4immlvuNxNL88pUjFTO+facHPRzeiSIS9t5BM3t1slm2+NiaMijc/9YDw+MhepWZ7/X71+nkfPrXENDocYHuGEYHwDeBfwSMAYcMU3zTPl9JMDXiyZz/PXZCeJNrhLJWwVen1vl+qLTXo+ucv9oiMN93m17tu0U4OVGQh7uHw1x32iInhr2J+3mAAA5v07W8AAvXrB8CjhbPBQA/tg0zU+W308CfLPFRJZPvThJMtf8cFxO5nhlOrY26afP5+LB8RD9gc3DKu0a4CUKCgfDXu4bDWEMBwl69nZJp5sDAOT8OllTeuCVkADf2sJqhk+dnSTVgmnmtm0zEU3z2kycdN4ZVjlUDEK/++5Fw3YP8HIKCof6fNw7Eqw4zLs5AEDOr5NJgDdBrW+g+XiGv3lpqiU9cXCGVcz5BNcWExRsUBU4NuDHGA7i0dWOCvByCgrjvV7uGQ5wz1CQgS3+uoDuDgCQ8+tkEuBNUI830OJqhk+/NEWihUGZyOR5fW6VyahTIaOrCieHAtw33oPdBVf5B/xuTg0HOTkY4GDYu1ZK2c0BAHJ+nUwCvAnq9QZaSmb5m5emWr7XZDSV4/XZ1bW1Sjy6yqmhAMcGfOhqd6y64NU1jg/4OTEY4I0nBsmnWrvXaCN1c8BBd5+fBHgT1PMNFEvn+MxLU0SSrQ+UxUSW12biLBUvdLo1hZODAY4P+nctPewkHo+Lfo/GsX4/xwb8HAj70NX2nOhUjW4OOOju85MAb4J6v4FSOYvPnZtiuomTfbZj2zaRVJ6LZRUrLlXh+KDTe/XonR/kG8f4XarKoT4fx/r9HOn3MRLytOU6MpXq5oCD7j4/CfAmaMQbKGcV+OKFWa4trtb1eavhdutkMjkWE1kuzyVYTDh/HWiKwpF+HycG/Xsu3Wsnu12k9eoah/t8a/86LdC7OeCgu89PArwJGvUGKtg23zIXODsRrftz78XGgIskspjziXXreY+GPJwc8jMYcHdUuMHeyyQ9usbBXi8H+3wcCvsY7/Ggt/GQUjcHHHT3+UmAN0Gj30Bn70T51pWFlq1Hsl3AraRyXF9MMhFNUdq/uNerc2Kws8aRay2T1BSFkZCHA2EfB3q9HOj10uurfmZovXVzwEF3n58EeBM04w10M5LkixdmSOebP+Fnt4DL5C1uRFLcjCTJFCcEuVRnQs2xAV9N09yboRF17gG3znivl/EeL+O9XkZ7PC1bUbGbAw66+/wkwJugWW+gaDLH374yzcJqc7ciqzTgrILNZDTFzaXUun05+/0ujvb7ORj2orVhr7xZE5XCPhejIQ+jPV5GQh5GezwE3I2/dtDNAQfdfX4S4E3QzDdQzirw1dfnuTgba8rrQXUBt5LKcXMpxcRyinxxfEVXFQ70ejnc52Mg4GqbsfJWzjQNeXSGgx5GQh6GQx6Gg276A+66rtfezQEH3X1+EuBN0Io30LnJKN8yF8kXGj9DspaAyxcKTEbT3IqkWC5bPtfvUjnY5+Nw2EfI29oKlnZbKkBTFAYCboaDHgaDboaCbgYDbsK+6n7pdXPAQXefnwR4E7TqDTQXz/ClV2fWyvoapV4BF0vnmYg6vfLSmuQAYZ/OgV5nrLgV5YjtFuDb0VWVgYCLAb+bgcDdf/1+144Tq7o54KC7z08CvAla+QbKWQW+dWWB85MrDXuNegecbdtEEjnuLKeYWkmvDbGAU8UyXgzznib1zDslwLejoBD06AwEXPT5XfT73PT5nc/7fC4GB4JdG3AgAQ4S4DVphzfQtcUE//TaHKsNCKJGBpxVsJmLZ5heSTMTy6wL85BHZ7zXw2jIQ5+/cWPmnR7gO1FQGOjx4lOdi6ilf70+F2GfTtCjt+3+qJVqh5+/RpEAb4J2eQOlchbfvLxQ9wuczQo4q2CzsJphaiXDTCxNzrr7nnRrylrlxnDIg7uOE2e6OcBh5/NTFYUer06P10WPV6e37PMer07Io+Nt8w2l2+XnrxEkwJug3d5AVxdW+dqlBeKZ+qxq2IqAK9g2C6tZZmMZZuMZktm79e8K0B9wSvKGgh7CPr2m3vl+DvCKHq+pBD1OoAc9TqiHPDpBj0bQoxN0O5+3ajZqu/381ZMEeBO04xsoky/w1LVFzk2u1DyDs9UBZ9s28YzFXDzDbCxDJJGl/IxcmsJgwM1w0M1Q0EPQo+0p0Ft9fo3WrPPz6hoBt0bAozsfi58H3Rr+4r+AS8fv1nDXcRG0dvz5qxcJ8CZo5zfQbCzN1y8vMLWSqvo52i3gclaB+XiWuXiGhdUsyQ3b0Xl1laHg3SqN0C6B3m7nV2/teH6aouB36/jdKj6Xht+l4XNp+NzO516Xc9xX9rlHV7ccr2/nn79aSYA3Qbu/gWzb5tWZON+9uljVRc52DIByiWyehXiWhdUsC4ns2nT+EpemFEvvXAwU66nLZ4S2+/nVqlvOT0HBpSlrYe51aXh1lf5eH1Y2j0dXnePFjx5dw60reDQVt67i1pzj7TKBrBIS4E3Q7gFeks0XOHNriRduR8ntYQJQJwWAbdvE0nkWElmWEjkiiezahs0lquJUZJTK7IZ7vbgVOuoHey866ftXjb2cX+mXgGst1BXcmhPuruLXLnXj58WPmoKuqbhU5/EuTUEvfq6rzueaqtT1fSQB3gSdEuAl8XSep29EuDAdq2h8vJMDwLZtkjmLSCLHUiJLJJkjlt58Li5VIex3Suz6fC7Cfh2/a29j6e2qk79/lWin81NwQrw80HVVQdcUNMX52jmmoqusfa2pZbcrCmrx4w+/9YgEeKN1WoCXRBJZnr4e4dLcKjbbf//b6QekHrJWgeVkjuVkjmjK+Vc+M7TEpSlrZXW9Xp0en1OJ0Wn7gnbb92+jbj6/P/j3b64owDt3OxVRtYGAmx9+cIx3rmZ45sYSl3cJ8m7h1lRGQs4CUuAEwEoiQzSZY7kY6MvJHFnLZjGRYzGxvhwz4Nbo9en0el1OSZ1XI+DW23KlRbE/SIDvY0NBD08+OEYkkeXMrSUuzsRbtnlEq/hcGr5ejbFeL+AMvaTzBWLpPCspZ9hlJZ0nns6TyFokshbTK+uX9Q24tbX66JD3br10PUvmhNiKBLhgIODm39w/yuMnBnnxzjLnp2JkWrCBRDtQFGWtjK3UUwcoFGziGSfMY8VAj2fuhnoia0F8/XO5NWV9bXRZrbS3w6oiRHuSABdrQl6d994zxGPHB7gwHePiQoLppe4cY9wrVVXoLa4lUs4q2CSyeeJpi3gmz2rGCfZ42iJr2WSLY+6bnk9xdu9xQr1YD+0u1UV3XtmbaA0JcLGJW1d5y+Ew73tgjPM3Fjk/ucLVhcS+G16phKYqxTVE1gd7aSgmkbVIbOipJzJ5spbTo49n8pt67uAEfOkvgdKkF+dzJ+C9uoZLq2/pmug8EuBiW4qicHwgwPGBAPF0ngszMS5MrazboEFsrXwoZjDg3nR7ziqsC/RUrkAya5HKWSRzFjnLvjs0k9j6NVTFmcbudTmTWEoTXDzFgA/5C2i2Lb35LiYBLioS8uq881g/7zzWz53lFK9Ox7g8v7pvx8pr5dJUwj6V8Da72OcLBVLZAsmcE+qprEUyV3ACPmuRyRfIF5xa941LCmykAB7dmcDi0dW1GYprxzTVmblYPObStp66LtqPBLjYs8N9Pg73+Xj/6SGuLSZ4fTbOtcVkU7Z52y90VSXkVXfcdi5vFUjni/9yVvFjgXTeIpMrkLEKpLLOWHzpfpVyaU6gr81WLM5cdGlKcbaiendm49pHVUoqm0wCXFRN11ROj4Q4PRIimy9wdTGBORfnRiRJ1pIwbzRdUwlqKkHP1reXJroUCk6AZ/NOqGfzBTIbPi+/LWvZ5CybnGUBe/sLS1VYF+qu4hR0feP08+J09bvHS7MUndvkL4DKSICLunDrKvePhrh/NETeKnBzKcnVhQRXFxIkunS2XKdQVWVtWddKFGybnFUKdpusVSBnFchZpc+d27NrH+8eK9jOksaZPOw1/MtpivMLqhTuWtkU9dLnbl1Dse0tbytNW9eK09adzxW0Llv7RgJc1J2uqZwaCnJqKIht28zGM1xfTHBjMcl0LC3VLG1OVRQ8uoZH39uOPLZtY9kUw/5uwOctm3zB6dXnC4X1n1s2uYJdvM/dry0brHyBzO4vu2fO2iNOsKvr1iFh3Xokpa/X7qcUH1O6vfQ8ZY9x7uv8H6pq8aNSdqzOv0AkwEVDKYrCWI+XsR4vjx0fIJ2zuL2U4tZSkltLSSLJbKubKOpEURR0BXTVqb6plm3bWIXy0HfCvXSs9BFFJZvLO8fsu7dtvF/5sYINlm3jjA61piOhwFqgK+uC/u7HSkmAi6byujSMkSDGSBCA1Uye28spJor/FhPZfbEui9ieohTHzDXwbl2kA1S3mJVtl4Lcxio4YV4o/gIoFL+21m4vBv7a7cW/DHa83Xnegm1j287thQ0fbZzfHZZtA3YtI021BbhhGN8H/AgwD9imaf6PWp5P7D9Bj742dg6QzllMrqSZiqaYWkkzvZKWC6KibhTFGVNvFdt2uieF4i+MTeFu28Vgr0zVAW4Yhh/4M+B+0zQzhmF83jCM95mm+c/VPqcQXpfGycEAJwcDgPOGjiSyzMScHetnYxnmVzMS6qIjKYpSHEJRoA5rndXSA387cNs0zdJ1hu8BHwDWBbjLpSlDQ6EaXqa9dfO5QXuc3zBwb6sbIUQbquV3wDDrV3GIFY8JIYRogloCfB4o7571FI8JIYRogloC/AxwxDCM0jywdwJfqb1JQgghKlHTnpiGYXw/8GPAApCTKhQhhGiehm1q3M0lhoZhjAIfAx4yTfORVren3gzDOIFzfueAg0DENM2PtrZV9WMYhgp8GXgecAMngP9ommaqpQ2rI8MwfDjn9w3TNH+u1e2pJ8MwngPSxS8t0zTf18r21JthGAbwQSAFPA78mmmaL2x134ZM5NkHJYaPAV8CHm5xOxqlH/isaZpfAjAM43XDML5imuZLLW5XPZ0xTfNjAIZhfAmns/Hp1japrj4GnG91Ixrka6Zp/lqrG9EIhmFowO8DP2SaZsEwjL8Ctp2t1KiZmBWVGHYq0zT/zjCM97S6HY1imuaLGw6pbLutQOcxTbOAE3AYhqHj/JVhtrRRdWQYxk/g/Mw9CARb3JxGeMAwjF8AfMCLpml207W3R3Bm23+42BGOAP9nuzs3attsKTHsEoZhPAl83TTNy61uS70ZhvF+4B+BfzRN82yr21MPhmHcB9xrmuYXWt2WBvpt0zR/G/h14JcNw3h3qxtUR0dwOsCfNE3zN4F3Az+13Z0bFeBSYtgFDMN4AngC+M+tbksjmKb5ddM0fxA4ZhjGf2p1e+rkSSBtGMYv4gz1PWoYxkda26T6Ko0Hm6ZpAU/jvEe7RQy4bJrmSvHrZ4D3bHfnRgW4lBh2OMMwPgC8H/hZYNQwjLe3uEl1YxjGfcXzK7kJHG9Ve+rJNM3fME3zo6Zp/hbOD/8Lpmn+YYubVTeGYZw2DONDZYdOAddb1Z4GeB4YKI6Fg9Mjv7LdnRtZhdK1JYaGYTwO/CTwg8CfAr/XZRUMbwaeAkrDCgHgj03T/GTLGlVHxSqb38GpsnHhzNT/GdM0Z1vasDoyDONHgZ/GqbL5Y9M0P9PiJtWFYRjjwB/hXKDtwfn+/ZfidY2uUBy2fC9Odh4GPrxdvjQswIUQQjRWo4ZQhBBCNJgEuBBCdCgJcCGE6FAS4EII0aEkwIUQokNJgAshRIeSABdCiA71/wET629m5poH6gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from matplotlib import pyplot\n", "\n", "pyplot.fill_between(coordinates, mean-std, mean+std, alpha=0.6)\n", "pyplot.plot(coordinates, mean)\n", "\n", "pyplot.axis([0, 6, 0, 10])\n", "\n", "pyplot.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generalized polynomial chaos\n", "\n", "Polynomial chaos expansions builds on the assumption of having an orthogonal\n", "polynomial expansion. However, the classical extension to the multivariate\n", "case assumes that the probability distribution consist of stochastic\n", "independent components. If the distribution has dependencies, the classical\n", "approach will not work.\n", "\n", "The recommended approach for addressing dependent distribution is to use\n", "*generalized polynomial chaos expansion* (g-pce). It assumes that there\n", "exists a smooth map $T$ between the dependent variables $Q$ and some other\n", "stochastic independent variables $R$, which we can build an expansion for. In\n", "other words:\n", "\n", "$$\n", "\\hat u(x, q) = \\hat u(x, T(r)) = \\sum_{n=0}^N c_n \\Phi_n(r)\n", "$$\n", "\n", "For multivariate normal distributions, the obvious choice is to select $R$ to\n", "be standard normal:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:17.518137Z", "iopub.status.busy": "2021-05-18T10:57:17.517822Z", "iopub.status.idle": "2021-05-18T10:57:17.524783Z", "shell.execute_reply": "2021-05-18T10:57:17.525043Z" } }, "outputs": [], "source": [ "distribution_q = distribution\n", "distribution_r = chaospy.J(chaospy.Normal(0, 1), chaospy.Normal(0, 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The $T$ is defined as a double Rosenblatt transformation:\n", "\n", "$$\n", "T(r) = F_Q^{-1}\\left( F_R(r) \\right)\n", "$$\n", "\n", "which in `chaospy` can be constructed as follows:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:17.527222Z", "iopub.status.busy": "2021-05-18T10:57:17.526910Z", "iopub.status.idle": "2021-05-18T10:57:17.533794Z", "shell.execute_reply": "2021-05-18T10:57:17.533529Z" } }, "outputs": [], "source": [ "def transform(samples):\n", " return distribution_q.inv(distribution_r.fwd(samples))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This formulation is general and can be used with any two distributions of the\n", "same size." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Point collocation method\n", "\n", "Implementing g-pce for point collocation require us to generate samples from\n", "$R$ and transform them using $T$:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:17.536135Z", "iopub.status.busy": "2021-05-18T10:57:17.535823Z", "iopub.status.idle": "2021-05-18T10:57:17.547092Z", "shell.execute_reply": "2021-05-18T10:57:17.546787Z" } }, "outputs": [], "source": [ "samples_r = distribution_r.sample(1000, rule=\"sobol\")\n", "samples_q = transform(samples_r)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The resulting samples can then be used to solve the equation above using\n", "regression-based method:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:17.549367Z", "iopub.status.busy": "2021-05-18T10:57:17.549048Z", "iopub.status.idle": "2021-05-18T10:57:18.215632Z", "shell.execute_reply": "2021-05-18T10:57:18.215909Z" } }, "outputs": [], "source": [ "expansion = chaospy.generate_expansion(7, distribution_r)\n", "evaluations = numpy.array([exponential_model(sample) for sample in samples_q.T])\n", "model_approx = chaospy.fit_regression(expansion, samples_r, evaluations)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that for generating the expansion and the model approximation, we use\n", "the distribution from $R$, while for the model evaluation we use the\n", "transformed samples from $Q$.\n", "\n", "The solution model can then be used to do analysis. Just remember that the\n", "model is defined with respect to $R$, not $Q$:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:18.218216Z", "iopub.status.busy": "2021-05-18T10:57:18.217900Z", "iopub.status.idle": "2021-05-18T10:57:18.754789Z", "shell.execute_reply": "2021-05-18T10:57:18.754501Z" } }, "outputs": [ { "data": { "text/plain": [ "(array([10. , 9.89956, 9.80022, 9.70198, 9.60482]),\n", " array([1. , 0.98159, 0.9644 , 0.94841, 0.93361]))" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mean = chaospy.E(model_approx, distribution_r)\n", "std = chaospy.Std(model_approx, distribution_r)\n", "\n", "(mean[:5].round(5), std[:5].round(5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting the final results:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:18.757297Z", "iopub.status.busy": "2021-05-18T10:57:18.756984Z", "iopub.status.idle": "2021-05-18T10:57:18.811062Z", "shell.execute_reply": "2021-05-18T10:57:18.811325Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD3CAYAAAAE2w/rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtRUlEQVR4nO3dZ3Ak6X3f8W+nicBgEBfA5nR9e7xEXmAQeccTJVE2i+U6SawyXSXaZb6zihJtS6VQsixTVEkqRVeZVCClomXxxExRDCIpUtTxlreXNlza3d7biEUGBhgMJoduv5gZ7AziYPIM/p+qLQA96enF4IcHT/+f51Ecx0EIIUTnUVvdACGEENWRABdCiA4lAS6EEB1KAlwIITqUBLgQQnQoCXAhhOhQ+k53ME1zFPg48IBlWY8UjnmAPwKmgJPA71uWdaWRDRVCCFGukh74O4GvAUrJsY8CE5Zl/R7wp8Bf179pQgghtrNjgFuW9SVgdd3h9wFnCre/Cjxgmmag/s0TQgixlR2HULYwQnmoRwrHIuvv+NDv/LPzE6dG7jyw180v/vjJKl+2vWiaSi5nt7oZDSPn19nk/DqXYWjKzveqPsDngd6SrwOFYxssxdLEEmkMLd/ZnwxlmZyL0OOu9qXbRzDoIxyOt7oZDSPn19nk/DrX8HDvznei+iqUbwJvBzBN8z7gZcuyNvS+ARxgMZYuO3ZrqTv/04UQopl2DHDTNB8Hfh4YM03zN03T9AL/GzhsmuZvAv8d+PB2z7EQLQ/wGxLgQghRsx3HMSzLehp4epObfqHSF5lfXRfgoUSlDxVCCLGFhk/k0VSF1VSWZCa3dmw1ldkwrCKEEGJ3Gh7gwz0uYJNhlFCs0S8thBBdreEBvi/gAWB+XYBfX5RxcCGEqEUTAtwNwEI0RenuP7fDCbK27AYkhBDVaniAB70GLk0hkbGJpe+Mg6dzNreX5WKmEEJUq+EBrigKwz3FXnj5MMq1RRkHF0KIajVlOdnihcz5aKrsuAS4EEJUr6kBvhBNl42Dh+JpwvFMM5oghBBdpykB7ndp+F0amZzDcqI8sK9KL1wIIarSlABXFIWR3vw4+PpZmRLgQghRnaZtqbavMIwyt1o+Dn5rKU46251LQgohRCM1LcCHelwoCizFM2WBnXMcrsusTCGE2LWmBbihqQz6Nq9GubIgAS6EELvV1F3p9/UWh1E21oPbjszKFEKI3WhygBcvZJZPq09kcjIrUwghdqmpAR7w6Hh0lWTWJpLMlt1mzUeb2RQhhOh4TQ3w0nLC9dUo1ryMgwshxG40NcBh63Hw1VSG6ZVks5sjhBAdq/HLyfZ6yr4eKSxsFYqnyeTK678vyzCKEEJUrOEBfvdob9nXLl2l32fgOBtXJ7w8t9ro5gghRNdoeICfGuvdcGy0MA4+GykfBw8nMsxEZBhFCCEq0fAA3x/00us2yo6NFXbpmV1XTghwaU6GUYQQohJN2dDhrmF/2bGAR8drqKSy9obVCS/NyjCKEEJUoilVKHeN9JR9rSgKo4XNjtcPo6wkM0yGZVKPEELspCkBfqjfi1vXyo6tDaOsC3CAi9ILF0KIHTUlwDVV4cSQr+zYkN+FpiqsJLPESzY7hvw4uKyNIoQQ22vaRJ6Tw+XDKJqqMFJYI3x23azMWDrLjVC8WU0TQoiO1LQAPz7kR1OUsmOja8MoG0sHX5+RYRQhhNhO0wLcrascHigfRinWgy9E02Tt8lmZVxaislOPEEJso6lrodw1Ul5O6DE0+n0GtrNxr8x0zpap9UIIsY3mBvhwDwrlwyhjhV74zCbVKK9OR5rSLiGE6ERNDfAet7427l001ndnHHx95cnEcoLwuok+Qggh8pq+nOz6ST29bp0el0Y65xCKlQ+jODjSCxdCiC00PcDNddPqFUVhvC8/K3N6ZeMwyivTkQ3rpQghhGhBgA/1uOkv7E5ftBbgkeSGsF5JZri5JFPrhRBivaYHOGzshQe9+cWtkpmNi1sBXJhaaVbThBCiY+i1PNg0zV8BjgCLwEngw5Zl7dhdvmukh+duLa99rSgK4wEP10JxpldSDKzroV+ZjxJP5/C5tPVPJYQQe1bVPXDTNEeBXwc+YlnW/wT8wM9U8tj9fR58Rvnvjjvj4BuHUXKOwyvT0gsXQohStfTA40AaCABhoAd4ff2dNE0lGPStP8wDh/s5O3GnFz7ar+HWVWLpHAk7P6xS6tJigvc+4EVZNx2/lbY6t24h59fZ5Py6X9UBbllWpDCE8nnTNGeASeDq+vvlcjbh8MaFqQ70GJxJZ8uOjQbc3FpKMLEYw7evvNxwJp3lwvUQRwfb5xsWDPo2PbduIefX2eT8Otfw8MatKDdTyxDKg8CvAO+zLOs/kR8H/61KH390wIehlr/8eGGTh6mVzffFPDcZrqapQgjRlWqpQtkPLFmWVexGzwCeSh+sayrH1q0RPtzjQlcVIsksq6nshse8sRAjkpSZmUIIAbUF+LeB103T/GPTNP8H8Ajwe7t5gs3WCB8vTK2fCm/shduOw/lJuZgphBBQ2xh4DviFWl785LAfVVHK1kDZ3+dlYjnJZDjJ3evGwQHOT0b4sWOD6Gr7XMwUQohWaMlEniKvoXEg6C07NtLrwtAUVlPZTYdL4pms7JkphBC0OMAB7lo3K1NVFPYXasInNxlGAXhpItzoZgkhRNtrgwDfOExyIHgnwDdbyGp2NcnEsqyPIoTY21oe4EGfwUhP+RrhQ37X2qSecGJjNQrACyVT8YUQYi9qeYDDxjXClZJhlK1qwt9YiLEUT296mxBC7AXtEeDrxsGhdBglsekwioPDi7fCjW6aEEK0rbYI8NGAhz6PUXZswGfgNVQSGZul+OaTd16ZjpDI5JrRRCGEaDttEeCQrwkvVTqMcju8+QXLjG1LRYoQYs9qmwA3RzZWoxzqz9eIT4aT5OzNt1V76XaYTM5uaNuEEKIdtU2AH+z34jXKN2zo8xr0eXQyOYe51Y37ZQIkMjnZsUcIsSe1TYCrirJhGAXywQ5sW/f93M3wlj10IYToVm0T4LD5pJ6DhWqU2dUUqezmQyWrqQyvzUQa2jYhhGg3bRXgxwZ9uLTyJnkMjZEeF46zdU04wLM3l8sWxRJCiG7XVgGuayrHNtlx51AFwyjL8TSXZJErIcQe0lYBDhtnZQKM9bnRVIXleIboJhs9FP3oxtKmk36EEKIbtV2AnxjKrxFeSldV9hc2etiuF74YS3NpLtrQ9gkhRLtouwD3GBpHBrYfRtmul336ekh64UKIPaHtAhw2n9Qz5Hfhc2kkMjbz0a0XsVqMpbkovXAhxB7QlgF+14gfhfJhFEVROFzohd9a2n4t8GeuhaQiRQjR9doywP0ufW0CT6ligE9HklvWhAMsxdO8NiMVKUKI7taWAQ5gjmyclel1aezrdeM421/MhHwvPCuzM4UQXayNA7xnwzAKwJGBO8Mo212sXElmeFnWSBFCdLG2DfCAx2C8sJxsqdGAG7eusprKbrlOeNHp60uktxlqEUKITta2AQ6bD6OoirJWUrjTxcxYOsuLsl64EKJLtXWA372vd9PjxYuZkyvJHdcCP3NziXhadu0RQnSftg7woNdgLLBxGKXXozPkd5GzHW7vcDEznbM5fT3UqCYKIUTLtHWAA9y9b+OkHmBt0avroe0vZgKcm1yRHeyFEF2n7QP81MjmwyhjfXcuZi7Gtg9n23H4wRuLjWieEEK0TNsHeNBnMNq7cRhFVRSOFkoKb4S2H0YBsOajO9aOCyFEJ2n7AAc4Nbr5MMqRQR8KML2SJJHZ+ULl968syEJXQoiu0RkBvsUwitfQGOtz47BzSSHATCTJK9Oy9ZoQojt0RIAHfZtXo8Cdi5k3QvGKFrD616uhbddREUKITtERAQ5wzxY14UN+Fz1ujWTWZmYltePzxNJZfiRlhUKILtAxAX5qdPO1URRFWeuFX12MVfRcL06ECe1QuSKEEO2uYwI84DHYv8naKJCfmWloCkvxTEX13jnH4Z+thXo3UQghmkqv5cGmaZrAB4EE8Djw25ZlvVCPhm3mntFeJlc2XqzUNZUjAz7eWIhxbSHOwGHXjs91PRTDmotibjFRSAgh2l3VPXDTNDXgT4CPWZb1B8CHgRv1athm7t7Xs2HD46JjhZLCqZVkxWuffNdakNUKhRAdq5YhlEcABfiIaZq/DrwfaOh0xx63vrYS4Xo+l8b+Pg8OcD0Ur+j5VlMZWSdFCNGxahlCOQy8HfigZVkrpmn+HZAGPlN6J01TCQY37jJfrbeeGGb6wtSmt909FmByJcnNpTj3HejD0Hb+/XRhNso77t7H6BZlitup97m1Gzm/zibn1/1qCfAIcNmyrOK2N6eBd7MuwHM5m3C4sh5xJfb7dHKZHLlNar57XSoDPoOleIarc6scH9q4nvhmPv/cLT70yAGULYZnthIM+up6bu1Gzq+zyfl1ruHhzcum16tlCOV5YLAwFg75HvmVGp6vIh5D48Tw1sFcvO3qYmUTewCmVhKcm5Tt14QQnaXqALcsawn4VeDPTNP8LWAY+NN6NWw7944GtrxtPOCmx6URT+eYCicrfs4fvLHIajJbj+YJIURT1FRGaFnWV4Gv1qktFTs+7MejaySzG6tNFEXh5Iif85MRrizEOBD0VDQ0ks7ZfPvyPB94cLwRTRZCiLrrmIk8pXRV2XKjB4BDQS9eQyWSzDIb2Xl6fdEbC1Euzq7Wo4lCCNFwHRngAPeObT2MoqoKJwoXMK2F2K6WkP3u5QXZQ1MI0RE6NsAPBj0EvcaWtx8Z9OLSFJbjmR137CkVz2T57uX5ejRRCCEaqmMDXFEU3rTNxUxdVdfKCK/MV7bIVdHFuVUuz8lQihCivXVsgAPcN759reSxQR+6qjAfTbMcz+zqub99aYFYWqpShBDtq6MDfMDn4kDf5lPrAVy6ytHCUrPWfHRXzx3PZPmnizKUIoRoXx0d4AD3jW89jAJwYsiHpsBMJLXrXviVhahswSaEaFsdH+D3jPaiq1ufhsfQ1nrhl+Z21wsH+O7lecKJ3QW/EEI0Q8cHuFtXuXtk+zW97xrpQVMV5lZTLO1yJ550zubrr81WPC1fCCGapeMDHOD+HYZR3LrK8aHqe+G3wwnO3Fiuqm1CCNEoXRHghwe829aEA5wc8q9VpOymLrzomeshJsMbdwMSQohW6YoAVxRlx164S1c5UeyFz+6+F247Dl97dZZkRmZpCiHaQ1cEOOSHUTbbtb7U8WE/hqawGEuzEK18jZSilWSGb12cq7aJQghRV10T4AGPwbHB7XfncGkqJwuzM1+fie5qjZSiy/NRzt4OV9NEIYSoq64JcIA3H+jb8T7Hh324dZXlRIaplcrXCy/1PWuB2Uh1jxVCiHrpqgA/Meyn1739Eue6qnKqsBTt67NRcvbue+E5x+Err8zIeLgQoqW6KsDVCi5mQr5qpdetE0/nuFHhDvbrhRMZvnxuqqphGCGEqIeuCnCAB/f37XgxU1UU3jSW74Vfno+SztpVvdal2QjP3pT6cCFEa3RdgPd5jYp2ox/tdTPkd5HJOVxZ2N1ys6WeuRaquhcvhBC16LoAB3jLwZ0vZiqKwr1j+eVory3Gqt6Fx3Yc/uHVGcK7XChLCCFq1ZUBfnzQt+PMTIB+n8GBoAfbgddmqt/AIZHJ8aWXp8nkqhuKEUKIanRlgCuKwlsqKCkEeNNoL5oCUyvJqib3FM1HU3z9tVm5qCmEaJquDHCAB/b3bbvMbJHPpXFXYTXDl6dWa1p18PJ8lNPXl6p+vBBC7EbXBrjX0HjT6PZbrhWdHPbjd2msprJcr/GC5OnrS1yS/TSFEE3QtQEO8PChYEX301SF+woXNC/PRkllq5+g4+Dw9dfmmK5ylqcQQlSqqwN8X6+bg8Gt98wsNRpws6/XRcZ2eH1m96sVlsraNl+8MC07+QghGqqrAxzg0UP9Fd1PURTuGw+gKHBrObHrnXvWi6WzfPH8tEy3F0I0TNcH+F0j/opKCgF63fraaoXnpyI1b6O2EEvx5ZdnyFax3ooQQuyk6wNcURQePhis+P7mvh78Lo1IMssbNczQLLq1HOcbUl4ohGiArg9wyJcUunWtovvqqsKD+/MLYl2eixJNZWt+/Ytzq3z/ymLNzyOEEKX2RIC7dZU37995lcKikV43B/vzMzTPT0bq0nt+YWKZMzelRlwIUT97IsABHjnUj6psv0phqfvGArgK269NLNdnM+MfvLHIy1MrdXkuIYTYMwHe69HXFq+qhFtX19YWf3VmtW7VJP90aZ7LMtFHCFEHeybAAd52uH/HtcJLHQh6GOnJLzl7fqo+QynF3e2vLtZ+gVQIsbftqQAf6nFzYnjntcKLFEXhzQf60FWF2UiKieX6zK7MOQ5feXmGm0uyjrgQonp7KsAB3nG0sok9RT6XtjaU8sp0pOp1w9fL2jZfujBdt/F1IcTeU3OAm6bpNU3zFdM0/6geDWq0/X1ejgz4dvWYQ/0exgJusrbDucmVutV0p3M2Xzg/xWRYQlwIsXv16IF/HDhfh+dpmnccGdjV/RUlXxvu0hQWoum6bqGWztl87pyEuBBi95RaepOmaf48EAPuB3osy/rl9fexbcfJteFONX/1zHUmdjkGfXs5zrPXltBUhZ86NULQ56p5un2RS1f50NsOc2Sw8jH6RtM0lXb83tWLnF9n6+bzMwytomoLvdoXME3zHuCUZVm/YZrm/VvdL5ezCYfb72LdQ2O9XJ2N7Oox+/wuDgY93A4nefZaiJ+8Zx+5GpaeLZVOw6f+9RofeHCcI4O7G+JplGDQ15bfu3qR8+ts3Xx+w8OVlTzXMoTyJJA0TfPXgHcCj5qm+dEanq+pjg/52d9X2VKzpR7YH8Dv0lhJZrkwGa5rmzK2zRcuTPPGQm3L2Qoh9oaqe+CWZf1u8XPTND3kh1D+rB6NapZ3Hhvg8+endvUYQ1N55FCQp6+FuDofY9BrMN7nqVubsrbNl1+e4f33jla8o5AQYm+qRxXKzwKPAW8zTfODtTepeY4P+TlQRS+832dwbyFcz02u1K20sMh2HP7x1VnO1bmHL4ToLlX3wIssy/oy8OU6tKUlHjsxyFNnJ3f9uONDPhbjGWZWkrw4EeZdxwd2tdbKThwcvn1pnng6xzuPDdbteYUQ3WPPTeRZ78iAb9d14ZAvLXz0SD8eXWUpnuG1mcasb/LDayG+fWle1hMXQmyw5wMc4N0nhqp6nMfQePRwEAW4thhvWC33uclwfmefLi2ZEkJURwIcGO/zYI70VPXYQb+L+8cL4+G3I6w0aCPjKwtRPnt2ili69g0mhBDdQQK84PETQ1WPYR8d9HEw6CHnODx/K0y6QT3lqZUE//eF2yzWuOGyEKI7SIAXDPlda4tW7ZaiKDx4oI8+j04sneOlifqtl7JeOJHhb1+4Xdfp/EKIziQBXuKx44O4tOr+S3RV4a1HghiawtxqqmEXNQGS2RyfPz/FSxPhhr2GEKL9SYCX6HHrvPXw7pabLeV36by1cFHz6mKcmw3sJduOw3eteb51cY6sLRUqQuxFEuDrvO1IP73u6svjh3vca7vaX5iKsBBN1atpm7owtcJTL00STcnFTSH2GgnwdQxNrbqssOjIoI8TQz4c4PlbYVYbHK6TKwn+5rkJbsvmEELsKRLgm7h3rJfxQG3rm9w71stYwE0m53DmxnLdNkXeSjSd5bNnJ3nh1nJDX0cI0T4kwDehKAo/dffIrjZA3uw5Hj7YR9Cbr0w5c3OZTIMn4tiOw/euLPCVl6dJZWXSjxDdTgJ8C+N9nqrLCot0TeXtR/rxuzTCiSzP3wqTa8IFx8vzUf7muVvMROqzCbMQoj1JgG/j3ScH8ehaTc/hMTTecbQft66yEE1z9nbjasRLLRfqxZ+/tSzrqAjRpSTAt+F36Tx+ovaVAHvcOu842o+uKkytJHllerUpoZpzHL5/ZYEvnJ+WKhUhupAE+A7ecqCPsRovaAIEvQZvOxJEVeB6KM5rM80JcYBroRifPjMhO/0I0WUkwHegKAr/5tRIXdb6Hu5x51cvVPITfS7ORpsW4vFMli9emOZbF+dIywVOIbqCBHgFRgMeHj0UrMtzjRWeSwGuLMS4PNfcXvGFqRU+/dwtbi3JWipCdDoJ8Aq96/ggQa9Rl+ca7/Pw8KE+AC7Px7CaHOLhRIanzk7xnUvz0hsXooNJgFfI0FT+7al9dXu+A0EvDx/Mh/jFuSiX5po3nAL5LdvOTob51JlbsrKhEB1KAnwXjgz6ePOBvro938F+Lw8VQvzyXLSpFzaLVpIZ/v7cJP/42mzdN2cWQjSWBPgu/fjJYfo89RlKATjU782PiRcubF6YirSkbvu1mQh/+exNXp5aafprCyGqIwG+S25d5X337Ktpmv16+4Me3na4H1WBm0sJXppYwW5BiCcyOb55cY7/9+Lthq+iKISonQR4FY4M+ni4TlUpRaMBNz92dABdVZhcSfLczXDLNjG+HU7w189N8K1XZxq+CJcQonoS4FV64sQgI73uuj7nUI+Ldx4bwFXY1eeZ60stC1DbcXj2eoi/+NEtXp5qzvR/IcTuSIBXSddUPvDQQbQ6TPAp1e8zePzE4NoCWP96NUQk2bpp8PFMlm9enONvnp+Q2nEh2owEeA3G+jy8+2Rtmz9spsedX4Ol32eQyNj88Gqo5WPSc6spPnt2ki9emGYxlm5pW4QQeRLgNXr0UJATQ/66P69bV3nXsYH8phC2w49uLLdFvfYbC1E+feYW37o4JwtkCdFiEuA1UhSF9987Sq+7fqWFRZqq8NbDQY4P+XCc/B6bF6Yi2C3exNh2HC5MrfDnp2/yL28skpALnUK0hAR4HXgNjSfvH63LglfrKYrC/eMB3nIggKrAjVCc0zeWSGVbH5oZ2+a5m0t88vRNnrkWkl2AhGgyCfA6ORD08kQDxsOLDg/4eNfxATy6SiiW4QdvhFiOZxr2eruRyuZ45nqITzxzg9PXQ1J6KESTSIDX0VsP93NqX2/Dnn/A5+KJk4MMFC9uXgtxbTHWNiV+yWyOH14L8cnTN/nhtZAMrQjRYBLgdfa+e/Yx7K9vfXgpj6HxzmMDHB30YTvwyvQqL9wKk27RpJ/NJLM5Thd65N+zFlhtYRmkEN1MArzOXLrKzz04hteobS/N7WiqwoP7Azx6KIiuKkxHUm01pFKUztm8MLHMJ0/f4Buvz7a8FFKIbiMB3gD9PhdP3jfWkIuapfYHPTxxcpCgVyeezvH0tRDWfLQl66hsJ+c4vDId4dNnJvjcuSmuh2KtbpIQXUECvEGODPr4qbuHG/46PW6dx44PcnwwX2p4cTbKD68ttWWNtoPD9VCMz52b4q+evcm5ybBsKCFEDZRqL4CZpnkc+DhwDjgAhCzL+tj6+2UyOSccbv0ElEYIBn3sdG7fsxZ4YWK5Ke2ZW01xbnKFZMZGUxTuHevl6KAXpcq/BFwunXS6sb8IPLrGfeMBHjrYx4DP1dDXWq+S718nk/PrXMPDvRX90NbSAx8APmdZ1h9alvVLwL83TfOhGp6vK73nriHMkZ6mvNa+XjfvOTnEwaCHnOPw8nSEZ28sE2twCNcimc3x4sQyf/mjWzx1dpJLc6vkWjxRSYhOoVf7QMuyXlx3SAVkcHMdRVH4d/eO8tTZKSZXEg1/PZeu8vChIGOBJOenVpiPpvm+FeLUaA/Hh3wNH5evloPDzaU4N5fi+Ayd+/cHeGA8wKC/ub1yITpJ1UMopUzTfBJ4d6EnXsa2HSfXRiVu9aRpKpWeWyKd41OnrzO/2rxKjEQmx/mJMLeX8784+n0GDx/uZ6DCUFQVpeUXRA8N+HjLoX7uGw/grnNlz26+f51Izq9zGYZWUU+r5gA3TfMJ4Engo5Zlbfjf3Otj4KVWk1n+9sXbrCSbW+43G0lyYSpCIpP/9hwf8nFqXw+Gtv0IWjPGwCulqyp3Dfu5dzzAscH6/CXRzWOoIOfXySodA68pwE3TfB/wLuDXgTHgsGVZZ0rvIwFebime5u9enCTa5GDM5mwuzkW5tphvr1tXedNoD4f6t77I2U4BXspn6Jza18ObxnrZ3+ep+iJtNwcAyPl1soYHeOGC5dPAS4VDfuATlmV9pvR+EuAbLUZT/N1LU8QzzQ/H5XiGV6YjLBUm/fR7De4f7910WKVdA7xUwGNwal8Pp/b1Mt7n2dVjuzkAQM6vkzWlB14JCfDNza2meOrsZEvWC3Ech9vhJK/PrJIs1GEfDHq4Z7QXn+vOOHMnBHipgMfAHOnh7pEeDgR37pl3cwCAnF8nkwBvglrfQK0MccgPq1gLMa4uxLAdUBU4OujDHPHj1rWOC/BSPkPnrhE/J4f9HB3woW8y3t/NAQByfp1MArwJ6vEGml9N8dTZ1gynFMVSWS7NRbkdTgKgqwonhnzcs78Ppwuu8huqypFBHyeG/JwY8tPryVfPdnMAgJxfJ5MAb4J6vYFCsTRPnZ1iNdXaxahWEhlen40yVyh1dOn5yo+jA95Ne7CdaqTHzfEhPw8eHSSg5hcH60bdHHDQ3ecnAd4E9XwDhRMZPnduiqV46zcMXoyleX1mde1Cp0tTOD7k59iQD1cXBbnLpUPO5vCAj6MDPo4M+hjqoolD3Rxw0N3nJwHeBPV+A8XSWb5wfpqZSLJuz1ktx3EIJbK8XlKxoqsKx4byQxFuvfODfLMx/l63zuEBH4f7vRwe8BH01n+v02bp5oCD7j4/CfAmaMQbKJ21+eorM1xrgyVXXS6dVCrDYiyNNR9jIZr/60BTFA4NeDk+5KPXXfVqDC1XyUXaPo/BwX4vh/u9HOz3Nn3BrVp0c8BBd5+fBHgTNOoNZDsO37k8z/nJlbo/926sD7hQIcjnSpYDGO11c2I4P/RQ7YSaVqmmysbv0jkQ9HAg6GV/n4fRgAe9TcfQuzngoLvPTwK8CRr9Bnr+1jL/cmURh9asR7JVwEWSGa4uxLkdTlBcOLDPo3NsyMeBoLdtA229epRJaorCaMDD/j4P44V/7TLs0s0BB919fhLgTdCMN9C1xRj/8OosqWzza8V3CrhUNseNUILroTipwoQgQ1U42O/lyICXvjYJsq00qs7dZ+iMBdyMBjyMBdyMBTxrpYvN1M0BB919fhLgTdCsN9BSPM2XLkyzGGtuhUqlAZezHSbDSW4sxcv25ez3GRwd8LK/TXvlzZyo5DN0RgNu9vXe+TfgMxo67NTNAQfdfX4S4E3QzDdQOmvzzYtzXJpbbcrrQXUBt5LIcHMpwcRygmxhfEVXFcb7PBzq97TVWHmrZ5oaqspQj4t9vW6Ge1wM+90M97rwu+rTW+/mgIPuPj8J8CZoxRvopYkw37+yQK4J63TXEnBZ22YqnOTGUqKsV+41VA4G8xUdgRYMK5RqdYBvxWfoDPW4GPa7GPS7GCp83O0wTDcHHHT3+UmAN0Gr3kCzkSRffXWW5QZP+qlXwK0ms9wOJ7i9nCCeuTM1v8+jsz+Yv/DXinLEdg3wrbg0lUG/i0GfiwG/wYDPxYDPoN/n2rQuv5sDDrr7/CTAm6CVb6B01uafrQVenm5cqWG9A85xHELxDLeXE0yFk2RK9r4MePS1So5m9cw7LcC343fp9HsNgj6Dfq9Bv8/g0L4AaiZLTwfX6m9HAlwCvCbt8Aay5qP808X5hiyG1ciAy9kO86spplaSzEZSZWHe49YY7/Mw2uALfd0U4Jspnp+hqvR5dfo8BkGvQaDwecCT/9jj1trmusRutMPPX6NIgDdBu7yBYuks37k0z+X5aF2ft1kBl7MdFqJppleSTEeSZHJ33pMuTSmr3HDVcQr/XgnwnaiKQq9bJ+DJ/+v1GATcOr1unR6Plv/o1ttuQ+x2+flrBAnwJmi3N9DluVW+c3mBWJ1CqRUBZzsOi9E0M5EUc6spYuny+vcBn8FowM1Ij5ugV6+p5ygBXjkFBZ+rGOYaPW4dv6vwuUvHX/joc2l1/SW7nXb7+asnCfAmaMc3UDKT4wdXF7kwGal5BmerA85xHKKpHLOrKWYjKUKxdNkZGZrCkN+VL8HrcdO7y6GAVp9fo7Xq/AxVxe/W8BkaPpeO36XhK/4zNLyFj8XPq10YrR1//upFArwJ2vkNNLWS4DuXFphdrX5lw3YLuEzOZn41zXw0xUI0vaF37tFVhnvyJXeDfoNe9/Y99HY7v3rrlPNTFQVvIdA9horX0Ar/VNyFjx49f1vpx31DPaxGEq1ufkNIgDdBOwc45Huw56dWePpqqKpt29o9AGLpLAvR9Nq/4nT+IkNTGPDlw3zQ76Lfa5Rt3tDu51ervXB+5Gzculr4p5V8ruLRVVy6iktTy+5jaArukuMuXW278X0J8CZo9wAvSmZyPHtjiRcnwruaANRJAeA4DqupfKCHYhlC8TTJTHmgKwoEvXfK7EYCHtwqHVmBUYlO+v5Vo57np6sqhqbg0vLB7tLzXxuailtT0Qu3FY8V76trKoZ655he8rlRuE1TlV2/xyTAm6BTArwonMjww6shXp9drWh8vJMDwHEcEhmbUCxNKJ4hFEsTSW48F11VCHoNgl6d/kINtc/VmWV163Xy968SnXR+xV8Quqqga2r+o3on8LWSr3VN4UOPnZAAb7ROC/Ci+dUUT18L8cbC9mWHnfQDUolMzmY5nmE5kWE5niGcyG46tGSoSr5G2luslc6X13XavqDd9v1br5vP70//w0MVBXh3TtES2xrpdfOBB8eZiSQ5fX1pxyDvFoamMtLrZqTXDeQDIBJLsZzIEI5nWE5kWY5nSOfsfK89Xr7JtN+lrQV6r0dfq4/u1k2RRfuTAN/DxgIePvDgOLORJGduLmPNR7GbsEhWO/EYGmOGxljAA+SHXlJZm5Vklkgyy0oiQySZJZLKEkvniKVzzERSZc/hd+XronsL9dHFcHdpux/7FGI3JMAFowEPT94/xnI8zQu3wrwyHSFj2zs/sAspioLH0PAYGvsKPXXITzCKprKsJLKsJLOsprJEk3dCPZbOsX6lX0NT1ia5+F3Ff/m6aI+hSriLmkmAizX9PhfvPTXCYycGOT+5wusLMRa6dIxxt1RFIeAxCHgMDpYct22HaDpHNJUP9Xyw51hNZcnknPx4eyKzyfOBryTQ/a78pBavoeIz8uVwEvBiJxLgYgOvofGOowP89AP7efGNBc5Ohrm1lGjZ3pztTC1c8Fy/gqLjOCSzdr53nsoRT5f31lNZm2gqRzS1eX2+olCY3HJnYouvEPDFrw0ZotnzJMDFllRVwdzXg7mvh6V4mpenIrw6HSEqvfIdKYXZhV5DY8i/8fZszi4L9Hg6RzyTI5HOkcjkSOec/LF0DtjYg4d8L96j52crenT1zkzFwuc9XhcaDu42nKgi6kMCXFRkwOfiiZNDPH5ikOuLcV6ZiXB1IUZ2j46V10rXVPq86pYbP2dtZy3M45kciYxd9nUyY5O1HeIZu2yTjK2UzlAsnYHo1pT8xzafmSg2JwEudkVVFE4M+zkx7CeZyWHNR3l9dpWJ5cSeq2BpJF1V8tUs22xukbVtkhmbZNYmlbFJZnNrXyczOVK5/C+BdM4mlbU3LDWwneJMQ7euYmgqrsLMQtfaTMTC53r5MU3p3pmt7UgCXFTNY2g8sL+PB/b3EUtnseajXJ6LSpg3ia6q9LhVetyb316c6GIXSiNTWZt04WMqd+fzYsDnv3ZI52wyOYdMLrdhwbCdqAplU82L08n1kmnmxdmGhqquHV9/m/wFUBkJcFEXfpfOWw4EecuBIIlMjqsLMd5YiHI9FCedk2GWVlJLxuMr4TgO6ZyTD/RC0GcKoX4n3G3ShY+ln9sOhV8WALtfQK1IU1ibcl6cgq6pKprK2ufFZWg1pXj7xo/rj6ld9heCBLioO6+hcd94gPvGA2Rth9vLCa4uxri2GGOpwRsxi9opioJbz6/Y17vLx+bsYsjfCfVsziFrO2RyDlk7/3XGXvf5uttyDuSyNqmdX3LX1sJdyV+o1xRl7aOm5n/haWXHKbk9/0tg7fPC7fnj6x6jKKiFxdI0pTG/PCTARUPpqsLRQR9HB338pDlMOJHhRijOzaU4t5YSDdnLU7SOpip41cp7+5txHIeckw/1TCH8s3b+om3OdtY+oiikMrmyY2ufO+X3LX60HQrP17ohvuIvAKUQ8Goh3Es/VkoCXDRV0Gvw5gN9vPlAH46T3wtzYjnBxHKcieWkBLpAURR0RUFXwbN5kQ5Q3WJWtnMn6PMhf+eYXfg65zjYhdvt4l8Da7fnv7ZLfknY654n5+R/UdiFx6997hTbQOEakVPLKBNQY4CbpvkTwM8A84BjWdb/qq05Yi9RFGVtcamHDwUBCMXSTIYTTIaTTK8kWYylZQKRqBtVUVA1hRr+QKiJXRLujuNgF39hOPm/PGyHXa3ZX3WAm6bpA/4CeJNlWSnTNL9smuZ7LMv6frXPKUR+OzQXD+zvA/IXxGYKu9XPRJLMRlKsJDef2CJEuysOk0B9xsJr6YG/HbhlWVbxOsOPgPcBZQFuGJoyPLzbSyGdo5vPDdrj/A6M9bW6CUK0pVpWqB8BStdfixSOCSGEaIJaAnweyqqMAoVjQgghmqCWAD8DHDZNszgP7MeAb9beJCGEEJWoaU9M0zR/Evg5YAHISBWKEEI0T8M2Ne7mEkPTNEeBjwMPWJb1SKvbU2+maR4nf37ngANAyLKsj7W2VfVjmqYKfB14HnABx4H/bFlWoqUNqyPTNL3kz++7lmX9cqvbU0+maT4HJAtf5izLek8r21NvpmmawAeBBPA48NuWZb2w2X0bMpFnD5QYvhP4GvBgi9vRKAPA5yzL+hqAaZoXTdP8pmVZZ1vcrno6Y1nWxwFM0/wa+c7GZ1vbpLr6OHC+1Y1okG9blvXbrW5EI5imqQF/ArzfsizbNM2/BbacrdSomZgVlRh2KsuyvmSa5rtb3Y5GsSzrxXWHVCDWirY0gmVZNvmAwzRNnfxfGVZLG1VHpmn+PPmfufuBnhY3pxHuM03zVwEv8KJlWd107e0R8kXiHyl0hEPAp7a6cy0XMbcjJYZdwjTNJ4HvWJZ1udVtqTfTNN8LfAP4hmVZL7W6PfVgmuY9wCnLsr7S6rY00B9YlvUHwO8Av2Ga5mOtblAdHSbfAf6MZVm/BzwG/Met7tyoAJcSwy5gmuYTwBPAf211WxrBsqzvWJb108BR0zT/S6vbUydPAknTNH+N/FDfo6ZpfrS1Taqv4niwZVk54Bny79FuEQEuW5a1Uvj6NPDure7cqACXEsMOZ5rm+4D3Ar8EjJqm+fYWN6luTNO8p3B+RTeAY61qTz1ZlvW7lmV9zLKs3yf/w/+CZVl/1uJm1Y1pmnebpvnhkkMngWutak8DPA8MFsbCId8jv7LVnRtZhdK1JYamaT4OfAj4aeDPgT/usgqGh4CngeKwgh/4hGVZn2lZo+qoUGXzh+SrbAzgFPCLlmXNtrRhdWSa5s8Cv0C+yuYTlmX9fYubVBemaY4D/4f8BdoA+e/ffytc1+gKhWHLHyefnYeAj2yVLw0LcCGEEI3VqCEUIYQQDSYBLoQQHUoCXAghOpQEuBBCdCgJcCGE6FAS4EII0aEkwIUQokP9f8QJb4GGzJhRAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pyplot.fill_between(coordinates, mean-std, mean+std, alpha=0.6)\n", "pyplot.plot(coordinates, mean)\n", "\n", "pyplot.axis([0, 6, 0, 10])\n", "\n", "pyplot.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Pseudo-spectral projection\n", "\n", "Implementing g-pce for pseudo-spectral projection require us to generate\n", "nodes and weights from $R$ and transform the nodes using $T$:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:18.813711Z", "iopub.status.busy": "2021-05-18T10:57:18.813389Z", "iopub.status.idle": "2021-05-18T10:57:19.008625Z", "shell.execute_reply": "2021-05-18T10:57:19.008280Z" } }, "outputs": [], "source": [ "nodes_r, weights_r = chaospy.generate_quadrature(10, distribution_r, rule=\"gaussian\")\n", "nodes_q = transform(nodes_r)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The resulting samples can then be used to solve the equation above using the\n", "quadrature-based method:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:19.010969Z", "iopub.status.busy": "2021-05-18T10:57:19.010656Z", "iopub.status.idle": "2021-05-18T10:57:19.639804Z", "shell.execute_reply": "2021-05-18T10:57:19.639438Z" } }, "outputs": [], "source": [ "expansion = chaospy.generate_expansion(7, distribution_r)\n", "evaluations = numpy.array([exponential_model(sample) for sample in nodes_q.T])\n", "model_approx = chaospy.fit_quadrature(expansion, nodes_r, weights_r, evaluations)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that for generating the expansion and the model approximation, we use\n", "the nodes and weights from $R$, while for the model evaluation we use the\n", "transformed samples from $Q$.\n", "\n", "The solution model, defined with respect to $R$ can then be used to do\n", "analysis:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:19.642216Z", "iopub.status.busy": "2021-05-18T10:57:19.641864Z", "iopub.status.idle": "2021-05-18T10:57:19.688007Z", "shell.execute_reply": "2021-05-18T10:57:19.687722Z" } }, "outputs": [ { "data": { "text/plain": [ "(array([10. , 9.89956, 9.80022, 9.70198, 9.60482]),\n", " array([1. , 0.98159, 0.9644 , 0.94841, 0.93361]))" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mean = chaospy.E(model_approx, distribution_r)\n", "std = chaospy.Std(model_approx, distribution_r)\n", "\n", "(mean[:5].round(5), std[:5].round(5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting the final results:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:19.690342Z", "iopub.status.busy": "2021-05-18T10:57:19.690028Z", "iopub.status.idle": "2021-05-18T10:57:19.743903Z", "shell.execute_reply": "2021-05-18T10:57:19.743593Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD3CAYAAAAE2w/rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtcElEQVR4nO3dZ3Aj6X3n8W8HZIAEM4ecHLZ3Vpu1uwqWtFrJtuxT6a7Wtu5OL+wLencu2bo7uxzK5/PJctkqx6s6OZxsl84lrZVlWZYsyZLl1a52Nk3YNDM9OTCTIEEQOXTfiwY4YAaRAf4/VVMkG2Dj6RnMjw+ffp7/o9i2jRBCiM6jtroBQgghqiMBLoQQHUoCXAghOpQEuBBCdCgJcCGE6FAS4EII0aH0nZ5gGMYo8HHgAdM0Hy0e8wJ/AEwCJ4DfM03zUiMbKoQQYq1KeuDvAL4GKGXHPgrcMk3zd4E/Bv6q/k0TQgixnR0D3DTNLwEr6w6/HzhVfPw14AHDMHrq3zwhhBBb2XEIZQvDrA31WPFYbP0T3/zb/2T/6MnhO98Y8vAL7zlR5cu2F01TKRSsVjejYeT6OptcX+dyuTRl52dVH+BzQKjs657isQ0WE1kSqSwuzensT0TyTMzGCHqqfen2EQ77iUaTrW5Gw8j1dTa5vs41NBTa+UlUPwvlG8DbAAzDuA94xTTNDb1vABtYSGTXHLu52J1/6UII0Uw7BrhhGI8DPwvsMwzjNwzD8AH/GzhkGMZvAP8d+PB255iPrw3w6xLgQghRsx3HMUzTfBp4epOHfr7SF5lbWRfgkVSl3yqEEGILDV/Io6kKK5k86Vxh9dhKJrdhWEUIIcTuNDzAh4JuYJNhlEii0S8thBBdreEBPtLjBWBuXYBfW5BxcCGEqEUTAtwDwHw8Q/nuP7ejKfKW7AYkhBDVaniAh30u3JpCKmeRyN4ZB88WLG4vyc1MIYSoVsMDXFEUhoKlXvjaYZSrCzIOLoQQ1WpKOdnSjcy5eGbNcQlwIYSoXlMDfD6eXTMOHklmiSZzzWiCEEJ0naYEeMCtEXBr5Ao2S6m1gX1FeuFCCFGVpgS4oigMh5xx8PWrMiXAhRCiOk3bUm2kOIwyu7J2HPzmYpJsvjtLQgohRCM1LcAHg24UBRaTuTWBXbBtrsmqTCGE2LWmBbhLUxnwbz4b5dK8BLgQQuxWU3elHwmVhlE2zge3bFmVKYQQu9HkAC/dyFy7rD6VK8iqTCGE2KWmBniPV8erq6TzFrF0fs1j5ly8mU0RQoiO19QAL59OuH42ijkn4+BCCLEbTQ1w2HocfCWTY2o53ezmCCFEx2p8OdmQd83Xw8XCVpFkllxh7fzvizKMIoQQFWt4gN89GlrztVtX6fO7sO2N1Qkvzq40ujlCCNE1Gh7gJ/eFNhwbLY6Dz8TWjoNHUzmmYzKMIoQQlWh4gI+HfYQ8rjXH9hV36ZlZN50Q4MKsDKMIIUQlmrKhw11DgTXHerw6PpdKJm9tqE54YUaGUYQQohJNmYVy13BwzdeKojBa3Ox4/TDKcjrHRFQW9QghxE6aEuAH+3x4dG3NsdVhlHUBDnBeeuFCCLGjpgS4piocH/SvOTYYcKOpCsvpPMmyzY7BGQeX2ihCCLG9pi3kOTG0dhhFUxWGizXCZ9atykxk81yPJJvVNCGE6EhNC/BjgwE0RVlzbHR1GGXj1ME3pmUYRQghttO0APfoKof61w6jlOaDz8ez5K21qzIvzcdlpx4hhNhGU2uh3DW8djqh16XR53dh2Rv3yswWLFlaL4QQ22hugA8FUVg7jLKv2Auf3mQ2ymtTsaa0SwghOlFTAzzo0VfHvUv29d4ZB18/8+TWUorouoU+QgghHE0vJ7t+UU/IoxN0a2QLNpHE2mEUG1t64UIIsYWmB7ixblm9oiiM9TqrMqeWNw6jvDoV21AvRQghRAsCfDDooa+4O33JaoDH0hvCejmd48aiLK0XQoj1mh7gsLEXHvY5xa3SuY3FrQDOTS43q2lCCNEx9Fq+2TCMXwYOAwvACeDDpmnu2F2+azjI8zeXVr9WFIWxHi9XI0mmljP0r+uhX5qLk8wW8Lu19acSQog9q+oeuGEYo8CvAR8xTfN/AgHgpyr53vFeL37X2p8dd8bBNw6jFGybV6ekFy6EEOVq6YEngSzQA0SBIPDG+idpmko47F9/mAcO9XH61p1e+GifhkdXSWQLpCxnWKXchYUU73vAh7JuOX4rbXVt3UKur7PJ9XW/qgPcNM1YcQjl84ZhTAMTwJX1zysULKLRjYWp9gddnMrm1xwb7fFwczHFrYUE/pG10w2ns3nOXYtwZKB9/sHCYf+m19Yt5Po6m1xf5xoa2rgV5WZqGUJ5EPhl4P2maf5HnHHw36z0+4/0+3Gpa19+rLjJw+Ty5vtinpmIVtNUIYToSrXMQhkHFk3TLHWjpwFvpd+saypH19UIHwq60VWFWDrPSia/4XsuzyeIpWVlphBCQG0B/i3gDcMw/tAwjP8BPAr87m5OsFmN8LHi0vrJ6MZeuGXbnJ2Qm5lCCAG1jYEXgJ+v5cVPDAVQFWVNDZTxXh+3ltJMRNPcvW4cHODsRIwfOTqArrbPzUwhhGiFlizkKfG5NPaHfWuODYfcuDSFlUx+0+GSZC4ve2YKIQQtDnCAu9atylQVhfHinPCJTYZRAF6+FW10s4QQou21QYBvHCbZH74T4JsVsppZSXNrSeqjCCH2tpYHeNjvYji4tkb4YMC9uqgnmto4GwXgxbKl+EIIsRe1PMBhY41wpWwYZas54ZfnEywms5s+JoQQe0F7BPi6cXAoH0ZJbTqMYmPz0s1oo5smhBBtqy0CfLTHS6/XteZYv9+Fz6WSylksJjdfvPPqVIxUrtCMJgohRNtpiwAHZ054ufJhlNvRzW9Y5ixLZqQIIfastglwY3jjbJSDfc4c8YlomoK1+bZqL9+OkitYDW2bEEK0o7YJ8AN9PnyutRs29Ppc9Hp1cgWb2ZWN+2UCpHIF2bFHCLEntU2Aq4qyYRgFnGAHtp33/fyN6JY9dCGE6FZtE+Cw+aKeA8XZKDMrGTL5zYdKVjI5Xp+ONbRtQgjRbtoqwI8O+HFra5vkdWkMB93Y9tZzwgGeu7G0piiWEEJ0u7YKcF1TObrJjjsHKxhGWUpmuSBFroQQe0hbBThsXJUJsK/Xg6YqLCVzxDfZ6KHkh9cXN130I4QQ3ajtAvz4oFMjvJyuqowXN3rYrhe+kMhyYTbe0PYJIUS7aLsA97o0DvdvP4yyXS/72WsR6YULIfaEtgtw2HxRz2DAjd+tkcpZzMW3LmK1kMhyXnrhQog9oC0D/K7hAAprh1EUReFQsRd+c3H7WuDPXI3IjBQhRNdrywAPuPXVBTzlSgE+FUtvOSccYDGZ5fVpmZEihOhubRngAMbwxlWZPrfGSMiDbW9/MxOcXnheVmcKIbpYGwd4cMMwCsDh/jvDKNvdrFxO53hFaqQIIbpY2wZ4j9fFWLGcbLnRHg8eXWUlk9+yTnjJs9cWyW4z1CKEEJ2sbQMcNh9GURVldUrhTjczE9k8L0m9cCFEl2rrAL97JLTp8dLNzInl9I61wE/dWCSZlV17hBDdp60DPOxzsa9n4zBKyKszGHBTsGxu73AzM1uwePZapFFNFEKIlmnrAAe4e2Tjoh5gtejVtcj2NzMBzkwsyw72Qoiu0/YBfnJ482GUfb13bmYuJLYPZ8u2+f7lhUY0TwghWqbtAzzsdzEa2jiMoioKR4pTCq9Hth9GATDn4jvOHRdCiE7S9gEOcHJ082GUwwN+FGBqOU0qt/ONyu9dmpdCV0KIrtEZAb7FMIrPpbGv14PNzlMKAaZjaV6dkq3XhBDdoSMCPOzffDYK3LmZeT2SrKiA1b9ciWxbR0UIITpFRwQ4wD1bzAkfDLgJejTSeYvp5cyO50lk8/xQphUKIbpAxwT4ydHNa6MoirLaC7+ykKjoXC/dihLZYeaKEEK0u44J8B6vi/FNaqOAszLTpSksJnMVzfcu2Db/ZM7Xu4lCCNFUei3fbBiGAXwISAGPA79lmuaL9WjYZu4ZDTGxvPFmpa6pHO73c3k+wdX5JP2H3Due61okgTkbx9hioZAQQrS7qnvghmFowB8BHzNN8xPAh4Hr9WrYZu4eCW7Y8LjkaHFK4eRyuuLaJ98x56VaoRCiY9UyhPIooAAfMQzj14APAA1d7hj06KuVCNfzuzXGe73YwLVIsqLzrWRyUidFCNGxahlCOQS8DfiQaZrLhmF8BsgCny5/kqaphMMbd5mv1luODzF1bnLTx+7e18PEcpobi0nu29+LS9v559O5mThvv3uE0S2mKW6n3tfWbuT6OptcX/erJcBjwEXTNEvb3jwLvJt1AV4oWESjlfWIKzHu1ynkChQ2mfMdcqv0+10sJnNcmV3h2ODGeuKb+fzzN/m5R/ejbDE8s5Vw2F/Xa2s3cn2dTa6vcw0NbT5ter1ahlBeAAaKY+Hg9Mgv1XC+inhdGseHtg7m0mNXFipb2AMwuZzizIRsvyaE6CxVB7hpmovArwB/YhjGbwJDwB/Xq2HbuXe0Z8vHxno8BN0ayWyByWi64nN+//ICK+l8PZonhBBNUdM0QtM0vwp8tU5tqdixoQBeXSOd3zjbRFEUTgwHODsR49J8gv1hb0VDI9mCxbcuzvHBB8ca0WQhhKi7jlnIU05XlS03egA4GPbhc6nE0nlmYjsvry+5PB/n/MxKPZoohBAN15EBDnDvvq2HUVRV4XjxBqY5n9hVCdnvXJyXPTSFEB2hYwP8QNhL2Ofa8vHDAz7cmsJSMrfjjj3lkrk837k4V48mCiFEQ3VsgCuKwpu2uZmpq+rqNMJLc5UVuSo5P7vCxVkZShFCtLeODXCA+8a2nyt5dMCPrirMxbMsJXO7Ove3LsyTyMqsFCFE++roAO/3u9nfu/nSegC3rnKkWGrWnIvv6tzJXJ5/PC9DKUKI9tXRAQ5w39jWwygAxwf9aApMxzK77oVfmo/LFmxCiLbV8QF+z2gIXd36MrwubbUXfmF2d71wgO9cnCOa2l3wCyFEM3R8gHt0lbuHt6/pfddwEE1VmF3JsLjLnXiyBYuvvz5T8bJ8IYRolo4PcID7dxhG8egqxwar74XfjqY4dX2pqrYJIUSjdEWAH+r3bTsnHODEYGB1Rspu5oWXPHMtwkR0425AQgjRKl0R4Iqi7NgLd+sqx0u98Jnd98It2+Zrr82QzskqTSFEe+iKAAdnGGWzXevLHRsK4NIUFhJZ5uOV10gpWU7n+Ob52WqbKIQQddU1Ad7jdXF0YPvdOdyayoni6sw3puO7qpFScnEuzunb0WqaKIQQddU1AQ7w0P7eHZ9zbMiPR1dZSuWYXK68Xni575rzzMSq+14hhKiXrgrw40MBQp7tS5zrqsrJYinaN2biFKzd98ILts1XXp2W8XAhREt1VYCrFdzMBGfWSsijk8wWuF7hDvbrRVM5vnxmsqphGCGEqIeuCnCAB8d7d7yZqSoKb9rn9MIvzsXJ5q2qXuvCTIznbsj8cCFEa3RdgPf6XBXtRj8a8jAYcJMr2Fya31252XLPXI1U3YsXQohadF2AAzx8YOebmYqicO8+pxzt1YVE1bvwWLbN3702TXSXhbKEEKJWXRngxwb8O67MBOjzu9gf9mLZ8Pp09Rs4pHIFvvTKFLlCdUMxQghRja4McEVReLiCKYUAbxoNoSkwuZyuanFPyVw8w9dfn5GbmkKIpunKAAd4YLx32zKzJX63xl3FaoavTK7UVHXw4lycZ68tVv39QgixG10b4D6XxptGt99yreTEUICAW2Mlk+dajTckn722yAXZT1MI0QRdG+AAjxwMV/Q8TVW4r3hD8+JMnEy++gU6NjZff32WqSpXeQohRKW6OsBHQh4OhLfeM7PcaI+HkZCbnGXzxvTuqxWWy1sWXzw3JTv5CCEaqqsDHOCxg30VPU9RFO4b60FR4OZSatc796yXyOb54tkpWW4vhGiYrg/wu4YDFU0pBAh59NVqhWcnYzVvozafyPDlV6bJV1FvRQghdtL1Aa4oCo8cCFf8fGMkSMCtEUvnuVzDCs2Sm0tJ/kGmFwohGqDrAxycKYUeXavoubqq8OC4UxDr4myceCZf8+ufn13he5cWaj6PEEKU2xMB7tFVHhrfuUphyXDIw4E+Z4Xm2YlYXXrPL95a4tQNmSMuhKifPRHgAI8e7ENVtq9SWO6+fT24i9uv3Vqqz2bG37+8wCuTy3U5lxBC7JkAD3n11eJVlfDo6mpt8demV+o2m+QfL8xxURb6CCHqYM8EOMBbD/XtWCu83P6wl+GgU3L27GR9hlJKu9tfWaj9BqkQYm/bUwE+GPRwfGjnWuEliqLw0P5edFVhJpbh1lJ9VlcWbJuvvDLNjUWpIy6EqN6eCnCAtx+pbGFPid+trQ6lvDoVq7pu+Hp5y+JL56bqNr4uhNh7ag5wwzB8hmG8ahjGH9SjQY023uvjcL9/V99zsM/Lvh4PecvmzMRy3eZ0ZwsWXzg7yURUQlwIsXv16IF/HDhbh/M0zdsP9+/q+YrizA13awrz8Wxdt1DLFiw+d0ZCXAixe0otvUnDMH4WSAD3A0HTNH9p/XMsy7YLbbhTzf995hq3djkGfXspyXNXF9FUhR8/OUzY7655uX2JW1f5ubce4vBA5WP0jaZpKu34b1cvcn2drZuvz+XSKpptoVf7AoZh3AOcNE3z1w3DuH+r5xUKFtFo+92se/O+EFdmYrv6npGAmwNhL7ejaZ67GuHH7hmhUEPp2XLZLHzqX67ywQfHODywuyGeRgmH/W35b1cvcn2drZuvb2iosinPtQyhPAmkDcP4VeAdwGOGYXy0hvM11bHBAOO9lZWaLffAeA8Bt8ZyOs+5iWhd25SzLL5wborL87WVsxVC7A1V98BN0/yd0ueGYXhxhlD+pB6NapZ3HO3n82cnd/U9Lk3l0YNhnr4a4cpcggGfi7Feb93alLcsvvzKNB+4d7TiHYWEEHtTPWah/DTwLuCthmF8qPYmNc+xwQD7q+iF9/ld3FsM1zMTy3WbWlhi2TZ//9oMZ+rcwxdCdJeqe+Alpml+GfhyHdrSEu86PsBTpyd2/X3HBv0sJHNML6d56VaUdx7r31WtlZ3Y2HzrwhzJbIF3HB2o23mFEN1jzy3kWe9wv3/X88LBmVr42OE+vLrKYjLH69ONqW/yg6sRvnVhTuqJCyE22PMBDvDu44NVfZ/XpfHYoTAKcHUh2bC53Gcmos7OPl06ZUoIUR0JcGCs14sxHKzqewcCbu4fK46H346x3KCNjC/Nx/ns6UkS2do3mBBCdAcJ8KLHjw9WPYZ9ZMDPgbCXgm3zws0o2Qb1lCeXU/y/F2+zUOOGy0KI7iABXjQYcK8WrdotRVF4cH8vvV6dRLbAy7fqVy9lvWgqx9+8eLuuy/mFEJ1JArzMu44N4Naq+yvRVYW3HA7j0hRmVzINu6kJkM4X+PzZSV6+FW3Yawgh2p8EeJmgR+cth3ZXbrZcwK3zluJNzSsLSW40sJds2TbfMef45vlZ8pbMUBFiL5IAX+eth/sIeaqfHj8U9Kzuan9uMsZ8PFOvpm3q3OQyT708QTwjNzeF2GskwNdxaWrV0wpLDg/4OT7oxwZeuBllpcHhOrGc4q+fv8Vt2RxCiD1FAnwT9+4LMdZTW32Te/eF2NfjIVewOXV9qW6bIm8lns3z2dMTvHhzqaGvI4RoHxLgm1AUhR+/e3hXGyBvdo5HDvQS9jkzU07dWCLX4IU4lm3z3UvzfOWVKTJ5WfQjRLeTAN/CWK+36mmFJbqm8rbDfQTcGtFUnhduRik04Ybjxbk4f/38TaZj9dmEWQjRniTAt/HuEwN4da2mc3hdGm8/0odHV5mPZzl9u3FzxMstFeeLv3BzSeqoCNGlJMC3EXDrPH689kqAQY/O24/0oasKk8tpXp1aaUqoFmyb712a5wtnp2SWihBdSAJ8Bw/v72VfjTc0AcI+F289HEZV4FokyevTzQlxgKuRBH956pbs9CNEl5EA34GiKPzkyeG61PoeCnqc6oWKs9Dn/Ey8aSGezOX54rkpvnl+lqzc4BSiK0iAV2C0x8tjB8N1Ode+4rkU4NJ8gouzze0Vn5tc5lOnbnJzUWqpCNHpJMAr9M5jA4R9rrqca6zXyyMHewG4OJfAbHKIL6dzPHV6km9fmJPeuBAdTAK8Qi5N5V+dHKnb+faHfTxywAnx87NxLsw2bzgFnC3bTk9E+dSpm1LZUIgOJQG+C4cH/Dy0v7du5zvQ5+PNxRC/OBtv6o3NkuV0jr89M8Hfvz5T982ZhRCNJQG+S+85MUSvtz5DKQAH+3zOmHjxxua5yVhL5m2/Ph3jL567wSuTzZmnLoSonQT4Lnl0lfffM1LTMvv1xsNe3nqoD1WBG4spXr61jNWCEE3lCnzj/CyfeXmi4VUUhRC1kwCvwuEBP4/UaVZKyWiPhx850o+uKkwsp3n+RrRlmxjfjqb4q+dv8c3XphtehEsIUT0J8Co9cXyA4ZCnruccDLp5x9F+3MVdfZ65ttiyALVsm+euRfjzH97k3IQMqwjRjiTAq6RrKh988wG0OizwKdfnd/H48YHVAlj/ciVCLN26ZfDJXJ5vXpjlr1+4JXPHhWgzEuA12Nfr5d0natv8YTNBj1ODpc/vIpWz+MGVSMvHpGdXMnz29ARfPDfFQiLb0rYIIRwS4DV67GCY44OBup/Xo6u882i/symEZfPD60ttMV/78nycvzx1k2+en2Wlhb8ZCCEkwGumKAofuHeUkKd+UwtLNFXhLYfCHBv0Y9vOHpvnJmNYLd7E2LJtzk0u82c/vME/X14gJTc6hWgJCfA68Lk0nrx/tC4Fr9ZTFIX7x3p4eH8PqgLXI0mevb5IJt/60MxbFs/fWOSTz1znB1cjMmNFiCaTAK+T/WEfTzRgPLzkUL+fdx7rx6urRBI5vn85wlIy17DX241sweLZaxH+9NkbPHtNglyIZpEAr6O3HOrj5EioYefv97t54sQA/aWbm1cjXF1ItM0Uv3S+wA+uRvjkszd4+kpEhlaEaDAJ8Dp7/z0jDAXqOz+8nNel8Y6j/RwZ8GPZ8OrUCi/ejJJt0aKfzWTyBX54PcInn7nOd815Yun2+E1BiG4jAV5nbl3lZx7ch89V216a29FUhQfHe3jsYBhdVZiKZdpqSKUkW7B48dYSf/rsDb7++kzLp0IK0W0kwBugz+/myfv2NeSmZrnxsJcnTgwQ9ukkswWevhrBnIu3pI7Kdizb5rXpGJ86dZPPnZlsq2EfITqZBHiDHB7w8+N3DzX8dYIenXcdG+DYgDPV8PxMnB9cXWzbTYyvRRJ8/uwknzp1kzMTUdlQQogaKNX2hAzDOAZ8HDgD7Acipml+bP3zcrmCHY22fgFKI4TDfna6tu+a87x4a6kp7ZldyXBmYpl0zkJTFO7dF+LIgA+lyt8E3G6dbLaxPwg8ulacJtnLQMDd0Ndar5J/v04m19e5hoZCFf2nraUH3g98zjTN3zdN8xeBf28YxptrOF9Xeu9dgxjDwaa81kjIw3tPDHIg7KVg27wyFeO560skGhzCtcjkC7x0a4m/eO4GT52e4PzMCvkWL1QSolPo1X6jaZovrTukAonamtN9FEXh39w7ylOnJ5lYTjX89dy6yiMHw+zrSXN2cpm5eJbvmRFOjgY5Nuhv+Lh8LW4sJrmxmMTv0rl/LMT9Yz0MBhs3o0eITlf1EEo5wzCeBN5d7ImvYVm2XWijKW71pGkqlV5bKlvgU89eY26leTMxUrkCZ29Fub3k/ODo87t45FAf/RUOVaiK0vIbogf6/Dx8MMy9Y7343PWd2bObf79OJNfXuVwuraKeVs0BbhjGE8CTwEdN09zwt7nXx8DLraTz/M1Lt1lu8rzomViac5MxUjnnn+fYoJ+TI0Fc2vYjaM0YA6+UrqqcGApw774QRwcCaGrtv0l08xgqyPV1skrHwGsKcMMw3g+8E/g1YB9wyDTNU+XPkQBfazGZ5TMvTRBvcjDmCxbnZ+NcXXDa69FV3jQa5GDf1jc52ynAy/lcGnePBHnTaA8Hwt6qb9J2cwCAXF8na3iAF29YPg28XDwUAD5pmuany58nAb7RQjzDZ16eJJlrfjguJXO8OhVjsbjop8/n4v6x0KbDKu0a4OVCHp27R0KcHAky3ru7MO/mAAC5vk7WlB54JSTANze7kuGp0xMtqRdi2za3o2nemF4hXZyHfSDs5Z7REP6yceZOCPByIY+OMRzkrmHnN4udbth2cwCAXF8nkwBvglrfQK0McXCGVcz5BFfmE1g2qAocGfBjDAfw6FrHBXg5n0vj+GCAu4aCHBnw49Y3jvd3cwCAXF8nkwBvgnq8geZWMjx1ujXDKSWJTJ4Ls3FuR9MA6KrC8UE/94z3YnfBXX5NUTjc7+f4UIBjgwHCPmfzjW4OAJDr62QS4E1QrzdQJJHlqdOTrGRaW4xqOZXjjZk4s8Wpjm5d5a6hAEf6feg7zFjpJAN+N8cGAzxwZIA+ja66tnLdHHDQ3dcnAd4E9XwDRVM5PndmksVk6zcMXkhkeWN6ZfVGp1tTODYY4OigH3cXhZ3brWPlLQ70+TjS7+Nwv5+RkKfqWS3tppsDDrr7+iTAm6Deb6BENs8Xzk4xHUvX7ZzVsm2bSCrPG2UzVnRV4eign+ODATybjCl3ms3G+H0ujUN9Pg71+znU5+volaDdHHDQ3dcnAd4EjXgDZfMWX311mquR1lclcLt1MpkcC4ks5lyC+bjz24GmKBzs93Fs0E/IU3U1hpar5Cat36VzsM/LgT4fB8I+hkOeti5HUK6bAw66+/okwJugUW8gy7b59sU5zk4s1/3cu7E+4CLFIJ8tKwcwGvJwfMjPYMDdcUMP1cyycWsqY71e9od97A97Gevx4m3g5h216OaAg+6+vkoDvHO7T11MVRR+8uQI/X43/3xpAZv2qM43EHDz9iNuYukcV+aT3I6mmFnJMLOSoderc3TQz/6wD70Oy9zbVbZgrRbdAlBQGAi4GOv1On96vAwFPXVZ6i/ETqQHXoNm9ACuLiT4u9dmyOSbP1d8px5qJl/geiTFtUiSTHFBkEtVONDn43C/j97idL121ah57pqiMBzyMNbjZbTHw0jI05JQ7+YeKnT39ckQShM06w20mMzypXNTLCSaO0Ol0oArWDYT0TTXF5Nr9uXs87s40u9jvE175c1cqKQpCoNBNyMhD8NBJ9SHQ56G7p3azQEH3X19EuBN0Mw3UDZv8Y3zs1yYXWnK60F1AbecynFjMcWtpdTqxgy6qjDW6+Vgn7etxsrbYaVpyKMzFPQwHHQzEHAzFPQwGHBvunJ0t7o54KC7r08CvAla8QZ6+VaU712ap9CEOt21BFzespiMprm+mFrTK/e5VA6EfRzo89Hjbe0tmHYI8M0oKAQ9OoNBNwN+F4NBN/1+NwN+N6Fd/J11c8BBd1+fBHgTtOoNNBNL89XXZlhq8KKfegXcSjrP7WiK20tpkmV1X3q9OuNh5+ZfK6YjtmuAb8etqfT5XfT53fT7XKufh306IY++5rebbg446O7rkwBvgla+gbJ5i38y53llqnFTDesdcLZtE0nmuL2UYjKaJle292WPV2e8OJOjWT3zTgzw7WiKQq/PRbj4Z3woiJYv0OvT6fW6CLi1thm+qgcJcAnwmrTDG8ici/OP5+caUgyrkQFXsGzmVjJMLqeZiWXWhHnQozHW62U05KHf72pY6HRbgK+3/vo0RaHHq9PjdRHy6vR4dEJep+de+uh3a7JQqQ1IgDdBu7yBEtk8374wx8W5eF3P26yAsyybuXiWqeU0U7E0ucKd96RbUxgJeVb/1OPm3uq591iAV0JVFAJujaBHX/Mx4NYJeDT8Ls352qPj1dWW9ujb5f9fI0iAN0G7vYEuzq7w7YvzJOoUSq0IOMu2WYhnmYk5C4QS2bXz3/v9LkZ7nKl4YZ9eU4BIgNdGVRS8uobfreJ36/hcKj6XVvzjfO51aXh1FW/xmEfXcGtKXYK/3f7/1ZMEeBO04xsonSvw/SsLnJuI1byCs9UBZ9s28UyBmZUMsysZFuLZNVfkUp251UNBZ/pdyLO7Md5WX1+jtev1qYqCW1Px6Bv/uHUVj6bi0pzP3ZqKS1OKH53PdVXBpakM9gdIxtPoqnOsm8b3JcCboB0DvGRyOcW3L8wzs1J9ZcN2C4BcwWIunmVuJcN8PLuhd+7RVYaCbgYDbgYCrg2zMtZrt+urt712fZqioKkKuqoWP4KmOsdURVl9XFUVVEBVFTQFFEVBLf+IgqJA6a2joLDZu8iG1U5SKUZt2/ktsvTRwumIWLYzVFhY93nBKv4pfp4vfv2Jf/eQ1ELZy8Z7ffyntxzg7OQyT1+JtGzbtnpyaSrjvV7Ge70AJLMF5uJOmM/Hs2TyFhPRNBPFnYVcmuLMnw64GAi46fO5pEZJFyvYNoWCTbYLdpGqlAR4F1MUhYf3h7lnJMRz1xd56Va0KQuAmsXv1jjc7+dwvx/btlnJFJiPZ1hM5Igks6RyFrPF4RdwelRhn4u+4vzp4R4vHpWu+tVb7C0S4HuA16XxnruGePhAmB9cifDGzErbVDisF2V1ipzOsUHnWDJbIJLIEknmiCSyxNJ5lpI5Z2VoBGAZXVUI+3Qn2P3O/Olumy8tupcE+B4S9rn41/eN8tbDfTx9NcLl+fpOO2w3freG3+0s2wdnDH0pmWMp5YT4cjpPMltgIZFjIXFnub+uOj8Men0ueos/FHq8Oq4u2k5OdAcJ8D1oOOThgw+OMR1L88zVRa4sdHeQl7g0leFiFUBwboLFEhmixUBfSuWJpnJk8haLydzqVnIlAbfmBLv3zsKXoEeXcXXRMhLge9i+Hi//9qExZmJpTt1YwpyLY3XRGHklvC6NUZfGaI939Vg6VyCWzrOczrOczhFL5VnJ5ElkCySyBaZjmTXn8Ls1Qh5n0Uv5qsZ6zXcWYisS4ILRHi9P3r+PpWSWF29GeXUqRs7aO3fy1/MWF6CUeurgTAmLZ/Isp/LE0k6gr2TyJDIFklnnz+zK2uJiLk25s4rRrRHwaKtf+1ytXcUouoMEuFjV53fzvpPDvOv4AGcnlnljPsF8F88j3g1VUejxuujxrt1lyLJs4tkC8cydUI+nC6xk8uQKNtFUnmhq49+hqjg991K4+93OMnWfS8Xn1lq+TF10BglwsYHPpfH2I/38xAPjvHR5ntMTUW4uprpu5ko9qOqd2S/lbNsmk7dWh10SZUMwiWyBTN4inikQz2w+P1/B+XfwudVisGv43Br+siXqMkQjJMDFllRVwRgJYowEWUxmeWUyxmtTMeLSK9+RoiirQzEDgY2P5wvWmkBPZQskcwVSuQLJrEW2YJHMOcci5DaeAGdeu1dX8eoaHpe6WnPE+agS8rnRsPHoasdUGBS7IwEuKtLvd/PEiUEePz7AtYUkr07HuDKfIL+Hx8proWsqvT51y42fC5ZdDHMn1FM5J9BTxa/TeYtcwSaVs0jlLEht/3p3ao8oq/VG3KX6I1pZHZLi1zKzpjNIgItdURWF40MBjg8FSOcKmHNx3phZ4dZSas/NYGkkTXW2VQtus1NRwbJJ5wukcxbpvEW6GOzpnEUmXyCTd34IZPJOjz5bsFjJbHm6NXRVWQ3z1WJSulNEqnTM+XztsW4rKtXuJMBF1bwujQfGe3lgvJdENo85F+fibFzCvEk0VSneBN388VKxJ8u2yeat1SDP5K3VrzOFzR/LWzb5bIEEu6uho0BZuBdDXXMKTLnU0ufO43r516q65rFScSmxPQlwURcBt87D+8M8vD9MKlfg8nyCy/NxrkeSe6q4UDtSy8bjK2HbNrliUahMcagmV+zBl46vP+Z87lTSyxZssoXdh385BdA1pRj6peqCTkVBvVhh0O3SUGx7teKgvu7jZt/jVCbsnh8OEuCi7nwujfvHerh/rIe8ZXNrKcnVhSRXFxIsNngjZlE7RVFwF8fKg56dn1/OsmxylkU274R6zrLJF5xj+YJTLjVXKPbwLecHRX7NY87Xlk3xB4MNufp2ABS4U2JWpfjRCXpVpfixVH72zuNOSVpWS9JqytbnUIuPqQprPxYfU6jPDxEJcNFQuqpwdCDA0YEAP2YMEU3luB5Jcn0xyc3FZFeUuRV3qKqCR9XYZui+Is4PgrXhXl4vO2/ZKKpCJldYPVb++Prnlh+zgbzlVPOu4ZeEmpVCXVkf8rsIdglw0VRhn4uH9vfy0P5ebNvZC/PmYpLb0RS3l9IN2ZxZdB7nB4GCh60LiFW7YUX5ZgqWTfHjFscsm0L5Zgzrvt7qHKsbOqz/aBU3fIDVY47q7hnVFOCGYfwo8FPAHGCbpvm/ajmf2FsU5c6GxY8d6gNgIZFlMppiIppmcjlFJJGTBUSirpwdeRQqvCXQEPa6cF/7deXv96oD3DAMP/DnwJtM08wYhvFlwzDea5rm96o9pxCDAWdLtAfGewGnsNR0LMNULM30cprpWIaVzOYLW4ToFEppPH3TzdoqV0sP/G3ATdM0SzNLfwi8H1gT4C6XpgwNhWp4mfbWzdcG7XF9B8Za3QIh2lMtFeqHgZWyr2PFY0IIIZqglgCfA8q7Zz3FY0IIIZqglgA/BRwyDKM0U/RHgG/U3iQhhBCVUOwaljwbhvFjwM8A80BOZqEIIUTz1BTg2+nmKYaGYYwCHwceME3z0Va3p94MwziGc31ngP1AxDTNj7W2VfVjGIYKfB14AXADx4D/bJrmDjX9OodhGD6c6/uOaZq/1Or21JNhGM8D6eKXBdM039vK9tSbYRgG8CGcGpOPA79lmuaLmz23IQt59sAUw3cAXwMebHE7GqUf+Jxpml8DMAzjvGEY3zBN83SL21VPp0zT/DiAYRhfw+lsfLa1TaqrjwNnW92IBvmWaZq/1epGNIJhGBrwR8AHTNO0DMP4G2DL1UqNWolZ0RTDTmWa5pcMw3h3q9vRKKZpvrTukAokWtGWRjBN08IJOAzD0HF+yzBb2qg6MgzjZ3H+z90PBFvcnEa4zzCMXwF8wEumaXbTvbdHccq1fKTYEY4An9rqybXcxNyOTDHsEoZhPAl82zTNi61uS70ZhvE+4B+AfzBN8+VWt6ceDMO4BzhpmuZXWt2WBvqEaZqfAH4b+HXDMN7V6gbV0SGcDvCnTdP8XeBdwH/Y6smNCnCZYtgFDMN4AngC+K+tbksjmKb5bdM0fwI4YhjGf2l1e+rkSSBtGMav4gz1PWYYxkdb26T6Ko0Hm6ZZAJ7BeY92ixhw0TTN5eLXzwLv3urJjQpwmWLY4QzDeD/wPuAXgVHDMN7W4ibVjWEY9xSvr+Q6cLRV7akn0zR/xzTNj5mm+Xs4//lfNE3zT1rcrLoxDONuwzA+XHboBHC1Ve1pgBeAgeJYODg98ktbPbmRs1C6doqhYRiPAz8H/ATwZ8AfdtkMhjcDTwOlYYUA8EnTND/dskbVUXGWze/jzLJxASeBXzBNc6alDasjwzB+Gvh5nFk2nzRN829b3KS6MAxjDPg/ODdoe3D+/f5b8b5GVygOW74HJzsPAh/ZKl8aFuBCCCEaq1FDKEIIIRpMAlwIITqUBLgQQnQoCXAhhOhQEuBCCNGhJMCFEKJDSYALIUSH+v+uRFi239cyuwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pyplot.fill_between(coordinates, mean-std, mean+std, alpha=0.6)\n", "pyplot.plot(coordinates, mean)\n", "\n", "pyplot.axis([0, 6, 0, 10])\n", "\n", "pyplot.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cholesky decomposition\n", "\n", "The assumption with generalized polynomial chaos expansion is that there\n", "exists a smooth mapping to a stochastic independent variable. However, such a\n", "mapping does not always exists. In those cases making an orthogonal expansion\n", "directly on the dependent variable using Cholesky decomposition. This can be\n", "done as follows:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:19.746111Z", "iopub.status.busy": "2021-05-18T10:57:19.745795Z", "iopub.status.idle": "2021-05-18T10:57:19.835661Z", "shell.execute_reply": "2021-05-18T10:57:19.836266Z" } }, "outputs": [ { "data": { "text/plain": [ "polynomial([1.0, q1-1.0, -0.9*q1+q0-9.1, q1**2-2.0*q1+0.9,\n", " -0.9*q1**2+q0*q1-8.2*q1-q0+9.1])" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "expansion = chaospy.generate_expansion(5, distribution_q, rule=\"cholesky\")\n", "expansion[:5].round(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The method is known to be numerical unstable, so it is important to verify\n", "that the expansion is indeed orthogonal:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:19.839586Z", "iopub.status.busy": "2021-05-18T10:57:19.839118Z", "iopub.status.idle": "2021-05-18T10:57:19.889550Z", "shell.execute_reply": "2021-05-18T10:57:19.890138Z" } }, "outputs": [ { "data": { "text/plain": [ "array([[ 1., -0., 0., -0., 0., -0., 0., 0., -0., 0.],\n", " [-0., 1., 0., -0., -0., -0., -0., -0., -0., -0.],\n", " [ 0., 0., 1., 0., 0., 0., 0., -0., -0., -0.],\n", " [-0., 0., 0., 1., 0., -0., -0., 0., 0., -0.],\n", " [-0., 0., -0., 0., 1., -0., 0., 0., -0., -0.],\n", " [ 0., 0., 0., 0., 0., 1., -0., -0., -0., -0.],\n", " [-0., -0., -0., -0., 0., -0., 1., -0., 0., 0.],\n", " [-0., -0., -0., 0., 0., -0., 0., 1., 0., 0.],\n", " [-0., 0., -0., 0., -0., -0., 0., -0., 1., 0.],\n", " [ 0., -0., 0., -0., 0., -0., 0., -0., 0., 1.]])" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "chaospy.Corr(expansion[-10:], distribution).round(5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This expansion can be used with point collocation method directly:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:19.893720Z", "iopub.status.busy": "2021-05-18T10:57:19.893054Z", "iopub.status.idle": "2021-05-18T10:57:20.115289Z", "shell.execute_reply": "2021-05-18T10:57:20.115023Z" } }, "outputs": [], "source": [ "samples_q = distribution_q.sample(1000, rule=\"sobol\")\n", "evaluations = numpy.array([exponential_model(sample) for sample in samples_q.T])\n", "model_approx = chaospy.fit_regression(expansion, samples_q, evaluations)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:20.117580Z", "iopub.status.busy": "2021-05-18T10:57:20.117306Z", "iopub.status.idle": "2021-05-18T10:57:20.143687Z", "shell.execute_reply": "2021-05-18T10:57:20.143354Z" } }, "outputs": [ { "data": { "text/plain": [ "(array([10. , 9.89956, 9.80022, 9.70198, 9.60482]),\n", " array([1. , 0.98159, 0.9644 , 0.94841, 0.93361]))" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mean = chaospy.E(model_approx, distribution_q)\n", "std = chaospy.Std(model_approx, distribution_q)\n", "\n", "(mean[:5].round(5), std[:5].round(5))" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "execution": { "iopub.execute_input": "2021-05-18T10:57:20.145917Z", "iopub.status.busy": "2021-05-18T10:57:20.145648Z", "iopub.status.idle": "2021-05-18T10:57:20.221982Z", "shell.execute_reply": "2021-05-18T10:57:20.222253Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD3CAYAAAAE2w/rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtcUlEQVR4nO3dZ5Bj2Xne8f8NyEA3OqfJ6c7MRoZd5l2uKFGyWSqbCmXzg2SX+c0qSrQtlULJskxRJbEUXSVKoimrKFqkSYpcikVSJFdiWO5w4+zOzO7OTN/ZydM5otHI6foDgG6gIxq4iP3+qrq6cRvh3Gn002fOfc85imVZCCGEaD9qsxsghBCiOhLgQgjRpiTAhRCiTUmACyFEm5IAF0KINiUBLoQQbUrf7Q6GYQwDHwceMk3zkcIxN/DHwCRwEvhD0zSv1bOhQgghylXSA3838DVAKTn2UeCuaZp/APwZ8H/sb5oQQoid7Brgpml+GVjdcPgDwHOF778GPGQYRpf9zRNCCLGdXYdQtjFIeaiHC8fCG+/4lt/7Z+vHzwyuPzDg4pd/7GSVL9taNE0lm801uxl1I+fX3uT82pfDoSm736v6AJ8DAiW3uwrHNlmKpojGUzi0fGd/YjHDxGwYv6val24dwaCXUCjW7GbUjZxfe5Pza18DA4Hd70T1VSjfBN4BYBjGA8Al0zQ39b4BLGAhmio7dmepM//RhRCikXYNcMMwHgd+ARgxDOO3DcPwAP8LOGwYxm8D/w348E7PMR8pD/BbEuBCCFGzXccxTNN8Gnh6i2/9UqUvMre6IcAX45U+VAghxDbqPpFHUxVWkxkS6ezasdVketOwihBCiL2pe4AP+J3AFsMoi9F6v7QQQnS0ugf4UJcbgLkNAX5zQcbBhRCiFg0IcBcA85Ekpbv/3AvFyeRkNyAhhKhW3QM86HHg1BTi6RzR1Po4eCqb496yXMwUQohq1T3AFUVhwF/shZcPo9xYkHFwIYSoVkOWky1eyJyLJMuOS4ALIUT1Ghrg85FU2Tj4YixFKJZuRBOEEKLjNCTAfU4Nn1MjnbVYjpcH9nXphQshRFUaEuCKojAYyI+Db5yVKQEuhBDVadiWakOFYZTZ1fJx8DtLMVKZzlwSUggh6qlhAd7vd6IosBRLlwV21rK4KbMyhRBizxoW4A5Npc+7dTXKtXkJcCGE2KuG7ko/FCgOo2yuB89ZMitTCCH2osEBXryQWT6tPp7OyqxMIYTYo4YGeJdbx62rJDI5wolM2ffMuUgjmyKEEG2voQFeWk64sRrFnJNxcCGE2IuGBjhsPw6+mkwztZJodHOEEKJt1X852YC77PZgYWGrxViKdLa8/ntchlGEEKJidQ/w08OBsttOXaXH68CyNq9OOD67Wu/mCCFEx6h7gJ8ZCWw6NlwYB58Jl4+Dh+JppsMyjCKEEJWoe4CPBT0EXI6yYyOFXXpmNpQTAlydlWEUIYSoREM2dDg14Cs71uXW8ThUkpncptUJr87IMIoQQlSiIVUopwb9ZbcVRWG4sNnxxmGUlUSaiZBM6hFCiN00JMAP9Xhw6VrZsbVhlA0BDnBFeuFCCLGrhgS4piqc6PeWHev3OdFUhZVEhljJZseQHweXtVGEEGJnDZvIc3KgfBhFUxUGC2uEz2yYlRlNZbi1GGtU04QQoi01LMCP9/vQFKXs2PDaMMrm0sHL0zKMIoQQO2lYgLt0lcO95cMoxXrw+UiKTK58Vua1+Yjs1COEEDto6FoopwbLywndDo0er4OctXmvzFQ2J1PrhRBiB40N8AE/CuXDKCOFXvj0FtUor02FG9IuIYRoRw0NcL9LXxv3LhrpXh8H31h5cnc5TmjDRB8hhBB5DV9OduOknoBLx+/USGUtFqPlwygWlvTChRBiGw0PcGPDtHpFURjtzs/KnFrZPIzy6lR403opQgghmhDg/X4XPYXd6YvWAjyc2BTWK4k0t5dkar0QQmzU8ACHzb3woCe/uFUivXlxK4CLkyuNapoQQrQNvZYHG4bxa8ARYAE4CXzYNM1du8unBv08f2d57baiKIx2ubmxGGNqJUnvhh76tbkIsVQWr1Pb+FRCCLFvVd0DNwxjGPhN4COmaf4PwAf8TCWPHet243WU/+1YHwffPIyStSxenZJeuBBClKqlBx4DUkAXEAL8wOWNd9I0lWDQu/EwDx3u4eW7673w4R4Nl64STWWJ5/LDKqWuLsT5yYc8KBum4zfTdufWKeT82pucX+erOsBN0wwXhlC+aBjGNDABXN94v2w2Ryi0eWGqA34Hz6UyZceGu1zcWYpzdyGKd6i83HA6leHizUWO9rXODywY9G55bp1Czq+9yfm1r4GBzVtRbqWWIZSHgV8DPmCa5n8kPw7+O5U+/mivF4da/vKjhU0eJle23hfzlYlQNU0VQoiOVEsVyhiwZJpmsRs9DbgrfbCuqRzbsEb4gN+JriqEExlWk5lNj3ljPko4ITMzhRACagvwbwOXDcP4E8Mw/jvwCPAHe3mCrdYIHy1MrZ8Mbe6F5yyLCxNyMVMIIaC2MfAs8Eu1vPjJAR+qopStgTLW7eHucoKJUILTG8bBAS5MhHnXsT50tXUuZgohRDM0ZSJPkcehcSDoKTs2GHDi0BRWk5kth0ti6YzsmSmEEDQ5wAFObZiVqSoKY4Wa8IkthlEAzt8N1btZQgjR8logwDcPkxwIrgf4VgtZzawmuLss66MIIfa3pgd40Otg0F++Rni/z7k2qScU31yNAvBiyVR8IYTYj5oe4LB5jXClZBhlu5rwN+ajLMVSW35PCCH2g9YI8A3j4FA6jBLfchjFwuKlO6F6N00IIVpWSwT4cJebbrej7Fiv14HHoRJP51iKbT1559WpMPF0thFNFEKIltMSAQ75mvBSpcMo90JbX7BM53JSkSKE2LdaJsCNwc3VKId68jXiE6EE2dzW26qdvxcinc3VtW1CCNGKWibAD/Z48DjKN2zo9jjoduuksxazq5v3ywSIp7OyY48QYl9qmQBXFWXTMArkgx3Yse77+duhbXvoQgjRqVomwGHrST0HC9UoM6tJkpmth0pWk2lenw7XtW1CCNFqWirAj/V5cWrlTXI7NAb9Tixr+5pwgGdvL5ctiiWEEJ2upQJc11SObbHjzqEKhlGWYymuyiJXQoh9pKUCHDbPygQY6XahqQrLsTSRLTZ6KPrRraUtJ/0IIUQnarkAP9GfXyO8lK6qjBU2etipF74QTXF1NlLX9gkhRKtouQB3OzSO9O48jLJTL/vczUXphQsh9oWWC3DYelJPv8+J16kRT+eYi2y/iNVCNMUV6YULIfaBlgzwU4M+FMqHURRF4XChF35naee1wJ+5sSgVKUKIjteSAe5z6msTeEoVA3wqnNi2JhxgKZbi9WmpSBFCdLaWDHAAY3DzrEyPU2Mo4MKydr6YCfleeEZmZwohOlgLB7h/0zAKwJHe9WGUnS5WriTSXJI1UoQQHaxlA7zL7WC0sJxsqeEuFy5dZTWZ2Xad8KJzN5dI7TDUIoQQ7axlAxy2HkZRFWWtpHC3i5nRVIaXZL1wIUSHaukAPz0U2PJ48WLmxEpi17XAn7u9RCwlu/YIITpPSwd40ONgpGvzMErArdPvc5LNWdzb5WJmKpvj3M3FejVRCCGapqUDHOD00OZJPcDaolc3F3e+mAnwysSK7GAvhOg4LR/gZwa3HkYZ6V6/mLkQ3Tmcc5bF999YqEfzhBCiaVo+wINeB8OBzcMoqqJwtFBSeGtx52EUAHMusmvtuBBCtJOWD3CAM8NbD6Mc6fOiAFMrCeLp3S9UfvfavCx0JYToGO0R4NsMo3gcGiPdLix2LykEmA4neHVKtl4TQnSGtgjwoHfrahRYv5h5azFW0QJWP7i+uOM6KkII0S7aIsABzm5TE97vc+J3aSQyOaZXkrs+TzSV4UdSViiE6ABtE+BnhrdeG0VRlLVe+PWFaEXP9dLdEIu7VK4IIUSra5sA73I7GNtibRTIz8x0aApLsXRF9d5Zy+KfzXm7myiEEA2l1/JgwzAM4ENAHHgc+F3TNF+0o2FbOTscYGJl88VKXVM50uvljfkoN+Zj9B527vpcNxejmLMRjG0mCgkhRKurugduGIYG/CnwMdM0PwF8GLhlV8O2cnrIv2nD46JjhZLCyZVExWufPGXOy2qFQoi2VcsQyiOAAnzEMIzfBH4aqOt0R79LX1uJcCOvU2Os240F3FyMVfR8q8m0rJMihGhbtQyhHAbeAXzINM0VwzD+HkgBnym9k6apBIObd5mv1ttODDB1cXLL750e6WJiJcHtpRgPHOjGoe3+9+niTIR3nh5ieJsyxZ3YfW6tRs6vvcn5db5aAjwMjJumWdz25hzwXjYEeDabIxSqrEdciTGvTjadJbtFzXfAqdLrdbAUS3N9dpXj/ZvXE9/KF5+/wy8+cgBlm+GZ7QSDXlvPrdXI+bU3Ob/2NTCwddn0RrUMobwA9BXGwiHfI79Ww/NVxO3QODGwfTAXv3d9obKJPQCTK3FemZDt14QQ7aXqADdNcwn4deDPDcP4HWAA+DO7GraT+4e7tv3eaJcLv1MjlsoyGUpU/Jzff2OB1UTGjuYJIURD1FRGaJrmV4Gv2tSWih0f8OHWNRKZzdUmiqJwctDHhYkw1+ajHAi6KxoaSWVzfHt8jp9/eLQeTRZCCNu1zUSeUrqqbLvRA8ChoAePQyWcyDAT3n16fdEb8xGuzKza0UQhhKi7tgxwgPtHth9GUVWFE4ULmOZ8dE9LyD41Pi97aAoh2kLbBvjBoJugx7Ht94/0eXBqCsux9K479pSKpTM8NT5nRxOFEKKu2jbAFUXhvh0uZuqqulZGeG2uskWuiq7MrjI+K0MpQojW1rYBDvDA6M61ksf6vOiqwlwkxXIsvafn/vbVeaIpqUoRQrSutg7wXq+TA91bT60HcOoqRwtLzZpzkT09dyyd4VtXZChFCNG62jrAAR4Y3X4YBeBEvxdNgelwcs+98GvzEdmCTQjRsto+wM8OB9DV7U/D7dDWeuFXZ/fWCwd4anyOUHxvwS+EEI3Q9gHu0lVOD+68pvepQT+aqjC7mmRpjzvxpLI5vv76TMXT8oUQolHaPsABHtxlGMWlqxzvr74Xfi8U57lby1W1TQgh6qUjAvxwr2fHmnCAk/2+tYqUvdSFFz1zc5GJ0ObdgIQQolk6IsAVRdm1F+7UVU4Ue+Eze++F5yyLr702QyItszSFEK2hIwIc8sMoW+1aX+r4gA+HprAQTTEfqXyNlKKVRJp/ujJbbROFEMJWHRPgXW4Hx/p23p3DqamcLMzOvDwd2dMaKUXjcxFevheqpolCCGGrjglwgDcd6N71PscHvLh0leV4msmVytcLL/Uv5jwz4eoeK4QQdumoAD8x4CPg2nmJc11VOVNYivbyTIRsbu+98Kxl8eSr0zIeLoRoqo4KcLWCi5mQr1oJuHRiqSy3KtzBfqNQPM1XXpmsahhGCCHs0FEBDvDwWPeuFzNVReG+kXwvfHwuQiqTq+q1rs6Eefa21IcLIZqj4wK82+OoaDf64YCLfp+TdNbi2vzelpst9cyNxap78UIIUYuOC3CANx/c/WKmoijcP5JfjvbGQrTqXXhylsU/vjZNaI8LZQkhRK06MsCP93l3nZkJ0ON1cCDoJmfB69PVb+AQT2f58qUp0tnqhmKEEKIaHRngiqLw5gpKCgHuGw6gKTC5kqhqck/RXCTJ1y/PykVNIUTDdGSAAzw01r3jMrNFXqfGqcJqhpcmV2tadXB8dpVzN5eqfrwQQuxFxwa4x6Fx3/DOW64VnRzw4XNqrCYz3KzxguS5m0tclf00hRAN0LEBDvDWQ8GK7qepCg8ULmiOz0RIZqqfoGNh8fXXZ5mqcpanEEJUqqMDfCjg4mBw+z0zSw13uRgKOEnnLC5P7321wlKZXI5/uDglO/kIIeqqowMc4NFDPRXdT1EUHhjtQlHgznJ8zzv3bBRNZfiHC1My3V4IUTcdH+CnBn0VlRQCBFz62mqFFybDNW+jNh9N8pVL02SqWG9FCCF20/EBrigKbz0YrPj+xpAfn1MjnMjwRg0zNIvuLMf4xuszUl4ohLBdxwc45EsKXbpW0X11VeHhsfyCWOOzESLJTM2vf2V2le9eW6j5eYQQotS+CHCXrvKmsd1XKSwaDLg42JOfoXlhImxL7/nFu8s8f1tqxIUQ9tkXAQ7wyKEeVGXnVQpLPTDShbOw/drdZXs2M/7eGwtcmlyx5bmEEGLfBHjAra8tXlUJl66urS3+2vSqbdUk37o6x7hM9BFC2GDfBDjA2w/37LpWeKkDQTeD/vySsxcm7RlKKe5uf32h9gukQoj9bV8FeL/fxYmB3dcKL1IUhTcd6EZXFWbCSe4u2zO7MmtZPHlpmttLso64EKJ6+yrAAd55tLKJPUVep7Y2lPLqVLjqdcM3yuRyfPniFPdsGl8XQuw/NQe4YRgewzBeNQzjj+1oUL2NdXs40uvd02MO9bgZ6XKRyVm8MrFiW013KpvjixcmmQhJiAsh9s6OHvjHgQs2PE/DvPNI757uryj52nCnpjAfSdm6hVoqm+MLr0iICyH2TqmlN2kYxi8AUeBBwG+a5q9uvE8uZ1nZFtyp5n8/c5O7exyDvrcc49kbS2iqwvvPDBL0Omuebl/k1FV+8e2HOdJX+Rh9vWmaSiv+7Owi59feOvn8HA6tomoLvdoXMAzjLHDGNM3fMgzjwe3ul83mCIVa72LdW0YCXJ8J7+kxQz4nB4Nu7oUSPHtjkZ84O0S2hqVnS6VS8Okf3ODnHx7lSN/ehnjqJRj0tuTPzi5yfu2tk89vYKCykudahlA+CCQMw/gN4N3Ao4ZhfLSG52uo4/0+xrorW2q21ENjXficGiuJDBcnQra2KZ3L8aWLU1y3YQ0WIUTnqzrATdP8fdM0P2aa5h8C54AXTdP8c9ta1gDvPra3sXAAh6byyKEgigLX56K2b9yQyeX48qUprszIZB8hxM7sqEL5WeAx4O2GYXyo9iY1zvF+Hweq6IX3eB3cX9iu7ZWJFdtKC4uKk31esbmHL4ToLFWPgReZpvkV4Cs2tKUpHjvRx+dfntjz4473e1mIpZleSfDS3RDvOd67p7VWdmNh8e2rc8RTOd5Vxf8UhBCdb99N5NnoSK93z3XhkC8tfPRID25dZSmW5vXp+gx5PH1jge9cnZP1xIUQm+z7AAd474n+qh7ndmg8ejiIAtxYiNWtlvvliRBPvjpNpkNLpoQQ1ZEAB0a73RiD/qoe2+dz8uBoYTz8XpiVOm1kbM5F+NzLk0RTtW8wIYToDBLgBY+f6K96DPton5eDQTdZy+KFOyFSdeopT67E+bsX77FQ44bLQojOIAFe0O9zri1atVeKovDwgW663TrRVJbzd+1bL2WjUDzNZ1+8Z+t0fiFEe5IAL/HY8T6cWnX/JLqq8LYjQRyawuxqsm4XNQESmSxfvDDJy/dCdXsNIUTrkwAv4XfpvO3w3pabLeVz6rytcFHz+kKM23XsJecsi++Mz/FPV2bJ5qRCRYj9SAJ8g7cf6SHgqr48fsDvWtvV/uJkmPlI0q6mbeni5AqfOz9BJCkXN4XYbyTAN3BoatVlhUVH+ryc6PdiAS/cCbFa53CdWInzt8/flSVphdhnJMC3cP9IgNEud83PMdLlIp21eO7Wsm2bIm8nksrw9+cnePHOcl1fRwjROiTAt6AoCu8/PbinDZC3eo63Huwm6MlXpjx3e5l0nSfi5CyLf7k2z5OvTpPMyKQfITqdBPg2RrvdVZcVFumayjuO9OBzaoTiGV64E2rIBcfx2VX+9oW7zITtXSlRCNFaJMB38N6Tfbh1rabncDs03nm0B5euMh9J8fK9+tWIl1qOpfi7F+/JkIoQHUwCfAc+p87jJ/pqfh6/S+edR3vQVYXJlQSvTq02JMSzhSGVL12QKfhCdCIJ8F28+UA3IzVe0AQIehy8/UgQVYGbizFen25MiANcX4jyN8/dlZ1+hOgwEuC7UBSFf3Vm0Ja1vgf8rvzqhUp+os+VmUjDQjyayvCli5N86+osKbnAKURHkACvwHCXm0cPBW15rpHCcynAtfko47MRW563UhcmVvib5+9wd1lqxoVodxLgFXrP8T6CHoctzzXa7eath7oBGJ+LYjY4xEPxNJ87P8FT43N1L20UQtSPBHiFHJrKvz4zZNvzHQh6eOvBfIhfmY1wdbZxwymQ37Lt/L0Qn37uTl3XbBFC1I8E+B4c6fPypgPdtj3fwR4PbymE+PhspKEXNotC8TSff2WCb16eJV7n2aJCCHtJgO/Rj50coNttz1AKwKEeT35MvHBh8+JkuCn7X16aWuFTP7rDa1Phhr+2EKI6EuB75NJVPnB2qKZp9huNBd28/XAPqgK3l+Kcv7tCrgkhHktn+PrlGT53foKFOq+iKISonQR4FY70eXmrTVUpRcNdLt51tBddVZhYSfD87VDTNjG+sxzjb56/y7cvz8iaKkK0MAnwKj1xoo/BgMvW5+z3O3n3sV6chV19nrm5VPdVDLeTsyzOXV/gUz+6zatTzRnWEULsTAK8Srqm8vNvOYhmwwSfUj1eB4+f6FtbAOsH1xcJJ5o3DT6SyvCNyzN85sV7st64EC1GArwGI91u3nuyts0ftuJ35ddg6fE6iKdz/PD6Yt139tnNdDjBZ1+6x5OXpliOpZraFiFEngR4jR49FOREv8/253XpKu851pvfFCJn8aNbyy2xE/34XIRPPXuHp8bnZIEsIZpMArxGiqLw0/cPE3DZV1pYpKkKbzsc5Hi/F8vK77F5cTJMrsmbGOes/CSgvzp3mx/eWJQLnUI0iQS4DTwOjQ8+OGzLglcbKYrCg6NdvPlAF6oCtxZjnLu1RDLT/Ek3qWyOczcX+ctzt3ju9pIskiVEg0mA2+RA0MMTdRgPLzrc6+U9x3tx6yqL0TTff2OR5Vi6bq+3F/F0lu+/scBfnrvN87eXZH0VIRpEAtxGbzvcw5mhQN2ev9fr5ImTffQWL27eWOTGQrRlSvxi6Qzfe2OBTz5zm2dvLcnQihB1JgFusw+cHWLAZ299eCm3Q+Pdx3o52uclZ8GrU6u8eCdEqoV6vbF0hh9cX+CTz9zi6esLxFLNH+4RohNJgNvMqav83MMjeBy17aW5E01VeHisi0cPBdFVhalwsqWGVIoSmSw/urXEXzxzi6fG5wi1WPuEaHcS4HXQ43XywQdG6nJRs9RY0M0TJ/sIenRiqSxP31jEnIs0ZR2VnWRyOc7fC/HXz97myUtTMiFICJtIgNfJkT4v7z89UPfX8bt0Hjvex/G+fKnhlZkIz9xYIpJsvRrtnGUxPhfhsy/d4zMv3OXyzCrZJpdECtHO9GofaBjGceDjwCvAAWDRNM2P2dWwTvDmA0GWomlevLtc19fRVIUHx7oY6nLxysQKS7E037u2yP2jAY72elDq/D+BakyFE3zttWm+69R5+EA3bxrrJuCu+u0oxL5USw+8F/iCaZp/ZJrmrwD/3jCMt9jUro7xvlP9GIP+hrzWUMDF+072czDoJmtZXJoM8+yt5Za+iBhJZTh3c5FPnrvFk5emuL0Ya5mqGiFaXdVdHtM0X9pwSAWitTWn8yiKwr+5f5jPvzzJxEr9x36duspbDwUZ6UpwYXKFuUiKfzEXODvs51i/t+7j8tUqDq+Mz0Xo8Tp5aLSLB0e78LukVy7EdhQ7ejuGYXwQeG+hJ14ml7OsbAuVuNlJ01QqPbd4Ksunz91kbrVxi1LF01ku3A1xr7ADfY/XwSNHeujxOit6vKooTb0gqioKxlCANx8KcmoogKba+8dnLz+/diTn174cDq2iN3vNAW4YxhPAB4GPmqa56V8znc5aoVDzF2Gqh2DQy17ObTWR4bMv3WMl0dhyuplwgouTYeLp/I/nRL+X00N+HNrOI2hOp06qRRas8jp0zg77eWC0i5Euty3PudefX7uR82tfAwOB+ge4YRgfAN4D/CYwAhw2TfO50vtIgJdbjqX4vy9NEGlwMGayOa7MRrixkG+vS1e5b9jPoZ7tL3K2UoCX6vc5uW+4i/uGAwS91S8i1skBAHJ+7azuAV64YPk0cL5wyAd80jTNz5TeTwJ8s4VIkr8/P0ks3fhwXI6luTQVXpv00+Nx8OBogF7f5mGVVg3wUqNdbs4MBzgz5Kdrj5tNd3IAgJxfO2tID7wSEuBbm11N8vmXJ4g3Ycs0y7K4F0pweXqVRGG9koNBN2eHA3id6zNI2yHAixQUxrrdnB7yYwz66fbsHuadHAAg59fOJMAboNY3UDNDHPLDKuZ8lOvzUXIWqAoc7fNiDPpx6WpbBfhGI11uTg36MQZ89Pu3XpumkwMA5PzamQR4A9jxBpqPJPlck4ZTiqLJDFdnI9wLJQDQVYUTAz7OjnZhdcBV/h6vk1MDPk70+zjY41krpezkAAA5v3YmAd4Adr2BFqMpPv/yJKvJ5i72tBJPc3kmwmyh1NGlq5wc8HG0z4OudsaqC25d42iflxP9Pt50vJ9MvHP39+zkgIPOPj8J8Aaw8w0Uiqf5wiuTLLXAhsEL0RSXp1dZKlzodGoKJ/p9HOv37lp62E5cLge9Lo1jfV6O9vkY63bbXmveTJ0ccNDZ5ycB3gB2v4GiqQxfujDFdDhh23NWy7IsFuMZXi+pWHGoCsf6vRzv9+HS2z/IN47xOzWVQz1ejvZ5ONLrZWCbsfN20ckBB519fhLgDVCPN1Aqk+Orr05zY7H5qxI4nTrJZJqFaIrx2SgL0fz/DjRF4XCvh+P93rae6r7bRVqfU+dQj4fDvR4O93jp26LUspV1csBBZ5+fBHgD1OsNlLMsvjM+x4WJFdufey82BtxiNIU5F10bIwcYDrg4MeCl3+dsyVUPd7LXKhufU+dgj4dDQQ8Hgm6GAq6WPudODjjo7POTAG+Aer+BXrizzPeuLWDRnPVItgu4lXiaGwsx7oXiFJfz7nbrHO/3Mhb0oLfJOHKtZZJOTWUs6OFAt5sDQQ+j3e6WGlrq5ICDzj4/CfAGaMQb6MZClH98bYZkpvG14rsFXDKT5eZinFuLsbUNjB2qwsEeD0f7PHueGdlodte5Kyj0+5yMBd2MdrkZ6XYz4Hc2bQXITg446OzzkwBvgEa9gZZiKb58cWptDLpRKg24bM5iIhTn1lK8bF/OXq+DI71exoLuluyVN2KikkNVGe5yMdzlZqTwuc/raMjQSycHHHT2+UmAN0Aj30CpTI5vXpnl6uxqQ14Pqgu4lXiaW0tx7i3HyRTGV3Q1P839UI+HPl9jwqsSzZpp6tRUBv0uhrtcDAXyH/1+l+1/5Do54KCzz08CvAGa8QY6fzfEd6/Nk23AOt21BFwml2MilOD2Ypzl+Hqv3ONQORj0cLDHQ1eTt1BrpaUCVEWhz+dkyO9iIOBkwO9i0O+saRiqkwMOOvv8JMAboFlvoJlwgq++NsNynSf92BVw4USGe6F8r7y4JjlA0KMz1u1mtNvdlHLEVgrw7bh1jX6fk36/kwG/kwGfiz6fs6L9Qzs54KCzz08CvAGa+QZKZXL8sznPpan6lRraHXCWZbEYTXN3Oc7kSmJtiAXyVSyjhTBvVM+8HQJ8O65CsPf5HPT5nPR6nfR5HfR4nWuzSTs54KCzz08CvAFa4Q1kzkX41pW5uiyGVc+Ay+YsZleTTK0kmA4ny8I84NIZ7XYxHHDRU8cLfu0c4NtRFYVut06v18mBAT8uy6LH66DX66Db42jZPVGr0Qq/f/UiAd4ArfIGiqYyfGd8nnGbL3A2KuCyOYv5SJLJlSTT4QTp7Pp70qkpDAXyF/wGAy6cNq7F0okBXmrj+RXDPeh1EvTo9HgcBEs+3A5th2drPa3y+1cPEuAN0GpvIHM2wnfG52zbrq0ZAZezLOYjKWbCSWZWk8RS6/XvCtDrczAccDHgdxH06DX1zvdbgO/GpWt0u3W6C4He5dYJuh10eXS63Do+Z2stm9Bqv392kgBvgFZ8AyXSWb5/fYGLE+GaZ3A2O+Asy2I1mWV2NclMOMliNFV2Rg4tP3Fm0J+v2vC7tD0FerPPr97sPj9dVQm4dbrdOgFXPtS73A4ChdsBl162o1O9teLvn10kwBugld9Akytxnhqfr2llw1YLuFQ2x9xqkrnVFPORFLENOxm5HSoDPif9Pie9PieBXQK91c7Pbs04P11V8bs0Ai4dfyHUi7d9hWN+p2bLcE0r//7VSgK8AVr9DWRZFhcnwzx9fbGqi5ytHnDRZIa5SD7M5yMpUht2D3JqCr2+fHVGn89J0OMoW++71c+vVq18fsWg9zn1tc8+p4av9Gunhtepb7u+TKv//tVCArwB2uUNlEhnefbWEi/dDe1pAlArB8BGlmURTmSYj6ZYiqZZjKbWNmwuUhUIehz0eB30eBwMdrtxKrTMzFC7tdPPbye6quItBrpDW/t6oMdLLpXB69TwFI57HRouXW37n6kEeAO0S4AXheJpfnhjkcvTqxWNj7dzAFiWRSydZbEQ5ovRNKvJzefi0JR8qHscBAvB7nG0fwBAe//8KrHd+amKglvX8DpVPA5t/cOp4XGoePT8bbdDxe0oHHNoLbXblAR4A7RbgBfNR5L88MYi5lxkx/t1WgCkMjmW42mWY/mPUCJNIr1502anptDldtBdqL7oducrMtptu7VO+/ltZPf5aYqyFuhuvRDwhc+utdv5Yy5dXfu+y6Hi0lVba+wlwBugXQO8aCac4NzNJd6Yj27ZI98PAbASTebDvBjs8XRZHXopv0uj2+3IV2G4CxfkXFrLTo7ZDz+/Vjo/p5YP8vyHhltXcRZuuwvHnLqCU1v/I+DU1x/j1PK3QQK8Ido9wItmV5M8d2uJ8bkIuZL3Q6v9gthtq/OzLItEOsdKIsNKIk04nmElkSGSzGw56KQAvrKqC20t3O2cdFSN/fjza3cKCg5N4RP/7k0VBXhrVeaLphgKuPi3D44QiqV58e4yr06FN1V07BeKouTHSp0aw13rmxpncxaryQwr8QzhRH48fTWZJZbKEknmPyBZ9lwuXcW/obLC69TwuzScWmeMswt7WViktvkf4FYkwMWaoNfB+08P8tjxPi5OrnB5PsZsh/VwqqWpytqUc/CsHc/mLCLJDKvJDJFkNh/shR57MpMjmcmxWLLJRZGuKmvVFPkPHY9TxVu44ObQFAl4sSsJcLGJ26Hx9iO9vP/BMc5fn+eViRVuLcaatjdnK9NUhW5PfqGoUpZlEU/niKYyRJNZoqniR4ZoKks6my97DCe2/gOpqQreQnWEp1AeV6yW8BYmwrTiLkeisSTAxbZUVeHUoJ9Tg35C8TSXJld4dWqV1eTmHqUopyj5HrbXqTHg3/z9VCZXFujRVJZ4Kks8nSOezpLJ5ZcRWE1uvxeqripllRJuPV8WV7xoFvBaaJYlvfkOJgEuKhL0OHj8RD+PHe/j1lKM16bCXJuLks7tz7HyWjkLFQg93s077liWRTprEU9niaULoZ7Krt9OZYlncmRyVsn4+/ZUhbVyt9JKifzXynolRKEKQlcl8NuFBLjYE0VRONbn41ifj1Qmx7X5CJenV7m1FCurYBHVUxQlX26mq5uGZoosK3+xK5nJkkjnSGRyJNLZwucciUyWZMYiUejNxwrhXwlVyf+BcRXC3qGpODUFRyH0HZqyfkxbPybB33gS4KJqTl3l/pEu7h/pIpHOYs5FGJ+NcHsp1pA9O/czRVFw6QouXaXLvfV9imV2mWw+4IsXVVPZwuctb1tkC6WUW01y2rFNUBLu60FfDHeHlu/d65qCQ1XR144r6IXbDvkjsCcS4MIWbofGQ2PdPDTWTTKT48ZClDfmI9xYiJHIVNbzE/Whayp+TcXv2v2+AJmcVRbuqUyOdDZHKmuVfS7/2so/LmuRymaJUv3PXFOUtTDXtfXg11Rl7bOmKrgcGlgWulL+vdL7lB7rxD8MEuDCdi5d5exwgLPDAXKWxUQowY2FKDcXo8ytpqSapcXpqoJeuAC7FzmrEPCZ8qDP5KxCwG/4OmuRzllkCreLfwSylkU2Y22oqq+dqlAW6KqyHvT5r/N/PFRVQVOKx/MVQevHKbl/8fjWj1OV/Od6/uGQABd1pSoKh3o8HOrx8MTJfiLJDLeXYoWPOOGEVLR0ClVRcOkarhpSxbIssrnSYF//I5DJ5b+XLXyNopBMZ9dur322rM3HchY5i/wkmT1MlLGDUhLm5Z8LX6ubj1VKAlw0lN+lr42bA4Riae4sx7gXinN3OU4oLoG+nymF4RNdA7a+frtmL1PpLcsia1EI9BzZHIVQL/xRsCxyOdbCP2eth3625I9C/jiF+29+XPGzZeWPWRZrX+cHlez941FTgBuG8ePAzwBzgGWa5v+0pVVi3wh6HQS9+bFzgEgyw0QozkQowdRKgpnVJBkpVRQ1UhQFvTCE4qJxa9RY1nqY56ziH5KSr3P5oaf1++S/V6mqA9wwDC/w18B9pmkmDcP4imEY7zNN87vVPqcQfpfO6aEAp4cCwPqO9VPhRH6j43CS+UhSqlxEW1AUJT+EQn3GwWvpgb8DuGOaZvFaw4+ADwBlAe5waMrAQKCGl2ltnXxu0BrnNzwEDzS7EUK0oFr+LzEIrJbcDheOCSGEaIBaAnwOKO2edRWOCSGEaIBaAvw54LBhGMXpAe8Cvll7k4QQQlSiph15DMP4CeDngHkgLVUoQgjROHXbUq2TSwwNwxgGPg48ZJrmI81uj90MwzhO/vxeAQ4Ai6Zpfqy5rbKPYRgq8HXgBcAJHAf+k2ma8aY2zEaGYXjIn99Tpmn+arPbYyfDMJ4HEoWbWdM039fM9tjNMAwD+BAQBx4Hftc0zRe3um9dJvLsgxLDdwNfAx5ucjvqpRf4gmmaXwMwDOOKYRjfNE3z5Sa3y07Pmab5cQDDML5GvrPxueY2yVYfBy40uxF18m3TNH+32Y2oB8MwNOBPgZ82TTNnGMZngW1nK9VrJmZFJYbtyjTNLxuG8d5mt6NeTNN8acMhFYg2oy31YJpmjnzAYRiGTv5/GWZTG2UjwzB+gfzv3IPAFttJtL0HDMP4dfJ7271kmmYnXXt7hPzCjh8pdIQXgU9vd+d6TUmSEsMOYRjGB4HvmKY53uy22M0wjJ8EvgF8wzTN881ujx0MwzgLnDFN88lmt6WOPmGa5ieA3wN+yzCMx5rdIBsdJt8B/oxpmn8APAb8h+3uXK8AlxLDDmAYxhPAE8B/aXZb6sE0ze+YpvlTwFHDMP5zs9tjkw8CCcMwfoP8UN+jhmF8tLlNsldxPNg0zSzwDPn3aKcIA+Omaa4Ubp8D3rvdnesV4FJi2OYMw/gA8JPArwDDhmG8o8lNso1hGGcL51d0CzjWrPbYyTTN3zdN82Omaf4h+V/+F03T/PMmN8s2hmGcNgzjwyWHTgI3mtWeOngB6CuMhUO+R35tuzvXswqlY0sMDcN4HPhF4KeAvwL+pMMqGN4CPA0UhxV8wCdN0/xM0xplo0KVzR+Rr7JxAGeAXzZNc6apDbORYRg/C/wS+SqbT5qm+f+a3CRbGIYxCvwF+Qu0XeR/fv+1cF2jIxSGLX+MfHYeAj6yXb7ULcCFEELUV+PWVRRCCGErCXAhhGhTEuBCCNGmJMCFEKJNSYALIUSbkgAXQog2JQEuhBBt6v8DmnibTvoSVIMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pyplot.fill_between(coordinates, mean-std, mean+std, alpha=0.6)\n", "pyplot.plot(coordinates, mean)\n", "\n", "pyplot.axis([0, 6, 0, 10])\n", "\n", "pyplot.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For more details on this methodology, you can read the journal article\n", "[Multivariate Polynomial Chaos Expansions with Dependent\n", "Variables](https://doi.org/10.1137/15M1020447)." ] } ], "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 }