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

# # Exploiting Active Subspaces to Quantify Uncertainty in the Numerical Simulation of the HyShot II Scramjet
# 
# __Paul Constantine__, Colorado School of Mines, <pconstan@mines.edu>
# 
# __Michael Emory__, Stanford University, <memory@stanford.edu>
# 
# __Johan Larsson__, University of Maryland, <jola@umd.edu>
# 
# __Gianluca Iaccarino__, Stanford University, <jops@stanford.edu>
# 
# __Ryan Howard__, Colorado School of Mines, <ryhoward@mines.edu>
# 
# 
# <br>
# 
# In this notebook, we'll be summarizing the active subspaces-based uncertainty quantification presented in [[1]][R1].
# 
# A scramjet is an engine designed to intake air and combust fuel at supersonic speeds. Because physically prototyping and experimenting is expensive, we must use numerical simulations to explore the behavior of these machines. The simulation we use is based on the HyShot II vehicle and has 7 input parameters:
# 
# Parameter | Lower bound | Upper bound
# --- | --- | --- 
# Angle of attack | 2.6 | 4.6 
# Turbulence intensity | 0.1 | 1.9 
# Turbulence length scale | 122.5 | 367.5 
# Stagnation pressure | 1.6448e7 | 1.9012e7 
# Stagnation enthalpy | 3.0551e6 | 3.4280e6 
# Cowl transition location | 0.03 | 0.07 
# Ramp transition location | 0.087 | 0.203 
# 
# These parameters are used in our complex multiphysics simulation modelling the fluid flow in the scramjet. From our simulation, we compute our quantity of interest, a normalized integral of pressure over the end of the engine:
# 
# $$
# \frac{1}{vol(V)}\iiint\limits_VP\ dx\ dy\ dz
# $$
# 
# where V is a 3D region at the end of the scramjet and vol(V) is the volume of this region. This quantity can be used as a metric for judging performance and determining whether _unstart_ has occured (unstart is a catastrophic failure resulting in loss of thrust). For each combination of parameters used in our simulation, we compute this quantity for fuel-air equivalence ratios of .30 and .35.
# 
# We'll use the notation $f(\mathbf{x})$ to denote the quantity of interest as a function of normalized input parameters ($\mathbf x$). For the purposes of uncertainty quantification, we'd like to estimate bounds of $f$, identify sets of inputs that result in $f$ remaining below a certain threshold, and estimate a cumulative distribution function of $f$ given a density on $\mathbf x$.
# 
# The most straightforward way to achieve our objective would be to exhaustively sample the space of input parameters and use the resulting data to calculate what we want. However, because our input space is high-dimensional and the simulation is very expensive, this is infeasible. To compute what we want, we reduce the dimensionality of the problem by identifying the 'most important' direction in the input space and examining our quantity of interest along this direction. We use an _active subspace_ approach (see [[2]][R2]) based on the matrix
# 
# $$
# \mathbf C = \int\nabla f(\mathbf x)\nabla f(\mathbf x)^T \rho(\mathbf x)\ d\mathbf x = \mathbf W \Lambda\mathbf W^T
# $$
# 
# where $\rho$ is a probability density on $\mathbf x$ and $\mathbf W \Lambda\mathbf W^T$ is the eigendecomposition of $\mathbf C$. The active subspace is the space spanned by the first several eigenvectors, and we can use it to construct an approximation of $f$:
# 
# $$
# f(\mathbf x)\approx g(\mathbf W_1^T\mathbf x)
# $$
# 
# where $\mathbf W_1$ contains the first $n$ columns of $\mathbf W$ and $g$ is a function from $R^n$ to $R$. Since we don't have access to the gradient of $f$, we use the following least squares-based algorithm to identify the single most important direction in the parameter space:
# 1. Using $M$ data points, determine the parameters $\hat u_0, \hat u_1,\dots,\hat u_m$ of the linear approximation $f(\mathbf x)\approx \hat u_0 + \hat u_1x_1 + \cdots + \hat u_mx_m$ such that $\hat{\mathbf u} = argmin_{\mathbf u}||\mathbf{Au}-\mathbf f||_2^2$ where
# $$
# \mathbf A = \left[\begin{matrix}1 & \mathbf x_1^T\\\vdots & \vdots\\ 1 & \mathbf x_M^T\end{matrix}\right],\ \mathbf f = \left[\begin{matrix}f_1\\\vdots\\f_M\end{matrix}\right],\ \hat{\mathbf u} = \left[\begin{matrix}\hat u_0\\\vdots\\\hat u_m\end{matrix}\right]
# $$<br>
# 2. Compute $\mathbf w = \hat{\mathbf u}'/||\hat{\mathbf u}'||$ where $\hat{\mathbf u}' = [\hat u_1, \cdots, \hat u_m]^T$. $\mathbf w$ is the vector pointing in the most important direction in the parameter space.
# 
# To examine the effects of variability in $\mathbf w$, we use a bootstrap algorithm:
# 1. Choose the number $N$ of bootstrap replicates
# 2. For $k$ from 1 to $N$:
#     1. Let $\boldsymbol{\pi}^k = [\pi_1^k,\cdots,\pi_M^k]$ be integers from 1 to $M$ sampled uniformly with replacement.
#     2. Let $\mathbf f^k$ be the vector whose $i^{th}$ element is $f_{\pi^k}$ and $\mathbf A^k$ be the matrix whose $i^{th}$ row is $[1,\ \mathbf x_{\pi^k}^T]$.
#     3. Compute $\mathbf w^k$ using $\mathbf A = \mathbf A^k$ and $\mathbf f = \mathbf f^k$ as above.
#     
# ### References:
# 
# [[1]][R1] P.G. Constantine, M. Emory, J. Larsson, and G. Iaccarino. _Exploiting active subspaces to quantify uncertainty in the numerical simulation of the HyShot II scramjet_. J. Comput. Phys., 302, 1-20, 2015
# 
# [[2]][R2] P.G. Constantine, E. Dow, and Q. Wang. _Active subspaces in theory and practice: applications to kriging surfaces_. SIAM J. Sci. Comput., 36(4), A1500–A1524, 2014
# 
# [R1]: http://dx.doi.org/10.1016/j.jcp.2015.09.001
# [R2]: http://dx.doi.org/10.1137/130916138
# 
# ### Acknowledgments
# 
# This material is based upon work supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Number DE-SC-0011077.
# 
# <br>
# 
# We'll now code up the analysis using the [Python Active-subspaces Utility Library](https://github.com/paulcon/active_subspaces)

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
import numpy as np
import pandas as pn
import active_subspaces as ac
import matplotlib.pyplot as plt
np.set_printoptions(formatter={'all' : '{:<15}'.format})


