STZ_Optimization/EGO_.py at master · dalixwill/STZ_Optimization · GitHub

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
import numpy as np
from Kriging_ import predictor
from gfunctions import *
import chaospy as cp
from scipy.stats import norm
from Grassmann_ import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import colormaps as cmaps


Tref = np.loadtxt('./tem.24', skiprows=1)
'''
#  Plot Temperature
fig = plt.figure()
ax = fig.add_subplot(111)
pos = ax.imshow(Tref, cmap=cmaps.parula, interpolation='bicubic', \
                origin='lower', )
ax.set_xlabel('x')
ax.set_ylabel('y')
fig.colorbar
#plt.show()
'''

u0, s0, v0 = svd(Tref, 0)

dim = 2
N = 21
# define the distribution and the boundaries of beta and U0
beta = cp.Uniform(lo=16, up=18)
U0 = cp.Uniform(lo=-3.6035, up=-3.6021)

# Generate N samples using latin hypercube sampling
X = np.zeros([N, dim])
X[:, 0] = beta.sample(N, rule="L")
X[:, 1] = U0.sample(N, rule="L")

# Calculate the response ffor each pair of (beta, U0)
Y = STZ_Darius(X, u0, 0).ravel()

# Set up the Kriging (Gaussian process) surrogate
gp = predictor(X, Y)

# Generate N_test random realizations of (beta, U0)
N_test = 200
X_test = np.zeros([N_test, dim])
X_test[:, 0] = beta.sample(N_test)
X_test[:, 1] = U0.sample(N_test)

# Use the surrogate in order to estimate the response at the N_test points
y_hat, s2 = gp.predict(X_test, return_std=True)

# Apply the EGO algorithm to find the next sampling point
fmin = min(y_hat)


EI= np.zeros(X_test.shape[0])
for i in range(N_test):
   EI[i] = (fmin-y_hat[i])*norm.cdf((fmin-y_hat[i])/np.sqrt(s2[i]), loc=0.0, scale=1.0) + np.sqrt(s2[i])*norm.pdf((fmin-y_hat[i])/np.sqrt(s2[i]), loc=0.0, scale=1.0)

# Find the position  of the maximum value of EI
J = np.argmax(EI)

'''
x0 = np.sort(X_test[:, 0])
index0 = np.argsort(x0)
x1 = np.sort(X_test[:, 1])
index1 = np.argsort(x1)

plt.figure()
plt.plot(x0, EI[index0], 'k', linewidth=2)

plt.figure()
plt.plot(x1, EI[index1], 'k', linewidth=2)
plt.show()
'''
Iter = 0
while Iter < 100:    # Generate 100 samples
    Iter = Iter +1
    # Estimate the continuum model for the realization that corresponds to the maximum value of EI
    Y_add = STZ_Darius(X_test[J, :], u0, 1).ravel()

    # Add this realization to the initial N points
    Y = np.concatenate((Y, Y_add), axis = 0)
    X = np.concatenate((X, X_test[J, :].reshape(1, dim)))

    # Set up the Kriging model again using also the added point
    gp = predictor(X, Y)

    N_test = 400
    X_test = np.zeros([N_test, dim])
    X_test[:, 0] = beta.sample(N_test, rule="L")
    X_test[:, 1] = U0.sample(N_test, rule="L")

    y_hat, s2 = gp.predict(X_test, return_std=True)
    fmin = min(y_hat)


    EI= np.zeros(X_test.shape[0])
    for i in range(N_test):
        EI[i] = (fmin-y_hat[i])*norm.cdf((fmin-y_hat[i])/np.sqrt(s2[i]), loc=0.0, scale=1.0) + np.sqrt(s2[i])*norm.pdf((fmin-y_hat[i])/np.sqrt(s2[i]), loc=0.0, scale=1.0)

    EIsort = np.sort(EI)
    J = np.argmax(EI)

    print(X_test[J, :])


np.savetxt('{0}.txt'.format('beta_U0_modeleval'), X)
np.savetxt('{0}.txt'.format('beta_U0_Krigeval'), X_test)
np.savetxt('{0}.txt'.format('Gr_dist_Krigeval'), y_hat)
np.savetxt('{0}.txt'.format('Gr_distance_modeleval'), Y)