Vinzenz Gregor Eck, Expert Analytics, Oslo
**Jacob Sturdy**, Department of Structural Engineering, NTNUDate: Jul 13, 2018
# ipython magic
%matplotlib notebook
%load_ext autoreload
%autoreload 2
# plot configuration
import matplotlib
import matplotlib.pyplot as plt
plt.style.use("ggplot")
# import seaborn as sns # sets another style
matplotlib.rcParams['lines.linewidth'] = 3
fig_width, fig_height = (7.0,5.0)
matplotlib.rcParams['figure.figsize'] = (fig_width, fig_height)
# font = {'family' : 'sans-serif',
# 'weight' : 'normal',
# 'size' : 18.0}
# matplotlib.rc('font', **font) # pass in the font dict as kwar
import chaospy as cp
import numpy as np
from monte_carlo import generate_sample_matrices_mc
from monte_carlo import calculate_sensitivity_indices_mc
unit_mmhg_pa = 133.3
unit_pa_mmhg = 1./unit_mmhg_pa
unit_cm2_m2 = 1. / 100. / 100.
unit_m2_cm2 = 1. / unit_cm2_m2
The elastic wall models are simplified algebraic functions $A(P)$ ([eck2015stochastic;@boileau_benchmark_2015]), which state the arterial lumen area $A$ as function of transmural pressure $P$.
For the calibration of the applied wall models, the wave speed in an arterial segment is required. The wave speed is given from fluid dynamics equations for one-dimensional arteries ([eck2015stochastic]):
with blood density $\rho= 1050\ [kg\ m^{-3}]$ and compliance $C(P) = dA / dP$.
The Quadratic area-pressure relationship ([sherwin2003]) is defined as:
where $\lambda$ is referred to as the stiffness coefficient and $A_s$ is the area at the reference pressure $P_s$.
The stiffness coefficient $\lambda$ is defined as:
The model is implemented in the following function:
# begin quadratic area model
def quadratic_area_model(pressure_range, samples):
"""
calculate the arterial vessel wall for a set pressure range
from 75 to 140 mmHg for a given set of reference wave speeds
area and pressure.
:param
pressure_range: np.array
pressure range for which to calculate the arterial area
samples: np.array with shape (4:n_samples)
sample matrix where
first indices correspond to
a_s: reference area
c_s: reference waves seed
p_s: reference pressure
rho: blood density
and n_samples to the number of samples
:return:
arterial pressure : np.array
of size n_samples
"""
pressure_range = pressure_range.reshape(pressure_range.shape[0], 1)
a_s, c_s, p_s, rho = samples
beta = 2*rho*c_s**2/np.sqrt(a_s)*a_s
#C_Laplace = (2. * ((P - Ps) * As / betaLaplace + np.sqrt(As))) * As / betaLaplace
return ((pressure_range - p_s)*a_s/beta + np.sqrt(a_s))**2.
The Logarithmic area-pressure relationship ([Hayashi1980]) is defined as:
where $\beta$ is called the stiffness coefficient and $A_s$ is the area at the reference pressure $P_s$.
The stiffness coefficient $\beta$ is defined as:
The model is implemented in the following function:
# begin exponential area model
def logarithmic_area_model(pressure_range, samples):
"""
calculate the arterial vessel wall for a set pressure range
from 75 to 140 mmHg for a given set of reference wave speeds
area and pressure.
:param
pressure_range: np.array
pressure range for which to calculate the arterial area
samples: np.array with shape (4:n_samples)
sample matrix where
first indices correspond to
a_s: reference area
c_s: reference waves seed
p_s: reference pressure
rho: blood density
and n_samples to the number of samples
:return:
arterial pressure : np.array
of size n_samples
"""
pressure_range = pressure_range.reshape(pressure_range.shape[0], 1)
a_s, c_s, p_s, rho = samples
beta = 2.0*rho*c_s**2./p_s
#C_hayashi = 2.0 * As / betaHayashi * (1.0 + np.log(P / Ps) / betaHayashi) / P
return a_s*(1.0 + np.log(pressure_range / p_s)/beta) ** 2.0
For a comparison of the wall models, we set the reference values to:
The two wall models give almost the same result around the reference area and pressure, for which the model was calibrated.
# start deterministic comparison
pressure_range = np.linspace(45, 180, 100) * unit_mmhg_pa
a_s = 5.12 * unit_cm2_m2
c_s = 6.25609258389
p_s = 100 * unit_mmhg_pa
rho = 1045.
plt.figure()
for model, name in [(quadratic_area_model, 'Quadratic model'), (logarithmic_area_model, 'Logarithmic model')]:
y_area = model(pressure_range, (a_s, c_s, p_s, rho))
plt.plot(pressure_range * unit_pa_mmhg, y_area * unit_m2_cm2, label=name)
plt.xlabel('Pressure [mmHg]')
plt.ylabel('Area [cm2]')
plt.legend()
plt.tight_layout()
# Create marginal and joint distributions
dev = 0.05
a_s = 5.12 * unit_cm2_m2
A_s = cp.Uniform(a_s * (1. - dev), a_s*(1. + dev))
c_s = 6.25609258389
C_s = cp.Uniform(c_s * (1. - dev), c_s*(1. + dev))
p_s = 100 * unit_mmhg_pa
P_s = cp.Uniform(p_s * (1. - dev), p_s*(1. + dev))
rho = 1045.
Rho = cp.Uniform(rho * (1. - dev), rho*(1. + dev))
jpdf = cp.J(A_s, C_s, P_s, Rho)
# scatter plots
pressure_range = np.linspace(45, 180, 100) * unit_mmhg_pa
Ns = 200
sample_scheme = 'R'
samples = jpdf.sample(Ns, sample_scheme)
for model, name, color in [(quadratic_area_model, 'Quadratic model', '#dd3c30'), (logarithmic_area_model, 'Logarithmic model', '#2775b5')]:
# evaluate the model for all samples
Y_area = model(pressure_range, samples)
plt.figure()
plt.title('{}: evaluations Y_area'.format(name))
plt.plot(pressure_range * unit_pa_mmhg, Y_area * unit_m2_cm2)
plt.xlabel('Pressure [mmHg]')
plt.ylabel('Area [cm2]')
plt.ylim([3.5,7.])
plt.figure()
plt.title('{}: scatter plots'.format(name))
for k in range(len(jpdf)):
plt.subplot(len(jpdf)/2, len(jpdf)/2, k+1)
plt.plot(samples[k], Y_area[0]*unit_m2_cm2, '.', color=color)
plt.ylabel('Area [cm2]')
plt.ylim([3.45, 4.58])
xlbl = 'Z' + str(k)
plt.xlabel(xlbl)
plt.tight_layout()
We apply the Monte Carlo method for both models as discussed before:
# start Monte Carlo
pressure_range = np.linspace(45, 180, 100) * unit_mmhg_pa
Ns = 50000
sample_method = 'R'
number_of_parameters = len(jpdf)
# 1. Generate sample matrices
A, B, C = generate_sample_matrices_mc(Ns, number_of_parameters, jpdf, sample_method)
fig_mean = plt.figure('mean')
ax_mean = fig_mean.add_subplot(1,2,1)
ax_dict = {}
figs = {}
for name in ['Quadratic model', 'Logarithmic model'] :
figs[name] = plt.figure(name)
ax = ax_dict[name]= figs[name].add_subplot(1,2,1)
print("\n Uncertainty measures (averaged)\n")
print('\n E(Y) | Std(Y) \n')
for model, name, color in [(quadratic_area_model, 'Quadratic model', '#dd3c30'),
(logarithmic_area_model, 'Logarithmic model', '#2775b5')]:
# evaluate the model for all samples
Y_A = model(pressure_range, A.T)
Y_B = model(pressure_range, B.T)
Y_C = np.empty((len(pressure_range), Ns, number_of_parameters))
for i in range(number_of_parameters):
Y_C[:, :, i] = model(pressure_range, C[i, :, :].T)
# calculate statistics
plotMeanConfidenceAlpha = 5.
Y = np.hstack([Y_A, Y_B])
expected_value = np.mean(Y, axis=1)
variance = np.var(Y, axis=1)
std = np.sqrt(np.var(Y, axis=1))
prediction_interval = np.percentile(Y, [plotMeanConfidenceAlpha / 2., 100. - plotMeanConfidenceAlpha / 2.], axis=1)
print('{:2.5f} | {:2.5f} : {}'.format(np.mean(expected_value)* unit_m2_cm2, np.mean(std)* unit_m2_cm2, name))
# 3. Approximate the sensitivity indices
# Better to centralize data before calculating sensitivities
expected_value_col = expected_value.copy()
expected_value_col.shape = (expected_value.shape[0],1)
Y_Ac = Y_A - expected_value_col
Y_Bc = Y_B - expected_value_col
Y_Cc = Y_C.transpose() - expected_value
Y_Cc = Y_Cc.transpose()
S, ST = calculate_sensitivity_indices_mc(Y_A, Y_B, Y_C)
Sc, STc = calculate_sensitivity_indices_mc(Y_Ac, Y_Bc, Y_Cc)
## Plots
fig = fig_mean
ax_mean = ax_mean
ax_mean.plot(pressure_range * unit_pa_mmhg, expected_value * unit_m2_cm2, label=name, color=color)
ax_mean.fill_between(pressure_range * unit_pa_mmhg, prediction_interval[0] * unit_m2_cm2,
prediction_interval[1] * unit_m2_cm2, alpha=0.3, color=color)
ax_mean.set_xlabel('Pressure [mmHg]')
ax_mean.set_ylabel('Area [cm2]')
ax_mean.legend()
fig.tight_layout()
fig = figs[name]
ax = ax_dict[name]
ax.set_title('{}: sensitvity indices'.format(name))
colorsPie = ['yellowgreen', 'gold', 'lightskyblue', 'lightcoral']
labels = ['Area', 'wave speed', 'pressure', 'density']
N = 4
ind = np.arange(N) # the x locations for the groups
width = 0.35 # the width of the bars
#firstOrderIndicesTimeAverage = S[:, 50]
#totalIndicesTimeAverage = ST[:, 50]
#rects1 = plt.bar(ind, firstOrderIndicesTimeAverage, width, color=colorsPie, alpha=0.5)
#rects2 = plt.bar(ind + width, totalIndicesTimeAverage, width, color=colorsPie, hatch='.')
firstOrderIndicesTimeAverage = np.mean(S*variance,axis=1)/np.mean(variance)
totalIndicesTimeAverage = np.mean(ST*variance,axis=1)/np.mean(variance)
rects1 = ax.bar(ind, firstOrderIndicesTimeAverage, width, color=colorsPie, alpha = 0.5)
rects2 = ax.bar(ind + width, totalIndicesTimeAverage, width, color=colorsPie, hatch = '.')
# add some text for labels, title and axes ticks
ax.set_ylabel('Si')
ax.set_xticks(ind + width)
ax.set_xticklabels(labels)
ax.set_xlim([-width, N + width])
ax.set_ylim([0, 1])
ax.spines['top'].set_visible(False)
ax.tick_params(axis='x', top='off')
ax.spines['right'].set_visible(False)
ax.tick_params(axis='y', right='off')
fig.tight_layout()
We apply the Polynomial Chaos method for both models as discussed before:
# start polynomial chaos
# orthogonal C_polynomial from marginals
polynomial_order = 3
Ns = 2*cp.bertran.terms(polynomial_order, len(jpdf))
print("\n Uncertainty measures (averaged)\n")
print('\n E(Y) | Std(Y) \n')
ax_mean = fig_mean.add_subplot(1,2,2)
ax_dict = {}
for name in ['Quadratic model', 'Logarithmic model'] :
figs[name] = plt.figure(name)
ax = ax_dict[name]= figs[name].add_subplot(1,2,2)
for model, name, color in [(quadratic_area_model, 'Quadratic model', '#dd3c30'),
(logarithmic_area_model, 'Logarithmic model', '#2775b5')]:
sample_scheme = 'R'
# create samples
samples = jpdf.sample(Ns, sample_scheme)
# create orthogonal polynomials
orthogonal_polynomials = cp.orth_ttr(polynomial_order, jpdf)
# evaluate the model for all samples
Y_area = model(pressure_range, samples)
# polynomial chaos expansion
polynomial_expansion = cp.fit_regression(orthogonal_polynomials, samples, Y_area.T)
# calculate statistics
plotMeanConfidenceAlpha = 5
expected_value = cp.E(polynomial_expansion, jpdf)
variance = cp.Var(polynomial_expansion, jpdf)
standard_deviation = cp.Std(polynomial_expansion, jpdf)
prediction_interval = cp.Perc(polynomial_expansion,
[plotMeanConfidenceAlpha / 2., 100 - plotMeanConfidenceAlpha / 2.],
jpdf)
print('{:2.5f} | {:2.5f} : {}'.format(np.mean(expected_value) * unit_m2_cm2, np.mean(std) * unit_m2_cm2, name))
# compute sensitivity indices
S = cp.Sens_m(polynomial_expansion, jpdf)
ST = cp.Sens_t(polynomial_expansion, jpdf)
fig = fig_mean
ax_mean = ax_mean
ax_mean.plot(pressure_range * unit_pa_mmhg, expected_value * unit_m2_cm2, label=name, color=color)
ax_mean.fill_between(pressure_range * unit_pa_mmhg, prediction_interval[0] * unit_m2_cm2,
prediction_interval[1] * unit_m2_cm2, alpha=0.3, color=color)
ax_mean.set_xlabel('Pressure [mmHg]')
ax_mean.set_ylabel('Area [cm2]')
ax_mean.legend()
fig.tight_layout()
fig = figs[name]
ax = ax_dict[name]
ax.set_title('{}: sensitvity indices'.format(name))
colorsPie = ['yellowgreen', 'gold', 'lightskyblue', 'lightcoral']
labels = ['Area', 'wave speed', 'pressure', 'density']
N = 4
ind = np.arange(N) # the x locations for the groups
width = 0.35 # the width of the bars
#firstOrderIndicesTimeAverage = S[:, 50]
#totalIndicesTimeAverage = ST[:, 50]
#rects1 = plt.bar(ind, firstOrderIndicesTimeAverage, width, color=colorsPie, alpha=0.5)
#rects2 = plt.bar(ind + width, totalIndicesTimeAverage, width, color=colorsPie, hatch='.')
firstOrderIndicesTimeAverage = np.mean(S*variance,axis=1)/np.mean(variance)
totalIndicesTimeAverage = np.mean(ST*variance,axis=1)/np.mean(variance)
rects1 = ax.bar(ind, firstOrderIndicesTimeAverage, width, color=colorsPie, alpha = 0.5)
rects2 = ax.bar(ind + width, totalIndicesTimeAverage, width, color=colorsPie, hatch = '.')
# add some text for labels, title and axes ticks
ax.set_ylabel('Si')
ax.set_xticks(ind + width)
ax.set_xticklabels(labels)
ax.set_xlim([-width, N + width])
ax.set_ylim([0, 1])
ax.spines['top'].set_visible(False)
ax.tick_params(axis='x', top='off')
ax.spines['right'].set_visible(False)
ax.tick_params(axis='y', right='off')
fig.tight_layout()