% gdama05
% W. Menke (2024), Geophysical Data Analysis and Inverse Theory
% with MATLAB(R) or Python, 5th Edition
% Code History:
% 2022-10-01, Imported from 4th Edition by W. Menke
% 2023-04-05, Edit and integration with text
% gdama05_01
% data Resolution Matrix example
clear all;
fprintf('gdama05_01\n');
% z-axis
N=101;
zmin=0;
zmax=10;
Dz = (zmax-zmin)/(N-1);
z = zmin + Dz*[0:N-1]';
% create synthetic data
% d = a + b*z + noise
a=2.0;
b=1.0;
sd=0.5;
dtrue = a+b*z;
dobs = dtrue+random('Normal',0,sd,N,1);
% least squares fit
M=2;
G=[ones(N,1), z]; % data kernel
GMG = (G'*G)\G'; % generalized inverse
mest = GMG * dobs;
dpre = G*mest;
% predicted data
dpre = G*mest;
e = dobs - dpre;
[emax, iemax] = max(abs(e));
% compute and plot data resolution matrix
Nres = G*GMG;
gda_draw(' ',dpre,'=',' ',Nres,' ',dobs);
fprintf('Caption: Graphical depiction of the equation dpre = N dobs\n');
% plot
figure();
clf;
% plot scale
pdmin=0;
pdmax=15;
% plot observed and predicted data
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [zmin, zmax, pdmin, pdmax ]' );
plot( z, dobs, 'ro', 'LineWidth', 2);
plot( z, dpre, 'b-', 'LineWidth', 2);
xlabel('z');
ylabel('d');
fprintf('Caption: Linear fit (blue) to observed data (red)\n');
% gdama05_02
% Model Resolution Matrix example
clear all;
% z-axis
M=101;
zmin=0;
zmax=10;
Dz=(zmax-zmin)/(M-1);
z=zmin+Dz*[0:M-1]';
% model, m(z), moztly zero but a few spikes
mtrue = zeros(M,1);
mtrue(5)=1;
mtrue(10)=1;
mtrue(20)=1;
mtrue(50)=1;
mtrue(90)=1;
% experiment: exponential smoothing of model
N=80;
cmin=0.00;
cmax=7.90;
Dc=(cmax-cmin)/(N-1);
c = cmin + Dc*[0:N-1]';
G = exp(-c*z'); % data kernel
% true data and synthetic observed data
sd=0.0;
dtrue = G*mtrue;
dobs = dtrue + random('Normal',0,sd,N,1);
% minimum length solution
epsilon=1e-12;
GMG = G'/(G*G'+epsilon*eye(N,N)); % generalized inverse
mest = GMG * dobs;
dpre = G * mest;
Rres = GMG*G; % model resolution matrix
Nres = G*GMG; % data resolution matrix
% plot data kernel
gda_draw(dobs,'caption dobs',' = ',G,'caption G',mtrue,'caption mtrue');
fprintf('Caption: Graphical depiction of the equation dobs = G mtrue\n');
% plot solution
gda_draw(mest,'caption mest',' = ',GMG,'caption GMG',dobs,'caption dobs');
fprintf('Caption: Graphical depiction of the equation mest = G^-g dobs\n');
% plot model resolution matrix
gda_draw(mest,'caption mest',' = ',Rres,'caption R',mtrue,'caption mtrue');
fprintf('Caption: Graphical depiction of the equation mest = R mtrue\n');
% plot data resolution matrix
gda_draw(dpre,'caption dpre',' = ',Nres,'caption N',dobs,'caption dobs');
fprintf('Caption: Graphical depiction of the equation dpre = N dobs\n');
% plot
figure();
clf;
hold on;
% plot scale
pmmin=-0.2;
pmmax=1;
% plot true and estimated model
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
axis( [zmin, zmax, pmmin, pmmax ]' );
plot( z, mtrue, 'r-', 'LineWidth', 3);
plot( z, mest, 'b-', 'LineWidth', 3);
xlabel('z');
ylabel('m');
fprintf('Caption: True model (red) and estimated model (blue)\n');
% gdama05_03
% least squares fit to two cases of synthetic data
% distinguished by the spacing of the z's
clear all;
fprintf('gdama05_03\n');
% CASE 1: date clumped in middle of interval
% z-axis
N=10;
zmin=0;
zmax=10;
z = sort(random('Uniform',zmin+4*(zmax-zmin)/10,zmin+6*(zmax-zmin)/10,N,1));
% create synthetic observed data
% d = a + b*z + noise
a=5.0;
b=0.5;
sd=1;
dobs = a+b*z+random('Normal',0,sd,N,1);
% least squares fit
M=2;
G=[ones(N,1), z];
mest = (G'*G)\(G'*dobs);
Cm = (sd^2)*inv(G'*G); % covariance of model parameters
sm = [sqrt(Cm(1,1)), sqrt(Cm(2,2)) ]; % sqrt(variance)
% predicted data at a new set of z's
No=10;
zeval = zmin + (zmax-zmin)*[0:No-1]'/(No-1);
Go = [ones(No,1), zeval];
deval = Go*mest; % predicted data
Cdeval = Go*Cm*Go'; % its covariance
sdeval = sqrt(diag(Cdeval)); % its sqrt(variance)
% plot
figure(1);
clf;
% plot scale
pdmin=0;
pdmax=15;
% plot observed data, predicted data and +/- one-sigma error bounds
subplot(1,2,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [zmin, zmax, pdmin, pdmax ]' );
plot( zeval, deval+sdeval, 'b-', 'LineWidth', 3);
plot( zeval, deval-sdeval, 'b-', 'LineWidth', 3);
plot( zeval, deval, 'r-', 'LineWidth', 3);
plot( z, dobs, 'ko', 'LineWidth', 3);
for i = (1:N)
plot( [z(i), z(i)], [dobs(i)-sd, dobs(i)+sd], 'k-', 'LineWidth', 2);
end
xlabel('z');
ylabel('d');
% CASE 2: data spread out over whole interval
% z-axis
N=10;
zmin=0;
zmax=10;
z = sort(random('Uniform',zmin,zmax,N,1));
% create synthetic observed data
% d = a + b*z + noise
a=5.0;
b=0.5;
sd=1;
dobs = a+b*z+random('Normal',0,sd,N,1);
% least squares fit
M=2;
G=[ones(N,1), z]; % data kernel
mest = (G'*G)\(G'*dobs); % estimated solution
Cm = (sd^2)*inv(G'*G); % model coverance
sm = [sqrt(Cm(1,1)), sqrt(Cm(2,2)) ]; % sqrt(variance)
% predicted data at a new set of z's
No=10;
zeval = zmin + (zmax-zmin)*[0:No-1]'/(No-1);
Go = [ones(No,1), zeval]; % data kernel
deval = Go*mest; % predicted data
Cdeval = Go*Cm*Go'; % its covariance
sdeval = sqrt(diag(Cdeval)); % sqrt(variance)
% plot observed data, predicted data and +/- one-sigma error bounds
subplot(1,2,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [zmin, zmax, pdmin, pdmax ]' );
plot( zeval, deval+sdeval, 'b-', 'LineWidth', 3);
plot( zeval, deval-sdeval, 'b-', 'LineWidth', 3);
plot( zeval, deval, 'r-', 'LineWidth', 3);
plot( z, dobs, 'ko', 'LineWidth', 3);
for i = (1:N)
plot( [z(i), z(i)], [dobs(i)-sd, dobs(i)+sd], 'k-', 'LineWidth', 2);
end
xlabel('z');
ylabel('d');
fprintf('Caption: Linear fit (red) to data (circles), showing covariance of\n');
fprintf('predicted data (blue). (Left) Data clumpted together. (Right) data spread apart.\n');
% gdama05_04
% verify bordering method solution
clear all;
fprintf('gdama05_04\n');
% random (M-1)x(M-1) symmetric matrix
M=6;
S = random('Normal',0.0,1.0,M-1,M-1);
S = 0.5*(S+S');
u = random('Normal',0.0,1.0,M-1,1);
z = 0.0;
% concateate and calculate inverse
Z = [S, u; u', z];
ZINV = inv(Z);
% Bordering method
SI = inv(S);
d = u'*SI*u;
b = SI*u/d;
c = -1.0/d;
% computationally, if a matrix A is known to be symmetric
% then its best to compute it with a manifestly symmetric formula
% A = (eye(M-1)-b*u')*SI; % form in book, is not manifestly symmetric
% A = SI-b*u'*SI % multiply out SI
% A = SI-(SI*u)*(u'*SI)/d % insert u
v = SI*u; % shorthand for SI*u
A = SI-v*v'/d; % this form is manifestly symmetric
ZZINV = [A, b; b', c];
disp('full inverse');
disp(ZINV);
disp(' ');
disp('bordering inverse');
disp(ZZINV);
disp(' ');
E = abs(ZINV-ZZINV); % difference
fprintf('maximum deviation in inverse: %e\n', max(E(:)) );
% gdama05_05
% Model Resolution Matrix example, Backus-Gilbert solution
clear all;
fprintf('gdama05_05\n');
% z-axis
M=101;
zmin=0.0;
zmax=10.0;
Dz=(zmax-zmin)/(M-1);
z=zmin+Dz*[0:M-1]';
% model, m(z), moztly zero but a few spikes
mtrue = zeros(M,1);
mtrue(5)=1;
mtrue(10)=1;
mtrue(20)=1;
mtrue(50)=1;
mtrue(90)=1;
% experiment: exponential smoothing of model
N=80;
cmin=0.00;
cmax=7.90;
Dc=(cmax-cmin)/(N-1);
c = cmin + Dc*[0:N-1]';
G = exp(-c*z'); % data kernel
% create synthetic observed data
sd=0.0;
dtrue = G*mtrue;
dobs = dtrue + random('Normal',0,sd,N,1);
% construct Backus-Gilbert solution row-wise
GMG = zeros(M,N);
u = G*ones(M,1);
CHECKS=0;
for k = (1:M)
% note that S is a symmetric matrix
S = G * diag(([1:M]-k).^2) * G';
if(CHECKS)
% elelemntwise calculatoion of S as a check
Sp=zeros(N,N);
for i=(1:N)
for j=(1:N)
tmp=0;
for el=(1:M)
tmp=tmp+((el-k)^2)*G(i,el)*G(j,el);
end
Sp(i,j)=tmp;
end
end
E = abs(S-Sp);
fprintf('Error in S: %e\n',max(E(:)));
end
epsilon=1e-6;
S = S+epsilon*eye(N);
uSinv = u'/S;
GMG(k,:) = uSinv / (uSinv*u); % generalized inverse
end
mest = GMG * dobs; % estimated model parameters
Rres = GMG*G; % model resolution matrix
% plot data kernel
gda_draw(dobs,'caption dobs',' = ',G,'caption G',mtrue,'caption mtrue');
fprintf('Caption: Graphical depiction of the equation dobs = G mtrue\n');
% plot generalized inverse
gda_draw(mest,'caption mest',' = ',GMG,'caption GMG',dobs,'caption dobs');
fprintf('Caption: Graphical depiction of the equation mest = G^-g dobs\n');
% plot model resolution matrix
gda_draw(mest,'caption mest',' = ',Rres,'caption R',mtrue,'caption mtrue');
fprintf('Caption: Graphical depiction of the equation mest = R mtrue\n');
% plot data resolution matrix
dpre = G*mest;
Nres = G*GMG;
gda_draw(dpre,'caption dpre',' = ',Nres,'caption N',dobs,'caption dobs');
fprintf('Caption: Graphical depiction of the equation dpre = N dobs\n');
% plot
figure();
clf;
% plot scale
pmmin=-0.2;
pmmax=1.0;
% plot true and estimated model
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [zmin, zmax, pmmin, pmmax ]' );
plot( z, mtrue, 'r-', 'LineWidth', 3);
plot( z, mest, 'b-', 'LineWidth', 3);
xlabel('z');
ylabel('m');
fprintf('Caption: True solution (red) and Backus-Gilbert estimated solution (blue)\n');
% gdmaa05_06
% trade off of resolution and variance
fprintf('gdmaa05_06\n');
% z-axis
M=101;
zmin=0.0;
zmax=10.0;
Dz=(zmax-zmin)/(M-1);
z=zmin+Dz*[0:M-1]';
% experiment: exponential smoothing of model
N=80;
cmin=0.00;
cmax=7.90;
Dc=(cmax-cmin);
c = cmin + Dc*[0:N-1]'/(N-1);
G = exp(-c*z'); % data kernel
% Part 1: Backus-Gilbert
% I hand-tuned the points on trade-off curve to make the plot look nice
alphavec = [0.999, 0.995, 0.99, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001]';
iahalf=find(alphavec==0.5);
Na=length(alphavec); % number of points on the trade-off curve
spreadR1=zeros(Na,1); % tabulate spread of resolution
sizeC1=zeros(Na,1); % tabulate size of covariance
for a=(1:Na)
alpha = alphavec(a);
% construct Backus-Gilbert solution row-wise
GMG = zeros(M,N);
u = G*ones(M,1);
for k = (1:M)
S = G * diag(([1:M]-k).^2) * G';
Sp = alpha*S + (1-alpha)*eye(N,N);
uSpinv = u'/Sp;
GMG(k,:) = uSpinv / (uSpinv*u); % generalized inverse
end
Cm1 = GMG*GMG'; % unit covariance matrix
R1 = GMG*G; % model resolution matrix
sizeC1(a) = sum(diag(Cm1)); % size of covariacne
% spread of resolution
spreadR1(a) = sum(sum( (([1:M]'*ones(1,M)-ones(M,1)*[1:M]).^2).*(R1.^2) ));
CHECKR=0;
if( CHECKR )
% element-wise, to check result
spreadR=0;
for i = (1:M)
for j = (1:M)
spreadR = spreadR + ((i-j)^2)*(R1(i,j)^2);
end
end
fprintf('Error in R: %e\n', abs(spreadR-spreadR1(a)))
end
end
% plot
figure(1);
clf;
subplot(1,2,1);
% plot trade-off curve
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [2, 3.5, -3, 5 ]' );
plot( log10(spreadR1), log10(sizeC1), 'k-', 'LineWidth', 4);
plot( log10([spreadR1(1), spreadR1(iahalf), spreadR1(Na)]') , log10([sizeC1(1), sizeC1(iahalf), sizeC1(Na)]'), 'ko', 'LineWidth', 2);
xlabel('log10 spread(R)');
ylabel('log10 size(Cm)');
% Part 2: Damped minimum length
% I hand-tuned the points on trade-off curve to make the plot look nice
alphavec = [0.9999, 0.999, 0.995, 0.99, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.01]';
iahalf=find(alphavec==0.5);
Na=length(alphavec); % number of points on the trade-of curve
spreadR2=zeros(Na,1); % tabulate spread of reolution
sizeC2=zeros(Na,1); % tabulate size of covariance
for a=(1:Na)
alpha = alphavec(a);
GMG = (alpha*G')/(alpha*(G*G')+ (1-alpha)*eye(N,N)); % generalized inverse
Cm2 = GMG*GMG'; % model covariance
R2 = GMG*G; % model resolution
sizeC2(a) = sum(diag(Cm2)); % size of covariance
spreadR2(a) = sum(sum((R2-eye(M)).^2)); % spread of resulution
end
% plot trade-off curve
subplot(1,2,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [90, 100, -3, 4 ]' );
plot( (spreadR2), log10(sizeC2), 'k-', 'LineWidth', 4);
% plot( (spreadR2), log10(sizeC2), 'ro', 'LineWidth', 2);
plot( ([spreadR2(1), spreadR2(iahalf), spreadR2(Na)]') , log10([sizeC2(1), sizeC2(iahalf), sizeC2(Na)]'), 'ko', 'LineWidth', 2);
xlabel('spread(R)');
ylabel('log10 size(Cm)');
fprintf('Caption: Trade-off curves, (left) L2, (rith) Backus-Gilbert\n');
% gdama05_07
% simple figures of rays crossing a box
clear all;
fprintf('gdama05_07\n');
% Part 1: coarse grid
figure();
clf;
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
% box has auxially variables (x,y)
xmin=0.0;
xmax=2.0;
ymin=0.0;
ymax=2.0;
plot( [xmin, xmax, xmax, xmin, xmin]', [ymin, ymin, ymax, ymax, ymin], 'k-', 'LineWidth', 2);
% box contains a circular object
Nth=101;
R=0.9;
x0 = 1.0;
y0 = 1.0;
th = 2*pi*[0:Nth]'/(Nth-1);
x = x0+R*sin(th);
y = y0+R*cos(th);
axis( [xmin, xmax, ymin, ymax]' );
axis off;
axis equal;
plot( x, y, 'm-', 'LineWidth', 4 );
% rays are generated randomly
Nr=21;
x1 = random('Uniform', xmin, xmax, Nr, 1);
x2 = random('Uniform', xmin, xmax, Nr, 1);
for i=(1:Nr)
plot( [x1(i), x2(i)]', [ymin, ymax]', 'b-', 'LineWidth', 2);
end
y1 = random('Uniform', ymin, ymax, Nr, 1);
y2 = random('Uniform', ymin, ymax, Nr, 1);
for i=(1:Nr)
plot( [xmin, xmax]', [y1(i), y2(i)]', 'b-', 'LineWidth', 2);
end
% (x,y) grid is superimosed on the box
Ng = 7;
xg = xmin + (xmax-xmin)*[0:Ng-1]'/(Ng-1)';
yg = ymin + (ymax-ymin)*[0:Ng-1]'/(Ng-1)';
for i=(2:Ng-1)
plot( [xg(1), xg(Ng)], [yg(i), yg(i)]', 'k-', 'LineWidth', 2 );
end
for i=(2:Ng-1)
plot( [xg(i), xg(i)]', [yg(1), yg(Ng)], 'k-', 'LineWidth', 2 );
end
fprintf('Caption: Rays crossing a coarse grid\n');
% Part 2: fine grid
figure();
clf;
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
% box has auxially variables (x,y)
xmin=0;
xmax=2;
ymin=0;
ymax=2;
plot( [xmin, xmax, xmax, xmin, xmin]', [ymin, ymin, ymax, ymax, ymin], 'k-', 'LineWidth', 2);
% box contains a circular object
Nth=101;
R=0.9;
x0 = 1.0;
y0 = 1.0;
th = 2*pi*[0:Nth]'/(Nth-1);
x = x0+R*sin(th);
y = y0+R*cos(th);
axis( [xmin, xmax, ymin, ymax]' );
axis off;
axis equal;
plot( x, y, 'm-', 'LineWidth', 4 );
% use same rays as in Part 1
for i=(1:Nr)
plot( [x1(i), x2(i)]', [ymin, ymax]', 'b-', 'LineWidth', 2);
end
for i=(1:Nr)
plot( [xmin, xmax]', [y1(i), y2(i)]', 'b-', 'LineWidth', 2);
end
% (x,y) grid is superimosed on the box
Ng = 21;
xg = xmin + (xmax-xmin)*[0:Ng-1]'/(Ng-1)';
yg = ymin + (ymax-ymin)*[0:Ng-1]'/(Ng-1)';
for i=(2:Ng-1)
plot( [xg(1), xg(Ng)], [yg(i), yg(i)]', 'k-', 'LineWidth', 2 );
end
for i=(2:Ng-1)
plot( [xg(i), xg(i)]', [yg(1), yg(Ng)], 'k-', 'LineWidth', 2 );
end
fprintf('Caption: Rays crossing a fine grid\n');
% gdama05_08
clear all;
fprintf('gdama05_08\n')
% Constructing just one column of the resolution matrix
% Note:
% A row of the resolution matrix R tells you the contributions
% that all the true model parameters make to a single estimated
% model parameter; that is, the estimated model parameter is
% a weighted average of the true model parameter, where the
% elements of the row gives the weights.
%
% On the other hand, a column of the resolution matrix tells
% you all the estimated model parameters that are influenced
% by a true model parameter. The column gives the "blur"
% of estimated model parameters that would be observed if
% the true model consisted of just a single spike.
%
% The distinction between row and column is important,
% since R is, in general, not symmetric.
%
% inverse problem considered here is an acoustic tomography
% problem, where the observations are along just rows and
% columns
% grid of unknowns is Lx by Ly
Lx = 20;
Ly = 20;
M = Lx*Ly;
% observations only along rows and columns
N=Lx+Ly;
% build backward index tables for convenience
iofk=zeros(M,1); % backward index table, i(k)
jofk=zeros(M,1); % backward index table, j(k)
kofij=zeros(Lx,Ly);
k=0;
for i=(1:Lx)
for j=(1:Ly)
k = k+1;
iofk(k)=i;
jofk(k)=j;
kofij(i,j)=k;
end
end
G=zeros(N,M);
% observations across rows
for i=(1:Lx)
for j=(1:Ly)
k = kofij(i,j);
G(i,k)=1;
end
end
% observations across columns
for j=(1:Ly)
for i=(1:Lx)
k = kofij(i,j);
G(j+Lx,k)=1;
end
end
% note that this problem is actually mix-determined
% since the sum of all the horizontal ray traveltimes
% equals the sum of all the vertical ray traveltimes
% so use the damped minimum-length solution when
% computing the solution
epsi = 0.0001;
GMG = G'/(G*G'+(epsi^2)*eye(N,N));
% compute the complete resolution matrix for comparison
% note that R is symmetric in this case: R = GMG*G =
% G' inv(G*G'+(epsi^2)*eye(N,N)) * G is symmetric, since
% the inverse of a symmetric matrix is symmetric
Rres = GMG*G;
% pull out just one column of the resolution matrix, one
% corresponding to a true model parameter near the middle of
% the model volume, at (ix, iy)
i=floor(Lx/2);
j=floor(Lx/2);
icolB=kofij(i,j);
RicolB=zeros(Lx,Ly);
for k=(1:M)
RicolB(iofk(k),jofk(k))=Rres(k,icolB);
end
% now construct just that column
% model parameter with unity in that col
mk = zeros(M,1);
mk(icolB) = 1;
% data it predicts
dk=G*mk;
% solve inverse problem, interpret the result as
% a column of the resolution matrix. In this case, I
% solve the inverse problem using the generalized
% inverse. But it could as well have been solved
% iteratively, using biconjugate gradients, using
% the following trick: First write
% rk = GMG dk = GT inv(GGT + e2 I) dk = GT x
% with x = inv(GGT + e2 I) dk or dk = (GGT + e2 I) x.
% Now solve dk = (GGT + e2 I) x
% with biconjugate gradients and then
% rk = GT x;
rk = GMG*dk;
% reorganize to 2D physical model space
Rk=zeros(Lx,Ly);
for k=(1:M)
Rk(iofk(k),jofk(k))=rk(k);
end
gda_draw(' ',RicolB,' ',Rk);
fprintf('Caption: One row of resolution matrix R, reorganized into an image\n');
fprintf('(Left) Just one row computed. (Right) All R computed.\n')
% gdama05_09
% Checkerboard resolution test
% inverse problem is an acoustic tomography problem, where
% the observations are along just rows and columns
clear all;
fprintf('gdama05_09\n');
% grid of unknowns is Lx by Ly
Lx = 20;
Ly = 20;
M = Lx*Ly;
% observations only along rows and columns
N=Lx+Ly;
% A pixel (i,j) of the Lx by Ly image S(i,j) maps to an element
% k of the model parameter vector m(k). There are
% a variet of ways to determine the mappings, k(i,j),
% and (i(k),j(k)). One is to use the formulas,
% k = Ly*(i-1)+j;
% i = floor((k-1)/Ly)+1;
% j = k - Ly*(i-1);
% these formulas work OK in most instances, but an alternative method
% is to build forwardb and backward index tables for convenience (which
% is what is done here. The advantage for the tables is that they are
% easy to generalize to higer dimensions and they are faster to evaluate
% than the formulas. A disadvantage is that they consume 2*M+Lx*Ly integers
% of memory
iofk=zeros(M,1); % backward index table, i(k)
jofk=zeros(M,1); % backward index table, j(k)
kofij=zeros(Lx,Ly); % forward index table, k(i,j)
k=0;
for i=(1:Lx)
for j=(1:Ly)
k=k+1;
iofk(k)=i;
jofk(k)=j;
kofij(i,j)=k;
end
end
% check consistency of formulas and index arrays
CHECK = 0;
if( CHECK )
E = 0;
for k=(1:M)
i = floor((k-1)/Ly)+1;
j = k - Ly*(i-1);
E = E + abs(i-iofk(k)) + abs(j-jofk(k));
end
fprintf('Error of k to (i,j): %d\n', E);
E = 0;
for i=(1:Lx)
for j=(1:Ly)
k=Ly*(i-1)+j;
E = E + abs(k-kofij(i,j));
end
end
fprintf('Error of (i,j) to k: %d\n', E);
end
G=zeros(N,M);
% observations across rows
for i=(1:Lx)
for j=(1:Ly)
k = kofij(i,j); % map model parameter at (i,j) into scalar index k
G(i,k)=1;
end
end
% observations across columns
for j=(1:Ly)
for i=(1:Lx)
k = kofij(i,j); % map model parameter at (i,j) into scalar index k
G(j+Lx,k)=1;
end
end
% note that this problem is actually mix-determined
% since the sum of all the horizontal ray traveltimes
% equals the sum of all the vertical ray traveltimes
% so use the damped minimum-length solution when
% computing the solution
epsi = 0.0001;
GMG = G'/(G*G'+(epsi^2)*eye(N,N));
% model parameter vector mk for crude checkerboard
mk = zeros(M,1);
for i=(1:4:Lx)
for j=(1:4:Ly)
k=kofij(i,j);
mk(k) = 1;
end
end
% data predicted by mk
dk=G*mk;
% solve inverse problem, interpret the result as
% a row of the resolution matrix.
mkest = GMG*dk;
% reorganize to 2D physical model space
Strue=zeros(Lx,Ly); % true checkerboard
Sest=zeros(Lx,Ly); % resulting checkerboard
for k=(1:M)
Strue(iofk(k),jofk(k))=mk(k);
Sest(iofk(k),jofk(k))=mkest(k);
end
gda_draw(Strue,'caption true', Sest,'caption est');
fprintf('Caption: (Left) True model with array of dots. (Right) Estimated model.\n');
% gdama05_10
% Backus-Gilbert matrix S for 2D weight function
clear all;
fprintf('gdama05_10\n');
% (x,y) grid
Lx=101;
xmin=0.0;
xmax=xmin+Lx-1;
x = xmin+(xmax-xmin)*[0:Lx-1]'/(Lx-1);
Ly=101;
ymin=0.0;
ymax=ymin+Ly-1;
y = ymin+(ymax-ymin)*[0:Ly-1]'/(Ly-1);
M = Lx*Ly;
% index tables for convenience
iofk=zeros(M,1); % backward table, i(j)
iofk=zeros(M,1); % backward table, j(k)
kofij=zeros(Lx,Ly); % forward table k(i,j)
k=0; % counter
for i=(1:Lx)
for j=(1:Ly)
k=k+1;
iofk(k)=i;
jofk(k)=j;
kofij(i,j)=k;
end
end
% coordinates of model parameters
xm = zeros(M,1);
ym = zeros(M,1);
for i = (1:Lx)
for j = (1:Ly)
k = kofij(i,j);
xm(k) = x(i);
ym(k) = y(j);
end
end
% random data kernel
N=10;
G = random('Normal',0.0,1.0,N,M);
k=kofij(50,30); % target row of S
wk=(xm-xm(k)).^2+(ym-ym(k)).^2; % weights
S = G*diag(wk)*G';
gda_draw(S,'caption S');
fprintf('Caption: The matrix S=G*diag(w)*GT\n');
% display weights on (x,y) grid
W = zeros(Lx,Ly);
for k =(1:M)
W(iofk(k),jofk(k)) = wk(k);
end
gda_draw(W,'caption W');
fprintf('Caption: The weight matrix W\n');