{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Orthogonality\n",
    "\n",
    "The core idea of polynomial chaos expansions is that the polynomials\n",
    "used as an expansion are all mutually orthogonal. The relation is\n",
    "typically written mathematically as:\n",
    "\n",
    "$$\\left\\langle \\Phi_n, \\Phi_m \\right\\rangle = 0 \\qquad n \\neq m$$\n",
    "\n",
    "In practice this relation is instead expressed by the equivalent\n",
    "notation using expected values:\n",
    "\n",
    "$$\\mbox E\\left(\\Phi_n \\Phi_m\\right) = 0 \\qquad n \\neq m$$\n",
    "\n",
    "In `chaospy` this property can be tested by taking the outer product of\n",
    "two expansions, and checking if the expected value of the resulting\n",
    "matrix is diagonal. For example, for a basic monomial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:04.334932Z",
     "iopub.status.busy": "2021-05-18T10:56:04.334479Z",
     "iopub.status.idle": "2021-05-18T10:56:04.344625Z",
     "shell.execute_reply": "2021-05-18T10:56:04.344267Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1, q0, q0**2, q0**3])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import chaospy\n",
    "\n",
    "monomial = chaospy.monomial(4)\n",
    "monomial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:04.346876Z",
     "iopub.status.busy": "2021-05-18T10:56:04.346564Z",
     "iopub.status.idle": "2021-05-18T10:56:04.362200Z",
     "shell.execute_reply": "2021-05-18T10:56:04.361856Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([[1, q0, q0**2, q0**3],\n",
       "            [q0, q0**2, q0**3, q0**4],\n",
       "            [q0**2, q0**3, q0**4, q0**5],\n",
       "            [q0**3, q0**4, q0**5, q0**6]])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "monomial2 = chaospy.outer(monomial, monomial)\n",
    "monomial2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:04.364138Z",
     "iopub.status.busy": "2021-05-18T10:56:04.363873Z",
     "iopub.status.idle": "2021-05-18T10:56:04.441459Z",
     "shell.execute_reply": "2021-05-18T10:56:04.441169Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  1.,  0.],\n",
       "       [ 0.,  1.,  0.,  3.],\n",
       "       [ 1.,  0.,  3.,  0.],\n",
       "       [ 0.,  3.,  0., 15.]])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "normal = chaospy.Normal(0, 1)\n",
    "chaospy.E(monomial2, normal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In other words, the basic monomial (beyond polynomial order 1) are not\n",
    "orthogonal.\n",
    "\n",
    "But if we replace the basic monomial with an explicit orthogonal\n",
    "polynomial constructor, we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:04.443434Z",
     "iopub.status.busy": "2021-05-18T10:56:04.443167Z",
     "iopub.status.idle": "2021-05-18T10:56:04.476480Z",
     "shell.execute_reply": "2021-05-18T10:56:04.476187Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1.0, q0, q0**2-1.0, q0**3-3.0*q0])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hermite = chaospy.generate_expansion(3, normal)\n",
    "hermite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:04.478557Z",
     "iopub.status.busy": "2021-05-18T10:56:04.478207Z",
     "iopub.status.idle": "2021-05-18T10:56:04.488859Z",
     "shell.execute_reply": "2021-05-18T10:56:04.488536Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0., 0.],\n",
       "       [0., 1., 0., 0.],\n",
       "       [0., 0., 2., 0.],\n",
       "       [0., 0., 0., 6.]])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hermite2 = chaospy.outer(hermite, hermite)\n",
    "chaospy.E(hermite2, normal).round(15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A fully diagonal matrix, which implies all the polynomials in the\n",
    "expansion are mutually orthogonal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate expansions\n",
    "\n",
    "Multivariate orthogonal expansion are (usually) created by doing a tensor product of univariate expansions together.\n",
    "To illustrate how this work, consider the distribution introduced in the [problem formulation](../main_usage/problem_formulation.ipynb):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:04.490644Z",
     "iopub.status.busy": "2021-05-18T10:56:04.490388Z",
     "iopub.status.idle": "2021-05-18T10:56:06.431756Z",
     "shell.execute_reply": "2021-05-18T10:56:06.431463Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "J(Normal(mu=1.5, sigma=0.2), Uniform(lower=0.1, upper=0.2))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from problem_formulation import joint\n",
    "\n",
    "joint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extracting the marginal density we can construct both one-dimensional expansions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:06.433888Z",
     "iopub.status.busy": "2021-05-18T10:56:06.433615Z",
     "iopub.status.idle": "2021-05-18T10:56:06.458070Z",
     "shell.execute_reply": "2021-05-18T10:56:06.458295Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1.0, q0-1.5, q0**2-3.0*q0+2.21])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expansion0 = chaospy.generate_expansion(2, joint[0])\n",
    "expansion0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:06.460295Z",
     "iopub.status.busy": "2021-05-18T10:56:06.459997Z",
     "iopub.status.idle": "2021-05-18T10:56:06.485078Z",
     "shell.execute_reply": "2021-05-18T10:56:06.485296Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1.0, q0-0.15, q0**2-0.3*q0+0.02167])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expansion1 = chaospy.generate_expansion(2, joint[1])\n",
    "expansion1.round(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When constructing a multivariate expansion, it is canonical to truncate the expansion at order and graded lexicographical sorting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:06.487285Z",
     "iopub.status.busy": "2021-05-18T10:56:06.487022Z",
     "iopub.status.idle": "2021-05-18T10:56:06.520245Z",
     "shell.execute_reply": "2021-05-18T10:56:06.519960Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1.0, q1-0.15, q0-1.5, q1**2-0.3*q1+0.02167,\n",
       "            q0*q1-1.5*q1-0.15*q0+0.225, q0**2-3.0*q0+2.21])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chaospy.generate_expansion(2, joint).round(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See [chaospy.generate_expansion()](../../api/chaospy.generate_expansion.rst) for variations in truncations and sorting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wiener-Askey scheme\n",
    "\n",
    "Polynomial chaos expansion often assume that the polynomial expansion used is of the Wiener-Askey scheme verity.\n",
    "The reason for this is that the expansion in the scheme correspond to orthogonality with respect to some standard probability distribution.\n",
    "These include:\n",
    "\n",
    "* Hermite polynomials which are orthogonal with a normal density weight function.\n",
    "* Legendre polynomials which are orthogonal with a uniform density weight function.\n",
    "* Laguerre polynomials which are orthogonal with a exponential density weight function.\n",
    "* Generalized Laguerre polynomials which are orthogonal with a gamma density weight function.\n",
    "* Jacobi polynomials which are orthogonal with a beta density weight function.\n",
    "\n",
    "In ``chaospy``, these can all be constructed using [chaospy.generate_expansion()](../../api/chaospy.generate_expansion.rst).\n",
    "Hermite and normal distribution is showed above.\n",
    "The others can be created in the same way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:06.522398Z",
     "iopub.status.busy": "2021-05-18T10:56:06.521973Z",
     "iopub.status.idle": "2021-05-18T10:56:06.556503Z",
     "shell.execute_reply": "2021-05-18T10:56:06.556222Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1.0, q0, q0**2-0.33333, q0**3-0.6*q0])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uniform = chaospy.Uniform(-1, 1)\n",
    "legendre = chaospy.generate_expansion(3, uniform)\n",
    "legendre.round(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:06.558381Z",
     "iopub.status.busy": "2021-05-18T10:56:06.558132Z",
     "iopub.status.idle": "2021-05-18T10:56:06.591896Z",
     "shell.execute_reply": "2021-05-18T10:56:06.592171Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1.0, q0-1.0, q0**2-4.0*q0+2.0, q0**3-9.0*q0**2+18.0*q0-6.0])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exponential = chaospy.Exponential()\n",
    "laguerre = chaospy.generate_expansion(3, exponential)\n",
    "laguerre"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:06.594205Z",
     "iopub.status.busy": "2021-05-18T10:56:06.593898Z",
     "iopub.status.idle": "2021-05-18T10:56:06.627675Z",
     "shell.execute_reply": "2021-05-18T10:56:06.627397Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1.0, q0-3.0, q0**2-8.0*q0+12.0, q0**3-15.0*q0**2+60.0*q0-60.0])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha = 2\n",
    "gamma = chaospy.Gamma(alpha+1)\n",
    "generalized_laguerre = chaospy.generate_expansion(3, gamma)\n",
    "generalized_laguerre"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-18T10:56:06.629826Z",
     "iopub.status.busy": "2021-05-18T10:56:06.629564Z",
     "iopub.status.idle": "2021-05-18T10:56:06.664782Z",
     "shell.execute_reply": "2021-05-18T10:56:06.664506Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "polynomial([1.0, q0+0.14286, q0**2+0.22222*q0-0.11111,\n",
       "            q0**3+0.27273*q0**2-0.27273*q0-0.0303])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha_, beta_ = 2, 3\n",
    "beta = chaospy.Beta(alpha_+1, beta_+1, lower=-1, upper=1)\n",
    "jacobi = chaospy.generate_expansion(3, beta)\n",
    "jacobi.round(5)"
   ]
  }
 ],
 "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
}