# First we import the data and separate inputs (X) from outputs (F)

# In[2]:


df = pn.DataFrame.from_csv('HyShotII35.txt')
data = df.as_matrix()
#The last 2 rows (for the extreme values) are already normalized
X = data[:-2,0:7]
F = data[:,7:9]
labels = df.keys()
in_labels = labels[:7]
out_labels = labels[7:9]

#The data for equivalence ratio .30 is already normalized
df30 = pn.DataFrame.from_csv('HyShotII30.txt', header=None, index_col=None)
data30 = df30.as_matrix()
XX30 = data30[:,0:7]
F30 = data30[:,7]


# Now we normalize each of the input parameters to the interval [-1, 1] and separate our outputs

# In[3]:


#normalize the ratio = .35 data
xl = np.array([2.6, 0.1, 122.5, 1.6448e7, 3.0551e6, 0.6*0.05, 0.6*0.145])
xu = np.array([4.6, 1.9, 367.5, 1.9012e7, 3.4280e6, 1.4*0.05, 1.4*0.145])
XX = ac.utils.misc.BoundedNormalizer(xl, xu).normalize(X)
XX = np.vstack((XX, data[-2:,:7]))

f30 = F30[:]
out_label30 = out_labels[0]
f35 = F[:,1]
out_label35 = out_labels[1]


# Set up our subspaces. The last two rows of each data set were from simulations done after initial subspace analysis was conducted, for the purpose of verifying predicted extrema; we exclude those data points initially to mirror the process of discovering the subspace.

# In[4]:


ss30 = ac.subspaces.Subspaces()
ss35 = ac.subspaces.Subspaces()

#Compute OLS subspaces
ss30.compute(X=XX30[:-2,:], f=f30[:-2,np.newaxis], sstype='OLS', nboot=500)
ss35.compute(X=XX[:-2,:], f=f35[:-2,np.newaxis], sstype='OLS', nboot=500)


