-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathborehole_cp.py
129 lines (99 loc) · 3.4 KB
/
borehole_cp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
"""
Filename: borehole_cp.py
Author(s): Takafumi Usui
E-mail: [email protected]
Description:
"""
import numpy as np
import chaospy as cp
import numpoly
import matplotlib.pyplot as plt
from matplotlib import rc
# Use TeX font
rc('font', **{'family': 'sans-serif', 'serif': ['Helvetica']})
rc('text', usetex=True)
# Font size
plt.rcParams["font.size"] = 12
plt.rcParams["axes.labelsize"] = 14
plt.rcParams["axes.titlesize"] = 14
plt.rcParams["legend.title_fontsize"] = 14
def borehole(x):
rw = x[0]
r = x[1]
Tu = x[2]
Hu = x[3]
Tl = x[4]
Hl = x[5]
L = x[6]
Kw = x[7]
r_rw_log = np.log(r / rw)
numerator = 2 * np.pi * Tu * (Hu - Hl)
denominator = r_rw_log * (
1 + (2 * L * Tu) / (r_rw_log * rw**2 * Kw) + (Tu / Tl))
return numerator / denominator
# rw_dist = cp.Normal(0.10, 0.0161812)
rw_dist = cp.Uniform(0.05, 0.15)
# r_dist = cp.LogNormal(7.71, 1.0056)
r_dist = cp.Uniform(100, 50000)
Tu_dist = cp.Uniform(63070, 115600)
Hu_dist = cp.Uniform(990, 1110)
Tl_dist = cp.Uniform(63.1, 116)
Hl_dist = cp.Uniform(700, 820)
L_dist = cp.Uniform(1120, 1680)
Kw_dist = cp.Uniform(9855, 12045)
# import ipdb; ipdb.set_trace()
joint_dist = cp.J(
rw_dist, r_dist, Tu_dist, Hu_dist, Tl_dist, Hl_dist, L_dist, Kw_dist)
poly_order = 4
phi, coeffs = cp.expansion.stieltjes(
poly_order, joint_dist, normed=True, graded=True, reverse=True,
retall=True, cross_truncation=1.0)
N_EDpoints = 3 * len(phi)
EDpoints = joint_dist.sample(N_EDpoints, 'latin_hypercube')
evals = np.array([borehole(x) for x in EDpoints.T])
M_hat, fourier_coeffs = cp.fit_regression(phi, EDpoints, evals, retall=1)
# Total variance using the Fourier coefficients
D_hat = np.sum(fourier_coeffs[1:]**2, axis=0)
# Information matrix
Ainfo = phi(*EDpoints).T
# Based on the least square minimization
fourier_coeffs_hat = np.linalg.inv(Ainfo.T @ Ainfo) @ Ainfo.T @ evals
M_hathat = numpoly.sum((phi * fourier_coeffs_hat.T), -1)
# --------------------------------------------------------------------------- #
# Compute the leave-one-out error (errLOO)
# errLOO <= 10^{-2}
# --------------------------------------------------------------------------- #
h_diag = np.diagonal(Ainfo @ np.linalg.inv(Ainfo.T @ Ainfo) @ Ainfo.T)
errLOO = 1 / N_EDpoints * np.sum(
((evals - M_hathat(*EDpoints)) / (1 - h_diag))**2)
print(r'errLOO error is {}'.format(errLOO))
import ipdb; ipdb.set_trace()
N_param = EDpoints.shape[0]
# Get a multi-index
alpha = cp.glexindex(
start=0, stop=poly_order+1, dimensions=N_param, cross_truncation=1.0,
graded=True, reverse=True).T
Sens_t_hat = np.empty(N_param)
for idx in range(N_param):
index = (alpha[idx, :] > 0)
Sens_t_hat[idx] = np.sum(fourier_coeffs[index]**2, axis=0) / D_hat
Sens_m_hat = np.empty(N_param)
for idx in range(N_param):
index = (alpha[idx, :] > 0) & (alpha.sum(0) == alpha[idx, :])
Sens_m_hat[idx] = np.sum(fourier_coeffs[index]**2, axis=0) / D_hat
xlabels = [r'$r_{w}$', r'$r$', r'$T_{u}$', r'$H_{u}$', r'$T_{l}$', r'$H_{l}$',
r'$L$', r'$K_{w}$']
fig, ax = plt.subplots(figsize=(9, 6))
ax.bar(xlabels, Sens_t_hat, color='royalblue', edgecolor='black')
plt.ylim([0, None])
plt.ylabel(r"$S^{T}_{i}$")
plt.show()
fig, ax = plt.subplots(figsize=(9, 6))
ax.bar(xlabels, Sens_m_hat, color='royalblue', edgecolor='black')
plt.ylim([0, None])
plt.ylabel(r"$S_{i}$")
plt.show()
import ipdb; ipdb.set_trace()