% gdama14
% 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
% gdama14_01
% Radon transform example. This version uses
% an image containing smooth features
% supports Section 11.6
% Note: The inversion of the Radom Transform
% requires working with 2D Fourier transforms,
% which is nitty-gritty stuff and requires
% extensive knowledge of their properties
% image mtrue(i,j) is N by N, in (x=i,y=j) space,
% range of both x, y is +/- 1
clear all;
fprintf('gdama14_01\n');
N = 256;
% independent variable x
Nx = N;
xmin = -1;
xmax = 1;
Dx = (xmax-xmin)/Nx;
x = xmin + Dx*[1:Nx]';
% independent variable y
Ny = N;
ymin = -1;
ymax = 1;
Dy = (ymax-ymin)/Ny;
y = ymin + Dy*[1:Ny]';
% true image m(x,y)
% test image is a sum of several Normal curves
xbar1=0.2;
ybar1=0;
s1=0.1;
mtrue = exp( -( ((x-xbar1).^2)*ones(1,Ny) + ones(Nx,1)*((y'-ybar1).^2))/(2*s1*s1));
xbar2=-0.3;
ybar2=0;
s2=0.05;
mtrue = mtrue+exp( -( ((x-xbar2).^2)*ones(1,Ny) + ones(Nx,1)*((y'-ybar2).^2))/(2*s2*s2));
xbar3=0;
ybar3=0.2;
s3=0.025;
mtrue = mtrue+exp( -( ((x-xbar3).^2)*ones(1,Ny) + ones(Nx,1)*((y'-ybar3).^2))/(2*s3*s3));
% some statistics, to compare with reconstructed image
fprintf('m has minimum value %f \n', min(min(mtrue)) );
fprintf('m has maximum value %f \n', max(max(mtrue)) );
fprintf('\n');
% plot true image m(x,y)
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
imagesc(mtrue);
%title('true image');
xlabel('y');
ylabel('x');
fprintf('Caption: True image\n');
% radon transform R(i,j) is in (u=i, q=theta=j) space
% range of u is +/- 1 (note corners of box cur off)
% range of q is 0 to pi
% independent variable u
Nu = N;
umin = -1;
umax = 1;
Du = (umax-umin)/Nu;
u = umin + Du*[1:Nu]';
% independent variable q
Nq = N+1;
qmin = -pi/2;
qmax = pi/2;
Dq = (qmax-qmin)/(Nq-1); % note includes both -pi/2 and pi/2
q = qmin + Dq*[0:Nq-1]'; % for ease of interpolation
% independent variable s
% integration over s, range +/- 1
Ns = 4*N;
smin = -1;
smax = 1;
Ds =(smax-smin)/Ns;
s = smin + Ds*[1:Ns]';
% brute force Radon transform via ray sums
R = zeros(Nu,Nq);
% loop over theta
for i = (1:Nq)
% object here is top build Nu by Ns array of image, m,
% evaluated at proper values of (u, s)
% (x,y) from (u,s); these are Nu by Ns matrices
xp = u*sin(q(i))*ones(1, Ns) - ones(Nu,1)*s'*cos(q(i))';
yp = u*cos(q(i))*ones(1, Ns) + ones(Nu,1)*s'*sin(q(i));
% indices corresponding to (x, y)
ixp = floor((xp-xmin)/Dx) + 1;
iyp = floor((yp-ymin)/Dy) + 1;
% handle out of range x indices
jxp1 = find(ixp<1);
ixp(jxp1)=1;
jxp2 = find(ixp>Nx);
ixp(jxp2)=Nx;
% handle out of range y indices
jyp1 = find(iyp<1);
iyp(jyp1)=1;
jyp2 = find(iyp>Ny);
iyp(jyp2)=Ny;
% nasty array manipulation to sample image, mtrue,
% at a series of (u, s) points, putting the values
% in an Nu by Ns matrix, S
% Step 1: convert Nu by Ns arrays to vectors
Lixp = ixp( [1:Nu*Ns]' );
Liyp = iyp( [1:Nu*Ns]' );
% Step 2: generate linear index into image
L = sub2ind( [Nx, Ny], Lixp, Liyp);
% Step 3: sample image, mtrue, and rearrange it back into Nu by Ns
% array
S = reshape(mtrue(L),Nu,Ns);
% now handle points that were off the edge of (x,y) space by
% setting values to zero
S( jxp1 ) = 0;
S( jxp2 ) = 0;
S( jyp1 ) = 0;
S( jyp2 ) = 0;
% integrate (sum) S(u,s) over s to get R(u) for fixed q
R(:,i) = Ds*sum(S,2);
end
% plot the radon transform
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
imagesc(R);
%title('Radon transform');
xlabel('theta');
ylabel('u');
fprintf('Caption: Radon transform of image.\n');
% statistics of the Radon transform
fprintf('R has minimum value %f \n', min(min(R)) );
fprintf('R has maximum value %f \n', max(max(R)) );
fprintf('\n');
% Inverse Radon Transform via Fourier transforms
% reorder array so that s=0, currenty at row N/2, is at row 1.
R = circshift(R, [-Nx/2+1,0] );
% now take 1D fourier transform over columns of R
Rt = fft(R);
% invert Radon transform via central slice theorem
% Fourier transform of image
mtt=zeros(Nx,Ny);
% fourier transform u -> ku
kumin=0;
kumax = 1/(2*Du);
Dku=kumax/(Nu/2);
ku=Dku*[0:Nu/2,-Nu/2+1:-1]';
Nku2 = Nu/2+1; % non-negative frequuencies
% fourier transform x -> kx
kxmin=0;
kxmax = 1/(2*Dx);
Dkx=kxmax/(Nx/2);
kx=Dkx*[0:Nx/2,-Nx/2+1:-1]';
Nkx2 = Nx/2+1; % non-negative frequuencies
% fourier transform y -> ky
kymin=0;
kymax = 1/(2*Dx);
Dky=kymax/(Ny/2);
ky=Dky*[0:Ny/2,-Ny/2+1:-1]';
Nky2 = Ny/2+1; % non-negative frequuencies
% Now use the Fourier Slice theorem. All the
% work is in sampling the radon transform at
% a series of points that correspond to wavenumber
% pairs on a rectangular grid. That requires
% the Radon tranform to be interpolated. I
% provide two ways of performing the interpolation.
% Note: Method 1 seems to work better. I think
% that's because my complex interpolator works
% better than using MatLab's TriScattered Interp
% on the real and imaginary parts, separately
if( 1 ) % Method 1: element by element
for i = (1:Nx)
for j = (1:Nky2) % only need left side, will rebuild right using symmetry
% theta calculation
theta = atan2( kx(i), ky(j) );
if( theta > (pi/2) ) % map large thetas onto negative thetas
theta=pi-theta;
thflag=1;
else
thflag=0;
end
% radial wavenumber
kr = sqrt(kx(i)^2 + ky(j)^2);
% index into randon transform array
if( kr > kumax )
continue; % ignore points that are out of bounds
end
% treat indices are real numbers
krindex = (kr-kumin)/Dku+1;
thindex = (theta-qmin)/Dq+1;
if( thindex<1 )
thindex=1;
elseif( thindex>Nq )
thindex=Nq;
end
% bracket real number indices with integers
% radial index
A = floor(krindex);
B = floor(krindex+1);
if( B >= Nku2 )
A = Nku2-1;
B= Nku2;
end
% theta index
C = floor(thindex);
D = floor(thindex+1);
if( D >= Nq )
C=Nq-1;
D=Nq;
end
% Lagrangian polynomial interpoaltion
% note that the indices are separated by unity
FC = (D-thindex);
CKR = ((B-krindex)*Rt(A,C)+(krindex-A)*Rt(B,C));
FD = (thindex-C);
DKR = ((B-krindex)*Rt(A,D)+(krindex-A)*Rt(B,D));
tmp = (FC*CKR+FD*DKR);
if (thflag==1)
mtt(i,j)=conj(tmp);
else
mtt(i,j)=tmp;
end
end
end
else % Method 2: all at once
% Rt is defined art these ku's and q's
kkuu = ku*ones(1,Nq);
qq = ones(Nu,1)*q';
% create an interpolant
% Rtinterpr = TriScatteredInterp( kkuu(:), qq(:), real(Rt(:)), 'natural' );
% Rtinterpi = TriScatteredInterp( kkuu(:), qq(:), imag(Rt(:)), 'natural' );
% MatLab's TriScatteredInterp() interpolator can't handle complex
% numbers, so interpolate real and imaginary parts separately. The
% Delaunay triangulation needs to be performed only once, so do it
% separately.
triangles = DelaunayTri( kkuu(:), qq(:) );
Rtinterpr = TriScatteredInterp( triangles, real(Rt(:)), 'natural' );
Rtinterpi = TriScatteredInterp( triangles, imag(Rt(:)), 'natural' );
% need all these elements of of mtt(i,j)
% ii = [1:Nx]'*ones(1,Nky2);
% jj = ones(Nx,1)*[1:Nky2];
% these indices have kr and theta
kkrr = sqrt( (kx(1:Nx)*ones(1,Nky2)).^2 + (ones(Nx,1)*ky(1:Nky2)').^2 );
qq = atan2( kx(1:Nx)*ones(1,Nky2), ones(Nx,1)*ky(1:Nky2)' );
% interpolate, set NaN's to zero, and populate mtt
tmp = complex( Rtinterpr( kkrr(:), qq(:) ), Rtinterpi( kkrr(:), qq(:) ) );
tmp(find(isnan(tmp)))=0;
mtt(1:Nx, 1:Nky2) = reshape( tmp, Nx, Nky2 );
end
% Fourier transform of image needs to be phase shifted
% else origin is in the wrong place
xshift = (xmax-xmin)/2;
yshift = (ymax-ymin)/2;
phase = 2*pi*(kx*ones(1,Ny)*xshift + ones(Nx,1)*(ky')*yshift);
mtt = mtt .* complex( cos(phase), sin(phase) );
% impose all necessary symmetries required for a real image
% these four elements must be real
mtt(1,1)=real(mtt(1,1));
mtt(Nx/2+1,1)=real(mtt(Nx/2+1,1));
mtt(1,Ny/2+1)=real(mtt(1,Ny/2+1,1));
mtt(Nx/2+1,Ny/2+1)=real(mtt(Nx/2+1,Ny/2+1));
% bottoms of these two columns must be conjugates of top
mtt(Nx:-1:Nx/2+2,1)=conj(mtt(2:Nx/2,1));
mtt(Nx:-1:Nx/2+2,Ny/2+1)=conj(mtt(2:Nx/2,Ny/2+1));
% right hand side flipped conjugate of left hand side
for m = (2:Ny/2)
mp=Ny-m+2;
mtt(1,mp) = conj(mtt(1,m));
mtt(2:Nx,mp) = flipud(conj(mtt(2:Nx,m)));
end
% now take inverse fourier transform
% Note: normalization not necessarily correct for non-square images
mest = ifft2(mtt)*(N/2);
mest = real(mest); % throw away imaginary part (which should be zero)
% plot reconstructed image
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
imagesc(mest);
% title('reconstructed image');
xlabel('y');
ylabel('x');
fprintf('Caption: Image reconstructed from Radon transform of image.\n');
% compare amplitudes of reconstructed image to original,
% to be sure that the normalization is correct
fprintf('mest has minimum value %f \n', min(min(mest)) );
fprintf('mest has maximum value %f \n', max(max(mest)) );
fprintf('\n')
% gdama14_02
% Radon transform example. This version uses
% an image from a file.
% Note: The inversion of the Radom Transform
% requires working with 2D Fourier transforms,
% which is nitty-gritty stuff and requires
% extensive knowledge of their properties
clear all;
fprintf('gdama14_02');
% image mtrue(i,j) is N by N, in (x=i,y=j) space,
% range of both x, y is +/- 1
N = 256;
% independent variable x
Nx = N;
xmin = -1;
xmax = 1;
Dx = (xmax-xmin)/Nx;
x = xmin + Dx*[1:Nx]';
% independent variable y
Ny = N;
ymin = -1;
ymax = 1;
Dy = (ymax-ymin)/Ny;
y = ymin + Dy*[1:Ny]';
% true image m(x,y)
% test image is a tiff file with a hypothetical depiction of a magma
% chamber; Image must be 256 by 256 in dimension to work.
mtrue=double(imread('../data/magma.tif','tif'));
mtrue = -(mtrue-mean(mean(mtrue)));
% some statistics, to compare with reconstructed image
fprintf('m has minimum value %f \n', min(min(mtrue)) );
fprintf('m has maximum value %f \n', max(max(mtrue)) );
fprintf('\n');
% plot true image m(x,y)
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [-1, 1, -1, 1] );
axis ij;
imagesc([-1, 1], [-1, 1], mtrue);
% title('true image');
xlabel('y');
ylabel('x');
fprintf('Caption: True image.\n');
% radon transform R(i,j) is in (u=i, q=theta=j) space
% range of u is +/- 1 (note corners of box cur off)
% range of q is 0 to pi
% independent variable u
Nu = N;
umin = -1;
umax = 1;
Du = (umax-umin)/Nu;
u = umin + Du*[1:Nu]';
% independent variable q
Nq = N+1;
qmin = -pi/2;
qmax = pi/2;
Dq = (qmax-qmin)/(Nq-1); % note includes both -pi/2 and pi/2
q = qmin + Dq*[0:Nq-1]'; % for ease of interpolation
% independent variable s
% integration over s, range +/- 1
Ns = 4*N;
smin = -1;
smax = 1;
Ds =(smax-smin)/Ns;
s = smin + Ds*[1:Ns]';
% radon transform R(u,q) tvia ray sums
R = zeros(Nu,Nq);
% loop over theta
for i = (1:Nq)
% object here is top build Nu by Ns array of image, m,
% evaluated at proper values of (u, s)
% (x,y) from (u,s); these are Nu by Ns matrices
xp = u*sin(q(i))*ones(1, Ns) - ones(Nu,1)*s'*cos(q(i))';
yp = u*cos(q(i))*ones(1, Ns) + ones(Nu,1)*s'*sin(q(i));
% indices corresponding to (x, y)
ixp = floor((xp-xmin)/Dx) + 1;
iyp = floor((yp-ymin)/Dy) + 1;
% handle out of range x indices
jxp1 = find(ixp<1);
ixp(jxp1)=1;
jxp2 = find(ixp>Nx);
ixp(jxp2)=Nx;
% handle out of range y indices
jyp1 = find(iyp<1);
iyp(jyp1)=1;
jyp2 = find(iyp>Ny);
iyp(jyp2)=Ny;
% nasty array manipulation to sample image, mtrue,
% at a series of (u, s) points, putting the values
% in an Nu by Ns matrix, S
% Step 1: convert Nu by Ns arrays to vectors
Lixp = ixp( [1:Nu*Ns]' );
Liyp = iyp( [1:Nu*Ns]' );
% Step 2: generate linear index into image
L = sub2ind( [Nx, Ny], Lixp, Liyp);
% Step 3: sample image, mtrue, and rearrange it back into Nu by Ns
% array
S = reshape(mtrue(L),Nu,Ns);
% now handle points that were off the edge of (x,y) space by
% setting values to zero
S( jxp1 ) = 0;
S( jxp2 ) = 0;
S( jyp1 ) = 0;
S( jyp2 ) = 0;
% integrate (sum) S(u,s) over s to get R(u) for fixed q
R(:,i) = Ds*sum(S,2);
end
% plot the radon transform
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [-1, 1, -1, 1] );
axis ij;
imagesc([-1, 1], [-1, 1], R);
% title('Radon transform');
xlabel('theta');
ylabel('u');
fprintf('Caption: Radon transform of true image.\n');
% statistics of the Radon transform
fprintf('R has minimum value %f \n', min(min(R)) );
fprintf('R has maximum value %f \n', max(max(R)) );
fprintf('\n');
% Inverse Radon Transform via Fourier transforms
% reorder array so that s=0, currenty at row N/2, is at row 1.
R = circshift(R, [-Nx/2+1,0] );
% now take 1D fourier transform over columns of R
Rt = fft(R);
% invert radon transform via central slice theorem
% fourier transform of image
mtt=zeros(Nx,Ny);
% fourier transform u -> ku
kumin=0;
kumax = 1/(2*Du);
Dku=kumax/(Nu/2);
ku=Dku*[0:Nu/2,-Nu/2+1:-1]';
Nku2 = Nu/2+1; % non-negative frequuencies
% fourier transform x -> kx
kxmin=0;
kxmax = 1/(2*Dx);
Dkx=kxmax/(Nx/2);
kx=Dkx*[0:Nx/2,-Nx/2+1:-1]';
Nkx2 = Nx/2+1; % non-negative frequuencies
% fourier transform y -> ky
kymin=0;
kymax = 1/(2*Dx);
Dky=kymax/(Ny/2);
ky=Dky*[0:Ny/2,-Ny/2+1:-1]';
Nky2 = Ny/2+1; % non-negative frequuencies
% Now use the Fourier Slice theorem. All the
% work is in sampling the radon transform at
% a series of points that correspond to wavenumber
% pairs on a rectangular grid. That requires
% the Radon tranform to be interpolated. I
% provide two ways of performing the interpolation.
% Note: Method 1 seems to work better. I think
% that's because my complex interpolator works
% better than using MatLab's TriScattered Interp
% on the real and imaginary parts, separately
if( 0 ) % Method 1: element by element
for i = (1:Nx)
for j = (1:Nky2) % only need left side, will rebuild right using symmetry
% theta calculation
theta = atan2( kx(i), ky(j) );
if( theta > (pi/2) ) % map large thetas onto negative thetas
theta=pi-theta;
thflag=1;
else
thflag=0;
end
% radial wavenumber
kr = sqrt(kx(i)^2 + ky(j)^2);
% index into randon transform array
if( kr > kumax )
continue; % ignore points that are out of bounds
end
% treat indices are real numbers
krindex = (kr-kumin)/Dku+1;
thindex = (theta-qmin)/Dq+1;
if( thindex<1 )
thindex=1;
elseif( thindex>Nq )
thindex=Nq;
end
% bracket real number indices with integers
% radial index
A = floor(krindex);
B = floor(krindex+1);
if( B >= Nku2 )
A = Nku2-1;
B= Nku2;
end
% theta index
C = floor(thindex);
D = floor(thindex+1);
if( D >= Nq )
C=Nq-1;
D=Nq;
end
% Lagrangian polynomial interpoaltion
% note that the indices are separated by unity
FC = (D-thindex);
CKR = ((B-krindex)*Rt(A,C)+(krindex-A)*Rt(B,C));
FD = (thindex-C);
DKR = ((B-krindex)*Rt(A,D)+(krindex-A)*Rt(B,D));
tmp = (FC*CKR+FD*DKR);
if (thflag==1)
mtt(i,j)=conj(tmp);
else
mtt(i,j)=tmp;
end
end
end
else % Method 2: all at once
% Rt is defined art these ku's and q's
kkuu = ku*ones(1,Nq);
qq = ones(Nu,1)*q';
% create an interpolant
% Rtinterpr = TriScatteredInterp( kkuu(:), qq(:), real(Rt(:)), 'natural' );
% Rtinterpi = TriScatteredInterp( kkuu(:), qq(:), imag(Rt(:)), 'natural' );
% MatLab's TriScatteredInterp() interpolator can't handle complex
% numbers, so interpolate real and imaginary parts separately. The
% Delaunay triangulation needs to be performed only once, so do it
% separately.
triangles = DelaunayTri( kkuu(:), qq(:) );
Rtinterpr = TriScatteredInterp( triangles, real(Rt(:)), 'natural' );
Rtinterpi = TriScatteredInterp( triangles, imag(Rt(:)), 'natural' );
% need all these elements of of mtt(i,j)
% ii = [1:Nx]'*ones(1,Nky2);
% jj = ones(Nx,1)*[1:Nky2];
% these indices have kr and theta
kkrr = sqrt( (kx(1:Nx)*ones(1,Nky2)).^2 + (ones(Nx,1)*ky(1:Nky2)').^2 );
qq = atan2( kx(1:Nx)*ones(1,Nky2), ones(Nx,1)*ky(1:Nky2)' );
% interpolate, set NaN's to zero, and populate mtt
tmp = complex( Rtinterpr( kkrr(:), qq(:) ), Rtinterpi( kkrr(:), qq(:) ) );
tmp(find(isnan(tmp)))=0;
mtt(1:Nx, 1:Nky2) = reshape( tmp, Nx, Nky2 );
end
% Fourier transform of image needs to be phase shifted
% else origin is in the wrong place
xshift = (xmax-xmin)/2;
yshift = (ymax-ymin)/2;
phase = 2*pi*(kx*ones(1,Ny)*xshift + ones(Nx,1)*(ky')*yshift);
mtt = mtt .* complex( cos(phase), sin(phase) );
% impose all necessary symmetries required for a real image
% these four elements must be real
mtt(1,1)=real(mtt(1,1));
mtt(Nx/2+1,1)=real(mtt(Nx/2+1,1));
mtt(1,Ny/2+1)=real(mtt(1,Ny/2+1,1));
mtt(Nx/2+1,Ny/2+1)=real(mtt(Nx/2+1,Ny/2+1));
% bottoms of these two columns must be conjugates of top
mtt(Nx:-1:Nx/2+2,1)=conj(mtt(2:Nx/2,1));
mtt(Nx:-1:Nx/2+2,Ny/2+1)=conj(mtt(2:Nx/2,Ny/2+1));
% right hand side flipped conjugate of left hand side
for m = (2:Ny/2)
mp=Ny-m+2;
mtt(1,mp) = conj(mtt(1,m));
mtt(2:Nx,mp) = flipud(conj(mtt(2:Nx,m)));
end
% now take inverse fourier transform
% Note: normalization not necessarily correct for non-square images
mest = ifft2(mtt)*(N/2);
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
colormap('jet');
hold on;
axis( [-1, 1, -1, 1] );
axis ij;
imagesc([-1, 1], [-1, 1], mest);
title('reconstructed image');
xlabel('y');
ylabel('x');
fprintf('Caption: Image reconstructed from Radon trsnsform.\n');
fprintf('mest has minimum value %f \n', min(min(mest)) );
fprintf('mest has maximum value %f \n', max(max(mest)) );
fprintf('\n');
% gdama14_03
% gradient descent inversion of d(x)=Lm(x),
% where L is a linear operator. The gradient of
% error dE/dm is calculated using the adjoint method
% The linear operator is 0.1*(D + 2*I) where D is a
% first derivative and I is an integral
clear all;
fprintf('gdama14_03');
% auxiliary variable x
M=101;
N=M;
Dx=1;
x = Dx*[0:M-1]'/(M-1);
xmax=max(x);
% a true model function m(x)that is zero at the boundaries
mtrue = sin(4*pi*x/xmax)+0.5*sin(7*pi*x/xmax)+0.25*sin(13*pi*x/xmax);
% derivative and integral operators
D = (1/Dx)*toeplitz( [1, -1, zeros(1,M-2)]', [1, zeros(1,M-1)] );
I = Dx*toeplitz( ones(M,1), [1, zeros(1,M-1)] );
% data is a linear operator L on the model parameters, d=Lm
% where the linear operator contains both a derivative and an intergal:
% d = 0.1*( dm/dx + 2 * (integral 0 x) m dm )
L = 0.1*(D + 2*I);
dtrue = L*mtrue;
% brute force dE/dm, to check
% mo = ones(M,1);
% do = L*mo;
% Eo = (dtrue-do)'*(dtrue-do);
% dEdm = zeros(M,1);
% Dm=0.00001;
% for i=[1:M]
% dm = zeros(M,1);
% dm(i)=Dm;
% m1 = mo+dm;
% d1 = L*m1;
% E1 = (dtrue-d1)'*(dtrue-d1);
% dEdm(i) = (E1-Eo)/Dm;
% end
% dEdmA = -2*(L'*(dtrue-do));
% fprintf('dEdm\n' );
% for i=[1:M]
% fprintf('%d %f %f\n', i, dEdm(i), dEdmA(i) );
% end
% trial solution
mgo = ones(M,1);
mgo1 = mgo;
% maximum iterations
Niter=100000;
% save error history
Ehist = zeros(Niter+1,1);
% error and its gradient at the trial solution
dgo = L*mgo;
Ego = (dtrue-dgo)'*(dtrue-dgo);
dEdmo = -2*Dx*(L'*(dtrue-dgo));
tau=0.5;
c1=0.0001;
alpha=0.9;
Ehist(1) = Ego;
% gradient descent method
for k = (1:Niter)
v = -dEdmo / sqrt(dEdmo'*dEdmo);
% backstep
for kk=(1:10)
mg = mgo+alpha*v;
dg = L*mg;
Eg = (dtrue-dg)'*(dtrue-dg);
dEdm= -2*(L'*(dtrue-dg));
if( (Eg<=(Ego + c1*alpha*v'*dEdmo)) )
break;
end
alpha = tau*alpha;
end
% change in solution
Dmg = sqrt( (mg-mgo)'*(mg-mgo) );
% update
mgo=mg;
dgo = dg;
Ego = Eg;
dEdmo = dEdm;
Ehist(k+1) = Ego;
% plot one intermediate result
if( k==40 )
mgo2=mgo;
end
% terminate on small change
if( Dmg < 1.0e-6 )
break;
end
end
fprintf('%d iterations\n',k);
% plot the true model parameters
figure();
clf;
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
plot(x,mtrue,'k-','LineWidth',6);
xlabel('x');
ylabel('m(x)');
plot(x,mgo1,'g:','LineWidth',3);
plot(x,mgo,'g-','LineWidth',3);
plot(x,mgo2,'g--','LineWidth',3);
fprintf('Caption: True solution (black), starting solution (green dotted line),\n');
fprintf('solution in progress (green dashed curve), estimated solution (green curve)\n');
% plot the data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
plot(x,dtrue,'r-','LineWidth',3);
plot(x,dgo,'k-','LineWidth',1);
xlabel('x');
ylabel('d(x)');
fprintf('Caption: True data (black curve), predicted data (red curve).\n');
% plot error history
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [1, 101, -4, 0 ] );
plot([1:101],log10(Ehist(1:101)/Ehist(1)),'k-','LineWidth',3);
xlabel('iteration i');
ylabel('log10 error, E(i)');
fprintf('Caption: Prediction error E Eas a function of iteration number i.\n');
% gdama14_04
% simple 1D backprojection example
% the relationship d(x) = L m(x) is "inveted"
% using the back-propagation approximation m=L'd
% (where L' is the adjoint of L). In this
% example, L is the integral from 0 to x, Linv
% is the first derivative, and its adjoint is the
% integral from infinity to x
clear all;
fprintf('gdama14_04\n');
% independent variable x
M=201;
Dx=1.0;
x = Dx*[0:M-1]';
xmax=max(x);
% true model
sigma=12;
sigma2=sigma^2;
xbar1 = 3*xmax/8;
xbar2 = 5*xmax/8;
f1=1.0;
f2=2.0;
mtrue = sin(f1*x).*exp(-((x-xbar1).^2)/(2*sigma2)) + sin(f2*x).*exp(-((x-xbar2).^2)/(2*sigma2));
mtrue=mtrue-mean(mtrue);
% data is integral of model
dobs = Dx*cumsum(mtrue);
% plot true model
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, M, min(mtrue), max(mtrue) ] );
plot( x, mtrue, 'k-', 'LineWidth', 2 );
xlabel('x');
ylabel('m');
%title('m true');
fprintf('Caption: True solution\n');
% exact solution is m = Linv d
% where Linv the derivative operator d/dx
mest1 = diff(dobs)/Dx;
% plot exact solution
figure();
clf;
set(gca,'LineWidth',3);
hold on;
axis( [0, M, min(mest1), max(mest1) ] );
plot( x(1:end-1), mest1, 'k-', 'LineWidth', 2 );
xlabel('x');
ylabel('m');
%title('m exact');
fprintf('Caption: Estimated solution using exact inversion.\n');
% solution via backprojection
mest2 = flipud(Dx*cumsum(flipud(dobs)));
% plot back-projected solution
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, M, min(mest2), max(mest2) ] );
plot( x, mest2, 'k-', 'LineWidth', 2 );
xlabel('x');
ylabel('m');
% title('m back-projected');
fprintf('Caption: Estimated solution using back-projection.\n');
% gdama14_05
% backprojection in a a simple tomography example
% The data are "travel times", that is integrals
% of m(x,y) along straight line rays between sources
% and receivers on the edges of a square region.
% Supports Figure 11.7
clear all;
fprintf('gdsma14_05');
clear gdaGsparse gdaepsilon;
global gdaGsparse gdaepsilon;
% shortest ray allowed
Rmin=5;
% true model m(x,y) is read 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 of 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
figure();
clf;
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' );
colorbar();
fprintf('Caption: True model.\n');
% define sources and receivers on the edges of the box
Nside = 25;
Nside = 50;
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;
% source receiver pairs
N = Ns*(Ns-1)/2;
ray = zeros(N,4);
k=0;
for i = (1:Ns)
for j = (i+1:Ns)
R = sqrt( (xs(i,1)-xs(j,1))^2 + (xs(i,2)-xs(j,2))^2 );
if (R<Rmin) % exclude really short rays
continue;
end
k = k+1;
ray(k,1) = xs(i,1);
ray(k,2) = xs(i,2);
ray(k,3) = xs(j,1);
ray(k,4) = xs(j,2);
end
end
N=k;
% build data kernel
% order of model parameters in vector: index = (ix-1)*Ny + iy
M = Nx*Ny;
% use row-column-value method, instead of gdaGsparse=spalloc(N,M,2*N*Nx);
NELEMENTS = 2*N*Nx;
Gr = zeros(NELEMENTS,1);
Gc = zeros(NELEMENTS,1);
Gv = zeros(NELEMENTS,1);
Nel=0;
Nr = 2000;
bins=[1:Nx*Ny];
for k = (1:N)
% also use build non-sparse row Gk first method
Gk = zeros(M,1);
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
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;
% gdaGsparse(k,q) = gdaGsparse(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 );
% gdaGsparse(k,icount) = gdaGsparse(k,icount) + count(icount)*dr;
Gk(icount)=Gk(icount)+count(icount)'*dr;
end
% now accumulate Gk into element arrays
for i=(1:M)
if( Gk(i)~=0.0 )
Nel=Nel+1;
Gr(Nel)=k;
Gc(Nel)=i;
Gv(Nel)=Gk(i);
end
end
end % next ray k
% delete unused elements
Gr=Gr(1:Nel,1);
Gc=Gc(1:Nel,1);
Gv=Gv(1:Nel,1);
gdaGsparse = sparse(Gr,Gc,Gv,N,M);
% dont clear element arrays yet
% compute true travel times
d = gdaGsparse*mtrue;
% inversion of data by damped least squares
gdaepsilon = 1.0e-2*max(max(abs(gdaGsparse)));
GTd = gda_dlsrhs(d);
mest1=bicg(@gda_dlsmul,GTd,1e-5,3*M);
% inversion of data by back projection
% normalize each row of the data kernel to unity
% with compensatory normalization of the data
rowsum = sum(gdaGsparse,2);
rowsumr = 1./rowsum;
rowsumr(find(~isfinite(rowsumr)))=0;
dp = d .* rowsumr;
% manipulaste row-col-value arrays instead of
% GP = gdaGsparse;
% for i=(1:M)
% GP(:,i) = gdaGsparse(:,i) .* rowsumr;
% end
for el=(1:Nel)
k = Gr(el);
Gv(el) = Gv(el) * rowsumr(k);
end
GP = sparse(Gr,Gc,Gv,N,M);
clear Gr Gc Gv; % now clear them
mest2=GP'*dp;
% rearrange m back into image
Sest1 = zeros(Nx,Ny);
Sest2 = zeros(Nx,Ny);
for k = (1:M)
ix=rowindex(k);
iy=colindex(k);
Sest1(ix,iy)=mest1(k);
Sest2(ix,iy)=mest2(k);
end
% plot estimated image (damped least squares)
figure();
clf;
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], Sest1, 'CDataMapping', 'scaled');
xlabel('y');
ylabel('x');
% title( 'LS estimate' );
colorbar();
fprintf('Caption: Estimated image, using damped least squares\n');
% plot estimated image (back projection)
figure();
clf;
set(gca, 'YDir', 'reverse' );
hold on;
set(gca, 'LineWidth', 3 );
set(gca, 'FontSize', 14 );
colormap('jet');
axis( [xmin, xmax, ymin, ymax] );
n=0; % exclues edges
Sest2 = Sest2-mean(Sest2(:));
zmin=-100;
zmax=10;
caxis( [zmin, zmax] ); % color scale fudged here
image( [xmin, xmax], [ymin, ymax], Sest2, 'CDataMapping', 'scaled');
xlabel('y');
ylabel('x');
% title( 'BP estimate' );
colorbar();
fprintf('Caption: Estimated image, using back projection\n');
% gdama14_06
% temperature / chemical reaction example
% adjoint method used to build data kernel
% and the problem is then solved usign least squares
% kernel for differential equation
clear all;
fprintf('gdama14_06');
% time variable t
M=501;
N=M;
Dt=1;
t = Dt*[0:M-1]';
tmax=max(t);
% differential equation (d/dt + c)u = f
% has green function G(t,t')=H(t-t')exp(-c(t-t'))
c = 0.04;
% Forward problem
% time series plots
figure();
clf;
% create a heat function m
% allow two options, complex and simple
if( 1 )
% complicated m built out of Gaussians
t1 = 30;
sigma1 = 6;
m = exp(-0.5*((t-t1).^2)/(sigma1^2));
t1 = 50;
sigma1 = 10;
m = m+2*exp(-0.5*((t-t1).^2)/(sigma1^2));
t1 = 80;
sigma1 = 14;
m = m+0.5*exp(-0.5*((t-t1).^2)/(sigma1^2));
t1 = 200;
sigma1 = 5;
m = m+exp(-0.5*((t-t1).^2)/(sigma1^2));
else
% very simple m for test purposes
m=zeros(M,1);
m(floor(1*M/4))=1;
m(floor(2*M/4))=0.5;
end
% perform the Green function integral to predict u
% use convolution function; faster than coding the integral
Ho = exp(-c*t);
u = Dt*conv(Ho,m); % Green function integral is a convolution
u=u(1:M);
% integrate u to get P
b=1;
P = Dt*b*cumsum(u);
% data is P with random noise
dtrue = P;
sigmad = 0.1;
dobs = dtrue+random('Normal',0,sigmad,N,1);
% plot Green function for t'=10;
subplot(4,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, tmax, 0, 1] );
tp=30;
H = (t>tp).*exp(-c*(t-tp));
plot( t, H, 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('H(t,10)');
% plot heat function
subplot(4,1,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, tmax, 0, 1.2*max(m)] );
plot( t, m, 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('m(t)');
subplot(4,1,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, tmax, 0, max(u)] );
plot( t, u, 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('u(t)');
% plot the observed data
subplot(4,1,4);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, tmax, 0, 1.1*max(dobs)] );
plot( t, dobs, 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('d(t)');
fprintf('Caption: (Top row) Green function, (2nd row) true solution,\n');
fprintf('(3rd row) field u, (bottom row) true data.\n');
% inverse problem
% adjoint differential equation (-d/dt + c)u = h
% has green function G(t,t')=H(t'-t)exp(c(t-t'))
% data kernel Gd
N=M;
Gd = zeros(N,M);
for i = (1:N)
for j = (1:M)
ti = i*Dt;
tee = j*Dt;
if( tee <= ti )
Gd(i,j)=Dt*(-b/c)*(exp(-c*(ti-tee))-1);
else
Gd(i,j)=0;
end
end
end
% solve with dampled least squares
e2 = 1.0e-1;
mest = (Gd'*Gd + e2*eye(M,M))\(Gd'*dobs);
% plot estimsted solution
subplot(4,1,2);
plot( t, mest, 'r--', 'LineWidth', 3 );
plot( t, m, 'k-', 'LineWidth', 2 );
% plot green function for t'=10;
figure();
clf;
subplot(4,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, tmax, 0, 1] );
tp=100;
H = Dt*(t<tp).*exp(c*(t-tp));
plot( t, H, 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('Q(t,tau=100)','FontSize',10);
% plot a few rowa of the data kernel
subplot(4,1,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
i=floor(1*N/4);
axis( [0, tmax, 0, 1.1*max(Gd(i,:)) ] );
plot( t, Gd(i,:)', 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('g(i=125,t)','FontSize',10);
subplot(4,1,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
i=floor(2*N/4);
axis( [0, tmax, 0, 1.1*max(Gd(i,:)) ] );
plot( t, Gd(i,:)', 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('g(i=250,t)','FontSize',10);
subplot(4,1,4);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
i=floor(3*N/4);
axis( [0, tmax, 0, 1.1*max(Gd(i,:)) ] );
plot( t, Gd(i,:)', 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('g(i=375,t)','FontSize',10);
fprintf('Caption: (top row) Adjoint Green function. (2nd-4th row)\n');
fprintf('selected rows of the data kernel.\n');
% draw the data kernel
gda_draw(Gd);
fprintf('Caption: Data kernel G\n');
% gdama14_07
% Uses heat flow problem to illustrate the use of
% the adjoint method to compute the data kernel dd/dm.
% This deriavtive is calculated two ways, using finite
% difference and using the adjoint method (achieving the
% same result, of course).
clear all;
fprintf('gdama14_07');
figure();
clf;
% time t setup
N=101;
Dt = 0.1;
i0 = floor(N/10);
t = Dt*[-i0+1:N-i0]';
% reference solution is analytic
c0 = 0.7;
u0 = (t>0) .* exp( -c0*t );
% plot the reference solution
subplot(3,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [t(1), t(end), -1.1, 1.1] );
xlabel('time t');
ylabel('u');
% the data is the b times the integral of the solution
% predicted data associated with the reference field
b = 0.5;
d0 = b*Dt*cumsum( (t>0) .* u0 );
% point scatterer at t0, that is
% c(t) = c0 + dm * delta(t-t0)
t0 = 0.5; % time of scatterer
it0 = find(t>=t0,1); % index of scatterer
dm = 0.3; % amplitude of scatterer
% field perturbation according the the Born approximation
du = -dm * u0(it0) * (t>t0) .* exp( -c0*(t-t0) );
u = u0 + du; % total field
d = b*Dt*cumsum( (t>0) .* u ); % predicted data
% plot the reference and total field
subplot(3,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
plot( t, u, 'k-', 'LineWidth', 3 );
plot( t, u0, 'r-', 'LineWidth', 2 );
% plot the predicted data associated with u and u0
% plot the data
subplot(3,1,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [t(1), t(end), -1.1, 1.1] );
xlabel('time t');
ylabel('d');
plot( t, d, 'k-', 'LineWidth', 3 );
plot( t, d0, 'r-', 'LineWidth', 2 );
% derivative by finite difference
dddm = (d-d0)/dm;
subplot(3,1,3);
set(gca,'LineWidth',3);
hold on;
yr = 1.1;
axis( [t(1), t(end), -yr, yr] );
plot( t, dddm, 'k-', 'LineWidth', 3 );
xlabel('time t');
ylabel('dd/dm');
% derivative by adjoint formula
% which should (and does) exactly match the adjoint formula
dddm2 = -(b/c0)*(t>0).*(t>t0).*(1-exp(-c0*(t-t0))).*exp(-c0*t0);
plot( t, dddm2, 'r-', 'LineWidth', 2 );
fprintf('Caption: (Top) Field u (black curve) and reference field u0 (red curve),\n');
fprintf('(middle) data predicted by the field (black) and by the reference field (red),\n');
fprintf('(Bottom) derivative dd/dm by finite differences (black) and the adjoint method\n');
% gdama14_08
% data kernel for heat flow problem with time-variable thermal constant c(t)
clear all;
fprintf('gdama14_08');
% data constant
b = 0.5;
% time t axis
N=101;
Dt = 0.1;
i0 = floor(N/10);
t = Dt*(-i0:(N-i0-1))';
tmax = t(N-1);
% thermal constant
c00 = 0.7;
c01 = 0.2;
% reference c
c0 = c00*ones(N,1);
% true c
ctrue = c0+c01*sin(pi*t/tmax).*(t>0.0);
% true solution and correponding data
utrue = zeros(N,1);
utrue(i0)=1.0;
for i = (i0+1:N)
utrue(i) = utrue(i-1) - ctrue(i)*Dt*utrue(i-1);
end
dtrue = b*Dt*cumsum(utrue);
% reference solution and corresponding data
u0 = zeros(N,1);
u0(i0)=1.0;
for i = (i0+1:N)
u0(i) = u0(i-1) - c0(i)*Dt*u0(i-1);
end
d0 = b*Dt*cumsum(u0);
figure();
clf;
hold on;
axis( [t(1), t(N), -1.1, 1.1] );
xlabel('time t');
ylabel('u');
plot( t, utrue, 'k-', 'LineWidth', 3);
plot( t, u0, 'r-', 'LineWidth', 2);
% title('u');
fprintf('Caption: True field u (black curve) and reference field (red curve)\n');
figure();
clf;
hold on;
axis( [t(1), t(N), -1.1, 1.1] );
xlabel('time t');
ylabel('d');
plot( t, dtrue, 'k-', 'LineWidth', 3);
plot( t, d0, 'r-', 'LineWidth', 2);
%title('d');
fprintf('Caption: True data d (black curve) and data predicted by the reference field (red curve)\n');
ddidmj = zeros(N,N);
for i = (1:N) % loop over data points
% solve adjoint equation
% ( -d/dt + c0 ) lami = hi where hi = b H(ti-t)
hi = b*(t(i)*ones(N,1)>=t);
lami = zeros(N,1);
% descrete version of adjoint equation
% (-lamj(k+1)+lamj(k))/Dt + c0 lamj(k+1) = hi
% lamj(k) = lamj(k+1) - c0 = Dt lamj(k+1) + Dt hi
for k = ((N-1):-1:1)
lami(k) = lami(k+1) - c0(k+1)*Dt*lami(k+1) + Dt*hi(k+1);
end
% perturbed equation is L u = 0 with
% L = ( d/dt + c0 + mj delta(t-tj) )
% so dLdmj = delta(t-tj)
% and xi = - dLdmj u0
% perform inner product
% ddidmj = ( lami, dLdmi u0 ) = -( lami, delta(t-tj) u0 ) =
% ddidmj = - lami(j) u0(j)
ddidmj(i,1:N) = - ( lami .* u0 )';
end
G = ddidmj; % data kernel
gda_draw(G,'caption G');
fprintf('Caption: Data kernel\n');
% gdama14_09
% solve inverse problem for the function c(t) in
% the cooling problem using adjoint methods
% the equation of cooling is
% ( d/dt + c(t) ) u = 0 with initial condition u(0)=1
% the true c(t) is a constant plus a sinusoid
% the problem is solved iteratively with Newton's
% methid using the linearized data kernel Gij = ddi/dmj
% calculated by adjoint methods.
% both the original equation and the adjoint equation
% are solved numerically with a very simple recursion based on
% Euler's method. A more accurate method, such as Runga
% Kutta could be substituted.
% At each iteration, the function c(t) is parameteried
% c = c0 + (sum over i) mi delta(t-ti)
% which implies that the update is c0 = c0 + m/Dt
% where the factor of Dt arises from the delta funcion
% that is, the integral of the mi delta(t-ti) over one
% pixel is mi, and the integral of mi/Dt over one pixel
% is also mi
% syntheic data are created by adding a small amount of
% random noise to the true data.
% The problem solved via generalized least squares, with
% both flatness and smallness prior information. A good
% solution is achieved in about five iterations.
clear all;
fprintf('gdama14_09');
% data constant
b = 0.5;
% time t axis
N=101;
Dt = 0.1;
i0 = floor(N/10);
t = Dt*(-i0:(N-i0-1))';
tmax = t(N-1);
% thermal constant
c00 = 0.7;
c01 = 0.2;
% reference c
c0 = c00*ones(N,1);
cstart = c0;
% true c
ctrue = c0+c01*sin(pi*t/tmax).*(t>0.0);
% true solution and correponding data
utrue = zeros(N,1);
utrue(i0)=1.0;
for i = (i0+1:N)
utrue(i) = utrue(i-1) - ctrue(i)*Dt*utrue(i-1);
end
dtrue = b*Dt*cumsum(utrue);
sd = 0.0001;
dobs = dtrue + random('Normal',0.0,sd,N,1);
% flatness matrix for generalized least squares
r = [-1.0; zeros(N-2,1)];
c = [-1.0, 1.0, zeros(1,N-2) ];
H = toeplitz(r,c);
h1 = zeros(N-1,1); % right hand side
h0=zeros(N,1);
maxiter = 10;
for iter =(1:maxiter)
% reference solution and corresponding data
u0 = zeros(N,1);
u0(i0)=1.0;
for i = (i0+1:N)
u0(i) = u0(i-1) - c0(i)*Dt*u0(i-1);
end
d0 = b*Dt*cumsum(u0);
ddidmj = zeros(N,N);
for i = (1:N) % loop over data points
% solve adjoint equation
% ( -d/dt + c0 ) lami = hi where hi = b H(ti-t)
hi = b*(t(i)*ones(N,1)>=t);
lami = zeros(N,1);
% descrete version of adjoint equation
% (-lamj(k+1)+lamj(k))/Dt + c0 lamj(k+1) = hi
% lamj(k) = lamj(k+1) - c0 = Dt lamj(k+1) + Dt hi
for k = (N-1:-1:1)
lami(k) = lami(k+1) - c0(k+1)*Dt*lami(k+1) + Dt*hi(k+1);
end
% perturbed equation is L u = 0 with
% L = ( d/dt + c0 + mj delta(t-tj) )
% so dLdmj = delta(t-tj)
% and xi = - dLdmj u0
% perform inner product
% ddidmj = ( lami, dLdmi u0 ) = -( lami, delta(t-tj) u0 ) =
% ddidmj = - lami(j) u0(j)
ddidmj(i,1:N) = - ( lami .* u0 )';
end
G = ddidmj; % data kernel
% inverse problem
% G Dm = dtrue-d0
% H Dm = h1 - H c0 or Hm = h1
% I Dm = 0
Dd = dobs-d0;
if(1)
epsi1 = 1e-2;
epsi2 = 0.4e-1;
Dh = h1 - H*c0; % apply prior info to solution
else
epsi1 = 1e-1;
epsi2 = 1e-2;
Dh = h1; % apply to perturbation, only
end
% generalized least squares with prior informaiton of flatness and smallness
% flatness can be applied to either solution or perturbation
% smallness is applied to perturbation
F = [G; epsi1*H; epsi2*eye(N) ];
f = [Dd; epsi1*Dh; epsi2*h0 ];
FTF = F'*F;
FTf = F'*f;
m = FTF\FTf;
c0 = c0 + m/Dt; % factor of Dt from delta function
end
figure();
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize',14);
axis( [t(1), t(N), 0.0, 1.5] );
xlabel('time t');
ylabel('c');
plot( t, ctrue, 'k-', 'LineWidth', 3);
plot( t, c0, 'r-', 'LineWidth', 2);
plot( t, cstart, 'g--', 'LineWidth', 2);
fprintf('Caprion: True model parameter c (black curve), reference value\n');
fprintf('(green curve) and estimated (red curve)\n');
% gdama14_10
% Data Assimilation example
% A insulated rod of length Lx cools because its ends
% are held at a temperature u(0)=u(Lx)=0. Heat moves
% according to the thermal diffusion equation
clear all;
fprintf('gdama14_10');
% x-axis
Dx = 1.0;
Lx = 128;
x = Dx*[0:Lx-1]';
xmax = x(Lx);
xmid = xmax/2.0;
M=Lx;
% t-axis
Dt = 1.0;
Lt = 128;
t = Dt*[0:Lt-1]';
tmax = t(Lt);
% thermal constant
kappa = 0.05/Dt;
% second derivative matrix
r = [-2.0; 1.0; zeros(Lx-2,1) ];
c = [-2.0, 1.0, zeros(1,Lx-2) ];
D2 = toeplitz( r,c) / (Dx^2);
D2(1,1)=1.0;
D2(1,2)=0.0;
D2(Lx,Lx-1)=0.0;
D2(Lx,Lx)=1.0;
% dynamics matrix
D = Dt*kappa*D2 + eye(Lx);
% source
sx = 5.0;
sx2 = sx^2;
m0 = exp( -0.5 * ((x-xmid).^2) / sx2 );
% true solution, by stepping the equation forward in time
utrue = zeros(Lx,Lt);
utrue(1:Lx,1) = m0;
for i = (2:Lt)
utrue(1:Lx,i) = D*utrue(1:Lx,i-1);
end
% data positions and time are randomly assigned
Ndt = 35; % data per time
N = (Lt-1)*Ndt; % total number of data
kd = zeros(Ndt,Lt-1);
for i = (2:Lt)
% no data at zero time and no data at boundaries
j = randperm(Lx-2)+1;
% index array kd gives locations of data
kd(1:Ndt,i-1) = sort(j(1:Ndt));
end
% although noise is defined everywhere, only values at
% data positions-times are used
sd = 0.1;
nse = random('Normal',0.0,sd,Lx,Lt);
% guessed source
sx = 2.0; % half-width of Gaussian source
sx2 = sx^2;
mg = exp( -0.5 * ((x-xmid).^2) / sx2 );
size(mg)
Dm = 0.05;
maxiter = 25;
Emax = zeros(maxiter,1);
for iter = (1:maxiter)
% guessed solution, treat source as initial condition
% and use recurion to propagate solution forward in time
u = zeros(Lx,Lt);
u(1:Lx,1) = mg;
for i = (2:Lt)
u(1:Lx,i) = D * u(1:Lx,i-1);
end
% error = observed temperature - predicted temperature
% observed temperature = true temperature + noise
e = zeros(Lx,Lt); % although the error is defined
% for all poitions-times, it is non-zero only at observation points
for i = (2:Lt)
k = kd(1:Ndt,i-1);
e(k,i) = (utrue(k,i)+nse(k,i))-u(k,i);
end
% total error
emax=sum(e(:).^2);
Emax(iter) = emax;
% simplified gradent method
if( (iter>1) && (emax>emaxlast) )
emaxlast = emax;
Dm = Dm/2.0;
continue;
end
emaxlast = emax;
% adjoint field, step backward in time, source is error
lam = zeros(Lx,Lt);
for i = (Lt:-1:2)
lam(1:Lx,i-1) = D*lam(1:Lx,i) + Dt*e(1:Lx,i);
end
% save starting u and e
if( iter==1 )
u0 = u;
e0 = e;
lam0 = lam;
end
% dE/dmi is adjoint field at t=0
dEdm = -2.0*lam(1:Lx,1);
% change in -dEdm direction
mg = mg - Dm * dEdm;
end
% plot data
figure();
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize',14);
axis( [0, maxiter, 0, 1.1*max(Emax) ] );
plot([0:maxiter-1]',Emax,'k-','LineWidth',2);
xlabel('iteration i');
ylabel('E');
fprintf('Caption: Error E as a function of iteration number i\n');
% plot images of solution and etc.
figure();
clf;
subplot(2,4,1);
hold on;
colormap('jet');
axis( [0, tmax, 0, xmax] );
imagesc(utrue);
colorbar();
xlabel('t');
ylabel('x');
title('u');
subplot(2,4,2);
colormap('jet');
axis( [0, tmax, 0, xmax] );
imagesc(u0);
colorbar();
xlabel('t');
ylabel('x');
title('u0');
subplot(2,4,3);
hold on;
colormap('jet');
axis( [0, tmax, 0, xmax] );
e0absmax = max(abs(e0(:)));
caxis( [-e0absmax, e0absmax ]);
imagesc(e0);
colorbar();
xlabel('t');
ylabel('x');
title('e0');
subplot(2,4,4);
hold on;
colormap('jet');
axis( [0, tmax, 0, xmax] );
lam0absmax = max(abs(lam0(:)));
caxis( [-lam0absmax, lam0absmax] );
imagesc(lam0);
colorbar();
xlabel('t');
ylabel('x');
title('lambda0');
subplot(2,4,5);
hold on;
colormap('jet');
axis( [0, tmax, 0, xmax] );
imagesc(utrue);
colorbar();
xlabel('t');
ylabel('x');
title('u');
subplot(2,4,6);
hold on;
colormap('jet');
axis( [0, tmax, 0, xmax] );
imagesc(u);
colorbar();
xlabel('t');
ylabel('x');
title('u pre');
subplot(2,4,7);
hold on;
colormap('jet');
axis( [0, tmax, 0, xmax] );
caxis( [-e0absmax, e0absmax ]);
imagesc(e);
colorbar();
xlabel('t');
ylabel('x');
title('e');
subplot(2,4,8);
hold on;
colormap('jet');
axis( [0, tmax, 0, xmax] );
caxis( [-lam0absmax, lam0absmax ]);
imagesc(lam);
colorbar();
xlabel('t');
ylabel('x');
title('lamda');
fprintf('Caption: Data assimilation problem. (Top row) (left) True field u,\n');
fprintf('(2nd) initial field u0, (3rd) initial error e0 (right) initial adjoint\n');
fprintf('field lambda0. (Bottom) (left) True field u, (2nd) final predicted u,\n');
fprintf('(3rd) final error e, (right) final adjoint field lambda\n');
% gdama11_11
% Supports Problem 11.6
clear all;
fprintf('gdama14_11');
% In order for this script to actually write the datafile
% you must set dowrite to 1
dowrite=0;
% known constants
c0 = 0.7;
b = 0.4;
figure();
clf;
% time setup
N=101;
Dt = 0.1;
i0 = floor(N/10);
t = Dt*[-i0+1:N-i0]';
% heterogeities
t01 = 1;
a1 = 0.1;
s1 = 0.25;
t02 = 4;
a2 = 0.05;
s2 = 0.5;
dm = a1*exp(-((t-t01).^2)/(2*s1*s1))+a2*exp(-((t-t02).^2)/(2*s2*s2));
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
mr = 0.11;
axis( [t(1), t(end), -mr, mr] );
plot( t, dm, 'k-', 'LineWidth', 3 );
xlabel('time t');
ylabel('dm');
% reference solution
u0 = (t>0) .* exp( -c0*t );
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
ur = 1.1;
axis( [t(1), t(end), -ur, ur] );
plot( t, u0, 'k-', 'LineWidth', 3 );
xlabel('time t');
ylabel('u');
% sum up point scatterers to get perturbation in solution
du = zeros(N,1);
for i=(1:N)
t0i = t(i);
dui = -dm(i) * u0(i) * (t>t0i) .* exp( -c0*(t-t0i) );
du = du + dui;
end
u = u0 + du;
plot( t, u, 'r-', 'LineWidth', 2 );
% data predicted by reference solution
d0 = b*Dt*cumsum( (t>0) .* u0 );
% data predicted by heterogenous solution
d = b*Dt*cumsum( (t>0) .* u );
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
dr = 1.1;
axis( [t(1), t(end), -dr, dr] );
plot( t, d0, 'k-', 'LineWidth', 3 );
plot( t, d, 'r-', 'LineWidth', 3 );
xlabel('time t');
ylabel('d');
fprintf('Caption: (Top) Heterogeneities in the model parsmeters. (Middle) True field u (red curve) and\n');
fprintf('reference field u0 (black curve), (bottom) true data (red curve) and reference data (black curve)\n');
% write it two places
if( dowrite )
dlmwrite('gda1106_data.txt', [t,dm,d], 'delimiter', '\t' );
dlmwrite('../data/gda1106_data.txt', [t,dm,d], 'delimiter', '\t' );
end