#!/usr/bin/env python
# coding: utf-8

# ## Orthogonality
# 
# The core idea of polynomial chaos expansions is that the polynomials
# used as an expansion are all mutually orthogonal. The relation is
# typically written mathematically as:
# 
# $$\left\langle \Phi_n, \Phi_m \right\rangle = 0 \qquad n \neq m$$
# 
# In practice this relation is instead expressed by the equivalent
# notation using expected values:
# 
# $$\mbox E\left(\Phi_n \Phi_m\right) = 0 \qquad n \neq m$$
# 
# In `chaospy` this property can be tested by taking the outer product of
# two expansions, and checking if the expected value of the resulting
# matrix is diagonal. For example, for a basic monomial:

# In[1]:


import chaospy

monomial = chaospy.monomial(4)
monomial


# In[2]:


monomial2 = chaospy.outer(monomial, monomial)
monomial2


# In[3]:


normal = chaospy.Normal(0, 1)
chaospy.E(monomial2, normal)


# In other words, the basic monomial (beyond polynomial order 1) are not
# orthogonal.
# 
# But if we replace the basic monomial with an explicit orthogonal
# polynomial constructor, we get:

# In[4]:


hermite = chaospy.generate_expansion(3, normal)
hermite


# In[5]:


hermite2 = chaospy.outer(hermite, hermite)
chaospy.E(hermite2, normal).round(15)


# A fully diagonal matrix, which implies all the polynomials in the
# expansion are mutually orthogonal.

# ### Multivariate expansions
# 
# Multivariate orthogonal expansion are (usually) created by doing a tensor product of univariate expansions together.
# To illustrate how this work, consider the distribution introduced in the [problem formulation](../main_usage/problem_formulation.ipynb):

# In[6]:


from problem_formulation import joint

joint


# Extracting the marginal density we can construct both one-dimensional expansions:

# In[7]:


expansion0 = chaospy.generate_expansion(2, joint[0])
expansion0


# In[8]:


expansion1 = chaospy.generate_expansion(2, joint[1])
expansion1.round(5)


# When constructing a multivariate expansion, it is canonical to truncate the expansion at order and graded lexicographical sorting:

# In[9]:


chaospy.generate_expansion(2, joint).round(5)


# See [chaospy.generate_expansion()](../../api/chaospy.generate_expansion.rst) for variations in truncations and sorting.

# ### Wiener-Askey scheme
# 
# Polynomial chaos expansion often assume that the polynomial expansion used is of the Wiener-Askey scheme verity.
# The reason for this is that the expansion in the scheme correspond to orthogonality with respect to some standard probability distribution.
# These include:
# 
# * Hermite polynomials which are orthogonal with a normal density weight function.
# * Legendre polynomials which are orthogonal with a uniform density weight function.
# * Laguerre polynomials which are orthogonal with a exponential density weight function.
# * Generalized Laguerre polynomials which are orthogonal with a gamma density weight function.
# * Jacobi polynomials which are orthogonal with a beta density weight function.
# 
# In ``chaospy``, these can all be constructed using [chaospy.generate_expansion()](../../api/chaospy.generate_expansion.rst).
# Hermite and normal distribution is showed above.
# The others can be created in the same way:

# In[10]:


uniform = chaospy.Uniform(-1, 1)
legendre = chaospy.generate_expansion(3, uniform)
legendre.round(5)


# In[11]:


exponential = chaospy.Exponential()
laguerre = chaospy.generate_expansion(3, exponential)
laguerre


# In[12]:


alpha = 2
gamma = chaospy.Gamma(alpha+1)
generalized_laguerre = chaospy.generate_expansion(3, gamma)
generalized_laguerre


# In[13]:


alpha_, beta_ = 2, 3
beta = chaospy.Beta(alpha_+1, beta_+1, lower=-1, upper=1)
jacobi = chaospy.generate_expansion(3, beta)
jacobi.round(5)