<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Pseudo-spectral projection
Pseudo-spectral projection method is one of two non-intrusive polynomial
chaos expansion methods. (The other being point collocation method.) In a
nutshell it can be performed as follows
Step1:  Evaluating model solver
Like in the case of problem formulation again,
evaluation is straight forward
Step2:  Expansion of orthogonal polynomials
Unlike point collocation method, the polynomials
in pseudo-spectral projection needs to be orthogonal for the method to work.
This can be achieved by using the
chaospy.generate_expansion()
function
Step3:  Fourier coefficients
Polynomial chaos expansion assumes that we can model function $u$ using an
expansion
Step4:  Note that if the Fourier coefficients are needed, then they are still
accessible by passing retall=True to
chaospy.fit_quadrature().
Descriptive statistics
The expected value and variance is calculated as follows
Step5:  Error analysis
It is hard to assess how well these models are doing from the final
estimation alone. They look about the same. So to compare results, we do
error analysis. To do so, we use the reference analytical solution and error
function as defined in  problem formulation.
Step6:  The analysis can be performed as follows
<ASSISTANT_TASK:>
Python Code:
import chaospy
from problem_formulation import joint
gauss_quads = [
    chaospy.generate_quadrature(order, joint, rule="gaussian")
    for order in range(1, 8)
]
sparse_quads = [
    chaospy.generate_quadrature(
        order, joint, rule=["genz_keister_24", "clenshaw_curtis"], sparse=True)
    for order in range(1, 5)
]
from matplotlib import pyplot
pyplot.rc("figure", figsize=[12, 4])
nodes, weights = gauss_quads[5]
pyplot.subplot(121)
pyplot.title("Gaussian")
pyplot.scatter(*nodes, s=weights*2e3)
nodes, weights = sparse_quads[3]
idx = weights > 0
pyplot.subplot(122)
pyplot.title("sparse-grid")
pyplot.scatter(*nodes[:, idx], s=weights[idx]*2e3)
pyplot.scatter(*nodes[:, ~idx], s=-weights[~idx]*2e3, color="grey")
pyplot.show()
Explanation: Pseudo-spectral projection
Pseudo-spectral projection method is one of two non-intrusive polynomial
chaos expansion methods. (The other being point collocation method.) In a
nutshell it can be performed as follows:
Generate nodes $Q_1, ..., Q_N$ and weights $W_1, ..., W_N$ from a
   quadrature integration scheme.
Use nodes to create model evaluations $U_1, ..., U_N$.
Select an expansion of orthogonal polynomials $P_1, ..., P_M$.
Estimate Fourier coefficients $c_1, ..., c_M$, creating model
   approximation $\hat u(q) = \sum_i c_i P_i$.
Perform model analysis on approximation $\hat u(q)$ as a proxy for the
   real model.
Let us go through the steps in more detail.
Generate quadrature nodes and weights
In pseudo spectral projection, the nodes used to evaluate the model have to
be taken from quadrature integration scheme.
For example, we can choose full tensor-grid with optimal Gaussian quadrature,
and Smolyak sparse-grid with Genz-Keister and Clenshaw-Curtis quadrature:
End of explanation
from problem_formulation import model_solver, coordinates
import numpy
gauss_evals = [
    numpy.array([model_solver(node) for node in nodes.T])
    for nodes, weights in gauss_quads
]
sparse_evals = [
    numpy.array([model_solver(node) for node in nodes.T])
    for nodes, weights in sparse_quads
]
pyplot.subplot(121)
pyplot.plot(coordinates, gauss_evals[6].T, alpha=0.3)
pyplot.title("Gaussian")
pyplot.subplot(122)
pyplot.plot(coordinates, sparse_evals[3].T, alpha=0.3)
pyplot.title("sparse-grid")
pyplot.show()
Explanation: Evaluating model solver
Like in the case of problem formulation again,
evaluation is straight forward:
End of explanation
expansions = [chaospy.generate_expansion(order, joint)
              for order in range(1, 10)]