# Plot the quantity of interest against $\mathbf w^T\mathbf x$ for each equivalence ratio

# In[5]:


#active variable values associated with XX, XX30
y30 = np.dot(XX30, ss30.W1)
y35 = np.dot(XX, ss35.W1)

#plots of response (pressure) vs. active variable
ac.utils.plotters.sufficient_summary(y30[:-2], f30[:-2], out_label30)
ac.utils.plotters.sufficient_summary(y35[:-2], f35[:-2], out_label35)


# We can see a very tight univariate trend in each case (see also Figure 7 of [[1]][R1]); we can use this low-dimensional structure to perform our uncertainty quantification.
# [R1]: http://dx.doi.org/10.1016/j.jcp.2015.09.001

# In[6]:


#Display the components of w
print "For ratio .30, the components of w are:"
print np.hstack((np.round(ss30.W1, 4), in_labels[:,np.newaxis]))
print "\nFor ratio .35, the components of w are:"
print np.hstack((np.round(ss35.W1, 4), in_labels[:,np.newaxis]))


# The components of $\mathbf w$ tell us how sensitive the quantity of interest is to each parameter. From the above output, we can see that Angle of Attack, Turbulence Intensity, and Stagnation conditions (pressure and enthalpy) dominate the active subspace, and (perhaps by coincidence) the components are roughly the same between the fuel-air equivalence ratios, except for an exchange of weight between stagnation pressure and stagnation enthalpy.
# 
# We'll now demonstrate the bootstrap for equivalence ratio .30.

# In[7]:


N = int(1e5) #Number of bootstrap replicates
M,m = XX30.shape
w_boot = np.empty(shape=(m, N)) #matrix containing replicates
for k in range(N):
    mask = np.random.randint(0, M, size=(M,))
    f = np.matrix(f30[mask].reshape(M,1))
    A = np.matrix(np.hstack((np.ones((M,1)), XX30[mask])))
    u = (A.T.dot(A)).I.dot(A.T).dot(f)
    w_boot[:,k] = np.squeeze(u[1:])/np.linalg.norm(np.squeeze(u[1:]))
w_br = np.empty(shape=(m, 2)) #min-to-max range of bootstrap values
w_br[:,0] = np.squeeze(np.apply_over_axes(np.amin, w_boot, axes=(1)))
w_br[:,1] = np.squeeze(np.apply_over_axes(np.amax, w_boot, axes=(1)))


# In[8]:


#plot eigenvector components with bootstrap ranges
ac.utils.plotters.eigenvectors(ss30.W1, w_br, in_labels, out_label30)
plt.figure(figsize=(14, 5*4))
for k in range(m): #plot histograms
    plt.subplot(4, 2, k+1)
    plt.grid(True)
    plt.xlabel(in_labels[k])
    plt.xlim([-1, 1])
    plt.vlines(ss30.W1[k], 0, .5, 'r', lw=2)
    plt.hist(w_boot[k,:], bins=10,\
             weights=1.0*np.ones_like(w_boot[k,:])/float(len(w_boot[k,:])))
plt.tight_layout()


# We can see the histograms are tightly concentrated around our estimated values (the red vertical lines), indicating that our estimates are close to their true values. Because we have only a few data points for the equivalence ratio of .35, the bootstrap ranges would be much wider in that case, but our confirmation of the results in the .30 case and the similarities between the two cases lead us to be confident for both ratios (see also Figures 8 and 9 of [[1]][R1]).
# [R1]: http://dx.doi.org/10.1016/j.jcp.2015.09.001
# 
# Now that we are confident in our subspace, we can begin the uncertainty quantification (estimating ranges on $f$, identifying safe operating conditions, and estimating a cdf on $f$). Range estimation is sraightforward given the monotonic nature of the subspace: we can maximize $f$ by maximizing the value of our active variable, $\mathbf w^T\mathbf x$, over our input space, and running our simulation at these values of $\mathbf x$ (four runs: 1 max and 1 min for each equivalence ratio). We plot the initial data (blue dots) and the acquired extreme values (red squares) below (see also Figure 7 of [[1]][R1]).
# [R1]: http://dx.doi.org/10.1016/j.jcp.2015.09.001

# In[9]:


