% gdama15
% 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-08, Edit and integration with text
% 2023-05-08, Some patches
% gdama15_01
% image deblurring example
clear all;
clear gdaGsparse gdaepsilon;
global gdaGsparse gdaepsilon;
fprintf('gdama15_01\n');
figure();
clf;
% black & white image
c = colormap( [[0:255]', [0:255]',[0:255]'] / 255 );
% read in true image and demean it
mtrue = double(imread('../Data/tern.tif', 'tif' ));
mtrue = mtrue-mean(mean(mtrue));
[I, J] = size( mtrue );
% plot true image
subplot(2,3,1);
set(gca,'LineWidth',2);
hold on;
axis( [1, I, 1, J] );
axis ij;
imagesc([1, I], [1, J], mtrue);
title('true image');
% data is blurred version of image
mblurred = zeros( I, J );
% blur is over rows of image with this filter
Lb = 100;
% try 2 kinds of blurs
if( 0 )
blur = [Lb:-1:1]'/Lb; % linear ramp-down
else
blur = ones(Lb,1)/Lb; % boxcar
end
% crease sparse data kernel
% in this case, one colud use the simple command
% G = spdiags( ones(J,1)*blur', [0:Lb-1], J, J );
% however, thst only works when the blur is parallel
% to the rows or columns of the image. So, instead
% I use the row-column-value method
NELEMENTS = Lb*J+20;
Gr = zeros(NELEMENTS,2);
Gc = zeros(NELEMENTS,2);
Gv = zeros(NELEMENTS,2);
Nel=0; % element counter
for i=(1:J) % loop over model pixel
for j=(i:min(i+Lb-1,J))
Nel=Nel+1;
Gr(Nel) = i;
Gc(Nel) = j;
Gv(Nel) = blur(j-i+1);
end
end
Gr=Gr(1:Nel); Gc=Gc(1:Nel); Gv=Gv(1:Nel); % truncate
gdaGsparse = sparse(Gr,Gc,Gv,J,J); % create sparse matrix
clear Gr Gc Gv; % clear index arrays
% data is blurred version of image
for i=(1:I)
dobs(i,:) = (gdaGsparse*(mtrue(i,:)'))';
end
% plot blurred image
subplot(2,3,2);
set(gca,'LineWidth',2);
hold on;
axis( [1, I, 1, J] );
axis ij;
imagesc([1, I], [1, J], dobs);
title('blurred image');
gdaepsilon=1.0e-6; % damping, just in case GGT is singular@
% damped minimum length solution
mest = zeros(I,J);
tol = 1e-4; % tolerance
maxit = 3*J; % maximum number of iterations
DOBICG = 0;
if(DOBICG) % I'm not happy with performance of bicg() ...
fprintf('This part slow, so show progress:\n',i,I);
% m = GT inv(G GT + eI) d so m = GT u with (G GT + eI) u = d
for i = (1:I) % process image row-wise
if( mod(i,100)==0 )
fprintf('%d out of %d\n',i,I);
end
dobsrow = dobs(i,:)';
% the flag supresses messages written to the screen
[u,flag] = bicg(@gda_dmlmul,dobsrow,tol,maxit);
mestrow = gdaGsparse'*u;
mest(i,:)=mestrow';
end
else
G=full(gdaGsparse);
GGT = G*G' + gdaepsilon*eye(J);
GMG = G'/GGT;
for i = (1:I) % process image row-wise
dobsrow = dobs(i,:)';
mestrow = GMG*dobsrow;
mest(i,:)=mestrow';
end
end
% plot reconstructed image
subplot(2,3,3);
set(gca,'LineWidth',2);
hold on;
axis( [1, I, 1, J] );
axis ij;
imagesc([1, I], [1, J], mest);
title('reconstructed image');
% plot closeup
I1=400;
I2=800;
J1=500;
J2=900;
subplot(2,3,4);
set(gca,'LineWidth',2);
hold on;
axis( [I1, I2, J1, J2] );
axis ij;
imagesc([I1, I2], [J1, J2], mtrue(I1:I2,J1:J2));
subplot(2,3,5);
set(gca,'LineWidth',2);
hold on;
axis( [I1, I2, J1, J2] );
axis ij;
imagesc([I1, I2], [J1, J2], dobs(I1:I2,J1:J2));
subplot(2,3,6);
set(gca,'LineWidth',2);
hold on;
axis( [I1, I2, J1, J2] );
axis ij;
imagesc([I1, I2], [J1, J2], mest(I1:I2,J1:J2));
fprintf('Caption: True image (left), blurred image (center) and reconstructed image (right)\n');
% rather brute force here
% plot central row of generalized inverse and resolution matrix
Gmg = gdaGsparse'/(gdaGsparse*gdaGsparse'+gdaepsilon*speye(J,J));
R = Gmg*gdaGsparse;
gda_draw(full(gdaGsparse),'caption G', Gmg, 'caption Gmg', R, 'caption R');
fprintf('Caption: Data kernel G (left), generlized inverse G^-g (middle),\n');
fprintf('Model resolution matrix (right)\n');
figure();
clf;
hold on;
subplot(2,1,1);
set(gca,'LineWidth',2);
tmp = Gmg( J/2, : )';
axis( [0, J, min(tmp), max(tmp)] );
xlabel('row');
ylabel('generalized inverse');
plot( [1:J]', tmp, 'k-', 'LineWidth', 2 );
fprintf('Caption: Central row of the generalized inverse\n');
figure();
clf;
hold on;
set(gca,'LineWidth',2);
hold on;
tmp = R( J/2, : )';
axis( [0, J, min(tmp), max(tmp)] );
xlabel('row i');
ylabel('R(i)');
plot( [1:J]', tmp, 'k-', 'LineWidth', 2 );
fprintf('Caption: Central row of the resolution matrix\n');
% some statistics of R
% Note: for some reason, max() returns Rmax as a 1x1 sparse
% matrix and fprint refuses to print it unless it is converted
% to a full matrix by a call to full()
[Rmax, jRmax] = max(R(J/2,:));
fprintf('R(i=%d, j) had maximum %f at j=%d\n', J/2, full(Rmax), jRmax );
Rpmax = max(abs([R(J/2,1:jRmax-1), R(J/2,jRmax+1:J)]));
fprintf('largest off diaginal element is %f as large\n', full(Rpmax)/Rmax );
% gdama15_02
% deconvolution example, d=g*m and m=ginv*d
% this uses generalized least squares
% the data equation is G m = rhs
% the constraint equation is H m = h
% the combined equations are Fm = f
% solution uses the bicg() solver set up to
% solve FTF = FT f, but in a clever way so that
% FTF, G and f are never explicitly
% however, the code also implements the brute
% force solution m = (FTF)\(FT f) inside if statements
% for test purposes
clear all;
fprintf('gdama15_02\n');
clear gdag gdaHsparse;
global gdag gdaHsparse;
D=load('../data/airgun.txt');
% digitized airgun pulse, data after:
% Smith, SG,
% Measurement of Airgun Waveforms,
% Geophys. J. R. astr. Soc. (1975) 42, 273-280.
t=D(:,1); % time
d=D(:,2); % this is the filter for which we seek an inverse filter
Mo = length(t);
tmin = t(1);
Dt = t(2)-t(1);
% pad with zeros
M=2*Mo;
t = tmin+Dt*[0:M-1];
gdag = [ d', zeros(1,M-Mo) ]';
tmax = t(M);
% pad with zeros
M=length(t);
N=M; % inverse filter same length as filter
figure();
clf;
hold on;
% plot g
set(gca,'LineWidth',2);
gmax = max(abs(gdag));
axis( [tmin, tmax, -gmax, gmax] );
plot( t, gdag, 'k-', 'LineWidth', 3 );
xlabel('time, t');
ylabel('g(t)');
fprintf('Caption: Airgun pulse\n');
% don't need G explicitly, but here it is for test purposes
test=1;
if( test )
G=zeros(N,M);
for j=(1:M)
G(j:N,j)=gdag(1:N-j+1);
end
end
% set up H and h
K=2*M;
L=N+K;
% use row-col-value method to make sparse matrix H
% an alternative would be the command
% gdaH=spalloc(K,M,3*L);
Hr = zeros(3*L,1);
Hc = zeros(3*L,1);
Hv = zeros(3*L,1);
Nel=0;
% two sets of constaint equations, curvature and length
% epsilon parameters adjusted empirically
% curvature
e=0.001*max(abs(gdag));
for j = (2:M)
% gdaH(j,j-1)=e;
Nel=Nel+1;
Hr(Nel) = j;
Hc(Nel) = j-1;
Hv(Nel) = e;
end
for j = (1:M)
% gdaH(j,j)=-2*e;
Nel=Nel+1;
Hr(Nel) = j;
Hc(Nel) = j;
Hv(Nel) = -2*e;
end
for j = (1:M-1)
% gdaH(j,j+1)=e;
Nel=Nel+1;
Hr(Nel) = j;
Hc(Nel) = j+1;
Hv(Nel) = e;
end
% damp size, especially towards the end
e2=0.1*max(abs(gdag));
for j = (1:M)
% gdaH(j+M,j)=e2*sqrt(j);
Nel=Nel+1;
Hr(Nel) = j+M;
Hc(Nel) = j;
Hv(Nel) = e2*sqrt(j);
end
% truncate to actual length used
Hr = Hr(1:Nel);
Hc = Hc(1:Nel);
Hv = Hv(1:Nel);
gdaHsparse = sparse(Hr,Hc,Hv,K,M);
clear Hr Hc Hv;
% rhs of data equaltion is a spike delayed by Nd samples
% (a little delay sometimes leads to a better result)
rhs=zeros(N,1);
Nd=20;
rhs(Nd+1)=1;
% rhs of constraint equation is zero
h=zeros(K,1);
if( test ) % don't need F, but here it is for test purposes
F=[ G; full(gdaHsparse) ];
end
if( test ) % don't need f, but here it is for test purposes
f=zeros(L,1);
f(1:N)=rhs;
f(N+1:N+K)=h;
end
if( test ) % don't need brute force solution, but here it is for test purposes
ginvbf = (F'*F)\(F'*f);
end
% don't need GT*rhs, but here it is for test purposes
if( test )
Gtrhs = G'*rhs;
end
% set up F'f = GT rhs + HT h
% note that GT v is similar to the convolution G v
% except that g is time-reversed and position of
% results in convolution is shifted to the bottom
FTf = gda_filterrhs(rhs,h);
% solve, filterfun does the (FT F v) multiplication cleverly
% so that G, GTG, HTH and f never computed
ginv=bicg(@gda_filtermul,FTf,1e-10,3*L);
% don't need G*ginv, but here it is for test purposes
if( test )
Gginv = G*ginvbf;
end
% plot inverse filter
figure();
clf;
hold on;
set(gca,'LineWidth',2);
ginvmax = max(abs(ginv));
axis( [tmin, tmax, -ginvmax, ginvmax] );
plot( t, ginv, 'k-', 'LineWidth', 3 );
if( test ) % plot brute force result
plot( t, ginvbf, 'r-', 'LineWidth', 2 );
end
xlabel('time, t');
ylabel('g^{inv}(t)');
fprintf('Caption: inverse filter to make airgun pulse more spiky\n');
% plot convolution of inverse filter and filter
c = conv( ginv, gdag );
c = c(1:N);
figure();
clf;
hold on;
set(gca,'LineWidth',2);
cmax = max(abs(c));
axis( [tmin, tmax, -cmax, cmax] );
plot( t, c, 'k-', 'LineWidth', 3 );
if( test ) % plot brute force result
plot( t, Gginv, 'r-', 'LineWidth', 2 );
end
xlabel('time, t');
ylabel('g^{inv}(t)*g(t)');
fprintf('Caption deconvolved (spiked) airgun pulse\n');
% error
E = (rhs-c)'*(rhs-c);
fprintf('damping factors %f %f prediction error %f\n', e, e2, E );
% as an example, deconvolve three pulses with additive noise
sd=0.1;
glong = [3*gdag', 2*gdag', gdag']' + random('Normal',0,sd, 3*M, 1);
clong = conv( ginv, glong );
clong = clong(1:3*N);
figure()
clf;
% plot data
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
glongmax = max(abs(glong));
axis( [tmin, tmin+3*tmax, -glongmax, glongmax] );
plot( tmin+Dt*[1:3*M]', glong, 'k-', 'LineWidth', 2 );
xlabel('time, t');
ylabel('d(t)');
% plot deconvolved data
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
clongmax = max(abs(clong));
axis( [tmin, tmin+3*tmax, -clongmax, clongmax] );
plot( tmin+Dt*[1:3*M]', clong, 'k-', 'LineWidth', 2 );
xlabel('time, t');
ylabel('g^{inv}(t)*d(t)');
fprintf('Caption: (top) Time series with three pulses (bottom) deconvolved timeseries\n');
% gdama15_03
% correcting gravity cross-over errors
clear all;
fprintf('gdama15_03\n');
% load gravity(lat,lon) data
D=load('../data/gravity.txt');
lon=D(:,1);
lat=D(:,2);
grav=D(:,3);
% setup axes
mingrav = min(grav);
maxgrav = max(grav);
meangrav = mean(grav);
Ng=length(lon);
lonmin=min(lon);
lonmax=max(lon);
latmin=min(lat);
latmax=max(lat);
Nlon=length(find(lat==lat(1)));
Nlat=length(find(lon==lon(1)));
Dlon = (lonmax-lonmin)/(Nlon-1);
Dlat = (latmax-latmin)/(Nlat-1);
% I tested the coordinate system by
% resetting a specific datum to a value
% that I could easily pick out, and then
% comparing its position on the plot
% with the position displayed below
if( 0 )
i=1500;
fprintf('%f %f\n', lon(i), lat(i) );
grav(i)=mingrav;
end
% put data into grid
% Note the messing about with trasnposes and
% flipup() to make it come out rightside up
LON=reshape(lon,Nlon,Nlat);
LAT=reshape(lat,Nlon,Nlat);
GRAV=reshape(grav,Nlon,Nlat);
LAT=flipud(LAT');
LON=flipud(LON');
GRAV=flipud(GRAV');
% examne order
debg = 0;
if( debg )
fprintf('LON(1,1)=%f LON(1,N)=%f\n', LON(1,1), LON(1,end) );
fprintf('LON(N,1)=%f LON(N,N)=%f\n', LON(end,1), LON(end,end) );
fprintf('\n');
fprintf('LAT(1,1)=%f LAT(1,N)=%f\n', LAT(1,1), LAT(1,end) );
fprintf('LAT(N,1)=%f LAT(N,N)=%f\n', LAT(end,1), LAT(end,end) );
fprintf('\n');
end
% at this point, ordering is
% LON(1,1)=0.041667 LON(1,N)=9.958330
% LON(N,1)=0.041667 LON(N,N)=9.958330
%
% LAT(1,1)=0.041667 LAT(1,N)=0.041667
% LAT(N,1)=9.958330 LAT(N,N)=9.958330
% set a small (i,1) block of Grav(1,j) to a low value
% and visually check where it plots on the map to
% confirm that the ordering is right
debg = 0;
if( debg )
GRAV(1:10,1:20)=mingrav;
end
% plot gravity data
figure(1);
clf;
colormap('jet');
% plot true gravity
subplot(2,2,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [lonmin, lonmax, latmin, latmax] );
axis xy;
imagesc([lonmin, lonmax], [latmin, latmax], GRAV);
title('true gravity');
% plot gravity and prepare to superimpose tracks
subplot(2,2,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [lonmin, lonmax, latmin, latmax] );
axis xy;
imagesc([lonmin, lonmax], [latmin, latmax], GRAV);
title('tracks');
% the tracks are synthetic; that is, constructed in this program
Ntracks = 35;
Ntlon = Nlon;
Dtlon = (lonmax-lonmin)/(Ntlon-1);
tlon = lonmin + Dtlon*[0:Ntlon-1]';
% construct ascending tracks
NA = Ntracks;
Atracklat = zeros(Ntlon, NA);
C = 1;
L = 2*(lonmax-lonmin);
for i = (1:NA)
Atracklat(:,i) = tlon + C*cos(2*pi*tlon/L) + i/2 - 9;
plot( tlon, Atracklat(:,i), 'k-', 'LineWidth', 2 );
end
% construct descending tracks
ND = Ntracks;
Dtracklat = zeros(Ntlon, NA);
C = 1;
L = 2*(lonmax-lonmin);
for i = (1:NA)
Dtracklat(:,i) = latmax - (tlon + C*cos(2*pi*tlon/L) + i/2 -9);
plot( tlon, Dtracklat(:,i), 'k-', 'LineWidth', 2 );
end
M = NA+ND; % number of model parameters are number of tracks
% model parameters are gravity shifts (offets) along each track.
% they are randomply chosen.
% store the shifts in vector with ascending tracks first
sm=0.1;
mtrue = random('Normal',0,sm*(maxgrav-mingrav)/2,M,1);
% create a dataset of all the data beneath the tracks,
% but offset by the shifts
Nsamples=0;
slonv = zeros(M*Ntlon,1);
slatv = zeros(M*Ntlon,1);
sgrav = zeros(M*Ntlon,1);
slon = tlon;
for i = (1:NA) % ascending tracks
slat = Atracklat(:,i);
ii = floor((slat-latmin)/Dlat)+1;
jj = floor((slon-lonmin)/Dlon)+1;
for kk = (1:Ntlon)
if( (ii(kk)>=1) && (ii(kk)<=Nlat) && (jj(kk)>=1) && (jj(kk)<=Nlon) )
Nsamples=Nsamples+1;
slonv(Nsamples) = slon(kk);
slatv(Nsamples) = slat(kk);
sgrav(Nsamples) = GRAV(ii(kk),jj(kk))+mtrue(i);
end
end
end
for j = (1:ND) % decending tracks
slat = Dtracklat(:,j);
ii = floor((slat-latmin)/Dlat)+1;
jj = floor((slon-lonmin)/Dlon)+1;
for kk = (1:Ntlon)
if( (ii(kk)>=1) && (ii(kk)<=Nlat) && (jj(kk)>=1) && (jj(kk)<=Nlon) )
Nsamples=Nsamples+1;
slonv(Nsamples) = slon(kk);
slatv(Nsamples) = slat(kk);
sgrav(Nsamples) = GRAV(ii(kk),jj(kk))+mtrue(NA+j);
end
end
end
% interpolate data along tracks to a 2D image
% Note: the code was originally
% ginterp = TriScatteredInterp( slonv, slatv, sgrav );
% however called complained that there were duplicate
% points. This patch explicitly removes duplicates
sorted = sortrows( [slonv, slatv, sgrav], [1,2] );
Nsort=length(sorted);
flag = ones(Nsort,1);
for i=(2:Nsort)
if( (sorted(i,1)==sorted(i-1,1)) && (sorted(i,2)==sorted(i-1,2)) )
flag(i)=0;
end
end
k = find( flag==1 );
ginterp = TriScatteredInterp( sorted(k,1), sorted(k,2), sorted(k,3) );
obsgrav = ginterp(LON(:),LAT(:));
OBSGRAV = reshape(obsgrav,Nlon,Nlat);
% plot gravity data (without cross-0ver corrections)
subplot(2,2,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [lonmin, lonmax, latmin, latmax] );
axis xy;
imagesc([lonmin, lonmax], [latmin, latmax], OBSGRAV);
title('obs gravity');
% find all track intersections (cross-overs)
subplot(2,2,2);
% Xlon(i,j) is lon where ascending track i and descending track j cross
Xlon = NaN(NA,ND);
% Xlat(i,j) is lat where ascending track i and descending track j cross
Xlat = NaN(NA,ND);
% Xgrav(i,j) is grav where ascending track i and descending track j cross
Xgrav = NaN(NA,ND);
% iiX, jjX are nearest indices in GRAV(iiX,iiJ) array
iiX = zeros(NA,ND);
jjX = zeros(NA,ND);
count=0;
for i = (1:NA)
for j = (1:ND)
signAD = sign(Atracklat(:,i)-Dtracklat(:,j));
sign1 = signAD(1);
k2 = find( signAD ~= sign1, 1 );
if( length(k2) == 0 ) % no intersection
continue;
end
if( k2 == 1 ) % intersects at left edge; ignore
continue;
end
% linearly interpolate
k1 = k2-1;
A = Atracklat(k1,i)-Dtracklat(k1,j);
B = Atracklat(k2,i)-Dtracklat(k2,j);
slope = (B-A)/Dtlon;
deltalon = -A/slope;
Xlon(i,j) = tlon(k1)+deltalon;
Xlat(i,j) = Xlon(i,j) + C*cos(2*pi*Xlon(i,j)/L) + (i/2-9);
jjXp=floor((Xlon(i,j)-lonmin)/Dlon)+1;
iiXp=floor((Xlat(i,j)-latmin)/Dlat)+1;
if( (iiXp>=1) && (iiXp<=Nlat) && (jjXp>=1) && (jjXp<=Nlon) )
count=count+1;
iiX(i,j) = iiXp;
jjX(i,j) = jjXp;
% check if these are right by making a dot on the
% image and comparing a displayed version of the lat, lon
% with the position of the dot
debg = 0;
if( debg && (count==20) )
GRAV( iiXp, jjXp ) = mingrav;
fprintf('lon %f lat %f\n', Xlon(i,j), Xlat(i,j) );
end
Xgrav(i,j) = GRAV( iiXp, jjXp );
% plot intersection point on map
debg = 0;
if( debg )
plot( Xlon(i,j), Xlat(i,j), 'ko', 'LineWidth', 2);
end
end
end
end
debg = 0;
if( debg ) % part of a check that indiced iXlon and iXlat are correct
% the position of the little spot in the image should match
% the (lon, lat) position displayed
figure(3);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [lonmin, lonmax, latmin, latmax] );
axis xy;
imagesc([lonmin, lonmax], [latmin, latmax], GRAV);
title('for debugging');
end
N = count; % number of data is number of intersections
% gravity values for each track at its intersection point
Agrav = zeros(N,1);
Dgrav = zeros(N,1);
% back indices to i,j
iA = zeros(N,1);
jB = zeros(N,1);
k=0;
for i = (1:NA)
for j = (1:ND)
if( iiX(i,j) == 0 )
continue;
end
k=k+1;
Agrav(k) = Xgrav(i,j)+mtrue(i);
Dgrav(k) = Xgrav(i,j)+mtrue(j+NA);
iA(k)=i;
jD(k)=j;
end
end
% data equations are that, at an intersection point,
% grav measured by ascending track - offset of that track =
% grav measured by decending track - offset of that track
G = spalloc( N, M, 2*N );
d = zeros(N,1);
k=0;
for k = (1:N) % loop row-wise
G(k,iA(k)) = 1;
G(k,NA+jD(k)) = -1;
d(k) = Agrav(k)-Dgrav(k);
debg = 0;
if( debg )
fprintf('%d %d %d %f\n', k, iA(j), jD(k), d(k) );
end
end
% solve with damped least squares
% damping is necessary, since the solution
% is only determined up to an additive constant
epsilon = 1.0e-5*max(max(G));
mest = (G'*G + epsilon*speye(M))\(G'*d);
% recreate the dataset of all the data beneath the tracks,
% but correct the shifts
Nsamples=0;
slonv = zeros(M*Ntlon,1);
slatv = zeros(M*Ntlon,1);
sgrav = zeros(M*Ntlon,1);
slon = tlon;
for i = (1:NA) % ascending tracks
slat = Atracklat(:,i);
ii = floor((slat-latmin)/Dlat)+1;
jj = floor((slon-lonmin)/Dlon)+1;
for kk = (1:Ntlon)
if( (ii(kk)>=1) && (ii(kk)<=Nlat) && (jj(kk)>=1) && (jj(kk)<=Nlon) )
Nsamples=Nsamples+1;
slonv(Nsamples) = slon(kk);
slatv(Nsamples) = slat(kk);
sgrav(Nsamples) = GRAV(ii(kk),jj(kk))+mtrue(i)-mest(i);
end
end
end
for j = (1:ND) % decending tracks
slat = Dtracklat(:,j);
ii = floor((slat-latmin)/Dlat)+1;
jj = floor((slon-lonmin)/Dlon)+1;
for kk = (1:Ntlon)
if( (ii(kk)>=1) && (ii(kk)<=Nlat) && (jj(kk)>=1) && (jj(kk)<=Nlon) )
Nsamples=Nsamples+1;
slonv(Nsamples) = slon(kk);
slatv(Nsamples) = slat(kk);
sgrav(Nsamples) = GRAV(ii(kk),jj(kk))+mtrue(NA+j)-mest(NA+j);
end
end
end
% interpolate data along tracks to make 2D image
% Note: same issue as above with duplicate points
% ginterp = TriScatteredInterp( slonv, slatv, sgrav );
sorted = sortrows( [slonv, slatv, sgrav], [1,2] );
Nsort=length(sorted);
flag = ones(Nsort,1);
for i=(2:Nsort)
if( (sorted(i,1)==sorted(i-1,1)) && (sorted(i,2)==sorted(i-1,2)) )
flag(i)=0;
end
end
k = find( flag==1 );
ginterp = TriScatteredInterp( sorted(k,1), sorted(k,2), sorted(k,3) );
pregrav = ginterp(LON(:),LAT(:));
PREGRAV = reshape(pregrav,Nlon,Nlat);
% plot predicted gravity data (with cross-over corrections)
subplot(2,2,4);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [lonmin, lonmax, latmin, latmax] );
axis xy;
imagesc([lonmin, lonmax], [latmin, latmax], PREGRAV);
title('pre gravity');
fprintf('Caption: Maps of (top left) True gravity, (top right) tracks\n');
fprintf('(bottom left) Observed gravity, (bottom right) predicted gravity\n');
% plot gravity data
figure(2);
clf;
colormap('jet');
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [lonmin, lonmax, latmin, latmax] );
axis xy;
imagesc([lonmin, lonmax], [latmin, latmax], PREGRAV);
colorbar();
fprintf('Caption: Predicted gravity\n');
% gdama15_04
% a simple tomography example with a 2D model m(x,y)
% the imaged region is a rectangle, the rays are straight lines
% sources and receivers are colocated on all four edges of the rectangle
clear all;
clear gdaGsparse gdaepsilon;
global gdaGsparse gdaepsilon;
fprintf('gdama15_04\n');
figure();
clf;
% read true model from a file
Sbig=double(imread('../data/magma.tif','tif'));
Sbig = -(Sbig-128);
[Nybig, Nxbig] = size(Sbig);
% decimate the model
idec=4;
Nx=floor(Nxbig/idec);
Ny=floor(Nybig/idec);
S=Sbig(1:idec:Nybig-1, 1:idec:Nxbig-1);
% independent variable, x
xmin = 0;
xmax = Nx;
dx=(xmax-xmin)/(Nx-1);
x = xmin+dx*[0:Nx-1]';
Dx = xmax-xmin;
% independent variable y
ymin = 0;
ymax = Ny;
dy=(ymax-ymin)/(Ny-1);
y = ymin+dy*[0:Ny-1]';
Dy = ymax-ymin;
% rearrage model into vector
% order or model parameters in vector: index = (ix-1)*Ny + iy
mtrue = zeros(Nx*Ny,1);
rowindex = zeros(Nx*Ny,1);
colindex = zeros(Nx*Ny,1);
for ix = [1:Nx]
for iy = [1:Ny]
q = (ix-1)*Ny + iy;
mtrue(q)=S(ix,iy);
rowindex(q)=ix;
colindex(q)=iy;
end
end
% plot true model
subplot(2,2,1);
set(gca, 'YDir', 'reverse' );
set(gca, 'LineWidth', 3 );
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [xmin, xmax, ymin, ymax] );
caxis( [-128, 128] );
image( [xmin, xmax], [ymin, ymax], S, 'CDataMapping', 'scaled');
hold on;
xlabel('y');
ylabel('x');
title( 'true model' );
% build sources and receivers
Nside = 25;
Ns = 4*Nside;
xs = zeros(Ns, 2);
% side 1
xs(1:Nside,1)=xmin+dx/10;
xs(1:Nside,2)=ymin+Dy*[1:Nside]'/(Nside+1);
% side 2
xs(Nside+1:2*Nside,1)=xmax-dx/10;
xs(Nside+1:2*Nside,2)=ymin+Dy*[1:Nside]'/(Nside+1);
% side 3
xs(2*Nside+1:3*Nside,1)=xmin+Dx*[1:Nside]'/(Nside+1);
xs(2*Nside+1:3*Nside,2)=ymin+dy/10;
% side 4
xs(3*Nside+1:4*Nside,1)=xmin+Dx*[1:Nside]'/(Nside+1);
xs(3*Nside+1:4*Nside,2)=ymax-dy/10;
% plot receivers
plot( xs(:,2), xs(:,1), 'k+' );
% rays from all source/receiver pairs
N = Ns*(Ns-1)/2;
ray = zeros(N,4);
k=1;
for i = (1:Ns)
for j = [i+1:Ns]
ray(k,1) = xs(i,1);
ray(k,2) = xs(i,2);
ray(k,3) = xs(j,1);
ray(k,4) = xs(j,2);
k=k+1;
end
end
% plot true model and rays
subplot(2,2,2)
set(gca, 'YDir', 'reverse' );
set(gca, 'LineWidth', 3 );
set(gca,'FontSize',14);
colormap('jet');
axis( [xmin, xmax, ymin, ymax] );
hold on;
caxis( [-128, 128] );
image( [xmin, xmax], [ymin, ymax], S, 'CDataMapping', 'scaled');
xlabel('y');
ylabel('x');
title( 'true model with rays' );
plot( xs(:,2), xs(:,1), 'k+' );
% plot every kRay ray, else image would be black
kRay=20;
for k = (1:kRay:N)
plot( [ray(k,2) ray(k,4)], [ray(k,1), ray(k,3)], 'k-', 'LineWidth', 2 );
end
% build data kernel
% order or model parameters in vector: index = (ix-1)*Ny + iy
M = Nx*Ny;
% I use the row-column-value method to build the sparse
% version of G, but first I set up one row (the kth row)
% as a non-sparse vector Gk
% gdaG=spalloc(N,M,2*N*Nx);
Gr = zeros(2*N*Nx,1);
Gc = zeros(2*N*Nx,1);
Gv = zeros(2*N*Nx,1);
Nel=0; % counts elements
Nr = 2000;
bins=[1:Nx*Ny];
for k = (1:N) % loop over rays, which is to say rows of the dat kernel
x1 = ray(k,1);
y1 = ray(k,2);
x2 = ray(k,3);
y2 = ray(k,4);
r = sqrt( (x2-x1)^2 + (y2-y1)^2 );
dr = r/Nr;
% I use a sloppy way of computing the length of the ray
% in each pixel. I subdivide the ray into Nr pieces, and
% assign each piece to exactly one pixel, the one its
% closest to
Gk = zeros(M,1); % one row of G (the kth row), here non-sparse
if( 0 ) % slow but sure way
for i = (1:Nr)
x = x1 + (x2-x1)*i/Nr;
y = y1 + (y2-y1)*i/Nr;
ix = 1+floor( (Nx-1)*(x-xmin)/Dx );
iy = 1+floor( (Ny-1)*(y-ymin)/Dy );
q = (ix-1)*Ny + iy;
% gdaG(k,q) = gdaG(k,q) + dr;
Gk(q)=Gk(q)+dr;
end
else % faster way, or so we hope
% calculate the array indices of all the ray pieces
xv = x1 + (x2-x1)*[1:Nr]'/Nr;
yv = y1 + (y2-y1)*[1:Nr]'/Nr;
ixv = 1+floor( (Nx-1)*(xv-xmin)/Dx );
iyv = 1+floor( (Ny-1)*(yv-ymin)/Dy );
qv = (ixv-1)*Ny + iyv;
% now count of the ray segments in each pixel of the
% image, and use the count to increment the appropriate
% element of G. The use of the hist() function to do
% the counting is a bit weird, but it seems to work
count=hist(qv,bins);
icount = find( count~=0 );
% gdaG(k,icount) = gdaG(k,icount) + count(icount)*dr;
Gk(icount) = Gk(icount) + count(icount)'*dr;
end
% build row-col-value arrays for this row of G
ii = find(Gk~=0.0)';
for i = ii
Nel=Nel+1;
Gr(Nel) = k;
Gc(Nel) = i;
Gv(Nel) = Gk(i);
end
end % next ray (next row)
% truncate to actual length
Gr = Gr(1:Nel);
Gc = Gc(1:Nel);
Gv = Gv(1:Nel);
% build sparse version of G
gdaGsparse = sparse(Gr,Gc,Gv,N,M);
% delete vectors
clear Gr Gc Gv;
% compute true travel times and plot them
d = gdaGsparse*mtrue;
% plot true travel times on a 2D diagram
% arranged by the rays's closest approach to the origin
% and its orientation
% want to plot traveltimes on NR by Ntheta grid
% independent variable R
NR = 100;
Rmin = 0.0;
Rmax = sqrt((Dx/2)^2+(Dy/2)^2);
dR = (Rmax-Rmin)/(NR-1);
% independent variable theta
Ntheta = 100;
thetamin = -180.0;
thetamax = 180.0;
dtheta = (thetamax-thetamin)/(Ntheta-1);
% build traveltime grid for plotting purposes
TT = zeros(NR,Ntheta);
s = zeros(N,1);
R = zeros(N,1);
theta = zeros(N,1);
for k = (1:N)
x1 = ray(k,1);
y1 = ray(k,2);
x2 = ray(k,3);
y2 = ray(k,4);
% compute closest approach to origin and angle
xd = (x2-x1);
yd = (y2-y1);
s(k) = -(x1*xd+y1*yd)/(xd^2+yd^2);
xs = x1+xd*s(k);
ys = y1+yd*s(k);
R(k) = sqrt( xs^2 + ys^2 );
theta(k) = (180/pi)*atan2(ys,xs);
iR = 1+round((R(k)-Rmin)/dR);
if( iR<1 )
iR=1;
elseif ( iR>NR )
iR=NR;
end
itheta = 1+round((theta(k)-thetamin)/dtheta);
if( itheta<1 )
itheta=1;
elseif ( itheta>Ntheta )
itheta=Ntheta;
end
TT(iR,itheta)=d(k);
% out in data for reverse orientation of ray, for symmetery
theta2 = (180/pi)*atan2(-ys,-xs);
itheta = 1+round((theta2-thetamin)/dtheta);
if( itheta<1 )
itheta=1;
elseif ( itheta>Ntheta )
itheta=Ntheta;
end
TT(iR,itheta)=d(k);
end
% plot TT(R,theta)
subplot(2,2,3);
set(gca, 'LineWidth', 3 );
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [thetamin, thetamax, Rmin, Rmax]);
imagesc( [thetamin, thetamax], [Rmin, Rmax], TT);
set(gca, 'YDir', 'reverse' );
xlabel('theta');
ylabel('R');
title( 'traveltimes' );
% inversion of data by damped least squares
gdaepsilon = 1.0e-2*max(max(abs(gdaGsparse)));
GTd = gda_dlsrhs(d);
mest=bicg(@gda_dlsmul,GTd,1e-5,3*M);
% rearrange m back into image
Sest = zeros(Nx,Ny);
for k = (1:M)
ix=rowindex(k);
iy=colindex(k);
Sest(ix,iy)=mest(k);
end
% plot estimated model
subplot(2,2,4);
set(gca, 'YDir', 'reverse' );
hold on;
set(gca, 'LineWidth', 3 );
set(gca,'FontSize',14);
colormap('jet');
axis( [xmin, xmax, ymin, ymax] );
caxis( [-128, 128] );
image( [xmin, xmax], [ymin, ymax], Sest, 'CDataMapping', 'scaled');
xlabel('y');
ylabel('x');
title( 'estimated model' );
fprintf('Caption: 2D tomography withh full data set (upper left) true model,\n');
fprintf(' (upper right) rays, (lower left) travel times, (lower right) estimated model\n');
% gdama15_05
% a simple tomography example with a 2D model m(x,y)
% the imaged region is a rectangle, the rays are straight lines
% sources and receivers are colocated on three of the four
% edges of the rectangle, a situation which is perhaps analagous
% to seismic tomography problems
clear all;
clear gdaGsparse gdaepsilon;
global gdaGsparse gdaepsilon;
fprintf('gdama15_05\n');
figure();
clf;
% read true model from a file
Sbig=double(imread('../data/magma.tif','tif'));
Sbig = -(Sbig-128);
[Nybig, Nxbig] = size(Sbig);
% decimate the model
idec=4;
Nx=floor(Nxbig/idec);
Ny=floor(Nybig/idec);
S=Sbig(1:idec:Nybig-1, 1:idec:Nxbig-1);
% independent varible x
xmin = 0;
xmax = Nx;
dx=(xmax-xmin)/(Nx-1);
x = xmin+dx*[0:Nx-1]';
Dx = xmax-xmin;
% independent variable y
ymin = 0;
ymax = Ny;
dy=(ymax-ymin)/(Ny-1);
y = ymin+dy*[0:Ny-1]';
Dy = ymax-ymin;
% rearrage model into vector
% order or model parameters in vector: index = (ix-1)*Ny + iy
mtrue = zeros(Nx*Ny,1);
rowindex = zeros(Nx*Ny,1);
colindex = zeros(Nx*Ny,1);
for ix = (1:Nx)
for iy = (1:Ny)
q = (ix-1)*Ny + iy;
mtrue(q)=S(ix,iy);
rowindex(q)=ix;
colindex(q)=iy;
end
end
% plot true model
subplot(2,2,1);
set(gca, 'YDir', 'reverse' );
set(gca, 'LineWidth', 3 );
set(gca, 'FontSize', 14 );
colormap('jet');
hold on;
axis( [xmin, xmax, ymin, ymax] );
caxis( [-128, 128] );
image( [xmin, xmax], [ymin, ymax], S, 'CDataMapping', 'scaled');
hold on;
xlabel('y');
ylabel('x');
title( 'true model' );
% construct sources and receivers on 3 sides of the model
Nside = 25;
Ns = 3*Nside;
xs = zeros(Ns, 2);
% side 1
xs(1:Nside,1)=xmin+dx/10;
xs(1:Nside,2)=ymin+Dy*[1:Nside]'/(Nside+1);
% side 2
xs(1*Nside+1:2*Nside,1)=xmin+Dx*[1:Nside]'/(Nside+1);
xs(1*Nside+1:2*Nside,2)=ymin+dy/10;
% side 3
xs(2*Nside+1:3*Nside,1)=xmin+Dx*[1:Nside]'/(Nside+1);
xs(2*Nside+1:3*Nside,2)=ymax-dy/10;
% plot receivers
plot( xs(:,2), xs(:,1), 'k+' );
% rays connect all source/receiver pairs
N = Ns*(Ns-1)/2;
ray = zeros(N,4);
k=1;
for i = (1:Ns)
for j = (i+1:Ns)
ray(k,1) = xs(i,1);
ray(k,2) = xs(i,2);
ray(k,3) = xs(j,1);
ray(k,4) = xs(j,2);
k=k+1;
end
end
% plot true model and rays
subplot(2,2,2)
set(gca, 'YDir', 'reverse' );
set(gca, 'LineWidth', 3 );
set(gca, 'FontSize', 14 );
colormap('jet');
axis( [xmin, xmax, ymin, ymax] );
hold on;
caxis( [-128, 128] );
image( [xmin, xmax], [ymin, ymax], S, 'CDataMapping', 'scaled');
xlabel('y');
ylabel('x');
title( 'true model with rays' );
plot( xs(:,2), xs(:,1), 'k+' );
% plot every kRay ray, else diagram would be black
kRay=20;
for k = (1:kRay:N)
plot( [ray(k,2) ray(k,4)], [ray(k,1), ray(k,3)], 'k-', 'LineWidth', 2 );
end
% build data kernel
% order or model parameters in vector: index = (ix-1)*Ny + iy
M = Nx*Ny;
% I use the row-column-value method to build the sparse
% version of G, but first I set up one row (the kth row)
% as a non-sparse vector Gk
% gdaG=spalloc(N,M,2*N*Nx);
Gr = zeros(2*N*Nx,1);
Gc = zeros(2*N*Nx,1);
Gv = zeros(2*N*Nx,1);
Nel=0; % counts elements
Nr = 2000;
bins=[1:Nx*Ny];
for k = (1:N) % loop over rays (rows of G)
x1 = ray(k,1);
y1 = ray(k,2);
x2 = ray(k,3);
y2 = ray(k,4);
r = sqrt( (x2-x1)^2 + (y2-y1)^2 );
dr = r/Nr;
% I use a sloppy way of computing the length of the ray
% in each pixel. I subdivide the ray into Nr pieces, and
% assign each piece to exactly one pixel, the one its
% closest to
Gk = zeros(M,1); % one row of G (the kth row), here non-sparse
if( 0 ) % slow but sure way
for i = (1:Nr)
x = x1 + (x2-x1)*i/Nr;
y = y1 + (y2-y1)*i/Nr;
ix = 1+floor( (Nx-1)*(x-xmin)/Dx );
iy = 1+floor( (Ny-1)*(y-ymin)/Dy );
q = (ix-1)*Ny + iy;
% gdaG(k,q) = gdaG(k,q) + dr;
Gk(q)=Gk(q)+dr;
end
else % faster way, or so we hope
% calculate the array indices of all the ray pieces
xv = x1 + (x2-x1)*[1:Nr]'/Nr;
yv = y1 + (y2-y1)*[1:Nr]'/Nr;
ixv = 1+floor( (Nx-1)*(xv-xmin)/Dx );
iyv = 1+floor( (Ny-1)*(yv-ymin)/Dy );
qv = (ixv-1)*Ny + iyv;
% now count of the ray segments in each pixel of the
% image, and use the count to increment the appropriate
% element of G. The use of the hist() function to do
% the counting is a bit weird, but it seems to work
count=hist(qv,bins);
icount = find( count~=0 );
% gdaG(k,icount) = gdaG(k,icount) + count(icount)*dr;
Gk(icount) = Gk(icount) + count(icount)'*dr;
end
% build row-col-value arrays for this row of G
ii = find(Gk~=0.0)';
for i = ii
Nel=Nel+1;
Gr(Nel) = k;
Gc(Nel) = i;
Gv(Nel) = Gk(i);
end
end % next ray (next row)
% truncate to actual length
Gr = Gr(1:Nel);
Gc = Gc(1:Nel);
Gv = Gv(1:Nel);
% build sparse version of G
gdaGsparse = sparse(Gr,Gc,Gv,N,M);
% delete vectors
clear Gr Gc Gv;
% compute true travel times and plot them
d = gdaGsparse*mtrue;
% plot true travel times on a 2D diagram
% arranged by the rays's closest approach to the origin
% and its orientation
% want to plot traveltimes on NR by Ntheta grid
% independent variable R
NR = 100;
Rmin = 0.0;
Rmax = sqrt((Dx/2)^2+(Dy/2)^2);
dR = (Rmax-Rmin)/(NR-1);
% independent variable theta
Ntheta = 100;
thetamin = -180.0;
thetamax = 180.0;
dtheta = (thetamax-thetamin)/(Ntheta-1);
% build traveltime grid for plotting purposes
TT = zeros(NR,Ntheta);
s = zeros(N,1);
R = zeros(N,1);
theta = zeros(N,1);
for k = (1:N)
x1 = ray(k,1);
y1 = ray(k,2);
x2 = ray(k,3);
y2 = ray(k,4);
% compute closest approach to origin and angle
xd = (x2-x1);
yd = (y2-y1);
s(k) = -(x1*xd+y1*yd)/(xd^2+yd^2);
xs = x1+xd*s(k);
ys = y1+yd*s(k);
R(k) = sqrt( xs^2 + ys^2 );
theta(k) = (180/pi)*atan2(ys,xs);
iR = 1+round((R(k)-Rmin)/dR);
if( iR<1 )
iR=1;
elseif ( iR>NR )
iR=NR;
end
itheta = 1+round((theta(k)-thetamin)/dtheta);
if( itheta<1 )
itheta=1;
elseif ( itheta>Ntheta )
itheta=Ntheta;
end
TT(iR,itheta)=d(k);
% out in data for reverse orientation of ray, for symmetery
theta2 = (180/pi)*atan2(-ys,-xs);
itheta = 1+round((theta2-thetamin)/dtheta);
if( itheta<1 )
itheta=1;
elseif ( itheta>Ntheta )
itheta=Ntheta;
end
TT(iR,itheta)=d(k);
end
% plot TT(R,theta)
subplot(2,2,3);
set(gca, 'LineWidth', 3 );
set(gca, 'FontSize', 14 );
colormap('jet');
hold on;
axis( [thetamin, thetamax, Rmin, Rmax]);
imagesc( [thetamin, thetamax], [Rmin, Rmax], TT);
set(gca, 'YDir', 'reverse' );
xlabel('theta');
ylabel('R');
title( 'traveltimes' );
% inversion of data by damped least squares
gdaepsilon = 1.0e-2*max(max(abs(gdaGsparse)));
GTd = gda_dlsrhs(d);
mest=bicg(@gda_dlsmul,GTd,1e-5,3*M);
% rearrange m back into image
Sest = zeros(Nx,Ny);
for k = (1:M)
ix=rowindex(k);
iy=colindex(k);
Sest(ix,iy)=mest(k);
end
% plot estimated image
subplot(2,2,4);
set(gca, 'YDir', 'reverse' );
hold on;
set(gca, 'LineWidth', 3 );
set(gca, 'FontSize', 14 );
colormap('jet');
axis( [xmin, xmax, ymin, ymax] );
caxis( [-128, 128] );
image( [xmin, xmax], [ymin, ymax], Sest, 'CDataMapping', 'scaled');
xlabel('y');
ylabel('x');
title( 'estimated model' );
fprintf('Caption: 2D tomography withh full data set (upper left) true model,\n');
fprintf(' (upper right) rays, (lower left) travel times, (lower right) estimated model\n');
% gdama15_06
% temperature inversion problem
%
% The data are a temperature profile made at
% a time tobs>=0. The model parameters are the
% initial temperature profile (at time t0=0).
% The problem is solved many times, for
% different tobs's, so that one can examine
% how the ability to reconstruct the initial
% temperature degrades as the data are measured
% later and later in time.
% Thus, time in the plot is the time of observation,
% tobs. In all cases, the quantity being determined is
% the inital temperature distributon at the inital time,
% t0=0.
clear all;
fprintf('gdama15_06\n');
% independent spartial variable x
Nx = 100;
xmin = -100;
xmax = 100;
Dx = (xmax-xmin)/(Nx-1);
x = xmin + Dx*[0:Nx-1]';
% independent temporal variable t
Nt = 100;
tmin = 0;
tmax = 200;
Dt = (tmax-tmin)/(Nt-1);
t = tmin + Dt*[0:Nt-1]';
% time/space matrices
tt = t*ones(1,Nx);
xx = ones(Nt,1)*x';
% model parameters, one at each distance
s = 1.0;
mtrue = zeros(Nx,1);
L=5;
mtrue(4*Nx/10:6*Nx/10)=(2+cos( 2*pi*([4*Nx/10:6*Nx/10]'-Nx/2)/L ));
% the true data as a function of observation time and distance
% are constructed by evaluating an analytic solution
% to a heat flow equation
dtxtrue = zeros(Nt,Nx);
for i=(1:Nx)
% data kernel, Gij, gives temperature at xi due to source
% at xj for a time, t0
t0=t(i);
if( t0==0 )
G=eye(Nx,Nx);
else
erfA = erf( (x*ones(1,Nx)-ones(Nx,1)*(x-Dx/2)') / sqrt(t0) );
erfB = erf( (x*ones(1,Nx)-ones(Nx,1)*(x+Dx/2)') / sqrt(t0) );
G = 0.5*(erfA-erfB);
end
dtxtrue(i,:) = (G*mtrue)';
end
% plot the true data
figure();
clf;
colormap('jet');
set( gca, 'LineWidth', 3);
set( gca, 'FontSize', 14);
hold on;
axis( [xmin, xmax, tmin, tmax] );
axis ij,
imagesc( [xmin, xmax], [tmin, tmax], dtxtrue);
ylabel('observation time');
xlabel('distance');
title('true temperature');
colorbar();
fprintf('Caption: True temperature as as function of time\n');
% true model, replicated at all times
mtxtrue = ones(Nt,1)*mtrue';
% plot the true model
figure();
clf;
subplot(1,3,1);
set( gca, 'LineWidth', 3);
set( gca, 'FontSize', 14);
hold on;
axis( [xmin, xmax, tmin, tmax] );
axis ij,
c1 = min(min(mtxtrue));
c2 = max(max(mtxtrue));
caxis( [c1, c2] );
image( [xmin, xmax], [tmin, tmax], mtxtrue, 'CDataMapping', 'scaled');
ylabel('observation time');
xlabel('distance');
title('true model');
R1at = 10;
R2at = 40;
% damped minimum length case
% estimate the model using the data each time
mtxest = zeros(Nt,Nx);
for i=(1:Nt)
% true data
dtrue=dtxtrue(i,:)';
% add noise
sd = 0.01;
dobs = dtrue + random('Normal',0,sd,Nx,1);
% compute data kernel
t0=t(i);
if( t0==0 )
G=eye(Nx,Nx);
else
erfA = erf( (x*ones(1,Nx)-ones(Nx,1)*(x-Dx/2)') / sqrt(t0) );
erfB = erf( (x*ones(1,Nx)-ones(Nx,1)*(x+Dx/2)') / sqrt(t0) );
G = 0.5*(erfA-erfB);
end
% damped minimum length solution
epsilon=1.0e-3;
mtxest(i,:) = (G'*((G*G'+epsilon*eye(Nx,Nx))\dobs))';
% save resolving kernels
if( i==R1at )
R1ML = G'*inv((G*G'+epsilon*eye(Nx,Nx)))*G;
elseif (i==R2at)
R2ML = G'*inv((G*G'+epsilon*eye(Nx,Nx)))*G;
end
end
% plot ML estimated solution
subplot(1,3,2);
set( gca, 'LineWidth', 3);
set( gca, 'FontSize', 14);
hold on;
axis( [xmin, xmax, tmin, tmax] );
axis ij,
caxis( [c1, c2] );
image( [xmin, xmax], [tmin, tmax], mtxest, 'CDataMapping', 'scaled');
ylabel('observation time');
xlabel('distance');
title('ML m^{est}');
% Backus-Gilbert case
% estimate the model using the data each time
mtxest = zeros(Nt,Nx);
for i=(1:Nt)
% true data
dtrue=dtxtrue(i,:)';
% add noise
sd = 0.01;
dobs = dtrue + random('Normal',0,sd,Nx,1);
% compute data kernel
t0=t(i);
if( t0==0 )
G=eye(Nx,Nx);
else
erfA = erf( (x*ones(1,Nx)-ones(Nx,1)*(x-Dx/2)') / sqrt(t0) );
erfB = erf( (x*ones(1,Nx)-ones(Nx,1)*(x+Dx/2)') / sqrt(t0) );
G = 0.5*(erfA-erfB);
end
% construct Backus-Gilbert generalized inverse
GMG = zeros(Nx,Nx);
u = G*ones(Nx,1);
for k = (1:Nx)
S = G * diag(([1:Nx]-k).^2) * G';
epsilon=1e-3;
S = S+epsilon*eye(Nx);
uSinv = u'/S;
GMG(k,:) = uSinv / (uSinv*u);
end
% estimate solution
mtxest(i,:) = (GMG*dobs)';
% save resolving kernels
if( i==R1at )
R1BG = GMG*G;
elseif (i==R2at)
R2BG = GMG*G;
end
end
% plot BG estimated solution
subplot(1,3,3);
set( gca, 'LineWidth', 3);
set( gca, 'FontSize', 14);
hold on;
axis( [xmin, xmax, tmin, tmax] );
axis ij,
caxis( [c1, c2] );
image( [xmin, xmax], [tmin, tmax], mtxest, 'CDataMapping', 'scaled');
ylabel('observation time');
xlabel('distance');
title('BG m^{est}');
fprintf('Caption: Comparison of results. (left) True temperature, (middle) Minumum Length (ML) estimate\n');
fprintf('(right) Backus-Gilbert (BG) estimate\n');
% plot resolution kernels
figure(3);
clf;
% plot ML resolution at time t1
subplot(2,2,1);
set( gca, 'LineWidth', 3);
set( gca, 'FontSize', 14);
hold on;
axis( [xmin, xmax, xmin, xmax] );
axis ij,
imagesc( [xmin, xmax], [xmin, xmax], R1ML);
ylabel('distance');
xlabel('distance');
title(sprintf('ML R at t=%.2f',t(R1at)));
% plot ML resolution at time t1
subplot(2,2,2);
set( gca, 'LineWidth', 3);
set( gca, 'FontSize', 14);
hold on;
axis( [xmin, xmax, xmin, xmax] );
axis ij,
imagesc( [xmin, xmax], [xmin, xmax], R2ML);
ylabel('distance');
xlabel('distance');
title(sprintf('ML R at t=%.2f',t(R2at)));
% plot BG resolution at time t1
subplot(2,2,3);
set( gca, 'LineWidth', 3);
set( gca, 'FontSize', 14);
hold on;
axis( [xmin, xmax, xmin, xmax] );
axis ij,
imagesc( [xmin, xmax], [xmin, xmax], R1BG);
ylabel('distance');
xlabel('distance');
title(sprintf('BG R at t=%.2f',t(R1at)));
% plot BG resolution at time t2
subplot(2,2,4);
set( gca, 'LineWidth', 3);
set( gca, 'FontSize', 14);
hold on;
axis( [xmin, xmax, xmin, xmax] );
axis ij,
imagesc( [xmin, xmax], [xmin, xmax], R2BG);
ylabel('distance');
xlabel('distance');
title(sprintf('BG R at t=%.2f',t(R2at)));
fprintf('Caption Minimum Length (ML) and Backus-Gilbert (BG) resolution\n');
fprintf('matrices using data collected at two different time\n');
% gdama15_07
%
% L1, L2 and Linf norm; fitting of a straight line
% L1 and Linf are via transformation to a linear
% programming problem; Ls is via least squares.
clear all;
fprintf('gdama15_07\n');
% auxillary variable, z
N=10;
zmin = 0;
zmax = 1;
Dz = (zmax-zmin)/(N-1);
z = zmin + Dz*[0:N-1]';
% set up for linear fit
M=2;
mtrue = [1, 3]';
G = [ones(N,1), z];
dtrue = G*mtrue;
% syntehtic data, with noise drawn from a two-sided exponential
% distribution with variance sd^2.
sd = 0.4 * ones(N,1);
dobs=zeros(N,1);
for i=(1:N)
r=random('exp',sd(i)/sqrt(2)).*(2*(random('unid',2)-1.5));
dobs(i) = dtrue(i) + r;
end
% Part 1: L2 Norm
% weighted least squares solution (sure the easiest!)
WI = diag(sd.^-2);
mls = (G'*WI*G)\(G'*WI*dobs);
dls = G*mls;
% Part 2: L1 norm
% set up for linear programming problem
% min f*x subject to A x <= b, Aeq x = beq
% variables
% m = mp - mpp
% x = [mp', mpp', alpha', x', xp']
% with mp, mpp length M and alpha, x, xp, length N
L = 2*M+3*N;
x = zeros(L,1);
% f is length L
% minimize sum aplha(i)/sd(i)
f = zeros(L,1);
f(2*M+1:2*M+N)=1./sd;
% equality constraints
Aeq = zeros(2*N,L);
beq = zeros(2*N,1);
% first equation G(mp-mpp)+x-alpha=d
Aeq(1:N,1:M) = G;
Aeq(1:N,M+1:2*M) = -G;
Aeq(1:N,2*M+1:2*M+N) = -eye(N,N);
Aeq(1:N,2*M+N+1:2*M+2*N) = eye(N,N);
beq(1:N) = dobs;
% second equation G(mp-mpp)-xp+alpha=d
Aeq(N+1:2*N,1:M) = G;
Aeq(N+1:2*N,M+1:2*M) = -G;
Aeq(N+1:2*N,2*M+1:2*M+N) = eye(N,N);
Aeq(N+1:2*N,2*M+2*N+1:2*M+3*N) = -eye(N,N);
beq(N+1:2*N) = dobs;
% inequality constraints A x <= b
% part 1: everything positive
A = zeros(L+2*M,L);
b = zeros(L+2*M,1);
A(1:L,:) = -eye(L,L);
b(1:L) = zeros(L,1);
% part 2; mp and mpp have an upper bound. Note:
% Without this constraint, a potential problem is that
% mp and mpp are individually unbounded, though their
% difference, m=mp-mpp, is not. Thus there might be cases
% where the algroithm strays to very large mp and mpp.
A(L+1:L+2*M,:) = eye(2*M,L);
mupperbound=10*max(abs(mls)); % might need to be adjusted for problem at hand
b(L+1:L+2*M) = mupperbound;
% solve linear programming problem
[x, fmin] = linprog(f,A,b,Aeq,beq);
fmin=-fmin;
mestL1 = x(1:M) - x(M+1:2*M);
dpreL1 = G*mestL1;
% Part 3: Linf Norm
% set up for linear programming problem
% min f*x subject to A x <= b, Aeq x = beq
% variables
% m = mp - mpp
% x = [mp', mpp', alpha', x', xp']
% with mp, mpp length M; alpha length 1, x, xp, length N
L = 2*M+1+2*N;
x = zeros(L,1);
% f is length L
% minimize alpha
f = zeros(L,1);
f(2*M+1:2*M+1)=1;
% equality constraints
Aeq = zeros(2*N,L);
beq = zeros(2*N,1);
% first equation G(mp-mpp)+x-alpha=d
Aeq(1:N,1:M) = G;
Aeq(1:N,M+1:2*M) = -G;
Aeq(1:N,2*M+1:2*M+1) = -1./sd;
Aeq(1:N,2*M+1+1:2*M+1+N) = eye(N,N);
beq(1:N) = dobs;
% second equation G(mp-mpp)-xp+alpha=d
Aeq(N+1:2*N,1:M) = G;
Aeq(N+1:2*N,M+1:2*M) = -G;
Aeq(N+1:2*N,2*M+1:2*M+1) = 1./sd;
Aeq(N+1:2*N,2*M+1+N+1:2*M+1+2*N) = -eye(N,N);
beq(N+1:2*N) = dobs;
% inequality constraints A x <= b
% part 1: everything positive
A = zeros(L+2*M,L);
b = zeros(L+2*M,1);
A(1:L,:) = -eye(L,L);
b(1:L) = zeros(L,1);
% part 2; mp and mpp have an upper bound. Note:
% A potential problem without this constraint is that
% mp and mpp are individually unbounded, though their
% difference, m=mp-mpp, is not. Thus there might be cases
% where the algroithm strays to very large mp and mpp.
A(L+1:L+2*M,:) = eye(2*M,L);
mupperbound=10*max(abs(mls)); % might need to be adjusted for problem at hand
b(L+1:L+2*M) = mupperbound;
% solve linear programming problem
[x, fmin] = linprog(f,A,b,Aeq,beq);
fmin=-fmin;
mestLinf = x(1:M) - x(M+1:2*M);
dpreLinf = G*mestLinf;
% plot the observed and presicted data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [zmin, zmax, 0, 5 ] );
plot( z, dtrue, 'k-', 'Linewidth', 3);
plot( z, dpreL1, 'g-', 'Linewidth', 3);
plot( z, dls, 'b-', 'Linewidth', 3);
plot( z, dpreLinf, 'r-', 'Linewidth', 3);
plot( z, dobs, 'ko', 'Linewidth', 3);
xlabel('z');
ylabel('d');
fprintf('Caption: Linear fits to data\n');
fprintf('black: d-true\n');
fprintf('circles: d-obs\n');
fprintf('green: d-L1\n');
fprintf('blue: d-L2\n');
fprintf('red: d-Linf\n');
fprintf(' \n');
fprintf('Solution:');
fprintf(' m1 m2\n');
fprintf('true %.4f %.4f\n', mtrue(1), mtrue(2));
fprintf('L1 %.4f %.4f\n', mestL1(1), mestL1(2));
fprintf('L2 %.4f %.4f\n', mls(1), mls(2));
fprintf('Linf %.4f %.4f\n', mestLinf(1), mestLinf(2));
% gdama15_08
% Depiction of a group of vectors in 3D space
% scattering about a central vector.
clear all;
fprintf('gdama15_08\n');
% grid
xmin = 0;
xmax = 1;
ymin = 0;
ymax = 1;
zmin = 0;
zmax = 1;
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [xmin, xmax, ymin, ymax, zmin, zmax]');
% improvise a 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 2 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 2 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 2 );
xlabel('x_1');
ylabel('x_2');
zlabel('x_3');
rbar=[1, 1, 1]'; % the central vector
rlen = sqrt( rbar'*rbar );
rbar = rbar/rlen;
% Note: this is not a Fisher distribution, but it is azimuthually uniform
for i=(1:20)
rvec = random('Uniform', -1, 1, 3, 1); % random vector
b = cross(rvec,rbar); % vector perpendicular to rvec
rnew = rbar + 0.15*b; % add small deviation to rvec
rnew = rnew/sqrt(rnew'*rnew); % normalize to unit length
gda_arrow3( rnew, 'k-', 2 );
end
gda_arrow3( rbar, 'r-', 4 )
daspect([1 1 1])
view(3)
camlight; lighting phong
fprintf('Caption: Observed vectors in a 3D space(black) and their mean (red)\n');
% gdama15_09
% depicts vector relationships in 3D space
clear all;
fprintf('gdama15_09\n');
% independent variable x
Nx = 51;
xmin = -1;
xmax = 1;
Dx = (xmax-xmin)/(Nx-1);
x = xmin + Dx*[0:Nx-1]';
% independent variable y
Ny = 51;
ymin = -1;
ymax = 1;
Dy = (ymax-ymin)/(Ny-1);
y = ymin + Dy*[0:Ny-1]';
% independent variable z
Nz = 51;
zmin = -1;
zmax = 1;
Dz = (zmax-zmin)/(Nz-1);
z = zmin + Dz*[0:Nz-1]';
figure();
clf;
set(gca,'LineWidth',2);
hold on;
axis( [xmin, xmax, ymin, ymax, zmin, zmax]');
axis off;
% central direction
q1=45;
f1=45;
b1=1;
a1= [b1*cos(pi*f1/180)*sin(pi*q1/180), b1*sin(pi*f1/180)*sin(pi*q1/180), b1*cos(pi*q1/180)]';
gda_arrow3( a1, 'k-', 3);
% another vector
q2=65;
f2=35;
b2=1;
a2= [b2*cos(pi*f2/180)*sin(pi*q2/180), b2*sin(pi*f2/180)*sin(pi*q2/180), b2*cos(pi*q2/180)]';
gda_arrow3( a2, 'b-', 3);
% co-latitudinal arc from central vector to the other vector
lena1 = sqrt(a1'*a1);
lena2 = sqrt(a2'*a2);
ua1 = a1/lena1;
ua2 = a2/lena2;
th=180*acos((ua1'*ua2))/pi;
Nth = floor(th)+1;
alpha = [0:Nth-1]'/(Nth-1);
arc1=zeros(Nth,3);
b3 = 0.85;
arc1 = alpha*ua1' + (1-alpha)*ua2';
lenarc1 = sqrt( arc1(:,1).^2 + arc1(:,2).^2 + arc1(:,3).^2 );
arc1 = b3 * arc1 ./ (lenarc1*ones(1,3));
plot3( arc1(:,1), arc1(:,2), arc1(:,3), 'k-', 'LineWidth', 2 );
% longitudinal arc
a3 = cross(a1/lena1,a2/lena2);
lena3 = sqrt(a3'*a3);
ua3 = a3/lena3;
Nth=20;
alpha = [-Nth:Nth-1]'/90;
arc2=zeros(Nth,3);
arc2 = alpha*ua3' + (1-alpha)*ua2';
lenarc2 = sqrt( arc2(:,1).^2 + arc2(:,2).^2 + arc2(:,3).^2 );
arc2 = b3 * arc2 ./ (lenarc2*ones(1,3));
plot3( arc2(:,1), arc2(:,2), arc2(:,3), 'k-', 'LineWidth', 2 );
% depict sphere in 3D space by contouring a Normal distribution
% 3D grid
[XX, YY, ZZ] = meshgrid( x, y, z );
% parameters for Normal distribution
rbar = [0, 0, 0]';
sd=1;
C = (sd^2)*eye(3,3);
CI = inv(C);
DC = det(C);
norm = ( ((2*pi)^(3/2)) * sqrt(DC) );
% tabulate Normal distribution
PP = zeros(Nx,Ny,Nz);
for i=(1:Nx)
for j=(1:Ny)
for k=(1:Nz)
r = [XX(i,j,k), YY(i,j,k), ZZ(i,j,k)]';
PP(i,j,k) = exp(-0.5*(r-rbar)'*CI*(r-rbar))/norm;
end
end
end
% contour Normal distribution in 3D space
maxP=max(max(max(PP)));
p=patch(isosurface( XX, YY, ZZ, PP, 0.7*maxP ));
isonormals(XX,YY,ZZ,PP, p)
set(p, 'FaceColor', 'red', 'FaceAlpha', 0.5, 'EdgeColor', 'none');
% set view angle, etc
daspect([1 1 1])
view(87,6)
camlight; lighting phong
fprintf('Caption: depiction of vector (blue) that has a deviation (dotted) from a central angle (black),\n');
fprintf('A longitudal arc (curve) is also shown.')
% gdama15_10
% plots Fisher distribution for two different
% values of the precision parameter
clear all;
fprintf('gdama15_10\n');
% independent variable x
Nx = 101;
xmin = -pi;
xmax = pi;
Dx = (xmax-xmin)/(Nx-1);
x = xmin + Dx*[0:Nx-1]';
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [-3.5, 3.5, 0, 1] );
xlabel('theta');
ylabel('p(theta)');
% precision parameters
k1=1;
k2=5;
% Fisher p.d.f.'s
p1 = (k1/(4*pi*sinh(k1)))*exp(k1*cos(x));
p2 = (k2/(4*pi*sinh(k2)))*exp(k2*cos(x));
% plot from -pi to pi so they look symmetrical,
% like a Normal distribution
plot( x, p1, 'r-', 'LineWidth', 3);
plot( x, p2, 'b-', 'LineWidth', 3);
plot( -pi, 0, 'k+', 'LineWidth', 2);
plot( pi, 0, 'k+', 'LineWidth', 2);
fprintf('Caption: Two Fisher pdfs on a sphere, with large (blue) and small (red) precision parameters\n');
% gdama15_11
% example of estimating central vector of Fisher distribution
% using principle of Maximum Liklihood. The data are P-axes of
% moment tensors of deep earthquakes on the Kurile-Kamchatka
% subduction zone. The boostrap method is used to compute
% confidence limits of the aximuth and plunge of the central vectorgdaFsparse
clear all;
fprintf('gdama15_11\n');
% plot set-up
figure();
clf;
set( gca, 'LineWidth', 1 );
set( gca, 'Xtick', [] );
set( gca, 'Ytick', [] );
axis( [-1.1 1.1 -1.1 1.1] );
hold on;
% draw stereonet
% steronet angles:
% theta, angle E of N, so theta=0 is north
% phi, vertical angle, 0 rhi=0 when down
theta = (2*pi/360)*[0:360];
plot( sin(theta), cos(theta), 'k', 'LineWidth', 3 );
theta = (pi/180.0)*0.0;
text( 1.05*sin(theta), 1.05*cos(theta), 'N' );
theta = (pi/180.0)*90.0;
text( 1.05*sin(theta), 1.05*cos(theta), 'E' );
plotgrid=1;
if( plotgrid == 1 )
phi = (pi/180.0)*[30, 60];
theta = (2*pi/360)*[0:360];
for p = phi
plot( gda_stereo(p)*sin(theta), gda_stereo(p)*cos(theta), 'k--', 'LineWidth', 2 );
end
theta = (pi/180.0)*linspace(0,330,12);
r = linspace(0,1,2);
for t = theta
plot( r*sin(t), r*cos(t), 'k--', 'LineWidth', 2 );
end
end
% load data from file
X=load('../data/moments.txt');
% file: lon lat depth mrr mtt mpp mrt mrp mtp iexp
% 1 2 3 4 5 6 7 8 9 10
% moment tensor elements as read from file:
RR=X(:,4); TT=X(:,5); PP=X(:,6); RT=X(:,7); RP=X(:,8); TP=X(:,9);
% restrict calculation to data in a specific depth range
dmin=300; dmax=600; % depth range
j=find( (X(:,3)>=dmin) & (X(:,3)<=dmax) );
N = length(j);
Pgs=zeros(N,3);
% plot all vectors
k=1;
for i = j' % loop over earthquakes
% I'm using an [up, east, north] coordinate system
% the CMT is in [R=up, T=south, P=east] coordinate system
M = zeros(3,3);
M(1,1)=RR(i); M(2,2)=PP(i); M(3,3)=TT(i);
M(1,2)=RP(i); M(2,1)=RP(i);
M(1,3)=-RT(i); M(3,1)=-RT(i);
M(2,3)=-TP(i); M(3,2)=-TP(i);
% extract eigenvectors
[V,D] = eig(M);
Pg = V(:,1); % P-axis in [up, east, north] coordinate system, is third eigenvector
if( Pg(1) > 0 ) % want axis to point down
Pg=-Pg;
end
% save in array
Pgs(k,:) = Pg';
k = k+1;
% angle east of north, north=0
theta = atan2( Pg(2), Pg(3) );
% vertical angle from down, down=0
phi = atan( sqrt(Pg(2)^2 + Pg(3)^2)/(-Pg(1)));
r = gda_stereo(phi);
plot( r*sin(theta), r*cos(theta), 'ko', 'LineWidth', 4 );
end
% estimate central vector
Pgsum = sum(Pgs)';
norm = sqrt(Pgsum'*Pgsum);
Pgsum = Pgsum/norm;
theta = atan2( Pgsum(2), Pgsum(3) );
phi = atan( sqrt(Pgsum(2)^2 + Pgsum(3)^2)/(-Pgsum(1)));
r = gda_stereo(phi);
plot( r*sin(theta), r*cos(theta), 'ro', 'LineWidth', 4 );
% bootstrap estimates saved in these arrays
Nboot=10000;
thetaboot=zeros(N,1);
phiboot=zeros(N,1);
% bootstrap calculation: resample data with duplications,
% and compute central vector estimate of each
for i=(1:N)
% resample
Pboot = Pgs(unidrnd(N,N,1),:);
% central vector calculation
Pbootsum = sum(Pboot)';
norm = sqrt(Pbootsum'*Pbootsum);
Pbootsum = Pbootsum/norm;
% convert to angles
thetaboot(i) = atan2( Pbootsum(2), Pbootsum(3) );
phiboot(i) = atan( sqrt(Pbootsum(2)^2 + Pbootsum(3)^2)/(-Pbootsum(1)));
% save, plot as cyan dot
r = gda_stereo(phiboot(i));
plot( r*sin(thetaboot(i)), r*cos(thetaboot(i)), 'c.', 'LineWidth', 2 );
end
% replot estimate, because dots from bootstrap obscure its symbol
r = gda_stereo(phi);
plot( r*sin(theta), r*cos(theta), 'ro', 'LineWidth', 4 );
fprintf('Caption: Earthquake P-axes (cirlces) for a subduction zone displayed on a\n');
fprintf('stereo net, with the Fisher estimate of the mean angle (red) and bootstrap\n');
fprintf('confidence region (blue)\n');
fprintf('\n');
% compute standard deviation as a quick estimate of confidence intervals
sigmatheta = std(thetaboot);
sigmaphi = std(phiboot);
% display results
fprintf('Maximum Liklihood Estimate\n');
fprintf('azimuthal angle (North=0), theta = %.2f +/- %.2f (2 sigma)\n', 180*theta/pi, 2*180*sigmatheta/pi);
fprintf('vertical angle (Down=0), phi = %.2f +/- %.2f (2 sigma)\n', 180*phi/pi, 2*180*sigmaphi/pi);
% gdama15_12
%
% model Mars rover Mosbauer spectra using both
% Gaussian and Lorentzian curves, fit via nonlinear
% least squares (Newton's Method).
%
% Note: The code requires you to click on the bottom of the
% peaks, and when you're done, to click to the left of the
% x-axis. This process defines the number of peaks and their
% starting positions and amplitudes.
clear all;
fprintf('gdama15_12\n');
% load data
D=load('../data/mars_soil.txt');
v=D(:,1);
d=D(:,2);
d=d/max(d); % normalize
N=length(d);
% delete negative velocities
i=find(v>=0,1);
v=v(i:N);
d=d(i:N);
N=length(v);
% plot data
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 12, min(d), max(d)] );
plot(v,d,'r-','LineWidth',2);
plot(v,d,'ro','LineWidth',3);
xlabel('velocity, mm/s');
ylabel('counts');
dolorentzian=0; % toggles between lorentzian curves and Normal curves
% lorentzian curve of peak amplitude a, center velocity v0 and width c
% f(v) = a c^2 / ( (v-v0)^2 + c^2) )
% df/da = c^2 / ( (v-v0)^2 + c^2) )
% df/dv0 = 2 a c^2 (v-v0) / ( (v-v0)^2 + c^2) )^2
% df/dc = 2 a c / ( (v-v0)^2 + c^2) ) - 2 a c^3 / ( (v-v0)^2 + c^2) )^2
% 3 model parameters per lorentzian, (a, v0, c)
% Normal curve of peak amplitude a, center velocity v0 and width c
% f(v) = ( a / ( sqrt(2 pi) c )) exp( -0.5 * (v-v0)^2 c^-2 )
% df/da = ( 1 / ( sqrt(2 pi) c )) exp( -0.5 * (v-v0)^2 c^-2 )
% df/dv0 = ( a / ( sqrt(2 pi) c )) ((v - v0 )/(c^2)) exp( -0.5 * (v-v0)^2 c^-2 )
% df/dc = ( a / ( sqrt(2 pi) c^2 )) (((v - v0)^2/(c^2))-1) exp( -0.5 * (v-v0)^2 c^-2 )
% 3 model parameters per Normal curve, (a, v0, c)
% estimate of background level
A = max(d);
rtp = sqrt(2*pi);
% number of peaks determined by clicking on graph
disp('click on the bottom each peak');
disp(' click to the left of zero when done');
disp(' ');
% input peaks
MAXPEAKS=100;
a = zeros(MAXPEAKS,1);
v0 = zeros(MAXPEAKS,1);
c = zeros(MAXPEAKS,1);
K=0;
for k = (1:20)
p = ginput(1);
if( p(1) < 0 )
break;
end
K=K+1;
a(K) = p(2)-A;
v0(K)=p(1);
c(K)=0.1;
end
% truncate vectors
a = a(1:K);
v0 = v0(1:K);
c = c(1:K);
% model parameters
M=K*3;
m = [a', v0', c']';
for iter=(1:10)
% predicted data
dpre = A*ones(N,1);
for i = [1:K]
temp = ((v-m(K+i)).^2+m(2*K+i)^2);
dpre = dpre + m(i)*(m(2*K+i)^2)./temp;
end
% data kernel
G=zeros(N,M);
for i = [1:K]
temp = ((v-m(K+i)).^2+m(2*K+i)^2);
G(:,i) = (m(2*K+i)^2)./temp; % d/da
G(:,K+i) = 2*m(i)*(m(2*K+i)^2)*(v-m(K+i))./(temp.^2); % d/dv0
G(:,2*K+i) = 2*m(i)*m(2*K+i)./temp - 2*m(i)*m(2*K+i)^3./(temp.^2); % d/dc
end
% deviations in data
dd = d - dpre;
E = dd'*dd;
% updated model
epsilon=0.001;
dm = (G'*G+epsilon*eye(M))\(G'*dd);
mold=m;
m = m+dm;
% output parameters, mostly for debugging purposes
verbose=0;
if( verbose )
fprintf('iteration %d A v0 c\n', iter);
for i=[1:K]
fprintf('%d %f %f %f\n', i, m(i), m(K+i), m(2*K+i));
end
fprintf('total error %f\n', E);
fprintf('\n');
end
end
% plot
dpre2 = A*ones(N,1);
for i = (1:K)
temp = ((v-m(K+i)).^2+m(2*K+i)^2);
dpre2 = dpre2 + m(i)*(m(2*K+i)^2)./temp;
end
dd = d - dpre2;
E = dd'*dd;
mlorentzian=m;
Elorentzian=E;
nulorentzian = (N-M);
plot(v,dpre2,'b-','LineWidth',3);
% output lorentzian results for debugginf prurposes
fprintf('Lorentzian parameters\n', iter);
for i=(1:K)
fprintf('%d %f %f %f\n', i, m(i), m(K+i), m(2*K+i));
end
fprintf('total error %f\n', E);
fprintf('\n');
% Normal Curve
% model parameters
M=K*3;
m = [a'.*(rtp.*c'), v0', c']';
for iter=(1:10)
% predicted data
dpre = A*ones(N,1);
for i = (1:K)
temp = exp(-0.5*(v-m(K+i)).^2/(m(2*K+i)^2));
dpre = dpre + (m(i)/(rtp*m(2*K+i)))*temp;
end
% data kernel
G=zeros(N,M);
for i = (1:K)
temp = exp(-0.5*((v-m(K+i)).^2)/(m(2*K+i)^2));
G(:,i) = (1/(rtp*m(K*2+i)))*temp; % d/da
G(:,K+i) = (m(i)/(rtp*m(2*K+i))) .* ((v-m(K+i))/(m(2*K+i)^2)) .* temp; % d/dv0
G(:,2*K+i) = (m(i)/(rtp*m(2*K+i)^2)) .* (((v-m(K+i)).^2/(m(2*K+i)^2))-1) .* temp; % d/dc
end
% deviations in data
dd = d - dpre;
E = dd'*dd;
% updated model via damped least squares
epsilon=0.001;
dm = (G'*G+epsilon*eye(M))\(G'*dd);
mold=m;
m = m+dm;
% output parameters, mostly for debugging purposes
verbose=0;
if( verbose )
fprintf('iteration %d A v0 c\n', iter);
for i=(1:K)
fprintf('%d %f %f %f\n', i, m(i), m(K+i), m(2*K+i));
end
fprintf('total error %f\n', E);
fprintf('\n');
end
end % end iterations
% plot results
dpre2 = A*ones(N,1);
for i = (1:K)
temp = exp(-0.5*(v-m(K+i)).^2/(m(2*K+i)^2));
dpre2 = dpre2 + m(i)/(rtp*m(2*K+i))*temp;
end
dd = d - dpre2;
E = dd'*dd;
mnormal=m;
Enormal=E;
nunormal = (N-M);
plot(v,dpre2,'g-','LineWidth',3);
fprintf('Caption: spectral data (red circles) with Lorentzian (blue) and Gaussian (green) fit\n');
fprintf('\n');
% output normal results
fprintf('Normal parameters\n', iter);
for i=(1:K)
fprintf('%d %f %f %f\n', i, m(i), m(K+i), m(2*K+i));
end
fprintf('total error %f\n', E);
fprintf('\n');
F = (Enormal/nunormal)/(Elorentzian/nulorentzian);
fprintf('Fest = E_normal/E_lorentzian: %f\n', F);
if( F<1 )
F=1/F;
end
p = 1-(fcdf(F,nunormal,nulorentzian)-fcdf(1/F,nunormal,nulorentzian));
fprintf('P(F<=1/Fest||F>=Fest) = %f\n', p);
% gdama15_13
% tidal periods using Fast Fourier Transform
% to compute amplitude spectrum of sea height time series
clear all;
fprintf('gdama15_13\n');
% load data
D = load('../data/tides.txt');
d = D(:,2); % sea surface in feet sampple every hour
N = length(d);
N = 2*floor(N/2); % truncate to an even number of points
d = d(1:N);
% t-axis in days
Dt = 1.0/24.0;
t = Dt*[0:N-1]';
tmin=0.0;
tmax=Dt*(N-1);
figure();
clf;
hold on;
axis( [0, 365, -10, 10]);
set(gca,'LineWidth',2);
plot(t,d,'k-','LineWidth',2);
xlabel('time (days)');
ylabel('sea height (ft)');
fprintf('Caption: Sea surface height time series (compete)\n');
figure();
clf;
hold on;
axis( [150, 210, -5, 5]);
set(gca,'LineWidth',2);
plot(t,d-mean(d),'k-','LineWidth',2);
xlabel('time (days)');
ylabel('sea height (ft)');
fprintf('Caption: Sea surface height time series (short)\n');
% frequency axis, cycles per day
fmax = 1.0/(2.0*Dt);
Df = fmax/(N/2);
M = N/2+1; % non-negative frequencies
f = Df*[0:M-1]';
p(2:M) = 1./f(2:M); % period
p(1)=0.0;
% hamming taper, a bell-shaped curve
a0 = 0.53836;
a1 = 0.46164;
ham = a0 - a1 * cos( 2*pi*[0:N-1]'/(N-1));
figure();
clf;
hold on;
axis( [0, 365, 0, 1]);
set(gca,'LineWidth',2);
plot(t,ham,'k-','LineWidth',2);
xlabel('time (days)');
ylabel('Hamming taper');
fprintf('Caption: Hamming taper\n');
% amplitude spectrim
dp = ham.*(d-mean(d)); % remove mean and hamming taper
m = fft(dp);
% normalization in a.s.d. appropriate for stationary time sereis
s = sqrt(2/N)*abs(m(1:M));
figure();
clf;
hold on;
set(gca,'LineWidth',2);
axis( [0, 2, 0, 1.1*max(s)]);
plot(p(2:M),s(2:M),'k-','LineWidth',2);
xlabel('period (days)');
ylabel('amplitude spectral density');
fprintf('Caption: Amplitude spectral amplitude as a function of period');
% gdama15_14
% earthquake location example
% rectangular earth volume, straight line rays
% data are travel times of P and S waves
% note: the problem is formulated here so that all
% the earthquakes are located simultaneously, with
% a single data kernel. Thus allows for the possibility
% of later adding constints thet involve several earthquakes
clear all;
fprintf('gdama15_14\n');
clear gdaGsparse gdaepsilon;
global gdaGsparse gdaepsilon;
gdaepsilon=1.0e-6;
% material constants
vpvs = 1.78;
vp=6.5;
vs=vp/vpvs;
% bounds of box
xmin=-10;
xmax=10;
ymin=-10;
ymax=10;
zmin=-10;
zmax=0;
% plot
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [xmin, xmax, ymin, ymax, zmin, zmax] );
% improvise outline of 3D box
plot3( [xmin,xmin], [ymin,ymin], [zmin,zmax], 'k-', 'LineWidth', 3 );
plot3( [xmin,xmin], [ymin,ymax], [zmin,zmin], 'k-', 'LineWidth', 3 );
plot3( [xmin,xmax], [ymin,ymin], [zmin,zmin], 'k-', 'LineWidth', 3 );
plot3( [xmax,xmax], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 3 );
plot3( [xmax,xmax], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 3 );
plot3( [xmax,xmin], [ymax,ymax], [zmax,zmax], 'k-', 'LineWidth', 3 );
plot3( [xmax,xmin], [ymin,ymin], [zmax,zmax], 'k-', 'LineWidth', 3 );
plot3( [xmax,xmax], [ymin,ymin], [zmax,zmin], 'k-', 'LineWidth', 3 );
plot3( [xmin,xmin], [ymax,ymin], [zmax,zmax], 'k-', 'LineWidth', 3 );
plot3( [xmin,xmin], [ymax,ymax], [zmax,zmin], 'k-', 'LineWidth', 3 );
plot3( [xmax,xmax], [ymax,ymin], [zmin,zmin], 'k-', 'LineWidth', 3 );
plot3( [xmax,xmin], [ymax,ymax], [zmin,zmin], 'k-', 'LineWidth', 3 );
xlabel('x_1');
ylabel('x_2');
zlabel('x_3');
% stations
sxm = [-8, -6, -4, -2, 0, 2, 4, 6, 8]'*ones(1,9);
sym = ones(9,1)*[-8, -6, -4, -2, 0, 2, 4, 6, 8];
sx = sxm(:);
sy = sym(:);
Ns = length(sx);
sz = zmax*ones(Ns,1);
plot3( sx, sy, sz, 'k+', 'LineWidth', 3 );
% earthquakes
Ne = 30; % number of earthquakes
M = 4*Ne; % 4 model parameters, x, y, z, and t0, per earthquake
extrue = random('uniform',xmin,xmax,Ne,1);
eytrue = random('uniform',ymin,ymax,Ne,1);
eztrue = random('uniform',zmin,zmax,Ne,1);
t0true = random('uniform',0,0.2,Ne,1);
mtrue = [extrue', eytrue', eztrue', t0true']';
Nr = Ne*Ns; % number of rays, that is, earthquake-stations pairs
N = 2*Ne*Ns; % data: P and S wave for each earthquake at each depth
% P times stored first in data array
% true data
% traveltime = length of ray divided by velocity
% ray is straight line connecting earthquake and station
dtrue=zeros(N,1);
for i = (1:Ns) % loop over stations
for j = (1:Ne) % loop over earthquakes
dx = mtrue(j)-sx(i);
dy = mtrue(Ne+j)-sy(i);
dz = mtrue(2*Ne+j)-sz(i);
r = sqrt( dx^2 + dy^2 + dz^2 );
k=(i-1)*Ne+j;
dtrue(k)=r/vp+mtrue(3*Ne+j);
dtrue(Nr+k)=r/vs+mtrue(3*Ne+j);
end
end
% observed data is true data plus random noise
sd = 0.1;
dobs=dtrue+random('normal',0,sd,N,1);
% inital guess of earthquake locations is random
mest = [random('uniform',xmin,xmax,1,Ne), random('uniform',ymin,ymax,1,Ne), ...
random('uniform',zmin+2,zmax-2,1,Ne), random('uniform',-0.1,0.1,1,Ne) ]';
% Geiger's method
for iter=(1:10)
% data kernel
% gdaG = spalloc(N,M,4*N);
% use row-col-value method of ceating sparse matrix
Gr = zeros(4*N,1);
Gc = zeros(4*N,1);
Gv = zeros(4*N,1);
Nel=0;
dpre = zeros(N,1);
for i = (1:Ns) % loop over stations
for j = (1:Ne) % loop over earthquakes
dx = mest(j)-sx(i);
dy = mest(Ne+j)-sy(i);
dz = mest(2*Ne+j)-sz(i);
r = sqrt( dx^2 + dy^2 + dz^2 );
k=(i-1)*Ne+j;
dpre(k)=r/vp+mest(3*Ne+j);
dpre(Nr+k)=r/vs+mest(3*Ne+j);
% gdaG(k,j) = dx/(r*vp);
Nel=Nel+1;
Gr(Nel) = k;
Gc(Nel) = j;
Gv(Nel) = dx/(r*vp);
% gdaG(k,Ne+j) = dy/(r*vp);
Nel=Nel+1;
Gr(Nel) = k;
Gc(Nel) = Ne+j;
Gv(Nel) = dy/(r*vp);
% gdaG(k,2*Ne+j) = dz/(r*vp);
Nel=Nel+1;
Gr(Nel) = k;
Gc(Nel) = 2*Ne+j;
Gv(Nel) = dz/(r*vp);
% gdaG(k,3*Ne+j) = 1;
Nel=Nel+1;
Gr(Nel) = k;
Gc(Nel) = 3*Ne+j;
Gv(Nel) = 1;
% gdaG(Nr+k,j) = dx/(r*vs);
Nel=Nel+1;
Gr(Nel) = Nr+k;
Gc(Nel) = j;
Gv(Nel) = dx/(r*vs);
% gdaG(Nr+k,Ne+j) = dy/(r*vs);
Nel=Nel+1;
Gr(Nel) = Nr+k;
Gc(Nel) = Ne+j;
Gv(Nel) = dy/(r*vs);
% gdaG(Nr+k,2*Ne+j) = dz/(r*vs);
Nel=Nel+1;
Gr(Nel) = Nr+k;
Gc(Nel) = 2*Ne+j;
Gv(Nel) = dz/(r*vs);
% gdaG(Nr+k,3*Ne+j) = 1;
Nel=Nel+1;
Gr(Nel) = Nr+k;
Gc(Nel) = 3*Ne+j;
Gv(Nel) = 1;
end
end
% truncate to actual length
Gr = Gr(1:Nel);
Gc = Gc(1:Nel);
Gv = Gv(1:Nel);
% build sparse matrix
gdaGsparse = sparse(Gr,Gc,Gv,N,M);
% delete lists
clear Gr Gc Gv;
% solve with dampled least squares
dd = dobs-dpre;
GTdd = gda_dlsrhs(dd);
dm=bicg(@gda_dlsmul,GTdd,1e-5,3*M);
mest = mest+dm;
end
% final predicted data
dpre=zeros(N,1);
for i = (1:Ns) % loop over stations
for j = (1:Ne) % loop over earthquakes
dx = mest(j)-sx(i);
dy = mest(Ne+j)-sy(i);
dz = mest(2*Ne+j)-sz(i);
r = sqrt( dx^2 + dy^2 + dz^2 );
k=(i-1)*Ne+j;
dpre(k)=r/vp+mest(3*Ne+j);
dpre(Nr+k)=r/vs+mest(3*Ne+j);
end
end
% display summary of results
expre = mest(1:Ne);
eypre = mest(Ne+1:2*Ne);
ezpre = mest(2*Ne+1:3*Ne);
t0pre = mest(3*Ne+1:4*Ne);
dd = dobs-dpre;
E = dd'*dd;
fprintf('RMS traveltime error: %f\n', sqrt(E/N) );
Emx = (extrue-expre)'*(extrue-expre);
Emy = (eytrue-eypre)'*(eytrue-eypre);
Emz = (eztrue-ezpre)'*(eztrue-ezpre);
Emt = (t0true-t0pre)'*(t0true-t0pre);
fprintf('RMS model misfit: x %f y %f z %f t0 %f\n', sqrt(Emx/Ne), sqrt(Emy/Ne), sqrt(Emz/Ne), sqrt(Emt/Ne) );
% plot results in 3D space
plot3( extrue, eytrue, eztrue, 'ro', 'LineWidth', 4 );
dx = 0.2;
dy = 0.2;
dz = 0.5;
for i = (1:Ne)
plot3( [expre(i)-dx, expre(i)+dx]', [eypre(i), eypre(i)]', [ezpre(i), ezpre(i)]', 'g-', 'LineWidth', 2 );
plot3( [expre(i), expre(i)]', [eypre(i)-dy, eypre(i)+dy]', [ezpre(i), ezpre(i)]', 'g-', 'LineWidth', 2 );
plot3( [expre(i), expre(i)]', [eypre(i), eypre(i)]', [ezpre(i)-dz, ezpre(i)+dz]', 'g-', 'LineWidth', 2 );
end
view(-12,52);
fprintf('Caption: Earthquake location problem with stations (black circles),\n');
fprintf('true eartquake locations (red) and estimated earthquake locations (green)\n')
% gdama12_15
% determine velocity structure v(x)=v0+v1(x)
% from frequencies of vibration
% unperturbed differetial equation -w^2 p = [v0+v1(x)]^2 d2p/dx2
% boundary conditions top: p=0 at x=0; bottom dp/dx=0 at x=h
% unperturbed solutions are analytic
% p(n,x) = A sin( (n-0.5) pi x / h ) n=1,2,3,...
% w(n) = (n-0.5)*pi*v0/h
% normalization so that (integral 0 to h) p(n,x) p(m,x) v0^(-2) dx = delta(n,m)
% A = sqrt( 2 v0^2 / h )
clear all;
fprintf('gdama12_15\n')
% independent variable x
h=2;
v0=5;
Nx = 101;
xmin=0;
xmax=h;
Dx = (xmax-xmin)/(Nx-1);
x = xmin + Dx*[0:Nx-1]';
% plot some wave p's and check normalization
figure(2);
clf;
verbose=0;
if( verbose )
fprintf('normalization (should be unity)\n');
end
for n =(1:5)
subplot(5,1,n);
set(gca,'LineWidth',3);
set(gca,'FontSize',12);
hold on;
p = sqrt(2*(v0^2)/h) * sin( (n-0.5) * pi * x / h );
plot(x,p,'k-','LineWidth',3);
I = Dx*sum(p.*p)/(v0^2);
if( verbose )
fprintf(' n %d area %f\n',n,I);
end
xlabel('x');
ylabel('p(x)');
end
fprintf('Caption, The first several eigenfunctions\n');
% build a true v1
v1true = zeros(Nx,1);
v1true(20:30)=v1true(20:30)-(v0/10);
v1true(60:90)=v1true(60:90)+(v0/10);
v1true(70:80)=v1true(70:80)+(v0/10);
figure(3);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [xmin, xmax, 0, 10] );
plot( x, v0+v1true, 'k-', 'LineWidth', 3 );
% unperturbed frequencies
N=200;
w0 = ([1:N]'-0.5)*pi*v0/h;
% perturbed frequencies
w1true = zeros(N,1);
for n =[1:N]
p = sqrt(2*(v0^2)/h) * sin( (n-0.5) * pi * x / h );
w1true(n)=w0(n)*Dx*sum(v1true.*p.*p)/(v0^3);
end
% observations are perturbed frequencies plus random noise
sigmad = 0.1*w0(1);
dobs = w1true + random('Normal',0.0,sigmad,N,1);
% ladder diagram; plot only a few frequencies
Nfew=10;
figure(4);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
axis( [0, max(w0(Nfew)), 0, 1] );
hold on
for n =(1:Nfew)
plot( [w0(n),w0(n)], [0, 1] ,'k-', 'LineWidth', 3);
plot( [w0(n)+w1true(n),w0(n)+w1true(n)], [0, 0.8] ,'r-', 'LineWidth', 3);
plot( [w0(n)+dobs(n),w0(n)+dobs(n)], [0, 0.6] ,'g-', 'LineWidth', 3);
end
xlabel('angular frequency w');
fprintf('Caption: Ladder diagram showing unperturbed eigenfrequencies (black),\n');
fprintf('true eigenfrequencies (red) and observed eigenfrequencies (green).\n');
% model parameters are velocity perturbations
M=Nx;
mest = zeros(M,1);
% data kernel is integral linking them
if( 0 ) % brute force way
G = zeros(N,M);
for i=(1:N)
for j=(1:M)
pij = sqrt(2*(v0^2)/h) * sin( (i-0.5) * pi * x(j) / h );
G(i,j) = Dx*w0(i)*(pij^2)/(v0^3);
end
end
else % tensor product way
pijm = sqrt(2*(v0^2)/h) * sin( pi * ([1:N]'-0.5)*x' / h );
G = Dx*(w0*ones(1,M)).*(pijm.*pijm)/(v0^3);
end
% damped least squares solution, but weight lower frequencies more
e2 = 1.0e-6;
wr = [1:M].^-1;
GMG = (G'*G+e2*diag(wr))\(G');
mest = GMG*dobs;
figure(3);
plot( x, v0+mest, 'r-', 'LineWidth', 3 );
xlabel('x');
ylabel('velocity m(x)');
fprintf('Caption: True model (black), estimated model (red).\n');
% resolution matrix
R = (GMG*G);
gda_draw(G,'caption G', ' ',R,'caption R');
fprintf('Caption: (left) Data kernel G and (right) model resolution matrix R.');