% kriging1.m generates field values at various locations, then predicts values at other locations using the known covariance

clear
sd       = 1;                                   % std. deviation of wavenumbers
% --------------------------------- Define points at which to observe the field
numpoints= 20;                                  % number of observations
r        = 5;
x        = r*rand(numpoints,1);
y        = r*rand(numpoints,1);

for i=1:numpoints,                              % calculate Sigma_tt matrix
  for j=i:numpoints,                            % just do the upper triangle
    Stt(i,j) = exp(-(sd*sd/2)*((x(i)-x(j))^2 + (y(i)-y(j))^2));
    Stt(j,i) = Stt(i,j);                        % make it symmetric
  end
end

Stti = inv(Stt);                                % inverse at known points

% --------------------------------- Generate variables for the field itself
m        = 100;                                 % number of Fourier modes
A        = randn(1,m)*sqrt(2/m);                % amplitudes of Fourier modes
phi      = 2 * pi * rand(m,1);                  % random phase
K        = sd * randn(m,2);                     % normally dist'd wavenumbers

for i=1:numpoints,                              % step through observations
  Zo(i,1) = A * cos(K*[x(i); y(i)] + phi);      % observed field value at x,y
end

% --------------------------------- Define points at which to predict field
gridsize = [30 30];                             % number of grid points in x, y
gridrange= [0 r 0 r];                           % [xmin xmax ymin ymax]
[xg,yg]  = easygrid(gridrange, gridsize);       % matrices of x and y values

for i=1:gridsize(1),                            % step through grid points
  for j=1:gridsize(2),
    Za(i,j) = A * cos(K*[xg(i,j); yg(i,j)] + phi); % actual field value at x,y
    for k=1:numpoints,
      Sst(1,k) = exp(-(sd*sd/2)*((x(k)-xg(i,j))^2 + (y(k)-yg(i,j))^2));
    end
    Zp(i,j) = Sst * Stti * Zo;                  % predicted field value
    Er(i,j) = 1 - Sst*Stti*Sst';
  end
end

figure(1)
clf

subplot(3,2,1)
contour(xg,yg,Er,20)
hold on
plot(x,y,'*');
axis(gridrange);
title('Locations at which the field is observed');
xlabel('with contours of expected error');

subplot(3,2,2)
tag3(x,y,Zo)
view(-10,30)
title('Observed data points in 3d')

rotate3d on

subplot(3,2,3)
plot3(x,y,Zo,'o')
hold on
mesh(xg,yg,Za)
view(-10,30)
xlabel('with observed values');
title('Actual field values at grid points');

subplot(3,2,4)
plot3(x,y,Zo,'o')
hold on
mesh(xg,yg,Zp)
view(-10,30)
title('Predicted field values at grid points');
xlabel('with observed values');

subplot(3,2,5)
mesh(xg,yg,abs(Zp-Za))
view(-10,30)
title('Absolute difference between predicted and actual');

subplot(3,2,6)
mesh(xg,yg,sqrt(abs(Er)))
view(-10,30)
title('Standard deviation of prediction error');