#Plot initial data (blue dots) and extrema (red squares)
plt.figure(figsize=(7, 7))
plt.plot(y30[:-2], f30[:-2], 'bo', y30[-2:], f30[-2:], 'rs', ms=12); plt.grid(True)
plt.xlabel('Active Variable'); plt.ylabel(out_label30)
plt.figure(figsize=(7, 7))
plt.plot(y35[:-2], f35[:-2], 'bo', y35[-2:], f35[-2:], 'rs', ms=12); plt.grid(True)
plt.xlabel('Active Variable'); plt.ylabel(out_label35)


# Suppose now that the engine remained safe so long as exit pressure is below 2.8 and becomes unsafe when pressure rises above this. We might reasonably be interested in what regions of our parameter space result in safe operating conditions and which don't. To constrain the exit pressure, we can construct response surfaces from the data points $\{(y_i,f_i)\}$, where $y_i$ are the active variable values and $f_i$ are the values of the quantity of interest, and use them to place a bound on how high the active variable can become before the engine becomes unsafe, and map the safe active-variable region to the original parameter space.
# 
# Below, we plot our response surfaces (obtained from OLS on a quadratic model) and an upper prediction limit that is 2.33 standard errors above our mean response (see also Figure 10 of [[1]][R1]). We also find the active variable value that marks the boundary between safe and unsafe operation. The shaded regions denote safe operation.
# [R1]: http://dx.doi.org/10.1016/j.jcp.2015.09.001

# In[10]:


#Quadratic response surfaces
surface30 = ac.utils.response_surfaces.PolynomialApproximation(2)
surface35 = ac.utils.response_surfaces.PolynomialApproximation(2)

#Training the surfaces with our data
surface30.train(y30, f30[:,None])
surface35.train(y35, f35[:,None])

M30 = len(y30); M35 = len(y35)#Number of data points for each ratio

#sums of squares required for computing standard error
sf30 = np.sqrt(sum((ac.utils.misc.atleast_2d_col(f30) -\
                    surface30.predict(y30)[0])**2)/(M30 - 3))
sy30 = sum((y30 - y30.mean())**2)
#predicted response using mean response + 2.33 std errors
predlim30 = lambda x: surface30.predict(x)[0] + 2.33*sf30*np.sqrt(1 +\
                      1/M30 + (x - y30.mean())**2/sy30)

sf35 = np.sqrt(sum((ac.utils.misc.atleast_2d_col(f35) -\
                    surface35.predict(y35)[0])**2)/(M35 - 3))
sy35 = sum((y35 - y35.mean())**2)
predlim35 = lambda x: surface35.predict(x)[0] + 2.33*sf35*np.sqrt(1 +\
                      1/M35 + (x - y35.mean())**2/sy35)

x = ac.utils.misc.atleast_2d_col(np.linspace(-2, 3, 800))
#maximum Active variable value s.t. f(x) <= 2.8 according to prediction interval
avmaxval30 = x[np.argmin(abs(predlim30(x) - 2.8))]
plt.figure(figsize=(7, 7))
plt.plot(x, surface30.predict(x)[0], 'b-', x, predlim30(x), 'r:', lw=2)
plt.plot(x, 2.8*np.ones_like(x), 'k-')
plt.vlines(avmaxval30, 2, 3, 'k')
plt.fill_between(x[x <= avmaxval30], 2, 2.8, alpha=.3)
plt.ylim([2, 3])
plt.xlabel('Active Variable'); plt.ylabel(out_label30)
plt.legend(['mean response', 'prediction limit'], loc=2)

avmaxval35 = x[np.argmin(abs(predlim35(x) - 2.8))]
plt.figure(figsize=(7, 7))
plt.plot(x, surface35.predict(x)[0], 'b-', x, predlim35(x), 'r:', lw=2)
plt.plot(x, 2.8*np.ones_like(x), 'k-')
plt.vlines(avmaxval35, 2, 3.6, 'k')
plt.fill_between(x[x <= avmaxval35], 2, 2.8, alpha=.3)
plt.ylim([2, 3.6]); plt.xlim([-2, 3])
plt.xlabel('Active Variable'); plt.ylabel(out_label35)
plt.legend(['mean response', 'prediction limit'], loc=2)


