{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sobol Sensitivity Indices\n",
    "\n",
    "Saleh Rezaeiravesh, salehr@kth.se <br/>\n",
    "SimEx/FLOW, Engineering Mechanics, KTH Royal Institute of Technology, Stockholm, Sweden\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import math as mt\n",
    "import matplotlib.pyplot as plt\n",
    "from UQit.sobol import sobol\n",
    "import UQit.analyticTestFuncs as analyticTestFuncs\n",
    "from UQit.pce import pce, pceEval\n",
    "import UQit.reshaper as reshaper\n",
    "import UQit.sampling as sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1: Sobol indices for 2 uniform parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following analytical function, \n",
    "$$\n",
    "\\begin{equation}\n",
    "f(\\mathbf{q}) = q_1^2+q_1 q_2\n",
    "\\label{eq:1dModel}\\tag{3}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{q}=\\{q_1,q_2\\}$, and both parameters are random with uniform distributions, i.e., $q_1\\sim \\mathcal{U}[\\mathbb{Q}_1]$ and $q_2\\sim \\mathcal{U}[\\mathbb{Q}_2]$. \n",
    "The aim is to compute the main Sobol indices using `UQit` and validate them with the analytical values..\n",
    "\n",
    "The above model function (simulator) is implemented in `analyticTestFuncs.fEx2D` with keyword `type3`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 1:** Set the number of samples for $q_1$ and $q_2$ and their associated admissible spaces $\\mathbb{Q}_1$ and $\\mathbb{Q}_2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=[100, 90]       #number of samples for q1 and q2\n",
    "qBound=[[-3,1],   #admissible space of q1 and q2\n",
    "        [-1,2]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 2:** Generate samples for parameters $q_1$ and $q_2$. We choose to have samples to be uniformly spaced over the parameter spaces. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "p=len(n)\n",
    "q=[]\n",
    "pdf=[]\n",
    "for i in range(p):\n",
    "    q.append(np.linspace(qBound[i][0],qBound[i][1],n[i]))\n",
    "    pdf.append(np.ones(n[i])/(qBound[i][1]-qBound[i][0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 3:** Run the simulator (3) at the parameter samples to generate training model outputs. We assume that a 2D sample parameter is obtained by tensor product of the samples for $q_1$ and $q_2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "fEx_=analyticTestFuncs.fEx2D(q[0],q[1],'type3','tensorProd')\n",
    "fEx=np.reshape(fEx_.val,n,'F')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 4:** Compute first- and second-order main and also total Sobol sensitivity indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sobol_=sobol(q,fEx,pdf)\n",
    "Si=sobol_.Si\n",
    "STi=sobol_.STi\n",
    "Sij=sobol_.Sij"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 5:** Validate the computed Sobol indices through comparison with the analytical values. To this end, we have derived expressions to compute the exact Sobol indices for the above model function. Note that, in general since the exact form of the model function (simulator) is not known, exact Sobol indices can not be evaluated. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "fEx_.sobol(qBound)\n",
    "Si_ex=fEx_.Si\n",
    "STi_ex=fEx_.STi\n",
    "Sij_ex=fEx_.Sij"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Results:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Main indices by UQit:\n",
      "\t S1=0.716474, S2=0.121511, S12=0.162015\n",
      "Analytical main Indices:\n",
      "\t S1=0.716472, S2=0.121512, S12=0.162016\n"
     ]
    }
   ],
   "source": [
    "print('Main indices by UQit:\\n\\t S1=%g, S2=%g, S12=%g' %(Si[0],Si[1],Sij[0]))\n",
    "print('Analytical main Indices:\\n\\t S1=%g, S2=%g, S12=%g' %(Si_ex[0],Si_ex[1],Sij_ex[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Extra:** In the above example, all the required samples to compute the Sobol indices were drawn directly from the simulator. This can be the case when running the simulator is computationally inexpensive. But in practice, we need to deal with the cases where the computational cost of running the simulator limits the number of training samples. \n",
    "In such situations, we need to construct a surrogate based on a limited number of training data and then use the surrogate to compute the Sobol indices. \n",
    "\n",
    "Following this procedure, we repeat the previous example. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 1:** Draw a few samples to construct the PCE surrogate. We choose to have the parameter samples to be Gauss-Legendre points. The simulator is run at these parameter samples. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "nQpce=[5,6]    #number of samples from q1, q2 to construct surrogate\n",
    "xi=[]\n",
    "qpce=[]\n",
    "for i in range(p):\n",
    "    samps=sampling.trainSample(sampleType='GQ',GQdistType='Unif',qInfo=qBound[i],nSamp=nQpce[i])\n",
    "    xi.append(samps.xi)\n",
    "    qpce.append(samps.q)\n",
    "fVal_pceCnstrct=analyticTestFuncs.fEx2D(qpce[0],qpce[1],'type3','tensorProd').val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 2:** Construct the PCE surrogate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "... A gPCE for a 2-D parameter space is constructed.\n",
      "...... Samples in each direction are Gauss Quadrature nodes (User should check this!).\n",
      "...... PCE truncation method: TP\n",
      "...... Method of computing PCE coefficients: Projection\n",
      "...... Number of terms in PCE, K=  30\n",
      "...... Number of Data point, n=  30\n"
     ]
    }
   ],
   "source": [
    "xiGrid=reshaper.vecs2grid(xi)\n",
    "pceDict={'p':2,'sampleType':'GQ','truncMethod':'TP','pceSolveMethod':'Projection',\n",
    "             'distType':['Unif']*p}\n",
    "pce_=pce(fVal=fVal_pceCnstrct,nQList=nQpce,xi=xiGrid,pceDict=pceDict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 3:** Run the surrogate at test samples which are used to compute the Sobol indices. The test samples are taken from $\\mathbb{Q}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "qpceTest=[]\n",
    "xiTest=[]\n",
    "for i in range(p):\n",
    "    testSamps=sampling.testSample('unifSpaced',GQdistType='Unif',qBound=qBound[i],nSamp=n[i])\n",
    "    xiTest.append(testSamps.xi)\n",
    "    qpceTest.append(testSamps.q)\n",
    "fPCETest_=pceEval(coefs=pce_.coefs,kSet=pce_.kSet,xi=xiTest,distType=['Unif']*2)\n",
    "fPCETest=fPCETest_.pceVal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 4:** Compute the Sobol indices using the surrogate values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sobolPCE_=sobol(qpceTest,fPCETest,pdf)\n",
    "Si_pce=sobolPCE_.Si\n",
    "Sij_pce=sobolPCE_.Sij"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can validate the computed indices by comparing them to the values resulted from the previous example (where actual simulator is used). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Indices by the surrogate:\n",
      "\t S1=0.716474, S2=0.121511, S12=0.162015\n"
     ]
    }
   ],
   "source": [
    "print('Indices by the surrogate:\\n\\t S1=%g, S2=%g, S12=%g' %(Si_pce[0],Si_pce[1],Sij_pce[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2: Sobol indices for 2 parameters with Gaussian distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do Example 15.7 in [Smith](https://rsmith.math.ncsu.edu/UQ_TIA/).\n",
    "\n",
    "$$\n",
    "f(\\mathbf{q})=c_1 q_1 +c_2 q_2\n",
    "$$\n",
    "\n",
    "where $c_1$ and $c_2$ are constant parameters, and $q_i\\sim\\mathcal{N}(0,\\sigma_i^2)$ for $i=1,2$. \n",
    "For different values of $c_1,\\,c_2,\\,\\sigma_1$ and $\\sigma_2$, we can compare the Sobol indices computed by `UQit` and associted exact values. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Settings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=[101, 100]        #number of samples for q1 and q2\n",
    "qBound=[[-20,20],   #admissible range of parameters\n",
    "        [-20,20]]\n",
    "sig=[1.,3.]\n",
    "c=[2,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate uniformly-spaced values for the parameters; construct and plot the PDFs (probability density functions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3iU1bX48e/K/UISAgQEAhIURFAO0IhW8UKlCtqKWq3aY8XLEW3ltD2ttvbYg9bW1lpr/Z3WVqmiqFXUVm1qUdQqWo83giByVUCUEC6R3CDXSWb9/njfCWOYJDOTmcybZH2eJ48z721WxmHNzt77XVtUFWOMMX1XUqIDMMYYE1+W6I0xpo+zRG+MMX2cJXpjjOnjLNEbY0wfl5LoANobMmSIjhkzJtFhGGNMr7Jq1arPVLUg1D7PJfoxY8ZQWlqa6DCMMaZXEZFPOtpnXTfGGNPHWaI3xpg+zhK9Mcb0cZ7rozfG9G0+n4+ysjIaGxsTHUqvlJGRQWFhIampqWGfY4neGNOjysrKyMnJYcyYMYhIosPpVVSVffv2UVZWRlFRUdjnhdV1IyKzRWSziGwRkRs7Oe4CEVERKQ7a9mP3vM0icmbYkRlj+qTGxkYGDx5sST4KIsLgwYMj/muoyxa9iCQD9wBfBsqAlSJSoqob2h2XA3wHeCdo20TgYmASMAJ4WUTGq2prRFEaY/oUS/LRi+a9C6dFPx3YoqrbVLUZWArMDXHcz4A7gOCvmrnAUlVtUtWPgS3u9YzpdVZs3su2igOJDsOYiIWT6EcCO4Kel7nb2ojIVGCUqj4X6bnu+fNFpFRESisqKsIK3Jie1NLq51uPvsdvXvww0aGYHrJixQq+8pWvANDU1MSsWbOYMmUKTzzxRFjnb9q0iSlTpjB16lS2bt0a1jmvv/4606ZNIyUlhb/85S9Rx95eOIk+1N8JbauViEgS8FvgB5Ge27ZBdZGqFqtqcUFByDt4jUmorRV1NPhaWbuzOtGhmARYvXo1Pp+PNWvWcNFFF4V1zrPPPsvcuXNZvXo1RxxxRFjnjB49moceeohvfOMb3Qn3EOEk+jJgVNDzQqA86HkOcAywQkS2AycAJe6AbFfnGtMrfLCzBoAdlQ1U1zcnOBrTHdu3b2fChAnMmzePyZMnc8EFF1BfXw/ACy+8wIQJE5gxYwZPP/00AHv37uXSSy9lzZo1TJky5ZDW+Zo1azjhhBOYPHky5513HlVVVSxbtoy7776b+++/n5kzZx4Sw4MPPsj48eM59dRTufrqq1mwYAHglICZPHkySUmxvcUpnOmVK4FxIlIE7MQZXG37ulHVGmBI4LmIrACuV9VSEWkAHhORu3AGY8cB78YufGN6xjo30TuPa5kxbkgnR5tw/fTv69lQXhvTa04ckcvNX53U6TGbN2/mgQce4KSTTuLKK6/kD3/4AwsWLODqq6/mlVde4cgjj2xruQ8dOpT777+fO++8k+eea987DZdddhm/+93vOPXUU1m4cCE//elPufvuu7n22msZMGAA119//eeO37VrFzfffDOrVq0iLy+PmTNnMnXq1Ni9ASF0+bWhqi3AAmA5sBF4UlXXi8itInJOF+euB54ENgAvANfZjBvTG60tq2b8sAHOY+u+6fVGjRrFSSedBMCll17KG2+8waZNmygqKmLcuHGICJdeemmX16mpqaG6uppTTz0VgHnz5vH66693es4777zDaaedRkFBAWlpaWF3BXVHWDdMqeoyYFm7bQs7OPa0ds9vA26LMj5jEq6l1c+GXbVcMn00jT7/51r3pnu6annHS/spioHnPTXts6enl1qtG2O6sLWijkafn2NH5nHsyLy2/nrTe3366ae89dZbADz++OPMmDGDCRMm8PHHH7f1wT/++ONdXicvL4/8/Hz+9a9/AfDII4+0te47cvzxx7NixQr27duHz+fjqaee6uZv0zVL9MZ0IZDYJxfmcczIPBuQ7QOOPvpolixZwuTJk6msrORb3/oWGRkZLFq0iLPPPpsZM2Zw+OGHh3WtJUuWcMMNNzB58mTWrFnDwoUhOzvaDB8+nFtuuYUvfvGLzJo1i2nTprXtW7lyJYWFhTz11FNcc801TJoUm794rNaNMV34oKyarLRkioYM4NiRTc62nTWcPM6mAvdWSUlJ3HvvvYdsnz17Nps2bTpk+2mnncZpp50W8lpTpkzh7bffPmT7Lbfc0uHrX3HFFVxxxRUAPPTQQ22LLR133HGUlZWF8RtExlr0xnThg501TBqRS3KScOzIvLZtxvQWluiN6URgIPYYN8HnZaUyelCWDcj2YmPGjGHdunWJDqPN5Zdfzu9///u4voYlemM6ERiInVyY17bNBmRNb2OJ3phOrC1z5swHumyAtgHZqjobkDW9gyV6YzqxbmdN20BsQCDpryu3Vr3pHSzRG9OJD3bWcMyIPJKTDt7gYgOyprexRG9MJz7cc4Cjh+d8blteViojB2ayeff+BEVlekIiyhTfddddTJw4kcmTJ3P66afzySefRB1/MEv0xnSg0dfKgaYWhuZmHLJvaG46ldZH32/0VJniqVOnUlpaytq1a7ngggv44Q9/2J2w21iiN6YDVe7dr4Oy0w7ZNygrzRJ9L+XlMsUzZ84kKysLgBNOOCFmN0/ZnbHGdCCQyPOzDk30+dlpbNwV2/K6/dLzN8LuD2J7zcOOhTm3d3pIbyhT/MADDzBnzpxuvBEHWYvemA5U1fmADlr02WlUWr2bXsvrZYofffRRSktLueGGG6L47Q5lLXpjOlDZ1nWTesi+/Kw0Gn1+GppbyUxL7unQ+o4uWt7x4uUyxS+//DK33XYbr732Gunp6TF5vbBa9CIyW0Q2i8gWEbkxxP5rReQDEVkjIm+IyER3+xgRaXC3rxGRQ6sIGeNRVZ103QSSv7XqeyevlilevXo111xzDSUlJQwdOjTaX+8QXbboRSQZuAf4Ms4asCtFpERVNwQd9piq3usefw5wFzDb3bdVVafELGJjekhlXTMikJcZukUPzpfByIGZPR2a6aZAmeJrrrmGcePGHVKmeMiQIcyYMSOsmjhLlizh2muvpb6+nrFjx/Lggw92enxwmeLhw4czbdo0WludhfduuOEGDhw4wIUXXgg4i4WXlJR0+/cNp+tmOrBFVbcBiMhSYC7O8oAAqGrwqFQ2oN2OzJgEq6pvJi8zlZTkQ//wDfTb28yb3smrZYpffvnlMKKPXDhdNyOBHUHPy9xtnyMi14nIVuAO4DtBu4pEZLWIvCYiJ4d6ARGZLyKlIlJaUVERQfjGxE9lXTODQnTbgDPrBg5OwTTGy8JJ9KFGDQ5psavqPap6BPAj4Cfu5l3AaFWdCnwfeExEckOcu0hVi1W1uKDAFnMw3lBV39yW0NsLfAFYi773sTLFoZUBo4KeFwLlnRy/FDgXQFWbVHWf+3gVsBUYH12oxvSsyjpfyIFYgNzMVJIEq2AZJVXr3Y1WNO9dOIl+JTBORIpEJA24GPjc6ICIjAt6ejbwkbu9wB3MRUTGAuOAbRFHaUwCVNU1h5xaCZCcJAzMsrn00cjIyGDfvn2W7KOgquzbt4+MjEPLcnSmy8FYVW0RkQXAciAZWKyq60XkVqBUVUuABSIyC/ABVcA89/RTgFtFpAVoBa5V1cqIIjQmAVSVyk66bgDys1Lbbqoy4SssLKSsrAwbj4tORkYGhYWFEZ0T1g1TqroMWNZu28Kgx9/t4Ly/An+NKCJjPKC+uZXmFn+Hg7Hg3h1rXTcRS01NpaioKNFh9CtWAsGYENrq3HTaok+zWTemV7BEb0wIbZUrrUVv+gBL9MaEEFaLPttp0dugovE6S/TGhNBZLfqAQVlp+FqVA00tPRWWMVGxRG9MCJWBEsWddN203R1rM2+Mx1miNyaEqrpmkpOEnIyOJ6ZZBUvTW1iiNyaEyvpm8rNSSUrquG54cAVLY7zMEr0xIVTVNXdY/iDAKlia3sISvTEhVNZ1flcsWAVL03tYojcmhKr6jksUB+Skp5CSJNaiN55nid6YECrrfF226EWkbS69MV5mid6YdlTVadF3ULky2KAsuzvWeJ8lemPaqW1sodWvXQ7GAuRnWwVL432W6I1pJzBdsrO7YgMGZVtNeuN9luiNaSeQuLvqowe3gqV13RiPs0RvTDttLfowum4GuYOxfr8VNjPeFVaiF5HZIrJZRLaIyI0h9l8rIh+IyBoReUNEJgbt+7F73mYROTOWwRsTD5URdN3kZ6XhV6httH56411dJnp3zdd7gDnAROCS4ETuekxVj1XVKcAdwF3uuRNx1pidBMwG/hBYQ9YYr6qKoOvG7o41vUE4LfrpwBZV3aaqzcBSYG7wAapaG/Q0Gwj8HTsXWKqqTar6MbDFvZ4xnlVZ5yMtOYnstK7bJHZ3rOkNwlkzdiSwI+h5GXB8+4NE5Drg+0Aa8KWgc99ud+7IEOfOB+YDjB49Opy4jYmbqrpm8rNTEem4oFlAoB+/0qZYGg8Lp0Uf6tN+yMiTqt6jqkcAPwJ+EuG5i1S1WFWLCwoKwgjJmPhxKld23W0Dzjx6sAqWxtvCSfRlwKig54VAeSfHLwXOjfJcYxKuqq45rIFYCOqjt64b42HhJPqVwDgRKRKRNJzB1ZLgA0RkXNDTs4GP3MclwMUiki4iRcA44N3uh21M/FTWd125MiAzNZn0lCRr0RtP67KPXlVbRGQBsBxIBhar6noRuRUoVdUSYIGIzAJ8QBUwzz13vYg8CWwAWoDrVLU1Tr+LMTFRVdd15coAEXHujrVEbzwsnMFYVHUZsKzdtoVBj7/bybm3AbdFG6AxPcnvV6obfORndV3QLGBgllWwNN5md8YaE2R/UwuqkJsZfqLPy0yhtqEljlEZ0z2W6I0JUtvgTJOMJNHnZqRS02DTK413WaI3JkggYedF1KJPtRIIxtMs0RsTJJCwczMiaNFnWoveeJslemOC1EbZoq9vbsXX6o9XWMZ0iyV6Y4IEBlVzM8OakOYcm5HinmuteuNNluiNCRJVH707FbO20WbeGG+yRG9MkNpGH0kC2WmRtOjdRG8teuNRluiNCVLT4CM3M5WkpK4rVwYEWv82IGu8yhK9MUFqG3wRzbiBg3PubYql8SpL9MYEqWnwRdQ/D9aiN95nid6YILWNLRHNuIHgPnobjDXeZInemCA1UXTdZKQmkZacZC1641mW6I0JUhtF142IkJuZYn30xrMs0RsTJDDrJlJWBsF4mSV6Y1yNvlaaWvwRt+jB6ae3efTGq8JK9CIyW0Q2i8gWEbkxxP7vi8gGEVkrIv8UkcOD9rWKyBr3p6T9ucZ4xX73ztZASYNI5GVaojfe1WWiF5Fk4B5gDjARuEREJrY7bDVQrKqTgb8AdwTta1DVKe7POTGK25iYq4miFn1AbmaqlUAwnhVOi346sEVVt6lqM7AUmBt8gKq+qqr17tO3gcLYhmlM/LWVKI4i0edlplgfvfGscBL9SGBH0PMyd1tHrgKeD3qeISKlIvK2iJwb6gQRme8eU1pRURFGSMbEXjQFzQICffSqGuuwjOm2cDojQxX9CPlpFpFLgWLg1KDNo1W1XETGAq+IyAequvVzF1NdBCwCKC4utn8pJiHalhGMcB49OF8OLX6lvrmV7PTI+/iNiadwWvRlwKig54VAefuDRGQWcBNwjqo2Bbararn7323ACmBqN+I1Jm6iWXQkwOrdGC8LJ9GvBMaJSJGIpAEXA5+bPSMiU4H7cJL83qDt+SKS7j4eApwEbIhV8MbEUmAwNdISCGD1boy3dfmJVtUWEVkALAeSgcWqul5EbgVKVbUE+DUwAHhKRAA+dWfYHA3cJyJ+nC+V21XVEr3xpJoGHxmpSaSnJEd8rtW7MV4WVtNFVZcBy9ptWxj0eFYH570JHNudAI3pKdGUKA6wFr3xMrsz1hhXtOUP4GB3j900ZbzIEr0xrtrGyAuaBQT+ErAWvfEiS/TGuJwSxdFNjcxxz7NZN8aLLNEb46ptaIm6RZ+SnMSAdLs71niTJXpjXLWN0ffRQ6Cwmc26Md5jid4YwO/XqBYdCZaTYS16402W6I0B6ppb8Gt05Q8C8jJTrY/eeJIlemPoXkGzgFyrSW88yhK9MRy8ozWa8gcBtviI8SpL9MbQvUVHAnIzbN1Y402W6I0haNGRbvbR1zW30tLqj1VYxsSEJXpjiFUfvdPts9+WFDQeY4neGIIWHenmPHqwMgjGeyzRG4OT6EUgpxurQ7WVKrYplsZjLNEbg7PoSE56CklJoVbODE9elrXojTdZojeG7pUoDrDFR4xXhZXoRWS2iGwWkS0icmOI/d8XkQ0islZE/ikihwftmyciH7k/82IZvDGx0t3yB2B99Ma7ukz0IpIM3APMASYCl4jIxHaHrQaKVXUy8BfgDvfcQcDNwPHAdOBmEcmPXfjGxEZtY/SrSwW0LT5iffTGY8Jp0U8HtqjqNlVtBpYCc4MPUNVXVbXeffo2UOg+PhN4SVUrVbUKeAmYHZvQjYmdmhi06DNTk0lNFmvRG88JJ9GPBHYEPS9zt3XkKuD5SM4VkfkiUioipRUVFWGEZExs1Ta0dKv8AYCIkJthZRCM94ST6ENNQ9CQB4pcChQDv47kXFVdpKrFqlpcUFAQRkjGxFYsWvTg9NNbi954TTiJvgwYFfS8EChvf5CIzAJuAs5R1aZIzjUmkZpb/DT4WrvdRw+QY4neeFA4iX4lME5EikQkDbgYKAk+QESmAvfhJPm9QbuWA2eISL47CHuGu80Yz6huaAZgYFb3E/1AS/TGg7rslFTVFhFZgJOgk4HFqrpeRG4FSlW1BKerZgDwlIgAfKqq56hqpYj8DOfLAuBWVa2My29iTJSq653EPDArrdvXys9KZdtnB7p9HWNiKazRJ1VdBixrt21h0ONZnZy7GFgcbYDGxFtVndOiz49Boh+YlUZ1nbXojbfYnbGm36tqa9F3v+smPyuN/U0t+KxUsfEQS/Sm36uud1v02THouslOda9prXrjHZboTb9X7Q6e5sdiMNbt/qlxB3iN8QJL9Kbfq6pvJi0liczU5G5fa6A7F7/KWvTGQyzRm36vus5HflYq7oyxbgkM6AYGeI3xAkv0pt+rqm9mYGb3++fh4ICu9dEbL7FEb/q96npfTGbcwMEB3ap6a9Eb77BEb/q9qvrmmMyhB8hOcypYWh+98RJL9Kbfq6r3tU2L7C4RcW6asha98RBL9KZfU1Wq65tjUv4gID8r1bpujKdYojf92oGmFlr8GpM59AEDs9Ks68Z4iiV606/FsqBZwMDMVGos0RsPsURv+rVAoo/VYGzgWtZ1Y7zEEr3p1wIJOaZdN9mpVNf7UA25EJsxPc4SvenXAok+VvPowWnRN7f6qW9ujdk1jemOsBK9iMwWkc0iskVEbgyx/xQReU9EWkTkgnb7WkVkjftT0v5cYxIpHn30gb8OrPvGeEWXC4+ISDJwD/BlnDVgV4pIiapuCDrsU+By4PoQl2hQ1SkxiNWYmGtr0cdgYfCAwJdGdb2PwvyYXdaYqIWzwtR0YIuqbgMQkaXAXKAt0avqdnefrbZgepXqeh85GSmkJMeuF7OtsJm16I1HhPPpHgnsCHpe5m4LV4aIlIrI2yJybkTRGRNnsSx/EHCw68amWBpvCKdFH6p2ayTTCUararmIjAVeEZEPVHXr515AZD4wH2D06NERXNqY7qmq98V0xg1AXlsFS2vRG28Ip0VfBowKel4IlIf7Aqpa7v53G7ACmBrimEWqWqyqxQUFBeFe2phui3X5A6Ct5HGVLRJuPCKcRL8SGCciRSKSBlwMhDV7RkTyRSTdfTwEOImgvn1jEq06Di36tJQkBqSnUG3LCRqP6DLRq2oLsABYDmwEnlTV9SJyq4icAyAix4lIGXAhcJ+IrHdPPxooFZH3gVeB29vN1jEmoari0KIHZ16+LT5ivCKcPnpUdRmwrN22hUGPV+J06bQ/703g2G7GaExctLT62d/YEvPBWLAyCMZb7M5Y029VNwRulopt103gmjbrxniFJXrTb1XHofxBQL4tPmI8JKyuG2P6oqpwKle2NMNHy2H7G6Du/YCpWTBxLoyYChJq9rG7+EidJXrjDZboTb8VSMQhE31tObx1D7y/FOo/g9RsSHGPazoA/3c3DDsGpl0GX7ji4D5XXlYatY0ttLT6Y3rXrTHRsERv+q2DBc3add18+CI8cw001cJRc2DqZXDElyDZ/efSUA3r/gqrH4Hnfwhrn4QLFkP+4W2XCEzZrGnwMXhAeo/8PsZ0xJoapt9qq0Wf7bbGW33w4v/AYxdC7kj49jtw0aMw/oyDSR4gcyAcdxXMXwEXLoHPPoT7ToaNz7UdcrDejQ3ImsSzRG/6reoGH6nJQnZaspPkn7gU3vxfpyvmP16CIUd2fZFJ58I1r0F+ETzx7/Dun4CDfyXU2E1TxgOs68b0W4HyB6IKz34LPnwBzroTpl8d2YUGjYWrXoSnLodl10NGHvn5ZwBWBsF4g7XoTb9VVecjPzMFnr8BPngKTl8YeZIPSEmHCx6EMSfDM9cyYs8K5zVsiqXxAEv0pt+qqm/mKv9fYOX9cOJ3YMb3u3fB1Ay45HEY/m8MeX4+0+RDK4NgPMESvem3impXcuGBR2HyxfDlWzucEx+R9By49K+QO5zfp/2OhpqK7l/TmG6yRG/6p/17uKHuTirSD4ev3BWbJB+QNQi58CEKpIZZH/4UNJLlG4yJPUv0pv/xt6JP/wdZ2kDJ+F9AWnbsX2PEVO5Lv5yJ+/8P3v5D7K9vTAQs0Zv+5193IR+/zsKWy2kdMiFuL7Mi73zezTgRXroZdq6K2+sY0xVL9KZ/2b0OXrudA+Pm8lTrqQyJ412rQ3IyuC3lOhgwFJ69zqmbY0wCWKI3/Ye/FUr+EzIGsmnq/wDC8LyMuL3cYXkZbN2fCl/5LVRshDfuittrGdOZsBK9iMwWkc0iskVEbgyx/xQReU9EWkTkgnb75onIR+7PvFgFbkzE3v4jlL8Hc35FWVMW4CTjeBmel8GBphb2j/4SHHshvH4n7N0Yt9czpiNdJnoRSQbuAeYAE4FLRGRiu8M+BS4HHmt37iDgZuB4YDpws4jkdz9sYyJUuQ1e+TmMnwPHfI1dNY0AHJYb3xY9wJ7aRph9uzP18m8LnL8sjOlB4bTopwNbVHWbqjYDS4G5wQeo6nZVXQv42517JvCSqlaqahXwEjA7BnEbEz5VeO6/ICkFzv4NiLCntpGcjBSy0+NXBSTwJbKrphGyh8CcX8HOUucGLWN6UDiJfiSwI+h5mbstHGGdKyLzRaRUREorKuwGExNjG56FbStg1s2Q53z8dtU0xLU1Dwdb9IG/Hjj2Qqfc8Ss/hwN74/raxgQLJ9GHupMk3DtAwjpXVReparGqFhcUFIR5aWPC0FwHy2+Cw46F4ivbNu+ubYpr/zzAMPeLZE8g0YvAnDvA1wAv3xLX1zYmWDiJvgwYFfS8ECgP8/rdOdeY7vvXb6B2p1OVMim5bfPumoa4zrgByEhNZlB2GrtqGw9uHDIOvngdrPkz7Hg3rq9vTEA4iX4lME5EikQkDbgYKAnz+suBM0Qk3x2EPcPdZkz87dsKb/7OqWUz+oS2zS2tfir2N8W96wacVn1biz7glBsgZ4RT0tgGZk0P6DLRq2oLsAAnQW8EnlTV9SJyq4icAyAix4lIGXAhcJ+IrHfPrQR+hvNlsRK41d1mTPw9/yNITncKlgWpONCEX+GwvMy4hzA8L+NgH31A+gA48+ew631Y9VDcYzAmrCkHqroMWNZu28KgxytxumVCnbsYWNyNGI2J3IcvwpaX4IyfQ86wz+0KJN54d92AMyC7Zkf1oTsmnQ+lDzoDs8ecD5k269jEj90Za/qeVh8s/28YfCRMv+aQ3YGulGE90HVzWG4GlXXNNPraddGIwOxfQmM1vHZH3OMw/ZsletP3vPsn2PcRnHEbpKQdsrunW/QAe2ubQuw8FqZdBu8ugooP4x6L6b8s0Zu+pW4fvHa7M199/JkhD9ld20haSlLbAt7xFBjw3V3bGPqAmT+B1Cx48aa4x2L6L0v0pm959TZoOgBn/rLDxUR21zQyPC8DieViIx0Y3nbTVEPoAwYUOLNwPnoRPnop7vGY/skSvek79qyHVQ/CcVfB0I7rzO+uaeyRqZXQrt5NR46/FgaNdcYVWm2NWRN7luhN36AKL9wIGXlw2o87PXRXbUPc74oNyMlIJTst+dAplsFS0uDMX8BnH1odHBMXluhN37DpH/Dx63Daf0PWoA4PU1X21MS//EGww/Iy2N1ZogcYPxvGzoQVv3TGGYyJIUv0pvdraYIXfwIFEz5XzyaUyrpmmlv9PdZ1A26i76zrBg5Ot2w6ACt+0TOBmX7DEr3p/d7+I1R97CTK5M7vAQwk3J6YWhlwWG5m1y16gKFHO+MLpYthz4b4B2b6DUv0pnfbv9tZuWn8HGdKZRcCCbcnyh8EDM/LYO/+Jlr9YRR9Pe3HkJ4LL/zIGXcwJgYs0Zve7aWF0NoEZ94W1uE9sbJUe8PyMmj1K58dCHHTVHtZg+BLP3HGGzY8G//gTL9gid70Xp+8CWufgBO/A4OPCOuUPbWNJCcJBTnpcQ7uoOG57RYg6Urxlc5ds8tvcvrsjekmS/Smd2ptgX9cD3mj4OQfhH3arppGCgakk5wU/5ulAgIzfMLqpwenbv5Zbh39f90Zx8hMf2GJ3vROK++Hveud+edpWWGftrumsUenVkJwou/g7thQRh8P//YNePP38NlHcYrM9BeW6E3vs38PvPoLZ9750V+N6NTdtY09OuMGYFBWGmnJSewOVdisM1/+KaRmwrIbbGDWdIsletP7PP9DaGmEs37dYT2bjuyuaeyR8sTBkpKEobnpkbXoAQYMhS/9D2x7FdY+GZ/gTL8QVqIXkdkisllEtojIjSH2p4vIE+7+d0RkjLt9jIg0iMga9+fe2IZv+p3NzzuzUU69wVl/NQL7G30caGrp8RY9dLDSVDiOuwoKj4PlP7Y7Zk3Uukz0IpIM3APMASYCl4jIxHaHXQVUqeqRwG+BXwXt26qqU9yfa2MUt2eeD00AAA/gSURBVOmPGmvhHz+AoRPhxO9GfPr2z+oBGDUo/D79WBmVn8XHn9VFfmJSMnz1f53fffl/xz4w0y+E06KfDmxR1W2q2gwsBea2O2YusMR9/BfgdOmJGrCmf3nlZ1BbDuf8LuSCIl1ZV14DwKQRubGOrEsTR+Syd38Te/dH0aofNhFm/BesXQpb/hn74EyfF06iHwnsCHpe5m4LeYy7mHgNMNjdVyQiq0XkNRE5OdQLiMh8ESkVkdKKioqIfgHTT3zyprNy1PHXQGFxVJdYt7OGnIwURiegRX/MyDwA1pfXRneBk38Ag8fB37/ntO6NiUA4iT5Uy7z9FICOjtkFjFbVqcD3gcdE5JDmlKouUtViVS0uKCgIIyTTrzTWwNPXQP4YZ3AySuvLa5k0IrdHFhxpb6L7V8SGaBN9agac+weoLYPnfxTDyEx/EE6iLwNGBT0vBMo7OkZEUoA8oFJVm1R1H4CqrgK2AuO7G7TpZ57/kXPz0Pl/gvQBUV2ipdXPxl21HDMiL8bBhSc3I5Uxg7NYt7Mm+ouMmg4nXw/vPwbrrTyCCV84iX4lME5EikQkDbgYKGl3TAkwz318AfCKqqqIFLiDuYjIWGAcsC02oZt+Yf0z8P7jcMr1MOq4qC+ztaKOphY/k0b2fP98wKQReW3jBFE79YcwYho89z1nvMKYMHSZ6N0+9wXAcmAj8KSqrheRW0XkHPewB4DBIrIFp4smMAXzFGCtiLyPM0h7rapWxvqXMH1UTZnTJz1imrOuajcEWtKJatEDTBqZy47KBmrqu7FcYHKq85dNSxM8+y3wt8YuQNNndV6826Wqy4Bl7bYtDHrcCFwY4ry/An/tZoymP/I1whPfBH+Lk9iSU7t1ufXltWSkJjG2ILqun1gIfMms31XDiUcMif5CQ46E2bfD37/jrEj1pZ/EKELTV9mdscabnr8Byt+D8+51Els3rSuvYeLw3B4tZtZeYFrn+p0xmDUz7TKYeim8/mvY+Fz3r2f6NEv0xntKH4T3HnamFEZYyyYUv1/ZUF7bNsUxUQYPSGdEXkb3++nBKf1w1m9gxFR45lorfGY6ZYneeMsnbzlFvI44HWbeFJtLVtZzoKklITdKtTdxRF73Zt4ES82Arz/i3Dy29BtQb8NfJjRL9MY7dq+Dxy6C/MPha/c7t//HwPq2O2IT26IHOGZkLts+q6O+uSU2Fxw4Cr7+MFRth8cvhub62FzX9CmW6I03VG2HR78GadnwzWecJfViZN3OWlKThfHDcmJ2zWgdMyIPVdi4K4Z3t46Z4QxY73gXnpoHrd2Y1WP6JEv0JvH274FHznNKD3/zaRg4OqaXX19ew1GH5ZCWkviPe2CcYF0sBmSDTToXvvJb+OhF+Nt1Nu3SfE5Y0yuNiZuq7U6S378bLvsbDD06ppdXVdbtrOGMiYfF9LrRGpabzuDstNj10wcrvgLq9znF31qb4bxFURV/M32PJXqTOHvWwyPnOy35y/7m3OIfY59W1lNV70voHbHBRIRJI/NYs6M6Pi9wyvWQnAYv/Y9TI+jrj0RdNsL0HYn/W9b0Tx//Cx6c40wTvPKFuCR5gCdLd5AkMPOooXG5fjRmHT2Uj/YeYPWnVfF5gZO+A3PvgW0r4OFznL+WTL9mid70LL8fXr/TSUDZQ+HK5THvrglo9LXy+Ls7OP3oYQlZbKQj508rZEB6Cg+/9Un8XmTqpXDRo7B3I9x7Mnz8evxey3ieJXrTc+r2wWNfd/qQJ50H8191plLGyT/W7qKyrpl5XxwTt9eIxoD0FC74QiHPrS2nYn+EC4ZHYsLZcPUrkDkQHp4Lr/3aBmn7KUv0Jv78fudO198Xw8evwdm/ga89AOnxm+6oqix5aztHFGRz0pGDuzy+p33zi4fja1Uef/fT+L7Q0KPh6lfhmK/Bqz+HP82EslXxfU3jOZboTXztfA8enA0l/wkFR8H81+C4/3D65uNozY5q1pbVMO/EMQlZaKQrRxQM4ORxQ/jzO5/ga/XH98XSBzjz7C9Y7Exlvf90eO6/4MDe+L6u8QxL9CY+PnnLmVHzp5mwbwvM/QNc8byz/mkPWPLmdgakp3D+tMIeeb1oXH7iGPbUNrF8fQ8Mloo4rfoFK+GEb8OqJXD3sbDsh045aNOn2fRKEzv1lbD+aVjzGOxcBVlD4PSFTgs+o+fKD2zZu59/fLCLfz/+cAake/cjftpRQxk9KIv7XtvG6ROGkZkWm5IPncrIhdm/gOIr4Y3fQukDULoYJpwF//YNOPL0bpeENt4jqu2Xf02s4uJiLS0tTXQYJlxV22HLP2HLy85PazMMnQRfmAdTvwlpPTvb5b1Pq7jyoZWkJAnPfPskT822CaXk/XK+u3Q100bn88C8YgZm9fANTtU74O0/wtonoP4zyC6Ao85yEn7Rqc5ArukVRGSVqhaH3BdOoheR2cD/A5KB+1X19nb704GHgS8A+4CLVHW7u+/HwFVAK/AdVV3e2WtZovcoVdi/yymHu+t9KF/ttNqr3SmCA0fDUWfDlEvgsMlx74MP5ZVNe/j2n99jWG4GD185ncMHZ/d4DNFY9sEuvrd0DaMHZ7HkyumMHJjZ80G0+uCjl5yEv/UVaKoFSYbDjnFW+BoxFYYd46wN0IN/nZnwdSvRu2u+fgh8GWcR8JXAJaq6IeiYbwOTVfVaEbkYOE9VLxKRicDjwHRgBPAyMF5VO5zjZYk+zlSdVZt8Dc5ydC0N0HQAmg84/7gbqp0umIZKJ7HX7nLWJq3cBr66g9fJG+X84z/8RDjyyzD4iB5N7n6/sr+xhR1V9bz+UQUrNlVQ+kklk0bk8eAVxzFkQHqPxRILb23dx/yHS/Grcsr4AmYeNZTjxw6iICedzNTknh1QbvVBWSls/adTKK18DTQFlWzILoD8MZAzHHJHwIChkDnIKUSXMdAZ/E3LcQrUpWQ45ZRTMmJWjdSE1t1E/0XgFlU9033+YwBV/WXQMcvdY94SkRRgN1CAu3Zs4Njg4zp6vWgTfc2+PVTdc3rE58VbpP88hc7+f2jb9QR1jlVnm+B3/6sk4ScJJYlWkvCTrH6SaSWZFtIIvzxuleSxTwbxWdJgdiUNZ0fSSHYmDWdbUhHVSfH5kz74t1dV57k6232tflr9SnOLn+oGH63+g0cfMzKXL00YxvxTxnq6X74zH+3Zz+L/+5hXN1Wwu7axbXt6ShK5mamkJgkpyUkkJ4nzORDn8xX8JRCPrwNRPyP8uzjcv4NC/04K/TsZ6q+gQPcxxL+PLBrCuk4rSbSQQgvJ+EmiVZz/+ts+sdLuB5CDzyHw+Qj9W2oUv723Oq5hX/Y4vvCDZ6I6t7NEH86/iJHAjqDnZcDxHR2jqi0iUgMMdre/3e7ckSECnA/MBxg9OrrKhZKcQmVWUVTnxlvkH8COj9fPfRUIKoH0HvQjQaleklBx0nyrpNAiqbRKMj5Jxydp+CSdpqRMGpMyaUrKoi4ph7rkXOqSc/FL6I9HgfsTLxL8+wclMyfRCanJSeRnpZGfnUZBTjonFA1iaG5GHCPqGeOG5fDL8yejqmzes5+1O2qorG+mqq6Z2kYfLa1Ki9/5Cf4SDNC4pq1cKjiKCmB1uz2p/iayW2vJ9teS1bqfDH896f4G0rWBVG0m1d9MijaTTAsp2kKy+khyGx+ifiSQ5tXfLs0rokpweu+4IRT57955oyoxmnNGxeW64ST6UFmn/TvU0THhnIuqLgIWgdOiDyOmQ+QOHMy06/8ezanGeIqIMOGwXCYc5o1CbKb3C2cefRkQ/DVTCJR3dIzbdZMHVIZ5rjHGmDgKJ9GvBMaJSJGIpAEXAyXtjikB5rmPLwBeUafzvwS4WETSRaQIGAe8G5vQjTHGhKPLrhu3z30BsBxneuViVV0vIrcCpapaAjwAPCIiW3Ba8he7564XkSeBDUALcF1nM26MMcbEnt0wZYwxfUBns26s1o0xxvRxluiNMaaPs0RvjDF9nCV6Y4zp4zw3GCsiFUB3FtMcAnwWo3BiyeKKjMUVGYsrMn0xrsNVNeRN655L9N0lIqUdjTwnksUVGYsrMhZXZPpbXNZ1Y4wxfZwlemOM6eP6YqJflOgAOmBxRcbiiozFFZl+FVef66M3xhjzeX2xRW+MMSaIJXpjjOnj+kSiF5Ffi8gmEVkrIs+IyMCgfT8WkS0isllEzuzhuC4UkfUi4heR4qDtY0SkQUTWuD/3eiEud1/C3q/2ROQWEdkZ9D6dlcBYZrvvyRYRuTFRcYQiIttF5AP3PUpYRUARWSwie0VkXdC2QSLykoh85P433yNxJfSzJSKjRORVEdno/lv8rrs9Pu+Xqvb6H+AMIMV9/CvgV+7jicD7QDpQBGwFknswrqOBo4AVQHHQ9jHAugS+Xx3FldD3K0SctwDXe+Dzley+F2OBNPc9mpjouILi2w4M8UAcpwDTgj/bwB3Aje7jGwP/Nj0QV0I/W8BwYJr7OAf40P33F5f3q0+06FX1RVUNrHr9Ns5KVgBzgaWq2qSqHwNbgOk9GNdGVd3cU68Xrk7iSuj75WHTgS2quk1Vm4GlOO+VCaKqr+OsRxFsLrDEfbwEOLdHg6LDuBJKVXep6nvu4/3ARpz1tOPyfvWJRN/OlcDz7uNQC5sfsjh5ghSJyGoReU1ETk50MC4vvl8L3C65xYn4s9/lxfclmAIvisgqEZmf6GDaGaaqu8BJbsDQBMcTzAufLURkDDAVeIc4vV/hLA7uCSLyMnBYiF03qerf3GNuwlnJ6s+B00IcH9P5pOHEFcIuYLSq7hORLwDPisgkVa1NcFxxf78OecFO4gT+CPzMjeFnwG9wvsh7Wo+/LxE6SVXLRWQo8JKIbHJbsaZjnvhsicgA4K/A91S1ViTUR637ek2iV9VZne0XkXnAV4DT1e3gogcWJ+8qrg7OaQKa3MerRGQrMB6I2UBaNHGRgMXcw41TRP4EPBfPWDrh6UXuVbXc/e9eEXkGp6vJK4l+j4gMV9VdIjIc2JvogABUdU/gcaI+WyKSipPk/6yqT7ub4/J+9YmuGxGZDfwIOEdV64N2eXJxchEpEJFk9/FYnLi2JTYqwGPvl/tBDzgPWNfRsXG2EhgnIkUikoazJnJJgmL5HBHJFpGcwGOciQmJep9CKQHmuY/nAR39NdmjEv3ZEqfp/gCwUVXvCtoVn/crUaPOMR7B3oLTh7rG/bk3aN9NODMmNgNzejiu83Bag03AHmC5u/1rwHqc2RvvAV/1QlyJfr9CxPkI8AGw1v0HMDyBsZyFMzNiK073V8Lel3ZxjXU/R++7n6mExQY8jtMt6XM/X1cBg4F/Ah+5/x3kkbgS+tkCZuB0G60Nyltnxev9shIIxhjTx/WJrhtjjDEds0RvjDF9nCV6Y4zp4yzRG2NMH2eJ3hhj+jhL9MYY08dZojfGmD7u/wM21MDT1Sm9XwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "q=[]\n",
    "pdf=[]\n",
    "for i in range(p):\n",
    "    q.append(np.linspace(qBound[i][0],qBound[i][1],n[i]))\n",
    "    pdf_=np.exp(-q[i]**2/(2*sig[i]**2))/(sig[i]*mt.sqrt(2*mt.pi))\n",
    "    pdf.append(pdf_)\n",
    "    plt.plot(q[i],pdf[i],label='pdf of q'+str(i+1))\n",
    "plt.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the function values at the samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "fEx=np.zeros(n)\n",
    "for j in range(n[1]):\n",
    "    for i in range(n[0]):\n",
    "        fEx[i,j]=c[0]*q[0][i]+c[1]*q[1][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the Sobol indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "sobol_=sobol(q,fEx,pdf=pdf)\n",
    "Si=sobol_.Si\n",
    "STi=sobol_.STi\n",
    "Sij=sobol_.Sij"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Validate the computed indices against corresponding exact values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " > Main Indices by UQit:\n",
      "\t S1=0.307692, S2=0.692308, S12=1.93492e-22\n",
      " > Main Analytical Reference:\n",
      "\t S1=0.307692, S2=0.692308, S12=0\n"
     ]
    }
   ],
   "source": [
    "Si_ex=[(c[0]*sig[0])**2./(c[0]**2*sig[0]**2+c[1]**2*sig[1]**2),\n",
    "           (c[1]*sig[1])**2./(c[0]**2*sig[0]**2+c[1]**2*sig[1]**2)]\n",
    "Sij_ex=[0]\n",
    "print(' > Main Indices by UQit:\\n\\t S1=%g, S2=%g, S12=%g' %(Si[0],Si[1],Sij[0]))\n",
    "print(' > Main Analytical Reference:\\n\\t S1=%g, S2=%g, S12=%g' %(Si_ex[0],Si_ex[1],Sij_ex[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 3: Sobol indices for 3 uniform parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following analytical function, \n",
    "$$\n",
    "\\begin{equation}\n",
    "f(\\mathbf{q}) = \\sin(q_1)+a\\sin^2(q_2)+b\\,q_3^4\\sin(q_1) \\,,\n",
    "\\label{eq:ishigami}\\tag{4}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "which is called [Ishigami](https://inis.iaea.org/search/search.aspx?orig_q=RN:21024954) function and widely used  to validate different methods in UQ. \n",
    "\n",
    "Here, $a$ and $b$ are two fixed parameters and the uncertain parameters $\\mathbf{q}=\\{q_1,q_2,q_3\\}$ have uniform distributions over the same range: $q_i\\sim \\mathcal{U}[-\\pi,\\pi]\\,,\\quad i=1,2,3$.\n",
    "The aim is to compute the main Sobol indices using `UQit` and validate them with the analytical values.\n",
    "\n",
    "The given model function (simulator) is implemented in `analyticTestFuncs.fEx3D` with keyword `Ishigami`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the number of samples for $q_1$, $q_2$ and $q_3$ and their associated admissible spaces $\\mathbb{Q}_1$, $\\mathbb{Q}_2$ and $\\mathbb{Q}_3$. Also set the model's fixed parameters $a$ and $b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=[100, 70, 80]        #number of samples for q1, q2, q3\n",
    "qBound=[[-mt.pi,mt.pi],      #admissible space of the parameters\n",
    "        [-mt.pi,mt.pi],\n",
    "        [-mt.pi,mt.pi]]\n",
    "a=7   \n",
    "b=0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate uniformly-spaced samples for parameters $q_1$ and $q_2$. Also define the PDFs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "q=[]\n",
    "pdf=[]\n",
    "for i in range(len(n)):\n",
    "    q.append(np.linspace(qBound[i][0],qBound[i][1],n[i]))\n",
    "    pdf.append(np.ones(n[i])/(qBound[i][1]-qBound[i][0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 4:** Run the simulator (4) at the parameter samples to generate training model outputs. We assume that a 3D sample parameter is obtained by tensor-product of the samples for $q_1$, $q_2$ and $q_3$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "fEx_=analyticTestFuncs.fEx3D(q[0],q[1],q[2],'Ishigami','tensorProd',{'a':a,'b':b})\n",
    "fEx=np.reshape(fEx_.val,n,'F')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute first- and second-order main Sobol sensitivity indices using `UQit`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "sobol_=sobol(q,fEx,pdf)\n",
    "Si=sobol_.Si\n",
    "Sij=sobol_.Sij\n",
    "SijName=sobol_.SijName\n",
    "STi=sobol_.STi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Validate the computed Sobol indices through comparison with the analytical values.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "fEx_.sobol(qBound)\n",
    "Si_ex=fEx_.Si\n",
    "STi_ex=fEx_.STi\n",
    "Sij_ex=fEx_.Sij"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Results:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " > Main Indices by UQit : S1=0.313905, S2=0.442312, S3=5.64351e-32\n",
      " >                        S12=9.90498e-32, S13=0.243783, S23=8.41268e-32\n",
      " > Total Indices        : ST1=0.557688, ST2=0.442312, ST3=0.243783\n",
      " > Main Analytical Reference: S1=0.313905, S2=0.442411, S3=0\n",
      " >                           S12=0, S13=0.243684, S23=0\n",
      " > Total Indices            : ST1=0.557589, ST2=0.442411, ST3=0.243684\n"
     ]
    }
   ],
   "source": [
    "#UQit \n",
    "print(' > Main Indices by UQit : S1=%g, S2=%g, S3=%g' %(Si[0],Si[1],Si[2]))\n",
    "print(' >                        S12=%g, S13=%g, S23=%g' %(Sij[0],Sij[1],Sij[2]))\n",
    "print(' > Total Indices        : ST1=%g, ST2=%g, ST3=%g' %(STi[0],STi[1],STi[2]))\n",
    "#Reference\n",
    "print(' > Main Analytical Reference: S1=%g, S2=%g, S3=%g' %(Si_ex[0],Si_ex[1],Si_ex[2]))\n",
    "print(' >                           S12=%g, S13=%g, S23=%g' %(Sij_ex[0],Sij_ex[1],Sij_ex[2]))\n",
    "print(' > Total Indices            : ST1=%g, ST2=%g, ST3=%g' %(STi_ex[0],STi_ex[1],STi_ex[2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 4: Sobol indices for 4 parameters with Gaussian distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do Example 15.8 in [Smith](https://rsmith.math.ncsu.edu/UQ_TIA/).\n",
    "\n",
    "$$\n",
    "f(\\mathbf{q})=q_1 q_3 +q_2 q_4\n",
    "$$\n",
    "\n",
    "where $q_i\\sim\\mathcal{N}(0,\\sigma_i^2)$ for $i=1,2,3,4$. \n",
    "For different values of $\\sigma_i$, we can compare the Sobol indices computed by `UQit` with associted exact values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU1f34/9d7lmSykJWwhLAEjbKUzQZwwSquaBW0YqVqXetW+XTx596KS2ur1Vpr1Y9StdJ+6lIXFJXWHZdvRQFZgyAQthCWkI1sM5nl/P6YmTCZTJJJMiGT8H4+HvPIzL3n3jknhPe859xzzxFjDEoppfouS09XQCmlVPfSQK+UUn2cBnqllOrjNNArpVQfp4FeKaX6OFtPVyBc//79zYgRI3q6Gkop1ausWLFivzEmJ9K+uAv0I0aMYPny5T1dDaWU6lVEZHtr+7TrRiml+jgN9Eop1cdpoFdKqT4u7vrolVJ9m9vtpqSkBKfT2dNV6ZUcDgd5eXnY7faoj9FAr5Q6pEpKSujXrx8jRoxARHq6Or2KMYby8nJKSkrIz8+P+rioum5EZIaIbBSRzSJye4T914vIWhFZJSKfi8iYwPYRItIQ2L5KRJ6KumZKqT7J6XSSnZ2tQb4TRITs7OwOfxtqN6MXESvwBHA6UAIsE5FFxpj1IcVeMMY8FSg/E3gEmBHYt8UYM7FDtVJK9Wka5DuvM7+7aLpupgCbjTHFgTd5CZgFNAV6Y8yBkPIpgM59rPqEPXV7WLhpIV7jbbY9LSGNS8dcikV0PIOKf9H8lQ4Bdoa8Lglsa0ZEbhSRLcAfgJ+F7MoXkZUi8omInNil2ip1iL25+U2eXP0k89fMb3o8veZpHlr+EJsqN/V09VQ3WrJkCeeccw4ALpeL0047jYkTJ/Lyyy9HdfyGDRuYOHEikyZNYsuWLVEd8+mnn3LMMcdgs9l49dVXO133cNFk9JG+J7TI2I0xTwBPiMjFwK+By4HdwDBjTLmIfBd4Q0TGhn0DQESuBa4FGDZsWAeboFT3qfPUkWBJYMWPVzRt+2/pf7nu/euo99T3YM3UobRy5UrcbjerVq2K+pg33niDWbNmce+990Z9zLBhw3j++ed5+OGHO1PNVkWT0ZcAQ0Ne5wGlbZR/CTgPwBjjMsaUB56vALYAR4UfYIyZb4wpNMYU5uREnKpBqR7h9Dhx2BzNtiXZkgBo8DT0RJVUF23bto1Ro0Zx+eWXM378eGbPnk19vf9D+z//+Q+jRo1i2rRpvP766wDs27ePSy+9lFWrVjFx4sQW2fmqVas49thjGT9+POeffz6VlZUsXryYRx99lGeeeYbp06e3qMPf/vY3jjrqKE466SSuueYa5s6dC/ingBk/fjwWS2y7BKPJ6JcBBSKSD+wC5gAXhxYQkQJjTPB77PeBTYHtOUCFMcYrIiOBAqA4VpVXqrs1eBpaBHqH1dG0T3XNvW8Vsb70QPsFO2BMbhp3nzu2zTIbN27k2Wef5YQTTuCqq67iySefZO7cuVxzzTV89NFHHHnkkVx00UUADBgwgGeeeYaHH36Yt99+u8W5LrvsMv7yl79w0kknMW/ePO69914effRRrr/+elJTU7n55publd+9ezd33303K1asID09nenTpzNp0qTY/QIiaPdjwxjjAeYC7wLfAP8yxhSJyH2BETYAc0WkSERWATfh77YB+B6wRkRWA68C1xtjKmLeCqW6idPjJNmW3GxbMKN3evSGn95q6NChnHDCCQBceumlfP7552zYsIH8/HwKCgoQES699NJ2z1NdXU1VVRUnnXQSAJdffjmffvppm8d8+eWXnHzyyeTk5JCQkND0gdKdorphyhizGFgctm1eyPOft3Lca8BrXamgUj0pUtdN8LUG+q5rL/PuLuFDFIOvD9Wwz0M9vFTHhinVhgZPQ1NXTZD20fd+O3bs4IsvvgDgxRdfZNq0aYwaNYqtW7c29cG/+OKL7Z4nPT2dzMxMPvvsMwD+8Y9/NGX3rZk6dSpLliyhvLwct9vNK6+80sXWtE8DvVJtaPA2NAX2oKauG69m9L3V6NGjWbBgAePHj6eiooIbbrgBh8PB/Pnz+f73v8+0adMYPnx4VOdasGABt9xyC+PHj2fVqlXMmzevzfKDBw/mnnvu4bjjjuO0007jmGOOadq3bNky8vLyeOWVV7juuusYOzY233h0rhul2tDgaSDbkd1sm91ixyIW6t06vLK3slgsPPVUyxlZZsyYwYYNG1psP/nkkzn55JMjnmvixIksXbq0xfZ77rmn1fe/8sorufLKKwF4/vnnmxZbmjx5MiUlJVG0oGM0o1eqDU6Ps0VGLyI4rA7N6FWvoRm9Um2IFOjB332jF2N7pxEjRrBu3bqerkaTK664giuuuKJb30MzeqXaEGkcPfhH3ujFWNVbaKBXqg2a0au+QAO9Uq1w+9x4jKfF8ErwB/oGr2b0qnfQQK9UK4JdM6123bg10KveQQO9Uq0Ids202nWjo276tJ6YpviRRx5hzJgxjB8/nlNPPZXt27d3uv6hNNAr1Yq2Ar3D6tA++sNI6DTF0c5NE5ymeOXKlRxxxBFRHTNp0iSWL1/OmjVrmD17NrfeemtXqt1EA71SrQh23UQM9DrqpteK52mKp0+fTnKyfxK9Y489NmY3T+k4eqVa0VYfvY66iZF/3w571sb2nIPGwVkPtFmkN0xT/Oyzz3LWWWd14RdxkGb0SrWiKdC3NupGM/peK96nKf6///s/li9fzi233NKJ1rWkGb1SrWjqo7dH7rpxep34jE8XCO+KdjLv7hLP0xR/8MEH3H///XzyySckJibG5P30L1SpVgRH1SRZI4+6AXB5XYe0Tio24nWa4pUrV3LdddexaNEiBgwY0NnmtaAZvVKtaHMcfchygpEu1qr4Fpym+LrrrqOgoKDFNMX9+/dn2rRpUc2Js2DBAq6//nrq6+sZOXIkf/vb39osHzpN8eDBgznmmGPwer0A3HLLLdTW1nLhhRcC/sXCFy1a1OX2aqBXqhVtjbrR5QR7t3idpviDDz6IovYdp103SrUiGMRbG3UTWkapeKYZvVKtCGb0idaWF8SCwV9H3vQ+Ok1xK0RkhohsFJHNInJ7hP3Xi8haEVklIp+LyJiQfXcEjtsoImfGsvJKdafgzJWRRtXourGqN2k30IuIFXgCOAsYA/woNJAHvGCMGWeMmQj8AXgkcOwYYA4wFpgBPBk4n1Jxz+l1RhxDDwczep3vRvUG0WT0U4DNxphiY0wj8BIwK7SAMeZAyMsUwASezwJeMsa4jDFbgc2B8ykV99oaURM66kapeBdNH/0QYGfI6xJganghEbkRuAlIAE4JOTb0cnRJYFv4sdcC14J/OJFS8aC11aUAkm3++Uj0YqzqDaLJ6CPdwmVabDDmCWPMEcBtwK87eOx8Y0yhMaYwJycniiop1f3aCvR6Mbbv64lpip966inGjRvHxIkTmTZtGuvXr+90/UNFE+hLgKEhr/OA0jbKvwSc18ljlYobrS0jCBroDzeHapriiy++mLVr17Jq1SpuvfVWbrrppq5Uu0k0gX4ZUCAi+SKSgP/iarNbtUSkIOTl94FNgeeLgDkikigi+UAB8FXXq61U93N6nO1m9Np10/vE8zTFaWlpTWXq6upiNvdOu330xhiPiMwF3gWswHPGmCIRuQ9YboxZBMwVkdMAN1AJXB44tkhE/gWsBzzAjcYYb0xqrlQ3a/A0MMg6KOI+u8WOzWLTjL6LHvzqQTZUtLwTtStGZY3itim3tVkmnqcpfuKJJ3jkkUdobGzko48+isFvJMpx9MaYxcaYo4wxRxhj7g9smxcI8hhjfm6MGWuMmWiMmW6MKQo59v7AcUcbY/4dk1ordQg4va133YAuJ9ibxfM0xTfeeCNbtmzhwQcf5Le//W0nW9ic3hmrVCvauhgL/lktteuma9rLvLtLPE9THDRnzhxuuOGGmLyfznWjVCvam5nSYXNQ76k/hDVSsRKv0xRv2rSp6fk777xDQUFBpFN0mGb0SkVgjGnzYizocoK9WbxOU/z444/zwQcfYLfbyczMZMGCBTFprwZ6pSJo9DViMO1m9Broe6d4nab4z3/+cxS17zjtulEqggZ363PRBzlsDh11o3oFzeiViiA4mqa1Sc3A/yFQ7ao+VFVSMaLTFCulAJousrY36kYzetUbaKBXKoJg37t23ai+QAO9UhG0tYxgkI66Ub2FBnqlIghm6sHpiCPRjF71FhrolYog2oze7XPj8XkOVbXUIdQT0xQHvfrqq4hI07DLrtJRN0pF0OD1Z+rtjboBcHld2Cz6X6kvC52mOFrBaYrvvffeDr1XTU0Njz32GFOntljfqdM0o1cqgmCXTJsXY3U5wV4pnqcpBrjrrru49dZbcThaTzI6StMQpSKIquvG7v8Q0EDfeXt+9ztc38R2muLE0aMYdOedbZaJ12mKV65cyc6dOznnnHN4+OGHY/Qb0YxeqYiiGl5p1cVHeqt4nKbY5/Pxy1/+kj/+8Y9dbF1LmtErFUGDpwGrWLFb7K2W0eUEu669zLu7xOM0xTU1Naxbt65pTp09e/Ywc+ZMFi1aRGFhYZfeTzN6pSIIzkXf1n/8YLavGX3vE4/TFKenp7N//362bdvGtm3bOPbYY2MS5EEzeqUiam8uejgY6DWj733idZri7qKBXqkInF5nm0MrIWTUjVcDfW8Tr9MUh1qyZEmrx3dUVF03IjJDRDaKyGYRuT3C/ptEZL2IrBGRD0VkeMg+r4isCjwWxazmSnWj9hYdgYOjbrTrRsW7djN6EbECTwCnAyXAMhFZZIxZH1JsJVBojKkXkRuAPwDBFW8bjDETY1xvpbpVg6ehzekPQMfR91Y6TXFkU4DNxphiY0wj8BIwK7SAMeZjY0xw8cylQF5sq6nUoRVVRq8XYzvNGNPTVei1OvO7iybQDwF2hrwuCWxrzdXAv0NeO0RkuYgsFZHzIh0gItcGyiwvKyuLokpKda/gqJu2BPdroO8Yh8NBeXm5BvtOMMZQXl7e4btmo7kYG2l8WcR/IRG5FCgEQscXDTPGlIrISOAjEVlrjGl2D7ExZj4wH6CwsFD/9VWPi2bUjUUsJFoTteumg/Ly8igpKUGTus5xOBzk5XWs0ySaQF8CDA15nQeUhhcSkdOAXwEnGWNcwe3GmNLAz2IRWQJMAjo2lZtSh1g0o27A332jgb5j7HY7+fn5PV2Nw0o0XTfLgAIRyReRBGAO0Gz0jIhMAp4GZhpj9oVszxSRxMDz/sAJQOhFXKXiUjR99ODvvgmuL6tUvGo3ozfGeERkLvAuYAWeM8YUich9wHJjzCLgISAVeCVwJ+EOY8xMYDTwtIj48H+oPBA2WkepuBTNqBvwj7zRjF7Fu6humDLGLAYWh22bF/L8tFaO+y8wrisVVOpQ8xkfLq8rqoxelxNUvYHOdaNUmGhmrgzSQK96Aw30SoUJdsVE20evXTcq3mmgVypM8OJq1KNudK4bFec00CsVpsEdWEbQ3n7XjcPmaCqvVLzSQK9UmGBGn2SNso9eh1eqOKeBXqkwHeqjtzr0YqyKexrolQoTDPTRjrpp8DTovC0qrmmgVypMMEOPdhy913jx+DzdXS2lOk0DvVJhmvroo8jomxYI15E3Ko5poFcqTNOom44Eeh15o+KYBnqlwnR0HH3oMUrFIw30SoXpyKib4BBMHXmj4pkGeqXCNHgasFvs2Cztz/nX1HWj0yCoOKaBXqkw0c5FDwe7bjTQq3imgV6pMNEsIxikGb3qDTTQKxXG6XFGHeibLsZqH72KYxrolQrT4I0+o9dRN6o30ECvVJgGT0NUQyvh4BBM7bpR8UwDvVJhOnQx1q4XY1X8iyrQi8gMEdkoIptF5PYI+28SkfUiskZEPhSR4SH7LheRTYHH5bGsvFLdoSN99AmWBATRPnoV19oN9CJiBZ4AzgLGAD8SkTFhxVYChcaY8cCrwB8Cx2YBdwNTgSnA3SKSGbvqKxV7DZ6GqDN6EdHlBFXciyajnwJsNsYUG2MagZeAWaEFjDEfG2PqAy+XAnmB52cC7xtjKowxlcD7wIzYVF2p7tGRjB50gXAV/6IJ9EOAnSGvSwLbWnM18O+OHCsi14rIchFZXlZWFkWVlOo+HRl1A7rKlIp/0QR6ibAt4ioLInIpUAg81JFjjTHzjTGFxpjCnJycKKqkVPfpyKgb8I+80a4bFc+iCfQlwNCQ13lAaXghETkN+BUw0xjj6sixSsULt8+Nx+fpcEavgV7Fs/ZnbYJlQIGI5AO7gDnAxaEFRGQS8DQwwxizL2TXu8DvQi7AngHc0eVaK9VNwleX+npHJX96/1t8YUsFpjnsPHzhBFISbXoxVsW9djN6Y4wHmIs/aH8D/MsYUyQi94nIzECxh4BU4BURWSUiiwLHVgC/wf9hsQy4L7BNqbgUDPTBjP6V5Tv5cmsFLrev6VFV7+bf6/awtLi8qaxejFXxLJqMHmPMYmBx2LZ5Ic9Pa+PY54DnOltBpQ6l8EC/emc1U/Oz+MfVU5vK1Dd6+M7d77K6pJpTRw/EYXNooFdxTe+MVSpEvcc/Sthhc+B0e9m4t4bxeenNyiQn2CgY0I81JVWA9tGr+KeBXqkQocsIFpUewOszTMjLaFFuwtB01pRUY4zR4ZUq7mmgVypEaNfN6p3+jH3C0JaBfnxeBhV1jZRUNujwShX3NNArFSIYsJNsSawpqWJgWiID01qOqQ9m+atLqpr66I2JeHuJUj1OA71SIUKHV64pqWZ8hG4bgKMH9SPBamFNSTVJtiQMBpfXFbGsUj1NA71SIYIZvcdjo3h/HRPCLsQGJdgsjMlNY/XOqqYx9zryRsUrDfRKhQgG+q1lbiBy/3zQhLx01u6qJtGii4+o+KaBXqkQwdEzG0r93TDjh7Qe6MfnZVDf6OVAg39KpwavBnoVnzTQKxUimJWv31XHiOxk0pPtrZadMNTfrVNa4W12rFLxRgO9UiGcHicOq4O1u2pavRAbNLJ/KqmJNkoCgV776FW80kCvVIgGTwOJVge7q50t7ogNZ7EI44akU7yvEdBAr+KXBnqlQjR4GrCQALR9ITZo/NB0tpU1Nh2rVDzSQK9UCKfHic9rx2oRxuamtVt+Ql4Gbo+/H18DvYpXGuiVCuH0Omn02CgYkEpyQvuTu47PSweT0HSsUvFIA71SIRrcDTS4rBEnMotkSEYSWUkpTccqFY800CsVotpVj8djY/zQti/EBokI44b41znWjF7FKw30SoWodtZhfPaoM3qAiXnZGGPhgKu+G2umVOdpoFcqRG1jPRYSOXpQv6iPmZCXAT47JVVV3VgzpTpPA71SIZxeJ5lJKdit0f/XGJ+XjvElsPvAgW6smVKdF9Vfs4jMEJGNIrJZRG6PsP97IvK1iHhEZHbYPm9gwfCmRcOVikden8FjXAxIjT6bB8hOTcQmiZTV1XZTzZTqmnbHj4mIFXgCOB0oAZaJyCJjzPqQYjuAK4CbI5yiwRgzMQZ1VapbbdpbA9JIblr74+fDJdmSqGzQQK/iUzQZ/RRgszGm2BjTCLwEzAotYIzZZoxZA/i6oY5KHRJf79yPiGFYZvQXYoNSE5NxepyU1+riIyr+RBPohwA7Q16XBLZFyyEiy0VkqYicF6mAiFwbKLO8rKysA6dWKnZWlewF6FRGn+lIAUsja3ZVx7paSnVZNIFeImzryOKYw4wxhcDFwKMickSLkxkz3xhTaIwpzMnJ6cCplYqdzWWVACTbkzt8bHZyKmJxs3mvdt+o+BNNoC8Bhoa8zgNKo30DY0xp4GcxsASY1IH6KXXIbA8MjwwuDdgR/RKTsVo9FO/XQK/iTzSBfhlQICL5IpIAzAGiGj0jIpkikhh43h84AVjf9lFKHXo1TjeV9f4g3ZlA77A5sFndbCmri3XVlOqydgO9McYDzAXeBb4B/mWMKRKR+0RkJoCITBaREuBC4GkRKQocPhpYLiKrgY+BB8JG6ygVF4rL6hCLf53YJFtSh49PsiUhFjfFGuhVHGp/ej7AGLMYWBy2bV7I82X4u3TCj/svMK6LdVSq2xXvr4UuBHqHzYFPGqmodXHA6SbN0foShEodanpnrFL4M3qr1R/oHdaOd90kWZPwmkbAp1m9ijsa6JXCH+izU/3PO9t1A4C42aoXZFWc0UCvFLClrJbswMwHnb0YC2C1aT+9ij8a6NVhz+czbCuvIzPVf8tIVzL63AyrBnoVd6K6GKtUX7b7gBOn20d6ioH6sEC//Qv45AHwepoflNIfvv8IpGQDBzP6vGwbW8q060bFF83o1WGvOBCYUxw+LGLBbgmMmHEegFevgn3ftDxowzvw7h1NL4MfDoPSLWwrr8Pn68jN40p1L83o1WEv2NWSlODzj4eXwKwfH9wDtXvg6g8g77vND/r49/5Mf9wPoeC0pkA/IN2C0+2jtLqBvMyOT6WgVHfQjF4d9orLaklNtGGxug8Ordz+X1j+LEy9oWWQBzjxJuh/NLz9C3DVNB3XP00C59R+ehU/NNCrw17x/jpG5qRQ76n397W7nbDoZ5AxDE75VeSDbIkw63GoLoGPftuU0aen+GfqLtZ+ehVHNNCrw15xWR0j+6dQ5aoiMzETPn0IyjfBOY9CQkrrBw6dAlOugS+fJqN8KwA+qSM10cbW/ZrRq/ihgV4d1hoaveyqamBkTioVDRVkig3+36Mw4WI48tT2T3DqPEgbQvq7/sy/0lnJyJwUijXQqziigV4d1oKZ98icFCqdlWTt/QYcGXDm/dGdILEfnPMn7GUbSbckUOGsYGT/FO2jV3FFA706rAXnj8/PTqGioYysmjI4+w+QnBX9SY46A8ZdSKarnorq7YzMSWVXVQMNjd5uqrVSHaOBXh3WtgYy7wHpPhqNl6y0oTD2Bx0/0Zm/J8vro7J8IyNz/P362k+v4oUGenVYK95fR266A+f2DwHIHHEySKTVM9uRmkNWygAq6ss4ItMeOLeOvFHxQQO9OqwVl9X6L8Sufw2ArJHTO32urJwxVFpgZMWngXNrRq/igwZ6ddgyxlBcVsfoLEPF9s8ByEoZ1OnzZeZ8h0qLBduafzIkI0nH0qu4oYFeHbbKal3UuDxMd39GJf4Lp1mODlyEDZOVlI0RoXrrEr6bWa999CpuaKBXh61g18p39r1FRZo/k890ZHb6fMEPiQqrhZl8QnFZHcbo5Gaq50UV6EVkhohsFJHNInJ7hP3fE5GvRcQjIrPD9l0uIpsCj8tjVXGluqq4rI4jpYS08tVUDBxFsi25U4uOBAUDfWXuBKZU/Zsal5uyWlesqqtUp7Ub6EXECjwBnAWMAX4kImPCiu0ArgBeCDs2C7gbmApMAe4Wkc6nTErFUHFZLT+yf4qx2KjoN7BL2Twc/DZQnn88aQ07mSwb9YKsigvRZPRTgM3GmGJjTCPwEjArtIAxZpsxZg3gCzv2TOB9Y0yFMaYSeB+YEYN6K9Vl28uq+YH1c+SoGVR668l2ZHfpfE0Zff8j8SWk8kPrEg30Ki5EMx/9EGBnyOsS/Bl6NCIdOyTKY5XqVgP2fEKmqYJJl1Lx7XMMTB7Yokzd0i+pfvstCFtIxJKSQvZVV2IfPLhpW0ZiBoJQ4alFxv6A73/9Mo/v2QsM6+6mKNWmaAJ9pLtHor3CFNWxInItcC3AsGH6n0J1v0aPj5Mb3qcuMYuUI0+nYvXDjM4a3bTfV1/Pvof/SOULL2BJS8OS0nwWS295OdULFzLwzjtJP/88RASrxUpGYgYVDRXIpEtJXvl3+u/4NzD5ELdOqeaiCfQlwNCQ13lAaZTnLwFODjt2SXghY8x8YD5AYWGhDlNQ3W5XyTamy9cUD7+cAouVCldFUx97/YoVlN5xJ+6dO8m6/DJyfvELLEnNFwxv3LGD0jvvZPedd1Lz3nsMuu9e7AMGkOnIpNJVCUOnsMc+lMmVi4F5PdBCpQ6Kpo9+GVAgIvkikgDMARZFef53gTNEJDNwEfaMwDalelTj1y9iEx/eCZdQ467B4/OQLf3Y+8CDbL/0x+DzMfzvCxh4xx0tgjxAwrBhDP/73xlw+23UffEFW8+dSfXb75CVmEl5QzmI8G3uLMb5vqFx78YeaKFSB7Ub6I0xHmAu/gD9DfAvY0yRiNwnIjMBRGSyiJQAFwJPi0hR4NgK4Df4PyyWAfcFtinVo/pveZ2vfUcypGAilc5KrF7D2AffpOL558m46IeMfPMNkie33eUiFgvZV1xB/sLXsQ8fTunNN/O9z6r8GT1Qe/RsvEao/fIfh6JJSrUqqsXBjTGLgcVh2+aFPF+Gv1sm0rHPAc91oY5Kxda+DWTXbeY529Uc47Cz5UAFl33oI3n1Zgbffz8ZF3Rs9srEkSMZ8cI/2fX/3czUhe+xNCkVzoPcofn81zeW7258E8xvOjdZmlIxoHfGqsNP0UJ8CJv6nwKAc+HbnLXCYOac2+EgHyQ2G7m//x11Q7O45tUa6os3MzInhbd9x5FctwN2r45lC5TqEA306vBiDKZoISvMaAYNGUH91ytJ+8vLrB4hZP7y5106tSU5mR2/uhSvBUpunEuK28nqlGl4sULR6zFqgFIdp4FeHV72rUf2b+RNz1QmOlyU/OxnuLJTefQ8C9mpOV0+ferwfB4534J3RwmlN9/C0MG5rLBOgKKFoPPeqB6igV4dXooWYrDwYeMkRj15P6a+ni9+cTKS1o8Ea0KXT5/lyGL9cAsNcy+m9pNPmLXyLV51FkLVDtj1dQwaoFTHaaBXhw9jYN3rbE/7Lj9a8yHy7QZyH36I7dm+Ls9zExScBmHvmRPJuPBCjnjvVSp2JuOz2LX7RvUYDfTq8LFnLVRsYeWukZy+fRnZ111Lv1NOocJZ0aV56EM1TVXsqmTQXb/GMmo0V61ezK6UKVD0BvjCp4NSqvtpoFeHj6LX8Xls5CxZS2XOEPr/9KcAVDorYxbo0xPTsYiFCmcFkpDAsN/fT1pjPTuW2uBACexaHpP3UaojNNCrw4MxULSQPduOIruuil1X/QJLgr9PPpYZvUUs/vlunP77ApNGj2Zp4Qyy1m6iriwF1mn3jTr0NCblYQ0AABqvSURBVNCrw8PuVTRs3kX1imreGXEcw046DgCf8VHlrIpZoAd/902ls7LpddkFl7I7JZvdK3PwrdHuG3XoaaBXhwWz+lV2L8vAlZbBc2PPZnRuGgA1jTV4jCdmF2PBH+iDGT3A6OEDeGzibNwVjexfWgs7l8bsvZSKhgZ61fcZQ/lLb+CqsvPBGZfRf2A2aQ47QFNA7s6MfkxuGqtyCjhw4imUb0jF+f6CmL2XUtHQQK/6vMavFrN/uZd+k4/mrbSjGRvI5uFgoI9lRp/pyKTcWd70emT/FBJtFpacdhnWZDu7F3yCcTfG7P2Uao8GetWnGWPYff+DiMWQMu9BtpfXNwv0wcy7q8sIhspyZFHTWIPb6wbAZrUwanAaK6t9DLzmBzjLhMrHfxez91OqPRroVZ924J13qP+2jJzTh7HJ5p/iYGxuetP+7sjom9aOdR3svhmbm8b60gP0u+JmUnLdlD3/Ku59+2L2nkq1RQO96rO8tXXs+/1vcWQ2knn5tRSVVgNE7rpJ7IZA72we6A84PeyqtzJozlSM28O+P/whZu+pVFs00Ks+a///PomnvJpBx3mRMedSVHqA/qkJ5PRLbCpT4aygX0I/7FZ7zN43+O0gtJ9+zGD/h0tR6QESTv0JWaNqOfD2O9Qv1xuoVPfTQK/6JNeWLVQsWED6ES6STr0A7A6KSg8wJjcdCVkApNJZGdP+eYic0Y8alIZFYH1pNeSfRP9j07Cl2dhz328wHk9M31+pcBroVZ9jjGHPb36LJcHGgHGVMPFiXB4vm/bWNOu2AX9GH8v+eQiZ7yZkLH1SgpUjclIpKj0AFguWwosZOL4M17ffUvnCCzF9f6XCaaBXfU7Nf/5D/dKl5ByXhG3o0ZB7DJv21uLxmYiBPpZj6AH6JfTDJrZmGT34++mLSg/4X0y8mH5DGkgZM4Syx/6Cp6wspnVQKpQGetWn+Orq2PvAgyQWjCQzez1MvBhEWB8IsKEjbqB7MnqLWMhwZDTL6IPvveeAk/JaF2SNRIYfz8BJ5fhcTvY9/MeY1kGpUFEFehGZISIbRWSziNweYX+iiLwc2P+liIwIbB8hIg0isirweCq21Vequf1PPYVn714GnTsCsVph/EUAFJVWk5JgZXhWclNZn/FR5YrtPDdB4dMggP8OWYD1uwNZ/aRLSPQWk33BDKrffJP6FStiXg+lIIpALyJW4AngLGAM8CMRGRNW7Gqg0hhzJPAn4MGQfVuMMRMDj+tjVG+lWnBt3kz58wtIP+88kms+hIIzoN9AwD/aZfTgNCyWgxdiq13V+IyvWwJ9piMzQkZ/cOQNAGPOA3sK/b9Th23wYPbcex/G7Y55XZSKJqOfAmw2xhQbYxqBl4BZYWVmAcEJPF4FTpXQoQ1KdTPj87H713dhTUlhwOzJULPb320D+HyGb3YfiNg/D7Gd5yYoUkafkZzAkIykg4E+MRXGzMLy7SIG3X4zrm+/pfy5v8W8LkpFE+iHADtDXpcEtkUsY4zxANVAcMxavoisFJFPROTESG8gIteKyHIRWV6mF6VUJ1S++CINq1Yx4PbbsG1dBElZcNQMALaV11HX6I3YPw+xvSs2KHxis6AxuWlNN24BMOkSaKyh36Aa+p1xBvufeALX1q0xr486vEUT6CNl5uHL2bdWZjcwzBgzCbgJeEFE0loUNGa+MabQGFOYk5MTRZWUOsi9Zw9lj/yJlOOPJ/30E2HDOzD+h2DzLywSzKDHHOKMvtZdS6O3+eRlY3PT2Lq/jjpXYOz8sOMhYzis+icDf/0rJDGRPfPuxuic9SqGogn0JcDQkNd5QGlrZUTEBqQDFcYYlzGmHMAYswLYAhzV1UorFWSM8fdte70MuvcepOh18DbCxEuayhSVHsBmEQoGpjY7Nphxd1cfPdDyguzgNIyBDXsC3TcWi7+uWz/FnuBkwC03U79sGVWvvRbzOqnDVzSBfhlQICL5IpIAzAEWhZVZBFweeD4b+MgYY0QkJ3AxFxEZCRQAxbGpulJQ8+671H78MTk/+xkJQ4fCqn/CoHEweHxTmfW7D1AwsB+JNmuzY4NBOCMxI+b1inTTFMDYIf7uo+BwTwAm/ggwsPolMmbPJnnyZPY99LBOeqZipt1AH+hznwu8C3wD/MsYUyQi94nIzECxZ4FsEdmMv4smOATze8AaEVmN/yLt9caY5n/5SnWSt7qaPb+9H8fYsWRd9mMoWQ6lK5tl8z6foWhXdYsLseAPwumJ6dgstpjXLdI0CAC56Q4yk+2sKQnpp88YBvknwYrnEZ+bQffdi3E62Xu/TmWsYiOqv3BjzGJgcdi2eSHPncCFEY57DdDvoKpb7H3oIbyVlQz763zEZoMlv4fkbJj046Yyy7dXUl7XyAlHtpzPpjvuig1qLaMXEY4/sj8fbtiH2+vDbg3kWtN+Af84H1b+g8TJP6H/T39K2aOPUvPhh/Q79dRuqaM6fOidsapXqv3sc6pffY3sq67EMXo07PwKNn8Ax//MP2wx4K3VpSTaLJw+ZlCLc1Q6K2M6PXGo1vroAc4dn0tFXSP/3XJwdktGToehx8Jnj4DHRfbVV5F41FHsufc+PBX6JVh1jQZ61eu49+yh9NZbSTzqKPrfeKN/YzCbn3JNUzmP18fitbs5bfRAUhNbfnmtcFaQnRTbmSuD+tn7YbPYIgb6k4/OoV+ijUWrQsY0iMD0O+DALvj674jdTu4Dv8dbVUXprbfpKBzVJRroVa9i3G52/fImjMvFkEcfxeJwwI4vYctHcMLPISGlqez/21JOeV0j507IjXiu7szoRYSsxMhj6R12K2d+ZxDvFe3B6fYe3JF/kn+45Wd/BLcTx5gxDLzzTuo+/5zy+fO7pZ7q8KCBXvUq+/70KA0rVzL4t78hcWS+f+OS30FKDkz+SbOyi1aV0i/RxslHt7w3w+vz+ue5SeqePnqArKSWd8cGzZyQS43Lw5KNISNrgll9zW742n+jecZFPyTtnHMoe+wv1C39stvqqvo2DfSq16j58EMqnnuOzIsvJu3ss/0bt/8Xipe0yOadbi/vFe3hjLGDcNitLc5V5arCYLotowf/8oSRMnqA44/IJjslgbdW726+I/97MHxaIKtvQEQYfO89JIwYwa6bb9Yhl6pTNNCrXqGxpITSO+7EMXYsA26/7eCOj38HKQOg8Opm5ZdsLKPG5WHmxMjdNk13xXZzRh+6nGAom9XC2eMG88E3e6l1ha0wNf0OqN0Ly/3z3lhSUsj786P4amsp/f9u1hWpVIdpoFdxz9fYyK6f/wKMYcifH8WS4J/agG2fw7bP/EMTE5KbHfPW6lKyUxI44YjIF1ub7opN7L5A31ZGDzBzYi4uj4/31+9pvmPENBhxInz+J2isByCxoIBBd99N/bJllD3+eLfVWfVNGuhVXDMeD6W33oazqIjcB35PQl7ewZ0f/x5SB0LhVc2OqXV5+OCbvZw9bjA2a+Q/8e6c5yYoOymbek89To8z4v7vDsskN93RfPRN0PQ7oW4fLH+uaVPG+eeRPvsCyp96WqdIUB2igV7FrWCQr/nPfxhw223Nbxz65i3Y/jlM+yXYk5od98H6vbg8vla7baB7Z64MCvb/t5bVWyzCuRNy+WzTfirrmk9+xvDj/aNwPn8Eag5m/IPuuouUE09k96/vomrhG91Wd9W3aKBXccl4vZTefgcHFi9mwC03k33lFQd37t8Mb/wUcie1yOYBFq0uZXC6g+8Oaz2IVzgrEKRb5rkJau3u2FDnTsjF4zP8e92eljvPehDcDfDKFeD1L0hiSUwk7y+PkXLccey+806qF4VPO6VUSxroVdwxXi+777yTA2+/Tc5NN5F9dciFVlctvHwpWGzww7+DLbHZsZV1jXz6bRnnTshttppUuEpnJRmJGVgtLUfkxEpbd8cGjc1NY2T/FBat3tVy54DRMPMvsOMLeO/XTZstDgd5TzxO8tSplN5+B9VvvR3zuqu+RQO9iivBlaKq31xEzs9/Rv9rrwnZaWDR/8D+jTD7Of9kYGH+vW4PHp9hZis3SQV15zw3QdmO7Kb3ao2Iv/vmy60V7KmO0Jc/bjYc+1P48ilY80rTZktSEkP/90mSCwspve02Dixe3PJYpQI00Ku44a2qYtfPf071woX0nzuX/jfc0LzA0ieh6HU45S44YnrEc7y1upSR/VMizlYZan/D/m7tn4eDGf3+hv1tlps5MRdj4J21uyMXOP0+GH6C/0Nuz7qmzZakJIY+9b8kH3MMu265lfJnn9OpElREGuhVXKj78iuKzzufmo+XMOC228iZe2PzAts+h/fuglHn+C/ARvD1jkqWbi3nnAm5tLVk8eqy1awqW8WkAZNi2YQWUuwpFGQW8Mq3r9DgaWi13BE5qYzNTeP/lm6nxhlhcXCrHWb/DZIy/N1WDVVNuyzJyQx9+in6nTKdfQ89xM6f/AT3Xr2pSjWngV71KNPYyL4/PsKOK67A4nAw4qWXml94Bajc7r8gmZUP5/2vf6qAMDvK67lmwXKGZiZzxfEjWn0/j8/Db5f+lgHJA/jJuJ+0Wi4WRIQ7ptzBrtpd/HXNX9sse+fZo9lZUc9P//k1bm+ErLzfQLhwAVSXwOvXgufgKB1LSgpDHnuMQffdS/3KVWydNYuaDz6IdXNUL6aBXvUYV3Ex2y6+hPK//pWM2ReQ//prJH1nbPNC616Hp08EtxMu+ic4WnbJVNU3csXzX+E1huevnExWSkKr7/nyxpfZULGB2ybfRrI9udVysTJ50GTOHXkufyv6G1urW1/0+4Qj+/O788fx2ab93PXGOowJX5YZGDYVznoANr0Lz5wKZRubdokImT/8IfmvvYY9N5eSuf/D7nl3462t645mqV5GA7065BqKiij55S8pPudcGnfuZMif/8zg3/wGS3JI4HUegIXXw6tXQnYBXPcJDBjV4lwuj5dr/7GCkooG5v+4kJE5qS3KBJXVl/H4ysc5IfcETh9+enc0LaKbCm8iyZrE7778XeQAHvDDyUP5n1OO5KVlO3lyyZbIhSb/BOa86J/O+OmT4Ku/+i9SBySOzGfESy+SdfVVVP3rX2w+5RT2/elRPPvbvk6g+jYN9OqQMMZQ98UX7LjqKrZdMJu6zz4n++qrOOLtt0g784zmhXd8CU9NgzUvw0m3wVX/gewjIp7z1lfX8NXWCh66cDxT8tseRfPw8odxeV3cMfWONvvwY61/Un/+55j/Yenupby77d02y950+lHMmpjLQ+9u5M1VEYZcAow6G274wn9T1eKb4YWLoLasabckJDDwllsY8corpBx7LOXz57P5lFPZfe+9NO7cGcumqV5C2sowekJhYaFZvnx5T1dDxYBpbKR+5SrqPv+M2iVLcG3ajDWnP1mXXUbmnDlY+/U7WNjrgW2fwtpXYfWLkJ4HP/grDDu21fM/8t5GHvtoM7eceTQ3Tj+yzbp8tfsrrn7vaq6fcD03TryxzbLdwevz8qN3fsT+hv0sOm8RqQmtf/Nwebz8+JmvWLWzin9eM5XJI1r5APP54Kv58P48f5fWd6+E71zQ4puPa+tWKp57juo33sR4vaQcO5WUE79H6onTSDjiiEP6oae6j4isMMYURtwXTaAXkRnAnwEr8Iwx5oGw/YnA34HvAuXARcaYbYF9dwBXA17gZ8aYNlMaDfS9kzEG7/79uLYU49q0ibqlS6n/4gt89fVgs5F8zDGknfN90mfNwpIYuMnJ54OdX8K612D9G1BXBgn9YPwP4bR7IvbH76pq4J01pby1ejdrd1VzUeFQHrhgXJvByu11c8FbF+D2ulk4ayEOm6N7fgntWFu2lksWX8Iloy/htim3tVm2qr6RH/zvf9lZUc9JRw3g3AmDOW30QFIirJTF3vXw7p2w9RMwPhgwFsZdAGN/4L+AHeDeu4/KF16g9qMPcW3aDIBt8GBSp00jufC7JIw8goT8fKypKS3fQ8W9LgV6EbEC3wKnAyXAMuBHxpj1IWV+Cow3xlwvInOA840xF4nIGOBFYAqQC3wAHGWM8Ya/T5AG+vhgfD6M04nP6cQ0NOBzOvHV1eGpqMBbUYm3ssL/fP9+XNu20Vi8FV9NTdPx9iFDSDl+CqmTx5P8nZFYLW7/iJHyzQcfFcXQWAs2Bxw1w5+NFpwO9iQaPT6q6hupqG+kss7Nt3treGt1Kcu3++eNmZCXzsyJQ7jsuOEHF9gOrb8x1LhrqHJWsXDzQp5Z+wxPnPoE38v73iH7HUZy3xf38fqm13nslMcoyCggw5FBki0pYtk91U6e/byYt1bvZs8BJw67hVNHD+TMsYMYkuEgIzmBrOQE0pPs/ruAa/b6PzDXveb/AAX/FM7ZR0L2yMDPIyF1IO6qRmpXfkPdVyupW7oMX21t0/vaBg4kYWQ+CUOHYc3KxJaVhTUzy/88MxNJSsKSlITF4UCSkpCEBP1WEAe6GuiPA+4xxpwZeH0HgDHm9yFl3g2U+UJEbMAeIAe4PbRsaLnW3q+zgX7Hxq/ZctUlHT4urkXZqxb6X0xMhH2m+U8xIQ/A4vM/rCE/rVHcd+O2Qm0y7M+EskzYlw1lWVCWZahJNc0r1tQkwYMVD3bcYsNFAvUk4cMSaK7BGPBG+LtMtFro57DTz2HDbrVgQn5BwefGGA40HqDaVY03JJ84ddipPDr90fYb1c2qXdXMfGNms7tlE62JZCT6A34wYErzf1Wcbi81Tje1Tg+esN+NABYRRA4eZ8NDCg3YcWM3Huy4sRI5vxKfkFUJORVCTgUMCDyyqiHZCZZ2/g59gM8KPgFvyE8j/j9hI/6HT2j6mzCBioef2oQ0O+pO5T70GVPV3875b67p1LFtBfoI3wNbGAKEXsEpAaa2VsYY4xGRaiA7sH1p2LFDIlTwWuBagGHDWt7WHg1rgoPqrGia07uYTvwRhx8T/E8Fgk8O/scLljUW8FrAWASvBXyBh9tuwW0Ht03w2KHRLjQkCXVJFhqSLTTaAbFgEIyI/yeCwUKSWPGKFR82vGLFiw23JRGXxYEJGQOQCDjEPzxQ8A+Rt4iQYLOQYLWQYLNgt1pIsltJTmw5L01oQPSfQ+iX0I+MxAzSE9PJdGSSkZjBsYNb7+s/lNIT03nl3FdYu38tVc4qqlwHH8GbqkKTLxMW7oyBGqebRo+PRq+PRo8Pt9dHo9dgjP9D0hj/ccb4g7Ar8LAaD4m+BmzGjRUPFuPFajxYrR5c2T5Ksg27MIgxCD4E/8mSnIbkekNygw+H02B3G+wesLkNdo//dTBZ8D+MP1EISSbE+PcFW0WEPCA8SYlGZ46JZ66M7uk2iyYyRgo14b/e1spEcyzGmPnAfPBn9FHUqYUh+WMY8tbazhyq1CE1IHkApw47tf2CSsVINMMrS4ChIa/zgPCVEprKBLpu0oGKKI9VSinVjaIJ9MuAAhHJF5EEYA4QPgn2IuDywPPZwEfG//1zETBHRBJFJB8oAL6KTdWVUkpFo92um0Cf+1zgXfzDK58zxhSJyH3AcmPMIuBZ4B8ishl/Jj8ncGyRiPwLWA94gBvbGnGjlFIq9vSGKaWU6gPaGnWjUyAopVQfp4FeKaX6OA30SinVx2mgV0qpPi7uLsaKSBmwvQun6A/0hcm3+0o7QNsSr/pKW/pKO6BrbRlujMmJtCPuAn1Xicjy1q489yZ9pR2gbYlXfaUtfaUd0H1t0a4bpZTq4zTQK6VUH9cXA/38nq5AjPSVdoC2JV71lbb0lXZAN7Wlz/XRK6WUaq4vZvRKKaVCaKBXSqk+rk8EehH5jYisEZFVIvKeiOQGtouIPCYimwP7j+npurZHRB4SkQ2B+i4UkYyQfXcE2rJRRM7syXpGQ0QuFJEiEfGJSGHYvt7WlhmBum4Wkdt7uj4dISLPicg+EVkXsi1LRN4XkU2Bn5k9WcdoichQEflYRL4J/G39PLC917VHRBwi8pWIrA605d7A9nwR+TLQlpcD08N3jX/5sd79ANJCnv8MeCrw/Gzg3/hXujoW+LKn6xpFW84AbIHnDwIPBp6PAVbjX30vH9gCWHu6vu20ZTRwNLAEKAzZ3qvagn967i3ASCAhUPcxPV2vDtT/e8AxwLqQbX8Abg88vz34dxbvD2AwcEzgeT/g28DfU69rTyAupQae24EvA3HqX8CcwPangBu6+l59IqM3xhwIeZnCweUKZwF/N35LgQwRGXzIK9gBxpj3jDGewMul+FflAn9bXjLGuIwxW4HNwJSeqGO0jDHfGGM2RtjV29oyBdhsjCk2xjQCL+FvQ69gjPkU/zoRoWYBCwLPFwDnHdJKdZIxZrcx5uvA8xrgG/zrUPe69gTiUm3gpT3wMMApwKuB7TFpS58I9AAicr+I7AQuAeYFNkda2LzF4uRx7Cr830ig97clVG9rS2+rbzQGGmN2gz94AgN6uD4dJiIjgEn4M+Fe2R4RsYrIKmAf8D7+b45VIcleTP7Wek2gF5EPRGRdhMcsAGPMr4wxQ4F/AnODh0U4VY+PJ22vLYEyv8K/Ktc/g5sinKpXtCXSYRG29Xhb2tDb6tvniUgq8Brwi7Bv9L2KMcZrjJmI/5v7FPzdnS2KdfV92l1KMF4YY06LsugLwDvA3cTp4uTttUVELgfOAU41gY46emlbWhGXbWlDb6tvNPaKyGBjzO5Ad+a+nq5QtETEjj/I/9MY83pgc69tD4AxpkpEluDvo88QEVsgq4/J31qvyejbIiIFIS9nAhsCzxcBlwVG3xwLVAe/3sUrEZkB3AbMNMbUh+zqSwut97a2LAMKAqMhEvCvibyoh+vUVYuAywPPLwfe7MG6RE1EBP8a1d8YYx4J2dXr2iMiOcFRdSKSBJyG/5rDx8DsQLHYtKWnrzzH6Or1a8A6YA3wFjAk5Kr2E/j7vdYSMvIjXh/4L0zuBFYFHk+F7PtVoC0bgbN6uq5RtOV8/NmwC9gLvNuL23I2/hEeW4Bf9XR9Olj3F4HdgDvw73E1kA18CGwK/Mzq6XpG2ZZp+Lsy1oT8Hzm7N7YHGA+sDLRlHTAvsH0k/sRnM/AKkNjV99IpEJRSqo/rE103SimlWqeBXiml+jgN9Eop1cdpoFdKqT5OA71SSvVxGuiVUqqP00CvlFJ93P8PHS6CgZ4hAQIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> Main Indices Si by UQit: [5.69332163e-35 4.46886262e-34 3.05804326e-35 1.63131702e-33]\n",
      "  Sij [1.24593888e-33 2.70270376e-02 1.20563049e-33 1.10254866e-33\n",
      " 9.72972962e-01 9.90293824e-34]\n",
      "> Reference Analytical Si: [0, 0, 0, 0]\n",
      "  Sij [0, 0.02702702702702703, 0, 0, 0.972972972972973, 0]\n"
     ]
    }
   ],
   "source": [
    "#Settings\n",
    "n=[60,60,60,60]       #number of samples for q1 and q2, Method1\n",
    "qBound=[[-30,30]]*4   #admissible range of parameters\n",
    "sig=[2.,3.,1,4]    #sdev of q's\n",
    "#(1) Samples from parameters space + the PDFs\n",
    "p=len(n)\n",
    "q=[]\n",
    "pdf=[]\n",
    "for i in range(p):\n",
    "    q.append(np.linspace(qBound[i][0],qBound[i][1],n[i]))\n",
    "    pdf_=np.exp(-q[i]**2/(2*sig[i]**2))/(sig[i]*mt.sqrt(2*mt.pi))\n",
    "    pdf.append(pdf_)\n",
    "#plot PDfs\n",
    "for i in range(p):\n",
    "   plt.plot(q[i],pdf[i],label='pdf of q'+str(i+1))\n",
    "plt.legend(loc='best')\n",
    "plt.show()\n",
    "#(2) Compute function value at the parameter samples\n",
    "fEx=np.zeros(n)\n",
    "for i3 in range(n[3]):\n",
    "    for i2 in range(n[2]):\n",
    "        for i1 in range(n[1]):\n",
    "            for i0 in range(n[0]):\n",
    "                fEx[i0,i1,i2,i3]=q[0][i0]*q[2][i2]+q[1][i1]*q[3][i3]\n",
    "#(3) Compute Sobol indices direct numerical integration\n",
    "sobol_=sobol(q,fEx,pdf=pdf)\n",
    "Si=sobol_.Si\n",
    "STi=sobol_.STi\n",
    "Sij=sobol_.Sij\n",
    "#(5) Exact Sobol indices (analytical expressions)\n",
    "Si_ex=[0]*p\n",
    "Sij_ex=[0,(sig[0]*sig[2])**2./((sig[0]*sig[2])**2.+(sig[1]*sig[3])**2.),0,0,\n",
    "        (sig[1]*sig[3])**2./((sig[1]*sig[3])**2.+(sig[0]*sig[2])**2.),0]\n",
    "#(6) Results\n",
    "print('> Main Indices Si by UQit:',Si)\n",
    "print('  Sij',Sij)\n",
    "print('> Reference Analytical Si:',Si_ex)\n",
    "print('  Sij',Sij_ex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 5: Sobol indices for arbitrary number of parameters with uniform distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we consider the test case which is considered in the original article by [Sobol](https://www.sciencedirect.com/science/article/abs/pii/S0378475400002706) as well as in Example 15.17 in [Smith](https://rsmith.math.ncsu.edu/UQ_TIA/):\n",
    "\n",
    "$$\n",
    "f(\\mathbf{q})=\\prod_{i=1}^p g_i(q_i)\\,,\\quad g_i(q_i)=\\frac{|4q_i-2|+a_i}{1+a_i}\\,, a_i\\geq 0\n",
    "$$\n",
    "\n",
    "Assume $q_i\\sim\\mathcal{U}[0,1]$ for $i=1,2,\\cdots,p$ for any $p>1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed Indices by UQit:\n",
      " Si: [0.13445718 0.21002473 0.06285799 0.15448376 0.1947089 ]\n",
      " Sij-Name: ['S12', 'S13', 'S14', 'S15', 'S23', 'S24', 'S25', 'S34', 'S35', 'S45']\n",
      " Sij: [0.0309405  0.00926014 0.0227583  0.0286842  0.01446451 0.0355489\n",
      " 0.04480528 0.01063938 0.0134097  0.03295654]\n",
      " STi: [0.22610031 0.33578393 0.11063172 0.25638687 0.31456462]\n",
      "Exact Sobol Indices:\n",
      " Si: [0.12932295 0.2020671  0.06011914 0.14845746 0.18622504]\n",
      " Sij: [0.02993587 0.00890654 0.0219937  0.0275889  0.01391647 0.03436515\n",
      " 0.04310765 0.01022434 0.01282542 0.03167093]\n"
     ]
    }
   ],
   "source": [
    "#---SETTINGS ---------------\n",
    "a=[0.5,0.2,1.2,0.4,0.25]   #coefficients in the model => p is determined from the length of this list\n",
    "nSamp=[20,20,21,22,25]   \n",
    "#---------------------------\n",
    "p=len(a)\n",
    "qBound=[[0,1]]*p\n",
    "#Generate samples, compute PDFs and evaluate model response at the samples\n",
    "q=[]\n",
    "pdf=[]\n",
    "for i in range(p):\n",
    "    q_=np.linspace(qBound[i][0],qBound[i][1],nSamp[i])\n",
    "    q.append(q_)\n",
    "    pdf.append(np.ones(nSamp[i])/(qBound[i][1]-qBound[i][0]))\n",
    "    fEx_=(abs(4*q_-2)+a[i])/(1+a[i])\n",
    "    if i==0:\n",
    "       fEx=fEx_\n",
    "    else:\n",
    "       fEx=np.tensordot(fEx,fEx_,0)\n",
    "#Computed Sobol indices\n",
    "sobol_=sobol(q,fEx,pdf)\n",
    "#Exact reference Sobol indices (Smith, p.336)\n",
    "Di=[]\n",
    "Dsum=1\n",
    "for i in range(p):\n",
    "    Di.append(1/(3*(1+a[i])**2.))\n",
    "    Dsum*=(1+Di[i])\n",
    "Dsum=-1+Dsum\n",
    "Di=np.asarray(Di)\n",
    "Si=Di/Dsum\n",
    "Dij=[]\n",
    "for i in range(p):\n",
    "    for j in range(p):\n",
    "        if i!=j and i<j:\n",
    "           Dij.append(Di[i]*Di[j])\n",
    "Dij=np.asarray(Dij)\n",
    "Sij=Dij/Dsum\n",
    "print('Computed Indices by UQit:')\n",
    "print(' Si:',sobol_.Si)\n",
    "print(' Sij-Name:',sobol_.SijName)\n",
    "print(' Sij:',sobol_.Sij)\n",
    "print(' STi:',sobol_.STi)\n",
    "print('Exact Sobol Indices:')\n",
    "print(' Si:',Si)\n",
    "print(' Sij:',Sij)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}