1D Example#
Todo
Write this section.
Load modules and start Fesslix engine#
import fesslix as flx
flx.load_engine()
import fesslix.gpr
import fesslix.tools
import matplotlib.pyplot as plt
import numpy as np
Random Number Generator: MT19937 - initialized with rand()=1307296760;
Random Number Generator: MT19937 - initialized with 1000 initial calls.
Create random samples#
N_smpls_max = 1000
eps_sd = 0.2 # »» 1.1 gives problems
## ------------------------------------------
## Uncertainties of model input
## ------------------------------------------
config_rv_X = { 'name':'X', 'type':'uniform', 'a':0., 'b':10. }
config_rv_noise = { 'name':'eps', 'type':'normal', 'mu':0., 'sd':eps_sd }
rv_set = flx.rv_set( {'name':'rv_set'}, [ config_rv_X, config_rv_noise ] )
## ------------------------------------------
## Model output
## ------------------------------------------
model = [ "sin(rbrv(rv_set::X))+1+rbrv(rv_set::eps)" ]
## ------------------------------------------
## Reference truth
## ------------------------------------------
plot_lim = [-5.,20]
X_vals = fesslix.tools.discretize_x(plot_lim[0], plot_lim[1], x_disc_N=200)
Y_truth = np.sin(X_vals)+1.
## ------------------------------------------
## Data box for storing samples
## ------------------------------------------
dbox = flx.dataBox(2,1)
dbox.write2mem( {
'N_reserve': N_smpls_max,
'cols': [0,1]
} )
## ------------------------------------------
## Generate random realizations of the model
## ------------------------------------------
sampler = flx.sampler(['rv_set'])
sampler.perform_MCS(N_smpls_max, model, dbox)
smpls_X = dbox.extract_col_from_mem(1)
smpls_Y = dbox.extract_col_from_mem(0)
print(np.std(smpls_X), np.std(smpls_Y))
2.8718147 0.6996426
fig, ax = plt.subplots(figsize=(10, 4))
plt.scatter(smpls_X,smpls_Y,alpha=0.3)
plt.grid()
plt.show()
Define Gaussian process#
config_gp = { 'name':'gp_1',
'Ndim':1,
#'mean_type':'zero',
'mean_type':'universal',
'mean_polyo':0,
'kernel_lst':['gauss'],
'useLSE':True
}
gp_1 = fesslix.gpr.gp(config_gp)
config_gp['name'] = 'gp_2'
config_gp['useLSE'] = False
gp_2 = fesslix.gpr.gp(config_gp)
def gp_summarize(gp):
gp_info = gp.info()
prior_intercept = None
if 'para_vec' in gp_info['mean']:
prior_intercept = gp_info['mean']['para_vec'][0]*gp_info['mean']['normalizef']
sd_noise = gp_info['noise']
sd_kernel = gp_info['kernel']['kernel_sd']
corr_l = gp_info['kernel']['para_vec'][1]*gp_info['kernel']['n_vec'][1]
logl = gp_info['logl_obsv']
print(f"{gp_info['name']}: {logl:10.2e} {sd_noise:8.4f} {sd_kernel:6.2f} {prior_intercept:6.2f} {corr_l:6.2f}")
Condition the Gaussian process on data#
N = 100
gp_1.condition_on(smpls_X[:N].reshape(-1, 1),smpls_Y[:N],opt_noise=False)
gp_2.condition_on(smpls_X[:N].reshape(-1, 1),smpls_Y[:N],opt_noise=False)
gp_summarize(gp_1)
gp_summarize(gp_2)
gp_1: -5.53e+01 0.1513 1513.47 65.70 2.90
gp_2: -2.44e+08 0.0001 0.69 1.18 2.90
if True:
gp_1.condition_on(smpls_X[:N].reshape(-1, 1),smpls_Y[:N],opt_noise=True)
gp_2.condition_on(smpls_X[:N].reshape(-1, 1),smpls_Y[:N],opt_noise=True)
else:
noise_val = 0.0
for i in range(20):
noise_val += 0.01
gp_1.noise_white(noise_val)
gp_2.noise_white(noise_val)
print( f"{noise_val:4.2f} » {gp_1.condition_on(smpls_X[:N].reshape(-1, 1),smpls_Y[:N],opt_noise=False):10.2f} » {gp_2.condition_on(smpls_X[:N].reshape(-1, 1),smpls_Y[:N],opt_noise=False):10.2f}" )
gp_summarize(gp_1)
gp_summarize(gp_2)
gp_1: 1.61e+01 0.1697 1.19 0.75 2.90
gp_2: 1.22e+01 0.1701 0.69 1.18 2.90
if False:
gp_1_info = gp_1.info()
print( "noise log:", gp_1_info['noise_log'] )
def plot_GP(gp, color):
Y_vals = np.empty(len(X_vals))
sigma_vals = np.empty(len(X_vals))
for i in range(len(X_vals)):
Y_vals[i], sigma_vals[i] = gp.predict([X_vals[i]],type="mean_sd",predict_noise=False)
plt.plot(X_vals,Y_vals, color=color)
plt.plot(X_vals,Y_vals+2*sigma_vals, color=color, ls=':')
plt.plot(X_vals,Y_vals-2*sigma_vals, color=color, ls=':')
for i in range(len(X_vals)):
Y_vals[i], sigma_vals[i] = gp.predict([X_vals[i]],type="mean_sd",predict_noise=True)
plt.plot(X_vals,Y_vals+2*sigma_vals, color=color, ls='--')
plt.plot(X_vals,Y_vals-2*sigma_vals, color=color, ls='--')
fig, ax = plt.subplots(figsize=(10, 4))
plt.plot(X_vals,Y_truth,color='black')
plot_GP(gp_1, color='blue')
plot_GP(gp_2, color='red')
ax.set_ylim([-1., 3.])
ax.set_xlim(plot_lim)
plt.grid()
plt.show()
optimize process parameters#
gp_1.optimize()
gp_1_info = gp_1.info()
print( gp_1_info['opt_log'] )
flxGPProj::likeli_f_89 16.0584 ( 0 ) no
initial point estimate: 16.05836 at ( 0 ) dim=1
flxGPProj::likeli_f_89 16.0584 ( 0 ) no
flxGPProj::likeli_f_89 -84.178 ( 1 ) yes
flxGPProj::likeli_f_89 8.91244 ( -1 ) yes
flxGPProj::likeli_f_89 12.8561 ( -0.5 ) yes
flxGPProj::likeli_f_89 -11.7054 ( 0.5 ) yes
flxGPProj::likeli_f_89 15.075 ( -0.25 ) yes
flxGPProj::likeli_f_89 11.8028 ( 0.25 ) yes
flxGPProj::likeli_f_89 15.8445 ( -0.125 ) yes
flxGPProj::likeli_f_89 15.2171 ( 0.125 ) yes
flxGPProj::likeli_f_89 16.0419 ( -0.0625 ) yes
flxGPProj::likeli_f_89 15.8268 ( 0.0625 ) yes
flxGPProj::likeli_f_89 16.076 ( -0.03125 ) yes
flxGPProj::likeli_f_89 16.0419 ( -0.0625 ) yes
flxGPProj::likeli_f_89 16.0742 ( -0.015625 ) yes
flxGPProj::likeli_f_89 16.0649 ( -0.046875 ) yes
flxGPProj::likeli_f_89 16.0767 ( -0.0234375 ) yes
MLE: 16.0767 at ( -0.0234375 ) :: niter = 16, dim = 1
flxGPProj::likeli_f_89 16.0767 ( -0.0234375 ) no
for i in range(10):
gp_1.optimize()
gp_summarize(gp_1)
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_2.optimize()
#gp_2_info = gp_2.info()
#print( gp_2_info['opt_log'] )
16.136169807651008
gp_summarize(gp_1)
gp_summarize(gp_2)
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
fig, ax = plt.subplots(figsize=(10, 4))
plt.plot(X_vals,Y_truth,color='black')
plot_GP(gp_1, color='blue')
plot_GP(gp_2, color='red')
ax.set_ylim([-1., 3.])
ax.set_xlim(plot_lim)
plt.grid()
plt.show()
optimize noise based on optimized parameters#
def fun_opt_noise(gp):
for i in range(10):
gp_summarize(gp)
print( i, gp.optimize() )
gp_summarize(gp)
gp_1.condition_on(smpls_X[:N].reshape(-1, 1),smpls_Y[:N], init_pvec=False)
gp_summarize(gp)
fun_opt_noise(gp_1)
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
0 16.07674584375785
gp_1: 1.61e+01 0.1692 1.19 0.77 2.83
gp_1: 1.61e+01 0.1697 1.14 0.77 2.83
gp_1: 1.61e+01 0.1697 1.14 0.77 2.83
1 16.10092661673511
gp_1: 1.61e+01 0.1692 1.14 0.79 2.76
gp_1: 1.61e+01 0.1697 1.10 0.79 2.76
gp_1: 1.61e+01 0.1697 1.10 0.79 2.76
2 16.118652736558545
gp_1: 1.61e+01 0.1697 1.10 0.80 2.73
gp_1: 1.61e+01 0.1697 1.08 0.80 2.73
gp_1: 1.61e+01 0.1697 1.08 0.80 2.73
3 16.12554965295074
gp_1: 1.61e+01 0.1694 1.07 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
4 16.129761310285943
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
5 16.1297613102863
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
6 16.1297613102863
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
7 16.1297613102863
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
8 16.1297613102863
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
9 16.1297613102863
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
fun_opt_noise(gp_2)
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
0 16.136169807651008
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
1 16.136169807651008
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
2 16.136169807651008
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
3 16.136169807651008
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
4 16.136169807651008
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
5 16.136169807651008
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
6 16.136169807651008
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
7 16.136169807651008
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
8 16.136169807651008
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
9 16.136169807651008
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
gp_summarize(gp_1)
gp_summarize(gp_2)
gp_1: 1.61e+01 0.1697 1.05 0.81 2.69
gp_2: 1.61e+01 0.1701 0.99 0.84 2.60
fig, ax = plt.subplots(figsize=(10, 4))
plt.plot(X_vals,Y_truth,color='black')
plot_GP(gp_1, color='blue')
plot_GP(gp_2, color='red')
ax.set_ylim([-1., 3.])
ax.set_xlim(plot_lim)
plt.grid()
plt.show()
automatic optimization of noise#
gp_1.noise_white(0.1)
gp_2.noise_white(0.1)
gp_1.condition_on(smpls_X[:N].reshape(-1, 1),smpls_Y[:N],opt_noise=False)
gp_2.condition_on(smpls_X[:N].reshape(-1, 1),smpls_Y[:N],opt_noise=False)
gp_summarize(gp_1)
gp_summarize(gp_2)
gp_1: -5.53e+01 0.1513 1513.47 65.70 2.90
gp_2: -2.44e+08 0.0001 0.69 1.18 2.90
gp_1.optimize(opt_noise=True)
16.13658536094384
gp_2.optimize(opt_noise=True)
-244017268.5156445
gp_summarize(gp_1)
gp_summarize(gp_2)
gp_1: 1.61e+01 255.5594 1513.47 0.84 2.61
gp_2: -2.44e+08 0.0001 0.69 1.18 2.90
print( gp_1.info()['opt_log'] )
flxGPProj::likeli_f_89 -55.2689 ( 0, -1.88818 ) yes
initial point estimate: -55.2689 at ( 0, -1.88818 ) dim=2
flxGPProj::likeli_f_89 -55.2689 ( 0, -1.88818 ) no
flxGPProj::likeli_f_89 -17.69 ( 1, -1.88818 ) yes
flxGPProj::likeli_f_89 -33.5659 ( 1, -3.77636 ) yes
flxGPProj::likeli_f_89 -33.0378 ( 2, -3.77636 ) yes
flxGPProj::likeli_f_89 -38.1138 ( 2, -1.88818 ) yes
flxGPProj::likeli_f_89 -23.0848 ( 1.25, -3.30432 ) yes
flxGPProj::likeli_f_89 -35.9559 ( 0.25, -1.41614 ) yes
flxGPProj::likeli_f_89 -17.5727 ( 1.5625, -3.18631 ) yes
flxGPProj::likeli_f_89 -11.8933 ( 1.3125, -1.77017 ) yes
flxGPProj::likeli_f_89 -11.4487 ( 1.34375, -1.0031 ) yes
flxGPProj::likeli_f_89 -32.1602 ( 1.90625, -2.30122 ) yes
flxGPProj::likeli_f_89 -14.0358 ( 1.226562, -1.99144 ) yes
flxGPProj::likeli_f_89 -3.55244 ( 1.007812, 0.191768 ) yes
flxGPProj::likeli_f_89 3.96256 ( 0.730469, 1.880806 ) yes
flxGPProj::likeli_f_89 -2.34901 ( 0.847656, 2.869151 ) yes
flxGPProj::likeli_f_89 6.05566 ( 0.234375, 5.753053 ) yes
flxGPProj::likeli_f_89 -92.1171 ( -0.320312, 9.131128 ) yes
flxGPProj::likeli_f_89 14.9965 ( 0.117188, 4.764708 ) yes
flxGPProj::likeli_f_89 15.9509 ( -0.248047, 5.712487 ) yes
flxGPProj::likeli_f_89 -97.4021 ( -0.744141, 9.584734 ) yes
flxGPProj::likeli_f_89 12.1498 ( 0.361816, 3.806788 ) yes
flxGPProj::likeli_f_89 5.28114 ( -0.120605, 3.766221 ) yes
flxGPProj::likeli_f_89 15.3793 ( 0.14563, 5.256345 ) yes
flxGPProj::likeli_f_89 -5.33916 ( -0.464233, 7.162044 ) yes
flxGPProj::likeli_f_89 14.7396 ( 0.155304, 4.645602 ) yes
flxGPProj::likeli_f_89 12.6393 ( -0.257721, 6.32323 ) yes
flxGPProj::likeli_f_89 15.6893 ( 0.0520477, 5.065009 ) yes
flxGPProj::likeli_f_89 14.9711 ( -0.341629, 5.52115 ) yes
flxGPProj::likeli_f_89 16.0188 ( 0.0238152, 5.322547 ) yes
flxGPProj::likeli_f_89 15.5577 ( -0.276279, 5.970024 ) yes
flxGPProj::likeli_f_89 15.9739 ( -0.0300341, 5.291263 ) yes
flxGPProj::likeli_f_89 14.9991 ( 0.241828, 4.901323 ) yes
flxGPProj::likeli_f_89 16.097 ( -0.125578, 5.509696 ) yes
flxGPProj::likeli_f_89 16.1144 ( -0.0717289, 5.540979 ) yes
flxGPProj::likeli_f_89 15.9911 ( -0.0925764, 5.665838 ) yes
flxGPProj::likeli_f_89 16.0187 ( -0.221122, 5.728129 ) yes
flxGPProj::likeli_f_89 16.1025 ( -0.0374192, 5.423942 ) yes
flxGPProj::likeli_f_89 15.9357 ( 0.01643, 5.455226 ) yes
flxGPProj::likeli_f_89 16.1289 ( -0.0900761, 5.496078 ) yes
flxGPProj::likeli_f_89 16.125 ( -0.124386, 5.613116 ) yes
flxGPProj::likeli_f_89 16.1145 ( -0.142733, 5.568214 ) yes
flxGPProj::likeli_f_89 16.131 ( -0.124982, 5.561406 ) yes
flxGPProj::likeli_f_89 16.0872 ( -0.0906723, 5.444368 ) yes
flxGPProj::likeli_f_89 16.1357 ( -0.115957, 5.570929 ) yes
flxGPProj::likeli_f_89 16.1168 ( -0.150863, 5.636256 ) yes
flxGPProj::likeli_f_89 16.1344 ( -0.105273, 5.531123 ) yes
flxGPProj::likeli_f_89 16.1363 ( -0.0962484, 5.540646 ) yes
flxGPProj::likeli_f_89 16.131 ( -0.0818816, 5.530266 ) yes
flxGPProj::likeli_f_89 16.1311 ( -0.106933, 5.580452 ) yes
flxGPProj::likeli_f_89 16.1366 ( -0.105688, 5.543455 ) yes
MLE: 16.1366 at ( -0.105688, 5.54345 ) :: niter = 50, dim = 2
flxGPProj::likeli_f_89 16.1366 ( -0.105688, 5.543455 ) no
print( gp_2.info()['opt_log'] )
flxGPProj::likeli_f_89 -2.44017e+08 ( 1, 0, 0, -9.58851 ) no
initial point estimate: -2.440173e8 at ( 1, 0, 0, -9.58851 ) dim=4
flxGPProj::likeli_f_89 -2.44017e+08 ( 1, 0, 0, -9.58851 ) no
flxGPProj::likeli_f_89 -2.44017e+08 ( 2, 0, 0, -9.58851 ) yes
flxGPProj::likeli_f_89 -2.40469e+08 ( 2, 1, 0, -9.58851 ) yes
flxGPProj::likeli_f_89 -2.81474e+08 ( 2, 1, 1, -9.58851 ) yes
flxGPProj::likeli_f_89 -inf ( 2, 1, 0, -19.177 ) yes
flxGPProj::likeli_f_89 -2.44017e+08 ( 1, 0, 0, -9.58851 ) yes