# The set of normalized parameters resulting in safe operation is $S = \{\mathbf x : \mathbf w^T\mathbf x \leq y_{max},\ -\mathbf 1 \leq \mathbf x\leq \mathbf 1\}$, where $y_{max}$ is the maximum safe active variable value. We will find ranges of input parameters that guarantee safety by finding the largest 'box' that can fit within the set $S$ by solving the optimization problem: maximize$_{\mathbf x}\prod_{i=1}^m|x_i-x_{i,min}|,\ s.t.\ \mathbf x\in S$, where the $x_{i,min}$ are the components of xmin above (the minimizer of the active variable, $\mathbf w^T\mathbf x$).

# In[11]:


import scipy.optimize as op

#Functions to minimize (-log(prod(x_i-x_{i,min})))
minfun30 = lambda x: -1.0*np.sum(np.log(abs(x - XX30[-2,:].reshape(x.shape))))
minfun35 = lambda x: -1.0*np.sum(np.log(abs(x - XX[-2,:].reshape(x.shape))))
#constraint functions of the form c(x) >= 0
cfun30 = lambda x: avmaxval30 - x.dot(ss30.W1)
cfun35 = lambda x: avmaxval35 - x.dot(ss35.W1)
#constraint dictionaries passed to minimizer
cdict30 = {'type': 'ineq', 'fun': cfun30}
cdict35 = {'type': 'ineq', 'fun': cfun35}

#Normalized parameter values defining the edge of the largest box assuring safe operation
xop30 = op.minimize(minfun30, XX30[-1,:],\
                    bounds=zip(-1.0*np.ones_like(xl), 1.0*np.ones_like(xl)),\
                    constraints=cdict30)
xop35 = op.minimize(minfun35, XX[-1,:],\
                    bounds=zip(-1.0*np.ones_like(xl), 1.0*np.ones_like(xl)),\
                    constraints=cdict35)

#un-normalizing inputs to original parameter space
xop30 = ac.utils.misc.BoundedNormalizer(xl, xu).unnormalize(xop30.x[None,:]).reshape(m, 1)
xop35 = ac.utils.misc.BoundedNormalizer(xl, xu).unnormalize(xop35.x[None,:]).reshape(m, 1)
min30 = ac.utils.misc.BoundedNormalizer(xl, xu).unnormalize(XX30[-2,:][None,:]).reshape(m, 1)
min35 = ac.utils.misc.BoundedNormalizer(xl, xu).unnormalize(XX[-2,:][None,:]).reshape(m, 1)

#sorting bounds for display convenience
A30 = np.round(np.sort(np.hstack((min30, xop30)), axis=1), 4)
A35 = np.round(np.sort(np.hstack((min35, xop35)), axis=1), 4)


# In[12]:


#Display the safe operating ranges
print "For ratio .30, the ranges providing safe operation are:"
for i in range(7):
    print "{:0.2e} \t {:0.2e}\t {}".format(A30[i,0],A30[i,1],in_labels[i])
print "\nFor ratio .35, the ranges providing safe operation are:"
for i in range(7):
    print "{:0.2e} \t {:0.2e}\t {}".format(A35[i,0],A35[i,1],in_labels[i])


# Lastly, we can construct a cumulative distribution of $f$ by randomly sampling the (normalized) parameter space, computing the active variable, and computing the response surface value at each point (see also Figure 11 of [[1]][R1]).
# [R1]: http://dx.doi.org/10.1016/j.jcp.2015.09.001

# In[13]:


N = int(1e5) #Number of random samples
x = np.random.uniform(-1, 1, (N, m))
a30 = x.dot(ss30.W1); a35 = x.dot(ss35.W1) #active variable values
s30 = surface30.predict(a30)[0]; s35 = surface35.predict(a35)[0] #response surface values
x = np.linspace(1.5, 3.5, 300)
c30 = np.vectorize(lambda x: sum(s30 <= x)/(1.0*N)) #CDF functions
c35 = np.vectorize(lambda x: sum(s35 <= x)/(1.0*N))

plt.plot(x, c30(x), 'k-')
plt.ylim([0, 1.1])
plt.xlabel(out_label30)
plt.vlines(f30[-2:], 0, 1.1, 'k', 'dashed')

plt.figure()
plt.plot(x, c35(x), 'k-')
plt.ylim([0, 1.1])
plt.xlabel(out_label35)
plt.vlines(f35[-2:], 0, 1.1, 'k', 'dashed')