-
Notifications
You must be signed in to change notification settings - Fork 0
/
emulator.m
83 lines (73 loc) · 2.32 KB
/
emulator.m
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
function em = emulator(sim_dat_input,sim_dat_output,...
pred_pts,omega,rho,lambda,Sigma_xx,inv_Sig_xx,verbose)
% This function outputs a covariance matrix for the input control and
% calibration vectors.
% sim_dat_input should be rescaled to [0,1] and sim_dat_output should be
% standardized. pred_pts on the same scales. Both should already
% have a dummy variable corresponding to output index as their first
% column.
% The covariance matrix of the training points may optionally be supplied.
% if not available, set argument Sigma_xx to be 0.
% Package things up for emulation
x = sim_dat_input;
sim_res = sim_dat_output;
xp = pred_pts;
num_cntrl = length(omega);
num_calib = length(rho);
nugsize = 1e-4;
% This will store the predicted response
ypred = zeros(size(pred_pts,1),1);
% Get Cov matrices
WN_xx = eye(size(x,1)) * nugsize; % For numerical stabilization
WN_xpxp = eye(size(pred_pts,1)) * nugsize; % For numerical stabilization
if verbose==true
fprintf('Getting three covariance matrices... ')
end
% Get Cov matrix of training points only if it is not supplied
if Sigma_xx == 0
Sigma_xx = gp_cov_old(omega,...
x(:,1:size(x,2)-num_calib),...
x(:,1:size(x,2)-num_calib),...
rho,...
x(:,size(x,2)-num_calib+1:end),...
x(:,size(x,2)-num_calib+1:end),...
lambda,verbose) +...
WN_xx;
inv_Sig_xx = inv(Sigma_xx);
rc=rcond(inv_Sig_xx);
end
if verbose==true
fprintf('1,... ')
end
Sigma_xpx = gp_cov_old(omega,...
xp(:,1:size(x,2)-num_calib),...
x(:,1:size(x,2)-num_calib),...
rho,...
xp(:,size(x,2)-num_calib+1:end),...
x(:,size(x,2)-num_calib+1:end),...
lambda,verbose);
if verbose==true
fprintf('2,... ')
end
Sigma_xxp = Sigma_xpx';
Sigma_xpxp = gp_cov_old(omega,...
xp(:,1:size(x,2)-num_calib),...
xp(:,1:size(x,2)-num_calib),...
rho,...
xp(:,size(x,2)-num_calib+1:end),...
xp(:,size(x,2)-num_calib+1:end),...
lambda,verbose) +...
WN_xpxp;
if verbose==true
fprintf('3.\n\n')
end
% fprintf('rcond(inv(Sigma_xx+WN))=%.3g\n\n',rc);
% Get mean and covariance at prediction points
mu = Sigma_xpx * inv_Sig_xx * sim_res(:);
cov_gp = Sigma_xpxp - Sigma_xpx * inv_Sig_xx * Sigma_xxp;
% Pack up and leave
em.Sigma_xx = Sigma_xx;
em.mu = mu;
em.cov_gp = cov_gp;
em.xp = xp;
end