Untitled

Live Script edama_05

This Live Script supports Chapter 5 of Environmental Data Analysis with MATLAB or PYTHON, Third Edition, by WIlliam Menke

The script naming system has the form edama_CC_SS, where CC is the chapter number and SS is the script number.

% edama_05_01: application of Bayes theorem to least squares

clear all; % clears all previousl-defined variabes

% set up vectors m1 and m2

L=40;

Dm = 1.0;

m1 = Dm*[0:L-1]';

m2 = Dm*[0:L-1]';

% make Normal p.d.f. of model parameters pp(m)

m1bar=10;

m2bar=20;

Cpm=[[6,0]',[0,10]'];

CpmI=inv(Cpm);

norm=1/(2*pi*sqrt(det(Cpm)));

Ppm=zeros(L,L);

for i = [1:L]

for j = [1:L]

dm = [m1(i)-m1bar, m2(j)-m2bar]';

pp(i,j)=norm*exp( -0.5 * dm' * CpmI * dm );

end

% search grid for maximum

[maxtemp, k] = max(pp);

[pmax, j] = max(maxtemp);

i=k(j);

m1bar = m1(i);

m2bar = m2(j);

fprintf('m1bar: %.2f m2bar: %.2f\n', m1bar, m2bar );

m1bar: 10.00   m2bar:  20.00

% make Normal p.d.f. pdgm(m)

% with a one-row data kernel that the

% difference of two model parameters is d1

G = [1, -1];

d1 = 0;

Cd=3;

CdI=1/Cd;

norm=1/sqrt(2*pi*Cd);

pdgm=zeros(L,L);

for i = [1:L]

for j = [1:L]

m = [m1(i), m2(j)]';

dd = d1-G*m;

pdgm(i,j)=norm*exp( -0.5 * dd' * CdI * dd );

end

% multiply the two to get Pmgd

pmgd = pp .* pdgm;

% search grid for maximum

[maxtemp, k] = max(pmgd);

[pmax, j] = max(maxtemp);

i=k(j);

m1est = m1(i);

m2est = m2(j);

fprintf('m1est: %.2f m2est: %.2f\n', m1est, m2est );

m1est: 13.00   m2est:  15.00

% use simple drawing function that encapsulates all the graphics

eda_draw(' ',pp,'caption pp(m)',pdgm,'caption p(d|m)',pmgd,'caption p(m|d)');

% edama_05_02: illustrate the product of twos Normal p.d.f.s

% set up vectors da and d2

L=40;

Dd = 1.0;

d1 = Dd*[0:L-1]';

d2 = Dd*[0:L-1]';

% first normal p.d.f.

d1bar=15;

d2bar=15;

C=[[10,8]',[8,10]'];

CI=inv(C);

norm=1/(2*pi*sqrt(det(C)));

Pa=zeros(L,L);

for i = [1:L]

for j = [1:L]

dd = [d1(i)-d1bar, d2(j)-d2bar]';

Pa(i,j)=norm*exp( -0.5 * dd' * CI * dd );

end

% second normal p.d.f.

d1bar=15;

d2bar=25;

C=[[10,-8]',[-8,10]'];

CI=inv(C);

norm=1/(2*pi*sqrt(det(C)));

Pb=zeros(L,L);

for i = [1:L]

for j = [1:L]

dd = [d1(i)-d1bar, d2(j)-d2bar]';

Pb(i,j)=norm*exp( -0.5 * dd' * CI * dd );

end

% product of two Normal p.d.f.'s, normalized to unit area

Pc=Pa.*Pb;

norm = (Dd^2)*sum(sum(Pc));

Pc = Pc/norm;

% plot results

eda_draw(Pa,'caption pA',Pb,'caption pB',Pc,'caption pC');

% edama_05_03: straight-line fit with unequal variance

% illustrate least squares fitting of a line

% when the data have unequal variance

% this version uses G and d

% standard set-up of x-axis

N=51;

N2=(N-1)/2;

Dx = 2;

xmin = 0.0;

xmax = xmin+Dx*(N-1);

x = xmin+Dx*[0:N-1]';

% prepare some test data hat consists of

% two sections, that differ slightly in

% interrcept and slope

a = 1.0;

b = 2.0;

dobs = a + b*x;

Da = 1;

Db = 0.3;

dobs(N2:N) = dobs(N2:N) + Da + Db*x(N2:N);

% Part A. high variance data at left

% prepare the variance, high for thr first section

% low for the second section

sd = zeros(N,1);

sd(1:N2-1) = 100;

sd(N2:N) = 5;

sd2i = sd.^(-2);

Cdi = diag(sd2i);

% set up the data kernel

G=zeros(N,2);

G(:,1)=ones(N,1);

G(:,2)=x;

% weighted least squares solution

FTF = G'*Cdi*G;

FTf = G'*Cdi*dobs;

mest = FTF\FTf;

% predict data

dpre = mest(1) + mest(2) .* x;

figure(1);

clf;

set(gca,'LineWidth',2);

hold on;

axis( [0, 100, -500, 500] );

plot( x, dobs, 'ko', 'LineWidth', 2 );

% plot 2-sigma error bars

for p = [1:N]

plot( [x(p), x(p)], [dobs(p)-2*sd(p), dobs(p)+2*sd(p)], 'k-', 'LineWidth', 2 );

hold on;

end

plot(x,dpre,'k-', 'LineWidth', 2);

% Part B. high variance data at right

% now swap the variances

sd = zeros(N,1);

sd(1:N2-1) = 5;

sd(N2:N) = 100;

sd2i = sd.^(-2);

Cdi = diag(sd2i);

% set up the data kernel

G(:,1)=ones(N,1);

G(:,2)=x;

% weighted least squares solution

FTF = G'*Cdi*G;

FTf = G'*Cdi*dobs;

mest = FTF\FTf;

% predict data

dpre = mest(1) + mest(2) .* x;

figure(2);

clf;

set(gca,'LineWidth',2);

hold on;

axis( [0, 100, -500, 500] );

plot( x, dobs, 'ko', 'LineWidth', 2 );

% plot 2-sigma error bars

for p = [1:N]

plot( [x(p), x(p)], [dobs(p)-2*sd(p), dobs(p)+2*sd(p)], 'k-', 'LineWidth', 2 );

hold on;

end

plot(x,dpre,'k-', 'LineWidth', 2);

hold on;

% edama_05_04: Straight-line fit with diagonal variance

% illustrate least squares fitting of a line

% when the data have unequal variance

% this version uses F and f

% standard set-up of x-axis

N=51;

N2=(N-1)/2;

Dx = 2;

xmin = 0.0;

xmax = xmin+Dx*(N-1);

x = xmin+Dx*[0:N-1]';

% prepare some test data hat consists of

% two sections, that differ slightly in

% interrcept and slope

M=2;

a = 1.0;

b = 2.0;

dobs = a + b*x;

Da = 1;

Db = 0.3;

dobs(N2:N) = dobs(N2:N) + Da + Db*x(N2:N);

% Part A, high variance on left

% prepare the variance, high for thr first section

% low for the second section

sd = zeros(N,1);

sd(1:N2-1) = 100;

sd(N2:N) = 5;

% set up the data kernel

G=zeros(N,2);

G(:,1)=ones(N,1);

G(:,2)=x;

% set up F and f

F=zeros(N,M);

for i=[1:N]

F(i,:)=G(i,:)./sd(i);

end

f=dobs./sd;

% generalized least squares

FTF = F'*F;

FTf = F'*f;

mest = FTF\FTf;

% predict data

dpre = mest(1) + mest(2) .* x;

figure(1);

clf;

set(gca,'LineWidth',2);

hold on;

axis( [0, 100, -500, 500] );

plot( x, dobs, 'ko', 'LineWidth', 2 );

% plot 2-sigma error bars

for p = [1:N]

plot( [x(p), x(p)], [dobs(p)-2*sd(p), dobs(p)+2*sd(p)], 'k-', 'LineWidth', 2 );

hold on;

end

plot(x,dpre,'k-','LineWidth',2);

% Part B, high variance on right

% now swap the variances

sd = zeros(N,1);

sd(1:N2-1) = 5;

sd(N2:N) = 100;

sd2i = sd.^(-2);

Cdi = diag(sd2i);

% set up the data kernel

G=zeros(N,2);

G(:,1)=ones(N,1);

G(:,2)=x;

% set up F and f

F=zeros(N,2);

for i=[1:N]

F(i,:)=G(i,:)./sd(i);

end

f=dobs./sd;

% generalized least squares

FTF = F'*F;

FTf = F'*f;

mest = FTF\FTf;

% predict data

dpre = mest(1) + mest(2) .* x;

figure(2);

clf;

set(gca,'LineWidth',2);

hold on;

axis( [0, 100, -500, 500] );

plot( x, dobs, 'ko', 'LineWidth', 2 );

% plot 2-sigma error bars

for p = [1:N]

plot( [x(p), x(p)], [dobs(p)-2*sd(p), dobs(p)+2*sd(p)], 'k-', 'LineWidth', 2 );

hold on;

end

plot(x,dpre,'k-','LineWidth',2);

hold on;

% edama_05_05: using smoothness to fill in data gaps

% simple example of smoothness applied

% to the problem of filliusing smoothness to fill in data gapsng in data gaps

% standard set-up of x-axis

M=101;

Dx=1.0;

Dx2 = Dx^2;

xmin=0.0;

xmax = xmin+Dx*(M-1);

x = xmin+Dx*[0:M-1]';

% the data are the vales the function measured on

% only a subset of points along the x axis

N=40;

% randomly choose N of the M points to have data

rowindex=sort(unidrnd(M,N,1));

% determine their x positions

xd=xmin+Dx*rowindex;

% evaluate a sinusoidal function at these points to

% simulate the data

freq=0.02;

dobs = sin(2*pi*freq*xd);

% plot the data

figure(1);

clf;

set(gca,'LineWidth',2);

hold on;

axis([xmin xmax -2 2]);

plot(xd,dobs,'ko','LineWidth',2);

xlabel('x');

ylabel('d');

% prior variances

% d = sin(2*pi*freq*x) has maximum value of 1

% dddx = 2*pi*freq*cos(2*pi*freq*x) has maximum value of 2*pi*freq

% d2ddx2 = - (2*pi*freq)^2 sin(2*pi*freq*x) has maximum value of (2*pi*freq)^2

% 1% error of signal level of 1 i 0.1

sigmad = 0.01; % observational error is 1% of data maximum

sigmadi = 1.0/sigmad;

sigmaD1 = 2*pi*freq;

sigmaD1i = 1.0/sigmaD1;

sigmaD2 = (2*pi*freq)^2;

sigmaD2i = 1.0/sigmaD2;

% build F and f

K = M;

F = zeros(N+K,M);

fobs=zeros(N+K,1);

nextrow=1;

% the model parameter equals the observed data

for p = [1:N]

F(nextrow,rowindex(p)) = sigmadi;

fobs(nextrow)=sigmadi*dobs(p);

nextrow=nextrow+1;

end

% smoothness (small second derivative) in the interior

for p = [1:M-2]

F(nextrow,p) = sigmaD2i/Dx2;

F(nextrow,p+1) = -2*sigmaD2i/Dx2;

F(nextrow,p+2) = sigmaD2i/Dx2;

fobs(nextrow)=0.0;

nextrow=nextrow+1;

end

% flatness (small first derivative) at the left end

F(nextrow,1)=-sigmaD1i/Dx;

F(nextrow,2)=sigmaD1i/Dx;

fobs(nextrow)=0.0;

nextrow=nextrow+1;

% flatness (small first derivative) at the right end

F(nextrow,M-1)=-sigmaD1i/Dx;

F(nextrow,M)=sigmaD1i/Dx;

fobs(nextrow)=0.0;

nextrow=nextrow+1;

fprintf('build %d rows out of %d\n', nextrow-1, N+K);

build 141 rows out of 141

% solve generalized least squares equation

FTF = F'*F;

FTf = F'*fobs;

mest= FTF\FTf;

plot(x,mest,'k-','LineWidth',2);

e = dobs-mest(rowindex(:));

E=e'*e;

% edama_05_06: sparse matrix from list of non-zero elements
% build lists
L=7;
M=6;
Frow = [1, 3, 4, 6]';
Fcol = [1, 3, 1, 5]';
Fval = [1.0, 2.0, 3.0, 4.0]';
KK=length(Fval);
% build sparse matrix
edaFsparse =sparse(Frow(1:KK),Fcol(1:KK),Fval(1:KK),L,M);
% clear no longer needed vectors
clear Frow
clear Fcol
clear Fval
F = full(edaFsparse);
F
F = 
     1     0     0     0     0     0
     0     0     0     0     0     0
     0     0     2     0     0     0
     3     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     0     4     0
     0     0     0     0     0     0

% eda05_07: sparse matrix implementation of Fm=f

% simple example of smoothness applied

% to the problem of filling in data gaps

% sparse matrix implementation

clear edaFsparse;

global edaFsparse;

% the model parameters are the values of a sinusoidal

% function evenly spaced along the x axis

M=101;

Dx=1.0;

Dx2 = Dx^2;

xmin=0.0;

xmax = xmin+Dx*(M-1);

x = xmin+Dx*[0:M-1]';

% the data are the vales the function measured on

% only a subset of points along the x axis

N=40;

% randomly choose N of the M points to have data

rowindex=sort(unidrnd(M,N,1));

% determine their x positions

xd=xmin+Dx*rowindex;

% evaluate a sinusoidal function at these points to

% simulate the data

freq=0.02;

dobs = sin(2*pi*freq*xd);

% plot the data

figure(1);

clf;

set(gca,'LineWidth',2);

hold on;

axis([xmin xmax -2 2]);

plot(xd,dobs,'ko','LineWidth',2);

xlabel('x');

ylabel('d');

% prior variances

% d = sin(2*pi*freq*x) has maximum value of 1

% dddx = 2*pi*freq*cos(2*pi*freq*x) has maximum value of 2*pi*freq

% d2ddx2 = - (2*pi*freq)^2 sin(2*pi*freq*x) has maximum value of (2*pi*freq)^2

% 1% error of signal level of 1 i 0.1

sigmad = 0.01; % observational error is 1% of data maximum

sigmadi = 1.0/sigmad;

sigmaD1 = 2*pi*freq;

sigmaD1i = 1.0/sigmaD1;

sigmaD2 = (2*pi*freq)^2;

sigmaD2i = 1.0/sigmaD2;

% arrays to hold non-zero elements of F

Frow = zeros((N+3*M)+100,1);

Fcol = zeros((N+3*M)+100,1);

Fval = zeros((N+3*M)+100,1);

% the first N rows of F are that the model parameter

% equals the observed data

fobs=zeros((N+M),1);

KK=1; % counts non-zero elements

nextrow=1;

% the model parameter equals the observed data

for p = [1:N]

Frow(KK) = nextrow;

Fcol(KK) = rowindex(p);

Fval(KK) = sigmadi;

KK=KK+1;

fobs(nextrow)=sigmadi*dobs(p);

nextrow=nextrow+1;

end

% smoothness (small second derivative) in the interior

for p = [1:M-2]

Frow(KK) = nextrow;

Fcol(KK) = p;

Fval(KK) = sigmaD2i/Dx2;

KK=KK+1;

Frow(KK) = nextrow;

Fcol(KK) = p+1;

Fval(KK) = -2*sigmaD2i/Dx2;

KK=KK+1;

Frow(KK) = nextrow;

Fcol(KK) = p+2;

Fval(KK) = sigmaD2i/Dx2;

KK=KK+1;

fobs(nextrow)=0.0;

nextrow=nextrow+1;

end

% flatness (small first derivative) at the left end

Frow(KK) = nextrow;

Fcol(KK) = 1;

Fval(KK) = -sigmaD1i*Dx;

KK=KK+1;

Frow(KK) = nextrow;

Fcol(KK) = 2;

Fval(KK) = sigmaD1i*Dx;

KK=KK+1;

fobs(nextrow)=0.0;

nextrow=nextrow+1;

% flatness (small first derivative) at the right end

Frow(KK) = nextrow;

Fcol(KK) = M-1;

Fval(KK) = -sigmaD1i/Dx;

KK=KK+1;

Frow(KK) = nextrow;

Fcol(KK) = M;

Fval(KK) = sigmaD1i/Dx;

KK=KK+1;

fobs(nextrow)=0.0;

nextrow=nextrow+1;

KK=KK-1;

fprintf('%d non-zero-elements, expected %d\n', KK, N+3*(M-2)+2*2 );

341 non-zero-elements, expected 341

fprintf('build %d rows out of %d\n', nextrow-1, N+M);

build 141 rows out of 141

% build sparse matrix

edaFsparse =sparse(Frow(1:KK),Fcol(1:KK),Fval(1:KK),N+M,M);

% clear no longer needed vectors

clear Frow

clear Fcol

clear Fval

% solve using conjugate gradient algorithm

mest=bicg(@FTFmul,edaFsparse'*fobs,1e-12,3*(M+N));

bicg converged at iteration 133 to a solution with relative residual 6.1e-13.

plot(x,mest,'k-','LineWidth',2);

e = dobs-mest(rowindex(:));

E=e'*e;

% edama_05_08: unfolding a matrix into a vector
I = 4;   % rows in matrix
J = 3;   % columns in matrix
M = I*J; % elements in matrix
A = [ [1, 2, 3, 4]', [5, 6, 7, 8]', [9, 10, 11 12]'];
A
% Part 1: with a formula
fprintf('PART 1: with a formula\n');
% unfold a matrix B into a vector m
m=zeros(M,1);
for i = [1:I]
for j = [1:J]
    k = (i-1)*J+j;  % row-wise unfold
    m(k) = A(i,j);
end
end
m
% fold a vector m into a matrix B
for k = [1:M]
    i = floor((k-1)/J)+1; % row-wise fold
    j = k-(i-1)*J
    B(i,j)=m(k);
end
B
% Part 2: with look-up tables
fprintf('\nPART 2: with look-up tables\n');
% I view the lookup table approach as being less
% error prone than a method that uses formulas to
% convert (i,j) into k (and vice versa).  The idea
% is to maintain look-up tables that do the conversion.
% Look-up tables are faster than formulas, but they
% do require some memory.
% build lookup tables
i_of_k=zeros(M,1);   % row i of B associated with m(k)
j_of_k=zeros(M,1);   % col j of B associated with m(k)
k_of_ij=zeros(I,J);  % row k of m associated with A(i,j)
k=1; % initialize counter
for i = [1:I]
for j = [1:J]
    i_of_k(k)=i;
    j_of_k(k)=j;
    k_of_ij(i,j)=k;
    k = k+1;  % increment counter
end
end
M=k-1;
% unfold a matrix A into a vector m
m=zeros(M,1);
for i = [1:I]
for j = [1:J]
    k = k_of_ij(i,j);
    m(k) = A(i,j);
end
end
m
% fold a vector m into a matrix B
B = zeros(I,J);
for k = [1:M]
    i = i_of_k(k);
    j = j_of_k(k);
    B(i,j)=m(k);
end
B

% edama_05_09: fill in data gaps in 2D using diff eqn

% use generealized least square to fill in

% data gaps for data that satisfy Laplace's

% equation in the in the interior of a grid,

% and the normal derivative is zero on the

% edges of the grid

% need this to communicate with edaFTFmul

clear edaFsparse;

global edaFsparse;

% load the data. These data are only 'simulated'

% and were created by the script eda05_08

D = load('../Data/pressure.txt');

% break-out individual columns

xobs=D(:,1);

yobs=D(:,2);

dobs=D(:,3);

N=length(dobs);

% prior variance of data

sigmad = 1.0;

rsigmad = 1/sigmad;

rsigmad2 = rsigmad^2;

% standard set-up of x-axis

I = 41;

Dx = 1;

xmin = 0;

xmax = 40;

x = xmin + (xmax-xmin)*[0:I-1]'/(I-1);

rDx = 1/Dx;

rDx2 = rDx^2;

% standard set-up of y-axis

J = 41;

Dy = 1;

ymin = 0;

ymax = 40;

y = ymin + (ymax-ymin)*[0:J-1]'/(J-1);

rDy = 1/Dy;

rDy2 = rDy^2;

% associate observations with points on the grid

% by scaling the range (xmin, xmax) to (1, N)

rowindex = 1+floor((I-1)*(xobs-xmin)/(xmax-xmin)+0.5);

colindex = 1+floor((J-1)*(yobs-ymin)/(ymax-ymin)+0.5);

% don't let any data fall outside the grid. In this

% case, we just put it on the edge of the grid. An

% alternative would be to throw it away.

for i = [1:N]

if( rowindex(i)<1 )

rowindex(i)=1;

elseif( rowindex(i)>I )

rowindex(i)=I;

end

if( colindex(i)<1 )

colindex(i)=1;

elseif( colindex(i)>J )

colindex(i)=J;

end

% copy the data to the grid A for later use

Aobs = zeros(I,J);

for i = [1:N]

Aobs(rowindex(i),colindex(i)) = dobs(i);

end

% one model parameter per grid point so I*J model parameters

M = I*J;

% Fm=f is LxM with L = (N rows of data) plus (M rows of prior information)

L = N+M;

fobs = zeros(L,1); % right hand side of equation Fm=f

% In this implementation, we first build a list of

% the non-zero elements of F consisting of three arrays,

% Frow, Fcol and Fval, with:

% Frow(k): the row index of a non-zero element of F(i,j)

% Fcol(k): the corresponding column index of F(i,j)

% Fval(k): the corresponding value of F(i,j)

% rows of G either have one non-zero elements (the G part)

% of 5 non-zero elements (the H part), so that

% F can have no more than 5L non-zero elements

Frow = zeros(5*L+100,1);

Fcol = zeros(5*L+100,1);

Fval = zeros(5*L+100,1);

% reminders on how to move between a grid of model parameters

% and a vector model parameter & vice-versa

% k = (i-1)*J+j;

% m(k) = A(i,j);

% i = floor((k-1)/J)+1;

% j = k-(i-1)*J;

% B(i,j)=m(k);

% we will also keep track of which row of F a particular datum is in

dindex = zeros(N,1);

% build F and f

% Part 1, Build the Gm=d part

% the variable nextrow is used to keep track of the current row of F and f

% the variable KK is used to keep track of the current row of Frow, Fcol and Fval

KK = 1; % initialize counter for non-zero elements of F

nextrow=1; % initialize counter for the rows of F

% first part of Fm=f: the data part Gm=d

for i = [1:N] % loop over data

% non-zero element of F

k = (rowindex(i)-1)*J+colindex(i);

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmad;

KK = KK+1; % Increment non-zero element counter

fobs(nextrow) = rsigmad*dobs(i);

dindex(i)=k;

nextrow = nextrow+1; % increment counter

end

% error of prior information

sigmah = 1000;

rsigmah = 1/sigmah;

% Part 2, Build the Hm=h part that consists of

% Laplaces equation in the interior of the grid

% second derivative is a (1,-2,1)/Dx^2 in x plus a (1,-2,1)/Dy^2 in y

for i = [2:I-1] % loop over x's of interior of grid

for j = [2:J-1] % loop over y's of interior of grid

% 1st non-zero element of F

k = (i-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*(-2*rDx2-2*rDy2);

KK = KK+1; % Increment non-zero element counter

% 2nd non-zero element of F

k = ((i-1)-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDx2;

KK = KK+1; % Increment non-zero element counter

% 3rd non-zero element of F

k = ((i+1)-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDx2;

KK = KK+1; % Increment non-zero element counter

% 4th non-zero element of F

k = (i-1)*J+(j-1);

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDy2;

KK = KK+1; % Increment non-zero element counter

% 5th non-zero element of F

k = (i-1)*J+(j+1);

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDy2;

KK = KK+1; % Increment non-zero element counter

fobs(nextrow) = 0.0;

nextrow = nextrow+1; % increment row counter

end

% third part: top edge haz zero slope

i = 1;

for j = [2:J-1]

% first non-zero element

k = (i-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = -rsigmah*rDx;

KK = KK+1; % Increment non-zero element counter

% second non-zero element

k = ((i+1)-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDx;

KK = KK+1; % Increment non-zero element counter

fobs(nextrow) = 0.0;

nextrow = nextrow+1; % Increment row counter

end

% fourth part: bottom edge

i = I;

for j = [2:J-1]

% first non-zero element

k = (i-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDx;

KK = KK+1; % Increment non-zero element counter

% second non-zero element

k = ((i-1)-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = -rsigmah*rDx;

KK = KK+1; % Increment non-zero element counter

fobs(nextrow) = 0.0;

nextrow = nextrow+1; % Increment row counter

end

% fifth part: left edge

j = 1;

for i = [2:I-1]

% first non-zero element

k = (i-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = -rsigmah*rDy;

KK = KK+1; % Increment non-zero element counter

% second non-zero element

k = (i-1)*J+(j+1);

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDy;

KK = KK+1; % Increment non-zero element counter

fobs(nextrow) = 0.0;

nextrow = nextrow+1; % Increment row counter

end

% sixth part: right edge

j = J;

for i = [2:I-1]

% first non-zero element

k = (i-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDy;

KK = KK+1; % Increment non-zero element counter

% second non-zero element

k = (i-1)*J+(j-1);

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = -rsigmah*rDy;

KK = KK+1; % Increment non-zero element counter

fobs(nextrow) = 0.0;

nextrow = nextrow+1;

end

% seventh part: four corners

% top left

rDz = 1/sqrt((Dx^2)+(Dy^2));

i=1;

j=1;

% first non-zero element

k = (i-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = -rsigmah*rDz;

KK = KK+1; % Increment non-zero element counter

% second non-zero element

k = ((i+1)-1)*J+(j+1);

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDz;

KK = KK+1; % Increment non-zero element counter

fobs(nextrow)=0.0;

nextrow=nextrow+1;

% top right

i=1;

j=J;

% first non-zero element

k = (i-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDz;

KK = KK+1; % Increment non-zero element counter

% second non-zero element

k = ((i+1)-1)*J+(j-1);

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = -rsigmah*rDz;

KK = KK+1; % Increment non-zero element counter

fobs(nextrow)=0.0;

nextrow=nextrow+1;

% bottom left

i=I;

j=1;

% first non-zero element

k = (i-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = -rsigmah*rDz;

KK = KK+1; % Increment non-zero element counter

% second non-zero element

k = ((i-1)-1)*J+(j+1);

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDz;

KK = KK+1; % Increment non-zero element counter

fobs(nextrow)=0.0;

nextrow=nextrow+1;

% bottom right

i=I;

j=1;

% first non-zero element

k = (i-1)*J+j;

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = rsigmah*rDz;

KK = KK+1; % Increment non-zero element counter

% second non-zero element

k = ((i-1)-1)*J+(j-1);

Frow(KK)=nextrow;

Fcol(KK)=k;

Fval(KK) = -rsigmah*rDz;

KK = KK+1; % Increment non-zero element counter

% dont increment element counter

fobs(nextrow)=0.0;

nextrow=nextrow+1;

KK = KK-1;

nextrow=nextrow-1;

% actual rows in F

LL = nextrow;

fprintf("%d non-zero elements\n", KK);

8093 non-zero elements

fprintf("%d rows, expected %d\n", LL, (N+M) );

1849 rows, expected 1849

% trim unused parts of arrays and build sparse matrix

Frow = Frow(1:KK);

Fcol = Fcol(1:KK);

Fval = Fval(1:KK);

edaFsparse=sparse(Frow,Fcol,Fval,L,M);

% clear no longer needed vectors

clear Frow

clear Fcol

clear Fval

% solve equation Fm=f in least squares sense

mest=bicg(@FTFmul,edaFsparse'*fobs,1e-12,3*(N+M));

bicg converged at iteration 581 to a solution with relative residual 1e-12.

% reorganize estimated model parameters back into grid

Apre = zeros(I,J);

for k = [1:M]

i = floor((k-1)/J)+1;

j = k-(i-1)*J;

Apre(i,j) = mest(k);

end

% plot results

eda_draw(' ',Aobs,'caption p-obs',Apre,'caption p-pre');

% compute weighted prediction error Ed

dpre=zeros(N,1);

for i = [1:N]

k=dindex(i);

dpre(i) = mest(k);

end

Ed = rsigmad2*(dobs-dpre)'*(dobs-dpre);

fprintf('weighted error %.4f\n', Ed);

weighted error 0.1419

% edama_05_10: make simulated pressure data
% Note that the file name has been changed from
% ../Data/pressure.txt to ../Data/mypressure.txt
% so as not to overwrite the original file that
% is used by edama_05_09.
% rows of grid are along x-axis
I = 41;
Dx = 1;
xmin = 0;
xmax = 40;
x = xmin + (xmax-xmin)*[0:I-1]'/(I-1);
% columns of grid are along y-axis
J = 41;
Dy = 1;
ymin = 0;
ymax = 40;
y = ymin + (ymax-ymin)*[0:J-1]'/(J-1);
% grid 10% populated with data
N = floor(I*J/10);         % 10 percent of data
rowindex=unidrnd(I,N,1);   % random rows
xobs = x(rowindex);        % x's corresponding to rows
colindex=unidrnd(J,N,1);   % random columns
yobs = y(colindex);        % y's corresponding to cols
% true data are an analytic solution to Laplace's equation
kappa = 0.1;
dtrue = 10*sin(kappa*xobs).*exp(-kappa*yobs);
% synthetic observed data are true data plus random noise
sigmad = 0.1;
dobs = dtrue + random('normal',0.0,sigmad,N,1);
% populate an array A with data for plotting purposed
A = zeros(I,J);
for i = [1:N]
    A(rowindex(i),colindex(i)) = dobs(i);
end
eda_draw(' ',A,'caption A');
% aggregate into array for output purposed
D=zeros(N,3);
D(:,1)=xobs;
D(:,2)=yobs;
D(:,3)=dobs;
dlmwrite('../Data/mypressure.txt',D,'\t');