{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lagrange Polynomials\n",
    "\n",
    "This tutorial uses the same example as the [problem formulation](./problem_formulation.ipynb).\n",
    "\n",
    "[Lagrange polynomials](https://en.wikipedia.org/wiki/Lagrange_polynomial) are not a method for creating orthogonal polynomials.\n",
    "Instead it is an interpolation method for creating an polynomial expansion that has the property that each polynomial interpolates exactly one point in space with the value 1 and has the value 0 for all other interpolation values.\n",
    "\n",
    "To summarize, we need to do the following:\n",
    "\n",
    "* Generate $Q_1, ..., Q_N = (\\alpha_1, \\beta_1), ..., (\\alpha_N, \\beta_N)$, using (pseudo-)random samples or otherwise.\n",
    "* Construct a Lagrange polynomial expansion $\\Phi_1, ..., \\Phi_M$ from the samples $Q_1, ..., Q_N$.\n",
    "* Evaluate model predictor $U_1, ..., U_N$ for each sample.\n",
    "* Analyze model approximation $u(t; \\alpha, \\beta) = \\sum_m U_m(t) \\Phi_n(\\alpha, \\beta)$ instead of actual model solver. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluation samples\n",
    "\n",
    "In chaospy, creating low-discrepancy sequences can be done using the `distribution.sample` method.\n",
    "Creating quadrature points can be done using `chaospy.generate_quadrature`.\n",
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:03.285670Z",
     "iopub.status.busy": "2021-05-18T10:57:03.285301Z",
     "iopub.status.idle": "2021-05-18T10:57:05.210102Z",
     "shell.execute_reply": "2021-05-18T10:57:05.209761Z"
    }
   },
   "outputs": [],
   "source": [
    "from problem_formulation import joint\n",
    "\n",
    "samples = joint.sample(5, rule=\"halton\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.212399Z",
     "iopub.status.busy": "2021-05-18T10:57:05.212130Z",
     "iopub.status.idle": "2021-05-18T10:57:05.377266Z",
     "shell.execute_reply": "2021-05-18T10:57:05.377525Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEICAYAAABxiqLiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcMUlEQVR4nO3df5xcdX3v8dfsZjfsssldNizFK5IQCp+2SojYcNlixRQwaOTRh1hNuAWUIvJD7EaL/P4Vyi9BwaRqk2JvU01rsM31oZdcim0e1FthwUAMeNV8uJAEQamBLCHZZNmdnZ37xznLmYy7mTnHmTMzu+/n48Ej+z1zvjOf/ZLdd875nvM9mXw+j4iISFxNtS5AREQakwJEREQSUYCIiEgiChAREUlEASIiIokoQEREJJFpaX2QmZ0BnAPsBPLuvrzo9auBI4GXgd8HbnL3reFr5wHvBHLA8+6+Oq26RURkfKkcgZhZO7AK+Iy73wLMM7PTi3brAD7r7ncD64F7wr5HAVcCV7r7VcAnzOy4NOoWEZGJpXUE0gO84O5DYftRYDGwcWwHd7+xYP8mYCD8ehHwlLuP3fHYB7wf+H9VrVhERA4qrQA5Athb0N4Tbvs1ZtYKfAz4VJy+o6Oj+VxOd9UDNDc3kcuN1rqMuqCxiGgsIhqLSEtLcyZp37QCZCcwo6A9M9x2gDA8/hq43t2fL+j720V9nyvum8vl2b17f8UKbmSdne0ai5DGIqKxiGgsIt3dM0rvNIG0rsLqA2ab2fSwfSqwwcy6zGwmvDlPshq4192fMrMPh/s+DLzLzMZSsgd4KKW6RURkAqkcgbj7fjO7DFhpZq8Az7j7RjO7G+gH7gLWAu8AjjEzgEOB9e7+kpl9AbjPzHLA19xd8x8iIjWWmSyr8WazubwOSQM6PI9oLCIai4jGItLdPSPxHIhuJBQRkUQUICIikogCREREElGAiIhIIgoQERFJRAEiIiKJKEBERCQRBYiIiCSiABERkUQUICIikogCREREElGAiIhIIgoQERFJRAEiIiKJKEBERCQRBYiIiCSS1jPRpYpaW5vJtEyjo62FgcFsrcsRkSlCAdLgWlubGco00bt2M5t29LNgThcrz51Pa2szw8O5WpcnIpOYTmE1uEzLNHrXbaFv2y5GRvP0bdvFn39zC5kW/dtARKpLAdLgOtpa2LSj/4Btm3b009HWUqOKRGSqUIA0uIHBLAvmdB2wbcGcLs2FiEjVKUAaXD47woql8+mZO4tpTRl65s5i5bnzyWdHal2aiExyOlHe4IaHc0xvhdXnnfTmVVhtLc3s3r2/1qWJyCSnAJkEhodzMJxjaN8QAG2d7TWuSESmAp3CEhGRRBQgIiKSiAJEREQSUYCIiEgiChAREUlEASIiIokoQEREJBEFiIiIJKIAERGRRBQgIiKSiAJEREQSUYCIiEgiqS2maGZnAOcAO4G8uy8fZ58lwB1Ar7s/WLD9c8Ac4FXgOOAidx9Mo24RERlfKkcgZtYOrAI+4+63APPM7PSifY4hCJcXi7YfCVwLfNrdbwYOJQgiERGpobSOQHqAF9x9KGw/CiwGNo7t4O7bge1mdnNR3/3AMDAT2A10AD8p/oDm5iY6tYw5oLEopLGIaCwiGovKSCtAjgD2FrT3hNtKcvc94SmsB8zsZeAl4Lni/XK5UT1EKdTZ2a6xCGksIhqLiMYi0t09I3HftCbRdwKFVc4Mt5VkZvOBzwGL3f3jBPMgN1W4PpGG0drazPRDpzPr8A6mHzqd1tbmWpckU1RaAdIHzDaz6WH7VGCDmXWZ2cwSfd8K9Lv72EO+XwYOqVKdInWttbWZoUwTl6zdzPHXP8QlazczlGlSiEhNpBIg7r4fuAxYaWa3Ac+4+0bgGuByADPLmNkNwGxgiZktCrv/C/ATM/uimd0ILADuTKNukXqTaZlG77ot9G3bxchonr5tu+hdt4VMi55OLenL5PP5WtdQEdlsLq9zmgGd341MtrGYdXgHx1//ECOj0c/ttKYMz97+fna9OnDQvpNtLH4TGotId/eMTNK+upFQpIEMDGZZMKfrgG0L5nQxMJitUUUylSlARBpIPjvCiqXz6Zk7i2lNGXrmzmLF0vnksyOlO4tUmE6cijSQ4eEc01th9Xkn0dHWwsBglnx2hOHhXK1LkylIASLSYIaHczCcY2jfUOmdRaoodoCY2W8BbcAv3X248iWJiEgjKCtAzKwJWA5cRDBvMgIcamb/TrDw4c+rVqGIiNSlcifR7wQ2A3Pd/Uh3P8rdDyMIlb80s85qFSgiIvWp5BFIePTxlfGOMtx9i5l9EphFsNChiIhMESUDxN1HgXFPUZlZt7u/Avyy0oWJiEh9izWJbmYdwJlECyOeDXyk0kWJiEj9i3sV1oPAFuC1sN018a4iIjKZxQ2Q59x92VjDzI6rbDkiItIo4gbIw2Z2IfB82D4fuLiyJYmISCOIGyAXAkNEV1ydUNFqRESkYcQNkFfd/YKxhpm9s8L1iIhIg4gbIE+b2UKiU1hnAz+qbEkiItII4gbIlcDWgvbRwK2VK0dERBpF3AC51t3XjDXM7IzKliMiIo0i7gOlLjWzk8Ya7v5vFa5HREQaRNwAedbdN481zGxWhesREZEGEftGQjM7C/gZkAeuAK6qeFUiIlL34gbIZcDCgvbRKEBERKYkTaKLiEgisQLE3deEoXEisEWT6CIiU1esSXQzuxH4LDAbuDJsi4jIFBT3FFaru39grGFmd1a4HhERaRBxL+MdLdEWEZEpoqwjEDM7DfgPYMTMvgtsA44FnqhibSIiUsfKPQI5N3w2+pPAl4GXgJXA2moVJiIi9a3cABk2szkEz0PfCnwLcODyKtUlIiJ1rtxJ9D7gb4DjgPnhtgy6kVBEZMoqK0Dc/ZvAN83sbHf/X2PbzewDB+kmIiKTWKyrsArDI2z/78qWIyIijaLcq7Ay7p43s6OLXrrU3a+rQl0iIlLnyp0DeQI4Gfg+sJ1g/gOCORAFiIjIFFTuHMjJ4Zd/XjQHckG5HxSuoXUOsBPIu/vycfZZAtwB9Lr7gwXbDTgXGAROA25x9x+W+9kiIlJ5cZcyecTMPgTMCNtnA18v1cnM2oFVwNvdfcjM1pvZ6e6+sWCfYwjC5cWivs3AvcDZ7j5qZl8HRmLWLSIiFRZ3KZMHgdOBY8L/usrs1wO84O5DYftRYHHhDu6+3d0fGafvAoJTZp82s2sJQuvVmHWLiEiFxT0C2e7uV4w1wqOGchwB7C1o7wm3lWM2QQCd6+6vm9laYBhYU7hTc3MTnZ3tZb7l5KaxiGgsIhqLiMaiMuIGyA4zOxN4luCRth8Dbimj306i014AM8Nt5dgDbHX318P2D4D3UhQgudwou3fvL/MtJ7fOznaNRUhjEdFYRDQWke7uGaV3mkDcU1iXElx1tQb4e+D8Mvv1AbPNbHrYPhXYYGZdZjazRN8ngFnhXAgERyTPxqpaREQqLpVH2rr7fjO7DFhpZq8Az7j7RjO7G+gH7jKzDHA9QUAsMbOsuz/s7v1mdjXwpbBvN3BrzLpFRKTCMvl8PnFnM2tz98EK1pNYNpvL65A0oMPziMYiorGIaCwi3d0zMqX3Gl+sIxAze0/RpvOBi5N+uIiINK64p7BWAD8iWok3V/GKRESkIcQNkE+6+6axhpldWuF6RGSKa21tJtMyjY62FgYGs+SzIwwP69+q9SjuaryF4dEB/LeKVyQiU1ZrazNDmSYuWbuZ469/iEvWbmYo00Rra3PpzpK6uHMgrwGvEZzC2kPweFsRkYrItEyjd+1m+rbtAqBv2y56121h9XkngY5C6k7JADGzJuAwd99FsHz7A0WvZ4Cu8HURkcQ62lrYtKP/gG2bdvTT0dbC0L6hCXpJrZQ8heXuo8ANZvbuccLjLQSLJMadSxER+TUDg1kWzDlwib0Fc7oYGMzWqCI5mHJ/8V9HcBPgAwQ3/o0As4CXCI5KflWl+kRkCslnR1ixdD6967awaUc/C+Z0sWLpfPJZLcBdj2LdSGhmhwLHAtOBF939P6tVWFy6kTCim6QiGotIo4xFGldhNcpYpCG1GwndfR/wTNIPExEpZXg4B8M5zXk0gLiLKYqIiAAKEBERSUgBIiIiiZQ9B2Jms4CPAm8APwF+XC8r8YqISPriHIF8m+BZHHcA9wCvm9nWqlQlIiJ1L06AzHD3W4FfuftpwLnAt6pTloiI1Ls4AfJG+OdQ+CCp9cD7qlCTiIg0gDj3gXzBzLqAB4D/YWaPAZ1VqUpEROpe2QESHnEA3Gtm5wMnAOdUpSoREal7iRZBdPdvVLoQERFpLLoPREREElGAiIhIIrECxMweN7OTqlWMiIg0jrhHIM+6++axRnh3uoiITEFxJ9GfM7OzgJ8BeeAK4KqKVyUiInUvboBcBiwsaB+NAkREZEqKGyDXuvuasYaZnVHZcqSepfGkOBFpHHGfSLjGzOYRLKrowMaqVCV1p7W1maFME71rNx/wrOrprShERKaouFdhfQ5YAVwAHAfcVY2ipP5kWqbRu24Lfdt2MTKap2/bLnrXbSHTkuheVBGZBOJehdXh7guBn7r7I8Duypck9aijrYVNO/oP2LZpRz8dbS01qkhEai1ugDSHf+bDPzsqWIvUsYHBLAvmdB2wbcGcLgYGszWqSERqLW6A5MzsX4APmtk/A/uqUJPUoXx2hBVL59MzdxbTmjL0zJ3FiqXzyWdHal2aiNRIJp/Pl96rgJm9D5gHPO3u/1qVqhLIZnP53bv317qMutDZ2U41xqIRr8Kq1lg0Io1FRGMR6e6ekUnaN+4k+kfd/Xvu/gXgV2Z2T9IPlsYzPJxjaN8Qu14dYGjfUN2Hh4hUV9xTWL8z9oW7PwMkTi4REWlsZV2DaWa9wDKg08w+ThAcI8CGqlUmIiJ1rawAcfcVwAoz+4i7/1OSDwrvWj8H2Ank3X35OPssAe4Aet39waLX2oAngO+5+5VJahARkcqJeyf6PxXdif4Ldy85C29m7cAq4O3uPmRm683sdHffWLDPMQTh8uIEb3Mb8KM49YqISPXECpDwTvQPAD8H1gCfBq4uo2sP8IK7D4XtR4HFFCyF4u7bge1mdvM4n3t+2GceE9x70tzcRGdne9nfy2SmsYhoLCIai4jGojLirkPR4e4Lzexqd3/EzE4ps98RwN6C9p5wW0lm9nvA77r7deHRz7hyuVFdlhfSJYoRjUVEYxHRWES6u2ck7pvWneg7gcIqZ4bbyvEh4A0zuwZ4N3CymS0rs6+IiFRJ3COQsTvR283sZGBzqQ6hPmC2mU0PT2OdCnzVzLqAEXffM1FHd7997GszO4TgKOhLMesWEZEKi3UE4u43A/cC3wVWA/eX2W8/wcOoVprZbcAz4QT6NcDlAGaWMbMbgNnAEjNbVPgeZvZh4D3AKWZ2bpy6RUSk8mItZWJmHcCZBKejMsAH3f0jVaotFi1lEtH53YjGIqKxiGgsIr/JUiZxT2E9CGwBXgvbXRPvKiIik1ncAHnO3ZeNNczsuMqWIyIijSJugDxsZhcCz4ft84GLK1uSiIg0grgBciEwRPQkwhMqWo2IiDSMuAHyqrtfMNYws3dWuB4REWkQcQPkaTNbSHQK62y0PpWIyJQUN0CuBLYWtI8Gbq1cOSIi0ijiBsi17r5mrBEu0S4iIlNQ3DvR1xRtSr4Kl4iINLS4y7mfDtwJHE5wJ/pM4NtVqEtEROpc3NV4zwUWEayDdTxwT8UrEhGRhhA3QNzdXwOmuXsWOKwKNYmISAOIO4l+mpk9BRxiZl8DtJSJiMgUFTdAlgA54HHgE+gUlojIlBU3QDYCl7v7ZmBlFeoREZEGEXcO5NkwPAAws1kVrkdERBpE7OXczews4GcEz0W/Ariq4lWJiEjdixsglwELC9pHowAREZmStJSJiIgkoqVMREQkES1lIiIiiWgpExERSURLmYiISCJaykRERBIpO0DC+z/uAp5ES5mIiEx5ZQWIma0E5gGzgNvdfR1aykREZEordw6kyd3fC7wT6KleOSIi0ijKDZCdAO4+Auwa22hmH6pGUSIiUv/KnQNZZGYd4dd/UPD1Keg+EBGRKancABkG9oVf/2vB9mxlyxERkUZRboBc5e6bijea2bsqXI9ITbS2NpNpmUZHWwsDg1ny2RGGh3O1LkukrpUVIOOFR7j9qcqWI5K+1tZmhjJN9K7dzKYd/SyY08WKpfOZ3opCROQg4t6JLjLpZFqm0btuC33bdjEymqdv2y56120h0xL3PluRqUUBIlNeR1sLm3b0H7Bt045+OtpaalSRSGNQgMiUNzCYZcGcrgO2LZjTxcCgrhEROZjUjtHDh0+dQ3BPSd7dl4+zzxLgDqDX3R8Mtx0L3AZsBo4Cdrn7rWnVLZNfPjvCiqXz6V235YA5kHx2pNalidS1VALEzNqBVcDb3X3IzNab2enuvrFgn2MIwuXFou5dwDp3/06430/NbIMm8KVShodzTG+F1eedpKuwRGJI6wikB3jB3YfC9qPAYuDNAHH37cB2M7u5sOM4V4A1Ed2TIlIRw8M5GM4xtG+o9M4iAqQXIEcAewvae8JtsYRLpzzs7luLX2tubqKzsz15hZOIxiKisYhoLCIai8pIK0B2cuDz02eG28pmZguBhcCy8V7P5UbZvXt/0vomlc7Odo1FSGMR0VhENBaR7u4ZpXeaQFpXYfUBs81setg+FdhgZl1mNrNUZzNbTPAo3V7gSDPTisAiIjWWSoC4+37gMmClmd0GPBNOoF8DXA5gZhkzuwGYDSwxs0Xh9ncBDxAs3PgI8B3A0qhbREQmlsnn87WuoSKy2Vxeh6QBHZ5HNBYRjUVEYxHp7p6RSdpXNxKKiEgiChAREUlEASIiIokoQEREJBEFiIiIJKIAERGRRBQgIiKSiAJEREQSUYCIiEgiChAREUlEASIiIokoQEREJBEFiIiIJKIAERGRRBQgIiKSiAJEREQSUYCIiEgiChAREUlEASIiIokoQEREJBEFiIiIJKIAERGRRBQgIiKSiAJEREQSUYCIiEgiChAREUlEASIiIokoQEREJBEFiIiIJKIAERGRRBQgIiKSiAJEREQSUYCIiEgiChAREUlEASIiIokoQEREJJFpaX2QmZ0BnAPsBPLuvnycfZYAdwC97v5gnL4iIpKuVI5AzKwdWAV8xt1vAeaZ2elF+xxDEBAvxu0rIiLpS+sIpAd4wd2HwvajwGJg49gO7r4d2G5mN8ftC9DS0pzp7p5RjdobksYiorGIaCwiGovfXFpzIEcAewvae8Jt1e4rIiJVklaA7AQK435muK3afUVEpErSCpA+YLaZTQ/bpwIbzKzLzGYm6VulOkVEpEyZfD6fygeZ2ZnAnwCvAFl3X25mdwP97n6XmWWA64GLgB8Aa9394Yn6plK0iIhMKLUAqQQzOxK4DTjR3ReM83oTwWXAe4E5wN+6++OpFpmSMsbi48ApwPPAScBfuftjqRaZklJjUbDfaQQXX8x39/+bVn1pKmcszOxjwKzwvxPd/YMplpiaMn5GDgNWA08DxwM/cPf7060yHWZ2LMFYbAaOAna5+61F+xwCfAH4BXAccJe7P3uw9220GwnfDXwHyEzw+keBme5+O3A18HUza06ruJSVGou3Asvc/R7gSwQ/KJNVqbHAzI4AlgAvpVVUjRx0LMzsD4HZ7n6vu18PXJdmcSkr9ffik8DL4e+LZcBfhf8InYy6gHXufo+79wJLzexdRfssA37u7ncC9wF/W+pNG2qw3P2fOfCKrGKLCeZMcPd+4A3g7SmUlrpSY+Hut7v7G2GzCRhIpbAaKDUWBUem16dWVI2U8TPyp0CTmfWa2R3AZP0HVjlj8SugO/y6G9ji7qNVL6wG3H2Tu3+nYFMTsK9ot8Lfnz8GTiw1R91QAVIGXfJbJJxb6gU+W+taauga4H53f63WhdSB2cDR7r6C4HTFt8NTOVPRWqDVzL4K/A3w1RrXkwoz+xDwsLtvLXop9u/PyRYguuS3QBge9wBr3L2v1vXUQnhe9x3AQjO7BvgvwEVTeDWDPcAT8OZR+n8CJ9a0otq5G3jK3S8H3g/cbmaT8ozFGDNbCCwEPjPOy7F/f6a2Fla1mNmhQLu7v0Jwee97gG+YWRdwCPCTWtaXpsKxCOd+7gPWu/v3zezD7r6+xiWmpujvxX8v2H4pwcUVk3ISfTxFY7ERODbc3gQcCWyrYXmpKhqLtwHPALj7oJm9Dkw/WP9GZmaLgT8kOCPxFjObDTgw4u57CH5/9gD/YWYnAE+H2yfUaFdhnQZcAJwF/DXwReDPgBPc/dLwB+JOYD9wNMFpi8l6FVapsbiX4NLnsV8Ox7r722pSbJWVGotwnxaCCyv+AvgGsMrdf1qbiqunjL8XrcDngV0EF1o84e5ralRuVZUxFr9LcGXSjwhO1fSH6+1NOuGE+feBJ8NNhwJfAX6P6FaKNoLTmi8Dvw3cUeoqrIYKEBERqR+TbQ5ERERSogAREZFEFCAiIpKIAkRERBJRgIiISCIKEBERSUQBIiIiiShARGrAzC4xswlXSK7HVaTrsSaprYZfykQmHzM7mWCdolbgewQrpY4CvbVeLdXMlrn7lyrwVicQLqNR9P7tBCsH307wALWaMrM/Bb7s7ocBR5nZp4Dr3H2kxqVJHdARiNQdd/8h8O/AY+5+i7t/CphHsCRFrS2r0PvMA348zvaVwP8M12qqOXf/B+D18OsXgMeBW2pZk9QPHYFI3TOzacDhwICZbQD+D2DAP7r7v5nZhQRroK0iWCjQgD8CHjjIvl8k+CV+OPB3wCKCp7B90N33mNmtBD8fOWCvu99tZh8FOs3sFmCru6+bYL9fq8fdTy76tt5B0RGImXUQPJPh4rA99j6rCZZhnztWX/j6TUALwQOTht391nE++6ywtlLfb8d44zXO/46HgK+Y2U21PhqU2tMRiNSzPwh/Wd9HcErnSeA+d/88cCXBL0rc/e+ArcAP3f184HKCU14H23dzuO8QMMPdLyJYVO9MM1sEnOLu17n7jcAHzGy+u38L2B0eFa07yH7j1fMmM3sbsM/ddxd9v8cSLGyXL6r1MXf/OMERy5nheywCTnb3G939BqDHzN43zme/v9T3G372uONVzN0HCR5C9VsT/U+TqUNHIFLPHitcHTVcivu9ZtYDZImeJjfmZwDu/mQZ+z4f/rm74OvXCJ6HMBdoD58fAvDiOP0h+Bf9wfZ7s56ifuPOfxAsJT7e3MLYiqivED2vYR4HLsP+HMFzPb5X/NlmBgf/fiE4ijnYeBXKAm0HeV2mCAWINJJPAP/V3f8sXJb80qLX8zH2PZingR53vwvAzP6I4Bc0QC58UNe8EvsV11NoovmPFwmeXV1svPd5GnhvQfs44LtlfPZEyhqv8HvvAH4Z8/1lEtJy7lJ3zOz3ia7Cum/sQVhm9jvA/QQTuf3AVQS/+PYQPJJ0PfD58IFapfb9e4JfuF8jOJXzeYJ5g9cITjldTPDMhL3AYcA17p4zs5UE/wLH3f/CzG4o3o9g/uWAeoq+v38ANrj7P47zvT8OfNjdf2FmZ4bv8w1gTVjra8Cl4fd4E8GRQAYYdPflBX3Wh9/T/DK/38MmGK9DgC8D17r7KjObDyx39z+e6P+fTB0KEJE6YmYnAlcQhESu1vUUCh8PfD9wo7vvqHE5UgcUICJ1xszeQvCY0bq4lHeMmb2V4Einv9a1SH1QgIiISCK6jFdERBJRgIiISCIKEBERSUQBIiIiiShAREQkEQWIiIgk8v8BLoWuw6m+MM0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot\n",
    "\n",
    "pyplot.scatter(*samples)\n",
    "    \n",
    "pyplot.rc(\"figure\", figsize=[6, 4])\n",
    "pyplot.xlabel(\"Parameter $I$ (normal)\")\n",
    "pyplot.ylabel(\"Parameter $a$ (uniform)\")\n",
    "pyplot.axis([1, 2, 0.1, 0.2])\n",
    "\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lagrange basis\n",
    "\n",
    "Creating an expansion of Lagrange polynomial terms can be done using the following constructor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.379804Z",
     "iopub.status.busy": "2021-05-18T10:57:05.379487Z",
     "iopub.status.idle": "2021-05-18T10:57:05.393605Z",
     "shell.execute_reply": "2021-05-18T10:57:05.393273Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial(-694.24*q1**2+28.97*q0*q1+170.1*q1-6.65*q0-5.96)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import chaospy\n",
    "\n",
    "polynomial_expansion = chaospy.expansion.lagrange(samples)\n",
    "polynomial_expansion[0].round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On can verify that the returned polynomials follows the property of evaluating 0 for all but one of the samples used in the construction as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.395549Z",
     "iopub.status.busy": "2021-05-18T10:57:05.395232Z",
     "iopub.status.idle": "2021-05-18T10:57:05.404406Z",
     "shell.execute_reply": "2021-05-18T10:57:05.404070Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1., -0., -0., -0., -0.],\n",
       "       [ 0.,  1.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  1.,  0.,  0.],\n",
       "       [-0., -0., -0.,  1.,  0.],\n",
       "       [-0., -0., -0., -0.,  1.]])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polynomial_expansion(*samples).round(8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model approximation\n",
    "\n",
    "Fitting the polynomials to the evaluations does not need to involve regression.\n",
    "Instead it is enough to just multiply the Lagrange polynomial expansion with the evaluations, and sum it up:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.406783Z",
     "iopub.status.busy": "2021-05-18T10:57:05.406469Z",
     "iopub.status.idle": "2021-05-18T10:57:05.419236Z",
     "shell.execute_reply": "2021-05-18T10:57:05.418903Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000,)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from problem_formulation import model_solver\n",
    "\n",
    "model_evaluations = numpy.array([\n",
    "    model_solver(sample) for sample in samples.T])\n",
    "model_approximation = chaospy.sum(\n",
    "    model_evaluations.T*polynomial_expansion, axis=-1).T\n",
    "\n",
    "model_approximation.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assess statistics\n",
    "\n",
    "The results is close to what we are used to from the other methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.421134Z",
     "iopub.status.busy": "2021-05-18T10:57:05.420871Z",
     "iopub.status.idle": "2021-05-18T10:57:05.488642Z",
     "shell.execute_reply": "2021-05-18T10:57:05.488340Z"
    }
   },
   "outputs": [],
   "source": [
    "expected = chaospy.E(model_approximation, joint)\n",
    "variance = chaospy.Var(model_approximation, joint)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.491006Z",
     "iopub.status.busy": "2021-05-18T10:57:05.490691Z",
     "iopub.status.idle": "2021-05-18T10:57:05.557041Z",
     "shell.execute_reply": "2021-05-18T10:57:05.556754Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEHCAYAAAC9TnFRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA96klEQVR4nO3deXRkV33o+2/Ng6o0S62eZ2/P7RFsbAzGZggG3gshlyEhkEC4JBcevmEIuQ8STLhhhkDIHO5jreQm4QaSADGGYIONMcaz3bbb/XPPg1qzVKUq1Ty8P06VTrWs7q5S61SVpN9nLS91nTpHZ/d21/nVnn7bVS6XUUoppc7F3eoCKKWUWhk0YCillKqLBgyllFJ10YChlFKqLhowlFJK1UUDhlJKqbp4nfzlxpidwKeAx4FNwJSIfHLBOUHgC8AwsBv4jIg8X3nv14ErgSJwSET+2snyKqWUOjOnWxi9wD+LyOdF5APAW4wxVy8453bguIh8Gvgy8HUAY8wm4EPAh0TkI8C7jTG7HS6vUkqpM3A0YIjIIyLynQX3m1tw2m3Ag5Xznwb2GGM6gVcDj4lIdWXhg8AvOVlepZRSZ+Zol1QtY8wvAz8Ukf0L3hoEEjWvZyvHznT8NKVSqVws6mp1AI/HTbFYanUx2oLWhU3rwqZ1YfP5PK5Gr2lKwDDG3AzcjNX9tNA4EK153Vk5Ng7sWnD84MKLi8UysVhq2cq6knV3h7UuKrQubFoXNq0L28BA9NwnLeD4LCljzG1Y3UsfAIaMMdcbY3or3U4AdwLXV869DHhKRGaBHwJXG2OqUfB64C6ny6uUUmpxjgaMygD3N4HrgJ8A3wEM8FHgdyunfQXYaoz5GPBB4F0AInISa/bUl40xXwT+TkQOOFlepZRSZ+Za6dlq8/liWZuYFm1u27QubFoXNq0L28BAtOExDF24p5RSqi4aMJRSStVFA4ZSSqm6aMBQSilVFw0YSiml6qIBQymlVF00YCillKqLBgyllFJ10YChlFKqLhowlFJK1UUDhlJKqbpowFBKKVUXDRhKKaXqogFDKaVUXTRgKKWUqosGDKWUUnXRgKGUUqouGjCUUkrVRQOGUkqpunid/OXGmCHgU8AeEbl2kfc/AbwSyFcOXQC8TUTuNcb8AshUjhdF5BYny6qUUursHA0YwI3Ad4ArzvD+o8DnRWTOGOMG7gTuq7z3AxH5hMPlU0opVSdHA4aIfMsY8/KzvP8fNS/fANwpIuXK68uMMb8PhIBHRORO50qqlFLqXJxuYTTiHcDba15/VkQeNsZ4gJ8aYxIi8tOFF3k8brq7w00rZDvTurBpXdi0LmxaF+enLQKGMeYK4KCIJKvHROThys+iMeZ+4GbgBQGjWCwRi6WaVdS21t0d1rqo0LqwaV3YtC5sAwPRhq9p+iwpY0yHMWZgweH3AV+rOedCY8y7at7fDRxqRvmUUkotzulZUi/D6mZab4z5GPBF4J3AZcB7K+esA4Iicqzm0lngNmPMBqATOAH8o5NlVUopdXaucrl87rPaWD5fLGsT06LNbZvWhU3rwqZ1YRsYiLoavWbFL9xLZgutLoJSSq0JKz5gPHkixngi2+piKKXUqrfiA0apXGb/eJIjUylWeveaUkq1sxUfMKpOxNLsG01QKJZaXRSllFqVVk3AAJhK5XlyeJZUrtjqoiil1KqzqgIGQCpf5InhOJNzuVYXRSmlVpVVFzAAiqUy+0YTHJ3W6XNKKbVcVmXAqDo+k+bZkVkd11BKqWWwqgMGWOMaTwzHdb2GUkqdp1UfMADS+RJPDc8ypus1lFJqydZEwAAolsvIeJIDE0lKul5DKaUatmYCRtXIbJanhmfJ5HXqrVJKNWLNBQyARLbA4yfjTOnUW6WUqtuaDBgAhVKZZ0cTHJ6a05QiSilVhzUbMKpOxjLsPaVdVEopdS5rPmAAxDPaRaWUUueiAaOi2kV1aHJOZ1EppdQiNGAsMBzP8ORwnLR2USml1Gk0YCwimS3y+Mm4LvRTSqkaGjDOoFiyFvrtH0tQKGkXlVJKeZ385caYIeBTwB4RuXaR918O/CkQqxy6U0Q+X3nvVuCNwDhQFpE7FruH01Nix5M5ZjMxLlwXoTPoc/ReSinVzhwNGMCNwHeAK85yzu0icm/tAWNMGPgr4BIRyRpjvm2MuUVE7ll48Vd/fIg3Xj5ENODcXyVTsHJRbekJsaUnhMvlcuxeSinVrhwNGCLyrUor4mzeboy5BugE/lZETgDXA8dEpDqI8ABwG/CCgLF3OM7RqTnecf1WLtnQtYylf6HJXIn8bI5L1ncS8nscvddSeDxuurvDrS5GW9C6sGld2LQuzo/TLYxz2Qf8sYgcNcZcAvzIGHMxMAgkas6brRx7AbMugowl+bOfHOKlO3q57ZJ1+DzODc3MzWUZm55jZ1+Yoc6gY/dZiu7uMLGYbhoFWhe1tC5sWhe2gYFow9e0dNBbRMZF5Gjlz88C3cBmrHGL2r9NZ+XYC3zolRfw2osHcbvg/sPTfOW+w4zOZhwtd7FU5vmJOZ4dTZAr6OZMSqm1oekBwxjTYYwZqPz5o8aY3sqfewE/MAY8CGw1xgQql90A3LnY73O7XdxywQDvv2kH/R1+RmazfPnewzxweNrxAfGpuRyPnYzp/uFKqTXB6VlSLwPeDqw3xnwM+CLwTuAy4L3AEeArxph9wMXAb4hIpnLt7wBfNcZMAHsXG/CutaUnxO/dvIN/3zvKw8dj/OveEfaPJ3nzlRuIODggni9a+4eviwTY2R/G62B3mFJKtZKrnm/hxphHgaeApyv/7RWRCYfLVpd79o2WYwu6oJ4ajvMvT54inS8RDXh569UbMYMRx8sS8Li5YLCDnrDf8XstRvtnbVoXNq0Lm9aFbWAg2vB0z3q/Dr8B+BesLqP/Chw1xhxr9GbNsmdjFx+8eRc7+sIksgX+5ufH+NenRsg6PN6QLZZ4eiTBgYmkLvZTSq06dfXViMgp4BTwAwBjzEXAmxws13nrCfv4nRu38ZMDk/zwuQkeODKNjCd529Ub2drr7LS6kdks06k8Fwy0rrWhlFLLra4WhjFma+1rEXkOuMCREi0jt8saEP/Ay7Yz1Blgci7Hn/30CHftG6NQcri1UbBaG8+Pa2tDKbU61Dsa/E/GmC1Yg9RPY6XyuNSpQi23jd0hbn/ZDn7w3Dj3HZzi7ucn2TeW5G1XbWR9l7NrKUYTWWZSeXYPdtCrrQ2l1ApWVwtDRF6CtT7iN4EfAYeA1ztYrmXn87h5/aVD/O6N2+gN+zgVz/Dl+w7zkwOTju9/kS2WeGYkgYwlyRd13YZSamWqa5ZUO1tsltS5ZPJFvvfsGL84OgPA9r4wb7lqI/0dzrcAfB4XO/s6GIwGzn1yg3QGiE3rwqZ1YdO6sDk5S2pVCfo8/OoVG3j3dVvoDHo5MpXiiz8+xM8OTzne2sgXy+wfT/LsiO4jrpRaWdZkwKi6aCjKh16xkys2dpIrlvi3vaP85c+OMpF0fuOkqVSex07GORlLO74iXSmllsOaDhgAHX4vb792M+940WaiAS+Hp1J84ceHuLcJYxvFUpnDUymeHJ4lmS04ei+llDpfdefMMMbswsogmwd+A/g3EWnbxXuNunxDJzv7w3zn6VEeOxHne8+O8dSpWd585QbHs9ImsgWeOBlnQ1eQrb1hvG7db0Mp1X4aaWF8HAgBXwI2AH/kSIlaqMPv5W1Xb+Jd122hK+jl+EyaL917mLtlgqLDaynKwHA8w2PHNZmhUqo9NRIwngCGgQtF5CPAfmeK1HoXD0X58C27uG5rD8VSmbueG+er9x3mVNzZtOlgTcHdN5rQQXGlVNtpJGBcDnwV+E9jTAjY4UyR2kPI5+FXr9zAf33JVnpCPk7GM3z53kPctW+sKWspplJ5HjsR5/hM2vGxFKWUqkcjAeMzWK2KTwPXYS3gW/UuGIzw4Vt2csP2XkpluPv5Sb7w40McmEg6fu9iuczR6RSPnYgxndJuKqVUazUSME4Bx4FfBbYAb3GkRG0o4PXwxj3red9Lt7MuauWk+qsHjvFPjw03ZXZTOm+tFH92NKHdVEqplmkkYPwH8ApgO7AN6HWiQI1q5nyi7X1hfu/mHfzSRYN43S4ePRHjs/cc5NHjsaaspZiay/HYiTjHplOOD8IrpdRCjWxFd0RE3l99YYzZ7kB5GnbZxi4eTedJNembt9ft5lYzwJ6NnXzryREOTs7xT48P8+iJGG/as57+yPKn/KhVLJc5NpNmLJFlR1/Y8fsppVRVIy2Mo8aYVxpjtlYy177DqUI1ojvs56rNXWztCdHM5QsDkQDvvWErb71qI2G/hwMTc3z+x4e4WyYcT50OkCmU2DeWZO+pWeZyuuhPKeW8upMPGmNGOH0q7RYR2elIqRqQzxfL1WRi6XyRAxNzxNL5ppYhmS3wvWdGefREHIB10QBvvHw9uwY6mnJ/F7C+M8jl2/tIJZ2f+rsSaJI5m9aFTevCtpTkg410Sf2BiHyj+sIYc+u5LjDGDAGfAvaIyLWLvP9OrBlXh4CrgD8TkZ9X3vsFUH36FUXklnPdL+TzcPmGTsYSWQ5PzZEvNqefPxLw8tarN3H1lm6+/eQIY4ksf/nAUa7a1MXrL11HZ9Dn6P3LwKnZDKkjU/T63GzoCuJ26WpxpdTyqjtgiMg3KkFiD/CkiNxdx2U3At8BrjjD+xuB20UkY4x5MfB3wGWV934gIp+ot3y11kUD9IV9HJ5KMZpwPpFg1QUDET70ip3ce3CKu2WCx0/GeXY0wWsuHOSGHb14HO4zyxfLHJ5NMTKbYXtfR1PStSul1o5GuqQ+DlwPHAR2Az8XkT+u47qXA18QkWvOcd71wJdE5PrK628DD2OlI3lERO5c7LpSqVwunmUhXSyVQ8aSTU/uN5nM8n8eO8nek1Y31cbuIG+5dgu7ByOO3dPjdp02e6on7GP3YISowy2cduTxuDnbv4u1ROvCpnVh8/k8jnZJ+UXktdUXxphPN3qzMzHGuIAPAL9Xc/izIvKwMcYD/NQYkxCRny68tlgsnbNP8oLuAMPxMsem0xSbtGo65IJ3XLOJfZu6+Le9IwzHMnzxR89zzeYuXnfJENFgI1Vfn46OAHNzdotqbi7LyYkk6yIBtvaGCPo8y37PdqV91TatC5vWhW1gINrwNY3MkloYlpclTFeCxeeBb4jIg9XjIvJw5WcRuB+4ean3cLlcbOoOcfXmrqZ301w8FOUjt+ziVWagsnYjzmfuPsD9h6aatpZiLJnl0RMxjkylKOi3K6XUEjUSMArGmO8aY/7UGPM9IL2UGxpjOowxA5U/e4CvAN8TkR8YY36lcvxCY8y7ai7bjTUwfl6CPg8XD0W5dH2UkK95W4H4PG5efdEgH75lFxeti5AplPj3p0f58r3NSTECUCrDiViaR07EOBnT/FRKqcY1tKe3MeZVWEkInxKRc+aSMsa8DGvvjNcAfwl8Efgt4DIRea8x5kvAm4DDlUt2ishmY8wG4GtYGXI7AR/weyLygq/HtdNqG1Eqlzkxk+ZkLNO0biqAcrnMs6MJvvP0KNMpa/rvZeujvP7SIfrOs/WzsEvqbIJeN9t6w47sLd4OtOvBpnVh07qwLWVa7TkDhjHGJSLlymK9Wu8Vkf/R6A2X21IDRlUmX+TQVIqpJu9BkS+WuO/gFPc8P0muWMLjdvGynX3cckH/kscaGgkY89f4PWzvC9MbXl0zqvTBYNO6sGld2Jxah/EQ8CLgPuBI5ZgLKwFhywPG+Qr6PFwyFGU6lePQ5BzpfHP6+H0eK8XItVu6+f6+MR49EefHByZ55HiM1148yDVbupuylmIuV+SZkQRdQS/besN0hdbejCqlVH0amVb7ehH5Xs3r14rI9x0rWZ3Ot4VRq1QuczKW4cRM82ZTVR2fSfHve0c5NmMNDW3uDvJ/Xbae7X3hun/HUloYC/WFfWztDRMJLP8srmbSb5I2rQub1oVtKS2MRkZ+Q9U/GGMu5zxmLbUrt8vFlp4Q12zpZiDS3C6aLT1h3nfTdt529UY6g15OxDJ87f4j/MOjJ5u6F8ZUKs/jJ+M8N5YgldNU6kopWyMB48LqH0RkrwNlaRsBr5uL1kW5fEMnHf7mrV1wu1xcvbmbj966i1sv6MfrdvHEyTifvfsg33tmlHQTH+ATyRyPnYgh40ndg0MpBdQ36P0B4HagG5jBGr8oAHeKyO3OFu/clrNLajHlcpmR2SzHZlJNy01VNT2X4/vPjfNEZbV42OfhlWaAl+zowet+Yaxfji6pxbhdMBgJsKVn5Sz+064Hm9aFTevC5sgsqSpjzK+KyL80XCqHOR0wqgrFEkdn0ozEMzR7BcOJmTTfe2aUQ1PW37Mv7OO1l6xjz4ZOXDUD404FjCq3C4aiQTb3hAh4m7eOZSn0wWDTurBpXdgcDRgLGWN+WUT+bUkXL6NmBYyquVyBw5MpZpqcQr1cLrNvNMl/PDvKeNIa09jSE+L1l65jR5+VRt3pgFG1EgKHPhhsWhc2rQub0y2MW4BPA/1Y3VKdItLX6A2XW7MDRtXUXI7DU82bhltVLJV5+NgMP9w/QaKSUPHS9VFuu3gd24c6mxIwqtwuKzPw5u7266rSB4NN68KmdWFzej+MtwKvBt4DfAn4YKM3W036Ovz0hH2cimc4PpOm0KS8UB63i+u393Llpi7uPTjFfQcneWYkwb7RBDfs7Ofmnb1NW0tRKsPIbJaxRHbFjXEopRrXSH+CiMgM4BWRPNDjUJlWDHclqeGLtnRXNi1q3r2DPg+vuWiQj966mxdv7aZchvsPTvInPzrA954Zbeq2raUyjCayPHI8howldTquUqtUI11S/4HVsrgZWA/sFpGXOVi2urSqS2oxqVyRI9PNTzMCMJbIcveBSR4/HgOsqcEv39XHTTv7mv6t3wX0R/xs7g61bAGgdj3YtC5sWhc2p8cwOrBSmruAdwM/FBFp9IbLrZ0CRlUslefw9BzJbHO/aXd0BNg/HOOufePIuJUFt8Pv4ZYL+nnJ9l58nuYPUPeGfWzuDjU95Yg+GGxaFzatC1uzZ0mFRGRJKc6XUzsGjKqxRJaj0ymyheYMjNfOkjo0Ocf3941xdNr6X9QV8vIqM8i1W7od3yp2MV1BL5u6Q+edkbde+mCwaV3YtC5sTrcwblpw6O0i8tuN3nC5tXPAACs/1XAsw4mY8wPjC6fVlstlnhtL8v19Y4zMWscHIn5eZQa4YlNXU5IbvqCMfg+bukMMRvynrSFZbvpgsGld2LQubE7PkvoK1v4U1Uy1OrJZB7fLxeaeEEOdAY7PpBmZzdCkCVW4XC4uHopy4boITw7P8sPnxplI5vjfjw3zI5nglRcOcsXGzqYGjrlcERlPcnTazcauIOs7gy1p8SilGtdIwHiPiDxSfWGMea8D5Vm1fB43O/s72NgV5Oh0an7xXTO4XS6u2tTFng2dPHo8xt3PTzCezPG/Hz3J3RLglWaAPU0OHNlCicNTKY7PpNnQGWRDVxB/my4CVEpZljSGYYyJAH8mIr+5/EVqTLt3SZ1JMlvg8FSK2DKuGK93pXexVJ4PHNVd/9ZFWxM4qqr5qjZ2B+nwn//MKu16sGld2LQubE6PYcxgJx+cBb4mIn/b6A2X20oNGFUzqRxHplPLMqOq0dQghVKJR4/HuVsm5lOdtDpwgDWzamNXkJ7z2AVQHww2rQub1oXN6YDxZhH5ZsOlcthKDxhVE0lrRtX5pBpZai4pK3DEuFsm5wPHUDTArS0OHB1+Dxu7ggxGAw2XQR8MNq0Lm9aFrdnTat8jIn9zjnOGgE8Be0Tk2kXedwN/AiSAbcDXReQXlfduBd4IjANlEbljsXusloAB1qym0USW49NpssXGA8f5Jh8slEo8cjzGPTWBo7/Dzysu6OfqzV2LplRvBr/HxVBnkA2d9Y9z6IPBpnVh07qwORIwjDHTQKx6PlCmzuSDxpg3AVngj0TkmkXefwtwk4j8rjGmF/gFcBEQAPYCl4hI1hjzbeAvROSehb9jNQWMqlK5zKm4NRW3kT04litbbTVw/OT5SaYqYxxdIS837+rnxVt7WjY47XbBQIc1znGuFeT6YLBpXdi0LmxOTat9n4j848KDxpi3netCEfmWMeblZznlNuA/K+dOG2MywCXAAHBMRKpPvwcq574gYKxG1RxVQ51BhmNpTsYzFJs1Fxfwut1cv62XF23p4cnhOD9+fpLRRJZ/f3qUu2WCm3b18ZLtvYSanHKkVIaxZJaxZJauoJcNXUH6Ovwt6zJTaq05Z8BYLFhULMcq70Gs7qiq2cqxgTMcfwGPx013d3gZitKe+ns7uLhY4thUipOx9FkDh8ftoqMjsKz3v+nCIDeaQfaejHPXM6Mcm07x/X3j/OTAFC+/YIBXXDhANNjctB9gbfl4PJlnLFNkY3eIDd2n782x2v9dNELrwqZ1cX7qnr+42H4YwPluoDQORGted1aOlc9w/AWKxdKaaGL2+9109oY4HkszeobFf05uoLS7N8Sul27j+Yk57nl+gkOTKe56dpR79o9x3bYebtrZT0+4+YFjDpiOp3n2uDXesr4zSFfIp10PNbQubFoXtoGB6LlPWqDp+2FUkhiGRWQCuBO4Cfj7yhhGEHgWawxjqzEmUOmWugH4i6XcbzXxe93s6u9gc3eIEzNpRhPNWzUO1spxMxjBDEY4MpXinucneG4syU8PTfOzw9NcubGLl+3qY2N3qHmFqiiVYTyZYzyZo8Pv4UJcBEtlXUWu1DJqJGCIiMwYY7wikjfGnHM/DGPMy4C3A+uNMR8Dvgi8E7gMeC/wf4ArjTF/hJVu5DdEpAikjDG/A3zVGDMB7F1swHutCnjd7BroYFN3kBOxNGOJbFMDB8D2vjDvvn4rw7E0Pz4wyd5Tszx2Ms5jJ+PsHujg5t39XDDQ4Wi+qDOZyxXZP5ogk86xLhpgfWdgWRYDKrXW6X4Yq0AmX+T4TJq5MiSSzduitdb0XI6fHprioWMxcpUpwes7A7x8Vz9XbOps+pTchd1zXUEv6zuD9EfW3iC5dsPYtC5suh/GGhcIB3j66BTjyea3OKpSuSIPHp3m/kPT83uOdwW9vHRnH9dt62nazKozjef4PK5KqyPY9FleraIPSZvWhc3pgPEGEfluw6VymAYMW/XDkMkXW9ZVVVUolnj8ZJx7D04xlrAe3AGvm+u29XDjjl56zyPtRz3qmQCwVlod+pC0aV3YnA4YDwP3Af8gIk81eiOnaMCwLfwwZAslTpxlVlUzlMplZCzJTw5OcmjSKpsLuGxDJzfu6GVHX9iRcY5GZox53VarY2iVjnXoQ9KmdWFzOmBcBBzFGsTeA/xERL7V6A2XmwYM25k+DLlCiZPxNCPxLMUlpoJZDidm0tx3aIqnhuPzAWxjV5Abd/Ry5aauZd1CdqlTjKMBL+s7A/RHAnhXyQwrfUjatC5sTgeMS4HnsabX/jdgVkRubfSGy00Dhu1cH4Z8scRwLMOp2Yzju/+dTTyd58GjMzx4ZJpkzsrSG/F7uH57Ly/Z3kPnMiwEPN81KR6Xi/6In3XRAN1N3o98uelD0qZ1YXM6YDyJtYDuLuDPRWRfozdzggYMW70fhkLJylV1Kp4m10CuquWWL5Z44mSc+w9PcyqeAawH9Z6Nnbx0Zx9bepa+nmM5FzGGfG7WRYMMRvwEV+BAuT4kbVoXNqcDxjeB3xaR2UZv4iQNGLZGPwzFUpnRRIaTsQzZwtLTqp+vcrnM4akU9x+a4pmRBNV/kVt7Q9ywvZc9GzrxNthd5cSqdxfQHfKxLhqgr8O/YhYF6kPSpnVhczpgeLATAwowLCKt+3paoQHDttQPQ7lcZjyZ42QszVyutVu1T8/leODINA8dm5nfG6TD7+FFW7u5flsvfR31za5yMk0KWHm7Bjr8DK6ALit9SNq0LmxOB4wPA68FjgPfAF4jIr/f6A2XmwYM23J8GKbmrMARzxSWqVRLky2UeOJkjAeOzMx3V7kAsy7CS7b1ctFQ5KxTYZ0OGLWCXjeD0QDrooG2XNuhD0mb1oXNqfTmVRERudkY8/si8hNjzHWN3ky1v74OP30dfuLpPCfjGabmci0ph7Vmo5cXb+3h+EyaB45M89TwLPvHkuwfS9IT8nHdth5evLWHaLC1U2EzhRLHZ9Icn0nTGfAyGA0wEPEv66wvpdpBI5+06lenapMkssxlUW2kK+SjK+QjlStyMpZu2epxl8vF1t4wW3vDvOHSAo8cj/Hg0Rmm5nLc9dw4/7l/gss2RHnJdufWdDRiNltgNlvg8NQcPSEfg9EAveGVM96h1Nk00iV1B/BiIIyVavxxEfkTB8tWF+2SsjnZ3M4VStbMqhZPyQVrMeDz43P8/Mg0+0btQfJ10QAv3trN1Zu7Wdfb0bQuqXPxuF30h/0MRv10h3xND2raDWPTurA5vqe3MeZVwOXAUyLyo0Zv5gQNGLZmfBiKpTJjiSzD8fT8oHQrzaTy/OLoNA8di83nrvK4XOzZ3MU1m7rYPdDRVmk/fB4XAx1Wl1VXkwbL9SFp07qwOR4w2pEGDFszPwzlcpmpVJ7hNhggByuQPTeW4KGjMzw3lpxvdfSEfLxoazfXbulpyQZPZxPwuhmI+BmMBM65R/n50IekTevCpgFjjWvVhyGZLTAcyzAx17pkh7Xi6TxPjib42YFJplN5wJ5h9eKtPVw8FGl6uvVzCfncVssj6l/2fFb6kLRpXdg0YKxxrf4wZAslRuIZRhIZ8i1cQQ7WtNpEMsPByTkeOjrD0yOJ+f3QI34P12zp5tot3Qx1BltazsWEfR76I34GIssTPFr976KdaF3YmhowjDHvEZG/WdLFy0gDhq1dPgylcpnxRJbheKZlCwEXrsOYyxV47ESch47OMJqwj2/qDnLN5m6u3NTlaLfQUnX4PfR1nF/waJd/F+1A68LmSMAwxkwDsYXXAZ0i0tfoDZebBgxbO34Y4uk8w5X1HM1sc5xp4V65XOb4TJqHj8V4cjhOppISxeNycdFQhGs2d3NRG3ZZQaXl0eGnP+JvKLi147+LVtG6sDm1cO99IvKPCw8aY97W6M3U2lNdz1HtrhpNZFqa8LB2Xcf/ffkQz4wkePREDBlL8sxIgmdGEoT9Hq7c2MW1W7rZ1B1s+dqOqlS+yPFYmuOxNCGfm/6OAP0d/pYvXFRrR6PTai+nwVxSxphbgTdird0oi8gdC97/OrCz5tBlwNUictQYcxRrDw4q9/u1hb9fWxi2lfDtqVQuM5nMcSqeYTbr3OyqRlODzGbyPHYizqMnYozO2tetiwa4ZnMXV2/ubto02EYFPG76In76w366Qt4XBLiV8O+iWbQubG2XS8oYEwb2ApeISNYY823gL0Tknppz3iwi36z8uRP4hoi8sfL6EyLyibPdQwOGbaV9GJLZAqfiGSaSuWXf2GmpuaTK5TLD8QyPHo/xxMn4/H4dLmBnfwdXburi8g2dhP3tlzMKrHUevWErvUtv2Ifb5Vpx/y6cpHVha8dcUtcDx0Sk+sl9ALgNmA8Y1WBR8VvA/6p5/VJjzEeAKHCXiPy8gfKqNhcJeLlgMMKOvhJjySwj8SypfGuz5bpcLjZ1h9jUHeL1lw6xfyzJoydiPDua4ODkHAcn5/jXp0a4cF2EqzZ3cfG6KH5v+4x35IvWwsqxRBaPy0V32McOlxtvsaS5rdR5czqX1CCQqHk9Wzn2AsYYN/Bq4Cs1h/9ARB6utFQeN8a8TkQOnlYoj5vu7nA95V/1VnJd9PdFuASYSeU4FUsznji/NR0et4uOjsB5l+tF0SAv2tVPKlfgiRMxHjk6g4wleHbU+i/gdXPFpm6u3dbDRes72y5nVAaQsQSlcpnukI/+SID+iJ/wKty7vB4r+TPSDhr5V1M0xvwACBtjXgQ8Xsc141itg6rOyrHFvAG4s3ZcREQervxMVXb8uwE4LWAUiyVtYlashua2C9gY9jEY8DCWyDIym1lSChIn0ptfMRTliqEos5k8Tw7P8sTJOMdn0jx0dJqHjk7T4fewZ2MnV27qYltvuG1SklTrIpnMcnIiaR3ze+gNW91WncEXjnusVqvhM7JcBgai5z5pgboDhoj8UW0uKaCe7qEHga3GmEClW+oG4C+MMb1AYcHufe8A3l59YYy5BfCJyA8qh3YBh+otr1rZfB73fNdQLJVnNJFhci7XFivJO4M+btrZx007+5hMZnlieJbHT8QYT+b4+ZEZfn5khp6Qj8s3dLJnYydbekJt90CeyxWZy6U5EUtb4x4hP30dPrrDfrxt1kpS7aOedRg3neGtt4vIb5/rBsaYVwJvAiaAvIjcYYz5HDAtIp+pnHMF8Gsi8uGa6y4DPgE8BmwATi2WHVcHvW2r/dtTvlia758/14LAZm6gBNZg+al4hsdPxnliOE48bc8A614QPJrd8mikLtwuKyD2hn30dfjbckOo87HaPyONcGrh3mNYM502AUHgMLADQERe2ngxl5cGDNta+jDMZvKMzmaZmMvNp/yo1eyAUatULnNsOs1Tw3H2npo9LTljV8jLng2dXL6hi629zQke51MXIZ+bnrCfvrC1nqZdutmWai19Rs7FqVlS/4+IPGCM+ZCIfKF60Bjzh43eTKnl0hn00Rn0sbNUZiJptTraIWsugNvlYntfmO19Yd5w2RDHptPsPRXnqVOzxNMFfnpomp8emqYr5OXy9Z3s2di84NGodL5EOp7hVDxjzboKeempjH0EV1nrQ53bOQOGiDxQ+eP2BW9tWv7iKNUYj9vFUGeQoc4g6XyRsVkreLSL2uDx+kuHOD6T5qnhWfaemiWWznP/4WnuPzxNV9DLpes7uXRDlJ19HW032wqgWElpP1XJABz2eegJ+yoD5762LLNaXo0s3PsK1sDzAeAC4JCIvN/BstVFu6Rs2ty2lMtlSn4fB07OMDWXX/ZFgcuhVC5zoiZ4zKTz8++FfB4uHopw6fpOzGCEwHmu82hG95zbBV1BHz1h67/lTtG+XPQzYmvGjnuvBS4B9onInY3ezAkaMGz6YbBV66JQKjPZZl1WC5XLZU7E0jwzkuDpU7OMJ3Pz73ndLsxghEvXR7l4KLqkjLqtGM8JeNzzwaM75GubRYP6GbE5vdIbIAeUsNYDKdX2vDVdVpl8kbFElvFkti22l61yuVxs6QmzpSfMay9ex1giy7Mjszw9kuD4THp+kaAL2N4X5rL1nVyyPkpfh7/VRT+jbLHEaCLLaCKLC2tVf0/ICh6dIW9bjteoc2ukS+rjWKk+DgK7gZ+LyB87WLa6aAvDpt+ebOeqi9lMnvFkjolktuWbPZ1NPJ3n2VEri+7BibnTutc2dAa5eCjCxUNRNp9lum4rZ4wtxuNy0Rn0zrc+mrkPiX5GbE63MPwi8trqC2PMpxu9mVLtYn6WVV+YmXSe8US2Lcc7ukI+XrK9l5ds7yWdL/LcWJJnRmbZP5bk1GyGU7MZ7n5+kojfw4Xrolw8FMEMRtp6BlOxXGYmnZ8ft/F5XHSHfPMtkHYu+1rXSMBY2IZvnza9UkvkcrkqKTL8FEtlpuasVsdMOt8Wq8prhXwertrUxVWbuigUSxycnGPfaJLnxhJMp/I8eiLGoydiuF2wo69jvvWxHDm1nJQvlplI5piojN0EvO75ANIV8p33oL9aPo0EjIIx5rtYC/d2Ag85UySlWsPjdjEYDTAYDVAolpicyzGezBFP55u6W2A9vB43F66LcuG6KOXyEGOJLPtGEzw3luTIVGo+s+53nxljXTTAhYMRLhqKsL0v3Ja7CdbKFuwV/WBN3+0KeekO+egK+toqO/Ba0+gsqflcUiLyI8dK1QAdw7Bp/6xtOesiWygxOZdlMplr25lWtVK5AvvHkuwbS7J/LEm6JmV80Otm90AHZp3VddUbbt+B8zPp8HvoCvroCnkbDiD6GbE5lRrkt0Tkfy1y/K0i8k+N3nC5acCw6YfB5lRdZPJFJudyTCZzju4YuFyKpTJj6QKPH53iudEkowsWNQ5E/Fw4GMGsi7Czr2NFfntvJIDoZ8TmVMAYBoYXeWu9iGxu9IbLTQOGTT8MtmbUxUoJHrWzpKZTOWQ8iYwlOTAxR6ZgD0V63S529IUxgxEuXBdhXTTQdll26xHyuSsBxEdX0HvaILp+RmxOzZL6B+Ba4M+BycoxF/Drjd5MqdUk6PPMp2CvBo+pufbutuoN+7l+Wy/Xb+ulWCpzbDplBZDxJCdiGZ6fmOP5iTm+9+wYXSGvFTwGI+we6Fgxmy6l8yXS+ex8ayrgddMV9NIV8uENrbwuuHZS1xiGMWYA+F2giLUn97Qxpl9EJs9xqeO0hWHTb0+2VtZFtlBiai7H5Fx7DJjXuw4jmS3w/HiS/eNJnh+fI1HTanIBG7uD7B7oYPdAhO294ZXZfdURIJvO0RX0VqZWe4kG1+ZCwmakBhkEPgjkROTjjd7MCRowbBowbO1SF9XZVlOpPLFUa9Z5LGXhXqlcZiSeYX+l9XF0On1aGnmP28W23hC7B6zWx+bu0IpIPrhYXbhdEA3YAaQz6G2bVCZOcnThnjHGB/wX4K3AfY3eSKm1yOtxz6cmKZbKzKRyTM3lmUrlKLTbQo8abpeLjd0hNnaHuOWCAXKFEkemUxyYSHJgfI7heIZDkykOTab4wXNWt8/OvrAVQAY7GFpB4x+lMsQzhdO6EsM+D9Ggl65KC6Rdkyk2Wz2D3l7gt4GPYm3L+kkRec4Y4xGRs2971gTawrC1y7fqdtDudVEul4mnC0ylrHGP2sHn5eZEapBUrsDBiRQHJq3B84mahIlgfWPfNdDBrv4OdvaH6e/wt0UAWWpdeN1WOhOrJeIlGvSt+K1snZoldRQYBT4JPAvzXbLvE5GPNHrD5aYBw9buD8lmWml1kcwWmE7lmZ5b/hlXzcglFUvnOTAxZ7VAJuaYXTDw3xn0sqMvzM7+Dnb2dzAYaU0AWc666PB7agKIl7DP0xZBsV5OdUkdAe7Fmil1Ddb4F8BVjd5MKbW4SMBLJOBlS0+IXKHEdKoy7pHOL7oFbbvpDvm4dks3127pplwuM57McWAiyaHJFIcnrQDy5PAsTw7PAhDxe9jR3zEfRIY6Aytu4HkuV2QuV5yfjeVxu6wAEvASCVo/V+LEgLOpJ2D8oYjcv/CgMebqem5gjLkVeCMwDpRF5I4F778TeC92yvSvi8jfV977deBKrNlZh0Tkr+u5p1Irmd9rj3uUymXi6bzV+kjl2iot+5m4XC7WRQOsiwa4cUcf5XKZsUSWw1MpDk3OcWgyRSJbYO8pa/MosMYMtveF2dlvBZANXcEVF0CKpTKxtBXkqwJe93wrpPqlYCV3ZTU0S6pRxpgwsBe4RESyxphvY03LvafmnHcC94rI0QXXbgL+A7hSRMrGmEeAt4nIgdrztEvKttK6YZy0WusinS8yNZdjJpUnnqkvQWK7pTcvl8tMzuXmWx+HplKnPWTBSmGytTfM9t4Q2/qsvUKWIwlhq+vCBYT9HiIBazwkGvDSEfC0JDg2YwOlRl0PHBOR6v+hB4DbgHsWnPc+Y8woEAa+JiLTwKuBx0Sk+pF4EPglrC1ilVqTQjWLBavfaGcqrQ8nB86Xk8vlYiASYCAS4LptPZTLZaZTeQ5Nzs23QqZT+fkFhWBNfd3YFWRbr7U/+rbeMF0hX4v/Jo0rY3dlVZMrul1WEKkGkEjAS4e/PcdDnA4Yg0Ci5vVs5Vit+4A7RWSisgXsvwC31HktHo+b7u7wshZ6pdK6sK2Vuuir+XMqV2BqLsd0pQVSHfvwuF1tn+I8EgmyZTDKzZXXM6kch8aTVhfWRJKTM2lOxDKciGW4//A0AH0dfnYOdLBzIMLOgQ42dIVwn6O7p13rogzMFmE2VYBUAbfLGtfqDPqIVmZldfg95/z7Oc3pgDEORGted1aOzRORIzUvfwx81xjjqZy3a8G1BxfeoFgsrcquh6VYrd0wS7FW6yLqgmjEz+YOH7PpAjPpPHmPm9GpuVYXrSF+4KKBDi4a6ICLBskWihybSXN0KsWR6RTHptNMVVKxPHx0Bji9G2tLb5gt3SFC/tM3Y2p1l1QjEsksIzWv3S7o8HuJBKwurYj//LqzBgai5z5pAacDxoPAVmNMoNItdQPwF8aYXqAgIrOVnfs+LiIFrK1fj4pI0RjzQ+D9xhhXpVvqeuDPHC6vUquC2+WiO+yz/usOMxH1M1OZdTWTzpNdId1XVQGvhwsGIlwwEAGsleijs1mOVALI0akUM+nTu7HAysa7tSfElp4wW3tD7FrBuaRKZUhkC5WULVbQqx0T6aj8dHJg3dFBbwBjzCuBNwETQF5E7jDGfA6YFpHPGGM+AFyKNX33MuArIvKLyrW/jjWVtwg8v9gsKR30tq3Vb9WL0bqwLVYX6XxxPoDE0vm2XnVer3g6X2l9pDg2k2Y4lnnB38vncbGxK2QFkd4QW3pC9IR8bTlecD5CPnelNWIHkoWTBhzPJdWONGDY9CFp07qwnasuyuUyc7ni/AD6bKbQdnubL0WhVGIknuXYTIrjM2mOTaeZnMu94LxoZf3LlkoQ2dwdIrQK9xX3eVx0+K0A0hn0cvH2fg0Ya5k+JG1aF7ZG66JULpOo5FaKpfPM1jl9d0Xwetg/HOPYTJrjM2mOT6dJ5V+Y4ai/w8+m7iCbuq0Asqk7eNq+GitdwOvmdddsabtptUqpFcbtclmbD4V8bOkJUSqXmc0UiKfzxDMFEiu4BdIR8M7vhQ5W62pqLsexSgvkRCzNqXjG2hhrLje/Mh2s8RArgFiBZGPX6goi9dCAoZQ6K7fLRXfIR3dl3UNtCySezjObLayI9CWLcblc9EcC9EcCXL25G7BWbI/OWlN4T8bSnIylOTWbZSKZYyKZ44mTcetaoD/in2+BbO4OsbE7SMC7eoOIBgylVENqWyD0hCiXyySzReIZqwUym8mTL67MAALWWo1qanfoAazxkNHZLCdj6flAMhK3g8jjNUGkr8PPxq4gG7qC8z87g95VMbCuAUMpdV5cLldlcZmXTZVjqVyR2UyeeLrAbDa/InJgnY3X7Z5fYX9d5VihWGIkkeXkTJqTsQwnYmlGZu3urKdO2d1ZEb+HDQuCyEAksCI2naqlAUMptezCfg9hv4ehTut1rlCyAkhlDCSZK6z4gXSvx83myqB4VaFYYiyRZTie4VQ8w6lZ62cyV5zfL33+ereL9Z12ALH+C7R1l5YGDKWU4/xe9/xYAdjjIIlsoRJE8uRWcDdWldfjrunOspTLZWbSeU7FM3YgiWeYTuU5EbMG2qtcQG/Yx1BnkPWdgfmf7dIa0YChlGq62nGQajdWJl9kthJEVksrBKwuu96wn96wn0vXd84fT+WKjMzaQWQ4nmFsNstUKs9UKs+zo3YqPY/LxWDUPx9A1ncGGYoG6Ak3d9GhBgylVFsI+jwEfR4Go3YrJJm1A0giW1jxYyG1wn7P/A6EVYVSiYlkjtHZLCOzmfmf06k8I7NZRmazPFHzOwJeN0PRSgDptH9GAs482jVgKKXaktvlojPoozPogy7rWKFYYrbaAqkEk9XQlVXldbtZ3xlkfWeQK6t/aSBbKDI6m7UDScL6mcxaSRmPzaRP+z2RgGd+E6t10QBD0SDrogEigfNLm64BQym1Yng97vnunapsoUQiWyBZaYUkc4UVPa13MQGvh629Ybb2np6yP5EtMFrTEhmdzTKayJLMFklmUxyaPH2Ff9jnYV1ngA2dQV53zZaGy6EBQym1ogW8bgJeP/0ddhDJ5IuVFoj1czUGEbDyYEUHIuyuZPEF5rf1HUtkGUtYAWQskWVsNksqX7Qy/E4tLW2OBgyl1KpTHQ/pt5+jZPJFPEE/wxMJ5rIFktki2eLqGROpcrtc9IT99IT98ylQwJqtNZspMJbILpqEsR4aMJRSa0LQ56E7GsBXtJMN5oslqwWSLTKXs36m80VWX1vEmq1VnZl22RL3R9eAoZRas3wed+XbuH2sWCqTytkBZC5XYC5XXBV7hpwvDRhKKVXD47ZTndTK5IvM5ar/FZhbxa2RM9GAoZRSdaiOi/TZyyYolautkeJ8q2QuV1xxW+DWSwOGUkotkdvlmt9Hu1ahVCaVK5wWTFK5lT/IrgFDKaWWmddds+iwRqFYIpWvBJC8HUgyK6RF4njAMMbcCrwRGAfKInLHgvd/HxgCRoBrgD8Ukf2V944CRyunDovIrzldXqWUcorX46bT435BICmWyqRrAkn1z+l8sa3yaTkaMIwxYeCvgEtEJGuM+bYx5hYRuafmtAjweyJSNsa8Gfg88PrKe98QkU84WUallGo1j3vxrq1yuUy2YLdK0vki6XyJdIu6t5xuYVwPHBORbOX1A8BtwHzAEJGP15zvBpI1r19qjPkIEAXuEpGfO1xepZRqGy6Xa36wfUFWEAqlMpmaQJLJl+ZbJ05NAXY6YAwCiZrXs5VjL2CM8QPvAP5bzeE/EJGHKy2Vx40xrxORg7XXeTxuursX1OQapXVh07qwaV3Y1kpd5Aolu1srVzitiytfLBP0tefCvXGs1kFVZ+XYaSrB4i+B/1dEDlWPi8jDlZ8pY8yTwA3AaQGjWCwRiy0tL8pq090d1rqo0LqwaV3Y1lpdhICQzw0+N4StcZNCsbTk7qylhZn6PQhsNcYEKq9vAO40xvQaYzphfpzjr4EvichjxphfqRy/xRjzmprftQs4hFJKqSXzetx0+JfWVnC0hVFpGfwO8FVjzASwV0TuMcZ8DpgGPgP8A3ApsN0YA9ABfBurJfIJY8xVwAbgX0XkZ06WVyml1Jm5yuU2mrO1BPl8sbyWmphns9aa22ejdWHTurBpXdgGBqIN76TkdJeUUkqpVUIDhlJKqbpowFBKKVUXDRhKKaXqogFDKaVUXTRgKKWUqosGDKWUUnXRgKGUUqouGjCUUkrVRQOGUkqpumjAUEopVRcNGEoppeqiAUMppVRdNGAopZSqiwYMpZRSddGAoZRSqi4aMJRSStVFA4ZSSqm6aMBQSilVFw0YSiml6uJ1+gbGmFuBNwLjQFlE7ljwfhD4AjAM7AY+IyLPV977deBKoAgcEpG/drq8SimlFudoC8MYEwb+CvjvIvIJ4HJjzC0LTrsdOC4inwa+DHy9cu0m4EPAh0TkI8C7jTG7nSyvUkqpM3O6hXE9cExEspXXDwC3AffUnHMb8D8ARORpY8weY0wn8GrgMREpV857EPgl4EDtDXw+j2tgIOrgX2Fl0bqwaV3YtC5sWhdL5/QYxiCQqHk9WzlWzzn1XKuUUqpJnA4Y40BtOO+sHKvnnHquVUop1SROB4wHga3GmEDl9Q3AncaY3kq3E8CdWF1XGGMuA54SkVngh8DVxhhX5bzrgbscLq9SSqkzcJXL5XOfdR6MMa8E3gRMAHkRucMY8zlgWkQ+Y4wJYc2SGgF2AX+yYJbUNVizpJ7XWVJKKdU6jgcMJ51ryu5aYYzZCXwKeBzYBEyJyCdbW6rWqXwJeQj4TxH5UKvL0yrGGAO8FUgDLwM+ISIPt7ZUrWGM+TCwDZjEmr7/LhFJt7RQTWSMGcJ6RuwRkWsrx864pOFMVuzCvTqn7K4VvcA/i8jnReQDwFuMMVe3ulAt9CngiVYXopWMMR7gS8AnReSzwLuAI60tVWtUHpZ/ALxfRP4I6MD6ormW3Ah8B3DVHLudRZY0nM2KDRicecrumiMij4jId2oOuYG5VpWnlYwxb8f6t7AmH441rsV6OLzfGPMHwOuxvl2vRSkghzVxBiACPNu64jSfiHyL02edgvW8fLDy/tPAnpqx5UWt5ICh024XYYz5ZeCHIrK/1WVpNmPMxcBFIvKvrS5LG9iK9aXqG5VvkDcB72htkVqjMonmw8A3jTHfAE4CB1taqPbQ8DN0JQcMnXa7gDHmZuBm4L+3uiwt8stAxhjzUawm+IuMMbe3tkgtMwvsF5F45fXPgJe3rjitY4y5Aitg3CYi78Rqaf1hK8vUJhp+hq7kgLHolN0WlqeljDG3Ya2O/wAwZIy5vsVFajoR+Z8i8kkR+QzWA/JhEfnTFherVR4C+ipjGWC1OM46oLmKbcSalVmovB4Bgi0sT7s405KGM1rps6ReMGW3xUVqicoA933Ao5VDHcCfi8g3WlaoFjLG/Arw3wA/Vj38U4uL1BKV7slXYH0+tmAN+q6ZmUFVlaD5VSADxIBLgdtFZKSV5WomY8zLgN8AXgP8JfDFyluLLmk4kxUdMJRSSjXPSu6SUkop1UQaMJRSStVFA4ZSSqm6aMBQSilVFw0YSiml6qIBQymlVF00YCillKqLBgylVhljzCZjzJtbXQ61+mjAUGqJjDE3GmMeN8a8vPL60ZpUHEv5fbcvU9FuAa5apt+l1DwNGEotkYj8DNhbc+haESmex6+8/fxKZAUxrH0w3mSMedIYs+N8f6dSVZoaRK1axpg/xMonlcXaaexNlWM+rL0ictWdCRc7boz5TeDTWBt17QSMiLzIGPPVyrmHsRI+fgor0+dXgZeLyNGaa/8aK/HfDuB1QAn4JvBTwAD/KCJ3G2P+C/A3wJ9iZZn9Z2PMJwEv1hbFCRH5XGUXvY8CTwNXAn+8MP+PMeYHwIdE5JnlrE+lNGCoVckY82rgAyLy2srrd2Htg/B+EXld5dhdWDuNuRY7LiL/aYy5F/iciHzfGHMNsA54n4j8UuXcnwEfE5F7K+e+U0SOVt67F/i0iPzQGPPnwI+Bu4CXVIJEL9beJdUtM4+KyLaa8n9QRF5V87tux0pRfhHwfqwsrJmFSfSMMUeA3TXZWZVaFtolpVary6nZJEdEvl45drjmnIPAnrMcr3qu8jseBS4BDtS8V3vdYqrf/iew9h5wAS83xnwceA8wcJbyh40xH63s73Gicu7fYu1ZcD9wB5CvvcgY0w/ENVgoJ2jAUKvVU1jdSAAYY35r4TGsje+fPMvxqtpm+D7ggprX5xojWNiEfzewQUT+GGusoVbRGOMyxuyplGlcRD5T2d/j/wMEeDHwGRF5MTCGlbK61jbg1DnKpNSSaJeUWrUq4xIhrH0QpkTkazXHXEC6uofKYscr+638DfBt4LMiMmGMcQFfAwJY3/pfBewHvos1/vBNEflozbV/D3wD+DtgpnLOnwC/AKaBjwDvFpFvV8ZG8gAi8kFjzMew9jZJAD1YYxe/DLwSq2VzEXCHiMzvX26MiWB1fYWB94jIz5etQtWapwFDKaVUXbRLSimlVF00YCillKqLBgyllFJ10YChlFKqLhowlFJK1UUDhlJKqbpowFBKKVWX/x9rIa6MMfe5mQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from problem_formulation import coordinates\n",
    "\n",
    "pyplot.fill_between(coordinates, expected-variance**0.5,\n",
    "                    expected+variance**0.5, alpha=0.3)\n",
    "pyplot.plot(coordinates, expected)\n",
    "\n",
    "pyplot.axis([0, 10, 0, 2])\n",
    "pyplot.xlabel(\"coordinates $t$\")\n",
    "pyplot.ylabel(\"Model evaluations $u$\")\n",
    "\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Avoiding matrix inversion issues\n",
    "\n",
    "It is worth noting that the construction of Lagrange polynomials are not always numerically stable.\n",
    "For example when using grids along , most often the expansion construction fails:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.559212Z",
     "iopub.status.busy": "2021-05-18T10:57:05.558951Z",
     "iopub.status.idle": "2021-05-18T10:57:05.567715Z",
     "shell.execute_reply": "2021-05-18T10:57:05.567970Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Lagrange abscissas resulted in invertible matrix'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes, _ = chaospy.generate_quadrature(1, joint)\n",
    "\n",
    "try:\n",
    "    chaospy.expansion.lagrange(nodes)\n",
    "except numpy.linalg.LinAlgError as err:\n",
    "    error = err.args[0]\n",
    "error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is impossible to avoid the issue entirely, but in the case of structured grid, it is often possible to get around the problem.\n",
    "To do so, the following steps can be taken:\n",
    "\n",
    "* Create multiple Lagrange polynomials, one for each marginal distribution.\n",
    "* Rotate dimensions so each expansion get a dimension corresponding to the marginal it was created for.\n",
    "* Use an outer product of the expansions to combine them into a single expansion.\n",
    "\n",
    "So we begin with making marginal Lagrange polynomials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.570362Z",
     "iopub.status.busy": "2021-05-18T10:57:05.570043Z",
     "iopub.status.idle": "2021-05-18T10:57:05.588691Z",
     "shell.execute_reply": "2021-05-18T10:57:05.588360Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[polynomial([-0.3041*q0+0.9562, 0.3041*q0+0.0438]),\n",
       " polynomial([-10.0*q0+2.0, 10.0*q0-1.0])]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes, _ = list(zip(*[chaospy.generate_quadrature(1, marginal)\n",
    "                      for marginal in joint]))\n",
    "expansions = [chaospy.expansion.lagrange(nodes_)\n",
    "              for nodes_ in nodes]\n",
    "\n",
    "[expansion_.round(4) for expansion_ in expansions]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we rotate the dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.591056Z",
     "iopub.status.busy": "2021-05-18T10:57:05.590739Z",
     "iopub.status.idle": "2021-05-18T10:57:05.631756Z",
     "shell.execute_reply": "2021-05-18T10:57:05.631417Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[polynomial([-0.3041*q0+0.9562, 0.3041*q0+0.0438]),\n",
       " polynomial([-10.0*q1+2.0, 10.0*q1-1.0])]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vars_ = chaospy.variable(len(joint))\n",
    "expansions = [\n",
    "    expans(q0=var)\n",
    "    for expans, var in zip(expansions, vars_)]\n",
    "\n",
    "[expans.round(4) for expans in expansions]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally construct the final expansion using an outer product:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.633900Z",
     "iopub.status.busy": "2021-05-18T10:57:05.633578Z",
     "iopub.status.idle": "2021-05-18T10:57:05.648001Z",
     "shell.execute_reply": "2021-05-18T10:57:05.647708Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([3.04*q0*q1-9.56*q1-0.61*q0+1.91,\n",
       "            -3.04*q0*q1+9.56*q1+0.3*q0-0.96,\n",
       "            -3.04*q0*q1-0.44*q1+0.61*q0+0.09,\n",
       "            3.04*q0*q1+0.44*q1-0.3*q0-0.04])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expansion = chaospy.outer(*expansions).flatten()\n",
    "expansion.round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The end results can be verified to have the properties we are looking for:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:57:05.650395Z",
     "iopub.status.busy": "2021-05-18T10:57:05.650028Z",
     "iopub.status.idle": "2021-05-18T10:57:05.660996Z",
     "shell.execute_reply": "2021-05-18T10:57:05.660655Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1., -0., -0., -0.],\n",
       "       [ 0.,  1.,  0.,  0.],\n",
       "       [ 0.,  0.,  1.,  0.],\n",
       "       [ 0.,  0.,  0.,  1.]])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes, _ = chaospy.generate_quadrature(1, joint)\n",
    "expansion(*nodes).round(8)"
   ]
  }
 ],
 "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
}