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

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

 ``` 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 %% %% 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 . % 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'); ```