Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

[aefc8d]: inst / test_optiminterp_err.m Maximize Restore History

Download this file

test_optiminterp_err.m    106 lines (79 with data), 2.5 kB

  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
%% Copyright (C) 2006 Alexander Barth <barth.alexander@gmail.com>
%%
%% This program is free software; you can redistribute it and/or modify it under
%% the terms of the GNU General Public License as published by the Free Software
%% Foundation; either version 3 of the License, or (at your option) any later
%% version.
%%
%% This program is distributed in the hope that it will be useful, but WITHOUT
%% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
%% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
%% details.
%%
%% You should have received a copy of the GNU General Public License along with
%% this program; if not, see <http://www.gnu.org/licenses/>.
% Tests 2D optimal interpolation and compares result
% (analysis and error field) to global OI solution
printf('Testing 2D-optimal interpolation (analysis and error field): ');
% grid of background field
[xi,yi] = ndgrid(linspace(0,1,30));
fi_ref = sin( xi*6 ) .* cos( yi*6);
% grid of observations
[x,y] = ndgrid(linspace(0,1,6));
x = x(:);
y = y(:);
on = numel(x);
var = 0.01 * ones(on,1);
f = sin( x*6 ) .* cos( y*6);
len = 0.1;
m = min(30,on);
% covariance function
% gaussian
%bcovar2 = @(d2) exp(-d2/len^2) ;
% diva
bcovar2 = @(d2) max(sqrt(d2)/len,eps) .* besselk(1,max(sqrt(d2)/len,eps));
% P: covariance between grid points (xi,yi) and grid points (xi,yi)
P = zeros(numel(xi),numel(xi));
for j=1:numel(xi)
for i=1:numel(xi)
P(i,j) = (xi(i) - xi(j))^2 + (yi(i) - yi(j))^2;
end
end
P = bcovar2(P);
% HPH: covariance between observation points (x,y) and observation points (x,y)
HPH = zeros(numel(x),numel(x));
for j=1:numel(x)
for i=1:numel(x)
HPH(i,j) = (x(i) - x(j))^2 + (y(i) - y(j))^2;
end
end
HPH = bcovar2(HPH);
% PH: covariance between grid points (xi,yi) and observation points (x,y)
PH = zeros(numel(xi),numel(x));
for j=1:numel(x)
for i=1:numel(xi)
PH(i,j) = (xi(i) - x(j))^2 + (yi(i) - y(j))^2;
end
end
PH = bcovar2(PH);
R = diag(var);
% call optiminterp
[fi,vari] = optiminterp2(x,y,f,var,len,len,m,xi,yi);
% Kalman gain
K = PH * inv(HPH + R);
% analysis
fi2 = K * f;
% error field
vari2 = diag(P - K * PH');
% transform vectors into 2d-arrays
fi2 = reshape(fi2,size(fi));
vari2 = reshape(vari2,size(fi));
rms = sqrt(mean((fi2(:) - fi(:)).^2));
if rms > 1e-4
error('unexpected large RMS difference (analysis)');
end
rms = sqrt(mean((vari2(:) - vari(:)).^2));
if rms > 1e-4
error('unexpected large RMS difference (error field)');
end
disp('OK');