expansions[0].round(10)
Explanation: Expansion of orthogonal polynomials
Unlike point collocation method, the polynomials
in pseudo-spectral projection needs to be orthogonal for the method to work.
This can be achieved by using the
chaospy.generate_expansion()
function:
End of explanation
gauss_model_approx = [
    chaospy.fit_quadrature(expansion, nodes, weights, evals)
    for expansion, (nodes, weights), evals in zip(expansions, gauss_quads, gauss_evals)
]
sparse_model_approx = [
    chaospy.fit_quadrature(expansion, nodes, weights, evals)
    for expansion, (nodes, weights), evals in zip(expansions, sparse_quads, sparse_evals)
]
model_approx = gauss_model_approx[4]
nodes, _ = gauss_quads[4]
evals = model_approx(*nodes)
pyplot.subplot(121)
pyplot.plot(coordinates, evals, alpha=0.3)
pyplot.title("Gaussian")
model_approx = sparse_model_approx[1]
nodes, _ = sparse_quads[1]
evals = model_approx(*nodes)
pyplot.subplot(122)
pyplot.plot(coordinates, evals, alpha=0.3)
pyplot.title("sparse-grid")
pyplot.show()
Explanation: Fourier coefficients
Polynomial chaos expansion assumes that we can model function $u$ using an
expansion:
$$
u(q) \approx \hat u(q) = \sum_{i=1}^N c_i P_i
$$
Since $P_1, \dots, P_N$ are orthogonal, the optimal selection for $c_1,
\dots, c_N$, the Fourier coefficients becomes:
$$
  c_i = \frac{\left\langle u, P_i\right\rangle}{
              \left\langle P_i, P_i\right\rangle}
$$
These coefficients can be estimated in pseudo-spectral projection using
quadrature integration. In chaospy we calculate the model approximation
directly through the
chaospy.fit_quadrature()
function:
End of explanation
expected = chaospy.E(gauss_model_approx[-2], joint)
std = chaospy.Std(gauss_model_approx[-2], joint)
expected[:4].round(4), std[:4].round(4)
pyplot.rc("figure", figsize=[6, 4])
pyplot.xlabel("coordinates")
pyplot.ylabel("model approximation")
pyplot.fill_between(
    coordinates, expected-2*std, expected+2*std, alpha=0.3)
pyplot.plot(coordinates, expected)
pyplot.show()
Explanation: Note that if the Fourier coefficients are needed, then they are still
accessible by passing retall=True to
chaospy.fit_quadrature().
Descriptive statistics
The expected value and variance is calculated as follows:
End of explanation
from problem_formulation import error_in_mean, error_in_variance
error_in_mean(expected), error_in_variance(std**2)
Explanation: Error analysis
It is hard to assess how well these models are doing from the final
estimation alone. They look about the same. So to compare results, we do
error analysis. To do so, we use the reference analytical solution and error
function as defined in  problem formulation.
End of explanation
gauss_sizes = [len(weights) for _, weights in gauss_quads]
eps_gauss_mean = [
    error_in_mean(chaospy.E(model, joint))
    for model in gauss_model_approx
]
eps_gauss_var = [
    error_in_variance(chaospy.Var(model, joint))
    for model in gauss_model_approx
]
sparse_sizes = [len(weights) for _, weights in sparse_quads]
eps_sparse_mean = [
    error_in_mean(chaospy.E(model, joint))
    for model in sparse_model_approx
]
eps_sparse_var = [
    error_in_variance(chaospy.Var(model, joint))
    for model in sparse_model_approx
]
pyplot.rc("figure", figsize=[12, 4])
pyplot.subplot(121)
pyplot.title("Error in mean")
pyplot.loglog(gauss_sizes, eps_gauss_mean, "o-", label="Gaussian")
pyplot.loglog(sparse_sizes, eps_sparse_mean, "o--", label="sparse")
pyplot.legend()
pyplot.subplot(122)
pyplot.title("Error in variance")
pyplot.loglog(gauss_sizes, eps_gauss_var, "o-", label="Gaussian")
pyplot.loglog(sparse_sizes, eps_sparse_var, "o--", label="sparse")
pyplot.show()
Explanation: The analysis can be performed as follows:
End of explanation
<END_TASK>