Live Script edama_10
This Live Script supports Chapter 10 of Environmental Data Analysis with MATLAB or PYTHON, Third Edition, by WIlliam Menke
The script naming system has the form edama_CC_SS, where CC is the chapter number and SS is the script number.
% edama_10_01: fit a high-order polynomial
clear all;
% unevenly samples times
t = [0.01, 0.11, 0.24, 0.38, 0.42, 0.57, 0.62, 0.77, 0.81, 0.91, 0.99]';
N=length(t);
% simple function
f0=2;
c=2;
d = t/2 + [0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1]';
% plot simple function
figure(1);
clf;
axis([0, 1, -5, 5]);
hold on;
set(gca,'LineWidth',2);
plot(t,d,'ko','LineWidth',2);
xlabel('t');
ylabel('d(t)');
% linear representation of polynomial relationship
G=zeros(N,N);
for i = [1:N]
G(:,i)=t.^(i-1);
end
% solve for coefficients
m = G\d;
% evaluate at equally-spaced values of t
Np=101;
tp = 0.01*[0:Np-1]';
Gp=zeros(Np,N);
for i = [1:N]
Gp(:,i)=tp.^(i-1);
end
dp=Gp*m;
% plot simple function
plot(tp,dp,'k-','LineWidth',2);
% edama_10_02: linear interpolation
% unevenly samples times
t = [0.01, 0.11, 0.24, 0.38, 0.42, 0.57, 0.62, 0.77, 0.81, 0.91, 0.99]';
N=length(t);
% simple function
f0=2;
c=2;
d = sin(2*pi*f0*t).*exp(-c*t);
d = t/2 + [0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1]';
% plot simple function
figure(1);
clf;
axis([0, 1, -5, 5]);
hold on;
set(gca,'LineWidth',2);
plot(t,d,'ko','LineWidth',2);
xlabel('t');
ylabel('d(t)');
% evaluate at equally-spaced values of t
Np=101;
tp = 0.01*[0:Np-1]';
% linear interpolation
dp=interp1(t,d,tp,'linear','extrap');
% plot simple function
plot(tp,dp,'k-','LineWidth',2);
% edama_10_03: cubic interpolation
% unevenly samples times
t = [0.01, 0.11, 0.24, 0.38, 0.42, 0.57, 0.62, 0.77, 0.81, 0.91, 0.99]';
N=length(t);
% simple function
f0=2;
c=2;
d = sin(2*pi*f0*t).*exp(-c*t);
d = t/2 + [0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1]';
% plot simple function
figure(1);
clf;
axis([0, 1, -5, 5]);
hold on;
set(gca,'LineWidth',2);
plot(t,d,'ko','LineWidth',2);
xlabel('t');
ylabel('d(t)');
% evaluate at equally-spaced values of t
Np=101;
tp = 0.01*[0:Np-1]';
% cubic interpolation
dp=interp1(t,d,tp,'spline','extrap');
% plot simple function
plot(tp,dp,'k-','LineWidth',2);
% edama_10_04: make correlated random time series
clear all;
% tunable parameters
CTYPE = 2; % oscillating exponential on 0, Gaussian on 1, Exponential on 2
MPLOT = 5; % scale of model parameter plot
CPLOT = 1; % scale of covariance plot
% standard time axis
N=256; % must be even, best if a power of 2
Dt=1.0;
tmax=Dt*(N-1);
t=Dt*[0:N-1]';
tmid=tmax/2.0;
% exact gamma*cos(w0T)exp(-sT) covariance function
n=6; % number of oscillations on (0,tmax) interal
w0 = n*2*pi/tmax; % angular frequency of oscilation
s = 20.0/tmax; % decay rate
gamma = 1.0; % sqrt variance
gamma2 = gamma^2; % variance
if( CTYPE==0 )
c = gamma2*exp(-0.5*(s*abs(t-tmid).^2));
elseif (CTYPE==1)
c = gamma2*cos(w0*(t-tmid)).*exp(-s*abs(t-tmid));
else
c = gamma2*exp(-s*abs(t-tmid));
end
fprintf('gamma: %.4f s: %.4f\n', gamma, s);
% realize a time series with a specified covariance
ctilde = fft(c); % Fourier transform of c
casd = sqrt(abs(ctilde)); % amp. spec. density of time series
d = random('Normal',0.0,1.0,N,1); % make random numbers
d = d - mean(d); % remove mean
dtilde = fft(d); % Fourier transform of d
dasd = abs(dtilde); % a.s.d. of d
dtilde = (casd./dasd).*dtilde; % adjust so has a.s.d. of c
d = real(ifft(dtilde)); % inverse Fourier transform
vard = var(d); % variance of d
d = (gamma/sqrt(vard))*d; % normalize mf to variance gamma^2
cest = xcorr(d,d)/N; % estimated covariance
% plot this covariance
figure(1);
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ -tmid, tmid, -CPLOT, CPLOT] );
xlabel('time lag t (s)');
ylabel('c(t)');
plot(t-tmid,c,'k-','LineWidth',3);
tc = Dt*[(-N+1):(N-1)]';
plot(tc,cest,'k:','LineWidth',3);
% time series figure
figure(2);
clf;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ 0, tmax, -MPLOT, MPLOT] );
xlabel('time t (s)');
ylabel('d(t)');
hold on;
plot(t,d,'k-','LineWidth',3);
% edama_10_05, Gaussian Process Regression example
% this example uses synthetic data
% and a gamma2*cos(w0T)*exp(-sT) covariance function,
% with T = |t(i)-t(j)|
clear all;
% tunable parameters
MPLOT = 5; % scale of model parameter plot
CPLOT = 1; % scale of covariance plot
gamma = 1; % sqrt(variance) of covariance
sigmad = 0.2; % noise level of observed data
wbias=1.0; % compute GPR with wrong w0 (e.g 1 and 3)
sbias=1.0; % compute GPR with wrong s
PLOTCONF=1; % plot confidence intervals on 1
USEFORRLOOP=0; % for loops on 1, outer products on 0
% axes definitons
% model parameter axis, tm, with sampling Dt
M=256; % must be even, best if a power of 2
Dtm=1.0;
tmax=Dtm*(M-1);
tm=Dtm*[0:M-1]';
% finely sampled model parameter vector, with sampling
% Dt/K, where K is a smallish power of 2, is used in
% the initialization process, as described below
K=8; % must be even, best if a power of 2
Mf=M*K;
Dtf=Dtm/K;
tf=Dtf*[0:Mf-1]';
% Part 1: realization of a finely sampled time
% series, mf, with the desired covariance
% PART 1: create true "fine" covariance function
n=6; % number of oscillations on (0,tmax) interal
w0 = n*2*pi/tmax; % angular frequency of oscilation
tmid=tmax/2.0;
s = 20.0/tmax; % decay rate
gamma2 = gamma^2;
cf = gamma2*cos(w0*(tf-tmid)).*exp(-s*abs(tf-tmid));
fprintf('gamma: %.4f w0: %.4f s: %.4f\n', gamma, w0, s);
% plot this covariance
figure(1);
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ -tmid, tmid, -CPLOT, CPLOT] );
xlabel('time t (s)');
ylabel('c(t)');
plot(tf-tmid,cf,'k-','LineWidth',3);
% PART 2: create "fine" time series with thus
% covariance using Fourier techniques
% amplitude spectral density of true time series
cftilde = fft(cf);
cfasd = sqrt(abs(cftilde));
% To create a time series with a specified covariannce,
% start with a times series of uniformly-distrbuted
% random noise, which has a uniform a.s.d., Fourier
% transorm it, multiply teh transform by the a.s.d
% of the specified covariance, trasnform back, and adjust
% the amplitude to the desired variance.
mf = random('Normal',0.0,1.0,Mf,1); % random numbers
mf = mf - mean(mf); % remove mean
mftilde = fft(mf); % Fourier transform of mf
mfasd = abs(mftilde); % a.s.d. of mf
mftilde = (cfasd./mfasd).*mftilde; % adjust so has a.s.d. of c
mf = real(ifft(mftilde)); % inverse Fourier transform
varmf = var(mf); % variance of mf
mf = (gamma/sqrt(varmf))*mf; % normalize mf to variance gamma^2
% PART 3: Regularly sub-sample mf to get mtrue, irregularly
% sub-sample mf to get true data dtrue, and add uncorrelated
% noise to dtrue to get observed data dobs
mtrue = mf(1:K:Mf);
Ntarget = 50; % fewer observed data than true data
% create list, j, of N unique, randomly chosen indices
% that are not in the sequence 1, K+1, 2K+1, ...
j=zeros(Ntarget,1);
Nj=0;
while (Nj<Ntarget)
i=unidrnd(Mf);
if( (mod(i-1,K)~=0) && (~ismember(i,j)) )
Nj=Nj+1;
j(Nj)=i;
end
end
j = unique(j); % sort for plotting purposes
j = sort(j,'ascend'); % sort for plotting purposes
N = length(j);
fprintf('target N: %d, actual N: %d\n',Ntarget, N)
td = tf(j); % times of observations
dtrue = mf(j); % irregularly sub-sample mf to get true data
% observed data is true data + uncorrelated random noise
sigmad2 = sigmad^2;
dobs = dtrue + random('Normal', 0, sigmad, N, 1 );
% the data kernel, G, for the equation d=Gm
G = [ zeros(N,M), zeros(N,N) ];
% time series figure
figure(2);
clf;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ 0, tmax, -MPLOT, MPLOT] );
xlabel('time t (s)');
ylabel('d(t) and m(t)');
hold on;
plot(tm,mtrue,'k-','LineWidth',3);
plot(td,dtrue,'ko','LineWidth',2);
% Since the true (w0, s) are in general unknown
% this code allows them to be assigned differently
% from the true values
w0est = wbias*w0; % w0est = 1.0*w0 is the true value
sest = sbias*s;
if( USEFORRLOOP )
% build covariances via for loop
Ccc = zeros(N,N); % covariance matrix, C(control,control)
for i=[1:N]
for j=[1:N]
T = abs(td(i)-td(j));
Ccc(i,j) = gamma2*cos(w0est*T) * exp(-sest*T);
end
end
Ctc = zeros(M,N); % covariance matrix, C(target,control)
for i=[1:M]
for j=[1:N]
T = abs(tm(i)-td(j));
Ctc(i,j) = gamma2*cos(w0est*T) * exp(-sest*T);
end
end
else % use outer product
% build covariances via outer product
Tdd = abs( td*ones(1,N) - ones(N,1)*(td') );
Ccc = gamma2*cos(w0est*Tdd) .* exp(-sest*Tdd);
Tmd = abs( tm*ones(1,N) - ones(M,1)*(td') );
Ctc = gamma2*cos(w0est*Tmd) .* exp(-sest*Tmd);
end
% solve GPE problem
A = (Ccc + sigmad2*eye(N,N) );
mest = Ctc * (A\dobs);
% plot mest
plot(tm,mest,'-','LineWidth',2,'Color',[0.8,0.8,0.8]);
% sqrt(variance) of mest
% plot 95% confidence intervals
if( PLOTCONF )
% as explained in the errate, the commented out lines are not correct
% B = Ctc/A;
% Cttest = sigmad2*B*B';
%
% Ctc (Ainv Ctc') = Ctc B with B = Ainv Ctc' or AB = Ctc'
Tmm = abs( tm*ones(1,M) - ones(M,1)*(tm') );
Ctt = gamma2*cos(w0est*Tmm) .* exp(-sest*Tmm);
Cttest = Ctt - Ctc*(A\(Ctc'));
sigmamest = sqrt(diag(Cttest));
plot(tm,mest-2*sigmamest,'k:','LineWidth',2);
plot(tm,mest+2*sigmamest,'k:','LineWidth',2);
end
% plot estimated covariance
figure(1);
cest = xcorr(mest,mest)/M;
figure(1);
tc = Dtm*[(-M+1):(M-1)]';
plot(tc,cest,'k:','LineWidth',3);
% edama_10_06, GPR example using bicg()
% this example uses synthetic data
% and a gamma2*cos(w0T)*exp(-sT) covariance function,
% with T = |t(i)-t(j)|
clear all;
% global variables used to pass info to CpImul()
global N
global td
global sigmad2
global gamma2
global w0est
global sest
% tunable parameters
MPLOT = 5; % scale of model parameter plot
CPLOT = 1; % scale of covariance plot
gamma = 1; % sqrt(variance) of covariance
sigmad = 0.2; % noise level of observed data
wbias=1.0; % compute GPR with wrong w0
sbias=1.0; % compute GPR with wrong s
PLOTCONF=1; % plot confidence intervals on 1
USEFORRLOOP=0; % for loops on 1, outer products on 0
% axes definitons
% model parameter axis, tm, with sampling Dt
M=256; % must be even, best if a power of 2
Dtm=1.0;
tmax=Dtm*(M-1);
tm=Dtm*[0:M-1]';
% finely sampled model parameter vector, with sampling
% Dt/K, where K is a smallish power of 2, is used in
% the initialization process, as described below
K=8; % must be even, best if a power of 2
Mf=M*K;
Dtf=Dtm/K;
tf=Dtf*[0:Mf-1]';
% Part 1: realization of a finely sampled time
% series, mf, with the desired covariance
% PART 1: create true "fine" covariance function
n=6; % number of oscillations on (0,tmax) interal
w0 = n*2*pi/tmax; % angular frequency of oscilation
tmid=tmax/2.0;
s = 20.0/tmax; % decay rate
gamma2 = gamma^2;
cf = gamma2*cos(w0*(tf-tmid)).*exp(-s*abs(tf-tmid));
fprintf('gamma: %.4f w0: %.4f s: %.4f\n', gamma, w0, s);
% plot this covariance
figure(1);
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ -tmid, tmid, -CPLOT, CPLOT] );
xlabel('time t (s)');
ylabel('c(t)');
plot(tf-tmid,cf,'k-','LineWidth',3);
% PART 2: create "fine" time series with thus
% covariance using Fourier techniques
% amplitude spectral density of true covariance
cftilde = fft(cf);
cfasd = sqrt(abs(cftilde));
% To create a time series with a specified covariannce,
% start with a times series of uniformly-distrbuted
% random noise, which has a uniform a.s.d., Fourier
% transorm it, multiply teh transform by the a.s.d
% of the specified covariance, trasnform back, and adjust
% the amplitude to the desired variance.
mf = random('Normal',0.0,1.0,Mf,1); % random numbers
mf = mf - mean(mf); % remove mean
mftilde = fft(mf); % Fourier transform of mf
mfasd = abs(mftilde); % a.s.d. of mf
mftilde = (cfasd./mfasd).*mftilde; % adjust so has a.s.d. of c
mf = real(ifft(mftilde)); % inverse Fourier transform
varmf = var(mf); % variance of mf
mf = (gamma/sqrt(varmf))*mf; % normalize mf to variance gamma^2
% PART 3: Regularly sub-sample mf to get mtrue, irregularly
% sub-sample mf to get true data dtrue, and add uncorrelated
% noise to dtrue to get observed data dobs
mtrue = mf(1:K:Mf);
Ntarget = 50; % fewer observed data than true data
% create list, j, of N unique, randomly chosen indices
% that are not in the sequence 1, K+1, 2K+1, ...
j=zeros(Ntarget,1);
Nj=0;
while (Nj<Ntarget)
i=unidrnd(Mf);
if( (mod(i-1,K)~=0) && (~ismember(i,j)) )
Nj=Nj+1;
j(Nj)=i;
end
end
j = unique(j);
j = sort(j,'ascend'); % sort for plotting purposes
N = length(j);
fprintf('target N: %d, actual N: %d\n',Ntarget,N);
td = tf(j); % times of observations
dtrue = mf(j); % irregularly sub-sample mf to get true data
% observed data is true data + uncorrelated random noise
sigmad2 = sigmad^2;
dobs = dtrue + random('Normal', 0, sigmad, N, 1 );
% time series figure
figure(2);
clf;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ 0, tmax, -MPLOT, MPLOT] );
xlabel('time t (s)');
ylabel('d(t) and m(t)');
hold on;
plot(tm,mtrue,'k-','LineWidth',3);
plot(td,dtrue,'ko','LineWidth',2);
% Since the true (w0, s) are in general unknown
% this code allows them to be assigned differently
% from the true values
w0est = wbias*w0; % w0est = 1.0*w0 is the true value
sest = sbias*s;
% solve GPE problem
% solve using conjugate gradient algorithm
u=bicg(@CpImul,dobs,1e-12,3*(M+N));
mest = zeros(M,1);
for i=[1:M]
mest(i)=0.0;
for j=[1:N]
Tabs = abs(tm(i)-td(j));
Cij = gamma2*cos(w0est*Tabs)*exp(-sest*Tabs);
mest(i)=mest(i)+Cij*u(j);
end
end
dpre = zeros(N,1);
for i=[1:N]
dpre(i)=0.0;
for j=[1:N]
Tabs = abs(td(i)-td(j));
Cij = gamma2*cos(w0est*Tabs)*exp(-sest*Tabs);
dpre(i)=dpre(i)+Cij*u(j);
end
end
% plot mest
plot(tm,mest,'-','LineWidth',2,'Color',[0.8,0.8,0.8]);
% plot estimated covariance
figure(1);
cest = xcorr(mest,mest)/M;
figure(1);
tc = Dtm*[(-M+1):(M-1)]';
plot(tc,cest,'k:','LineWidth',3);
% edama_10_07, GPR Tuning Example
% this script tunes angular frequency, q, in the
% covariance c(t) = gamma^2 cos(qt) exp(-st)
clear all;
% tunable parameters
MPLOT = 5; % scale of model parameter plot
CPLOT = 1; % scale of covariance plot
gamma = 1; % sqrt(variance) of covariance
sigmad = 0.2; % noise level of observed data
wbias=3.0; % compute GPR with wrong w0
sbias=1.0; % compute GPR with wrong s
% axes definitons
% model parameter axis, tm, with sampling Dt
M=512; % must be even, best if a power of 2
Dtm=1.0;
tmax=Dtm*(M-1);
tm=Dtm*[0:M-1]';
% finely sampled model parameter vector, with sampling
% Dt/K, where K is a smallish power of 2, is used in
% the initialization process, as described below
K=8; % must be even, best if a power of 2
Mf=M*K;
Dtf=Dtm/K;
tf=Dtf*[0:Mf-1]';
% Part 1: realization of a finely sampled time
% series, mf, with the desired covariance
% PART 1: create true "fine" covariance function
n=6; % number of oscillations on (0,tmax) interal
w0 = n*2*pi/tmax; % angular frequency of oscilation
tmid=tmax/2.0; % idpoint
s = Dtm/60.0; % decay rate
gamma2 = gamma^2; % variance
cf = gamma2*cos(w0*(tf-tmid)).*exp(-s*abs(tf-tmid));
% plot this covariance
figure(1);
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ -tmid, tmid, -CPLOT, CPLOT] );
xlabel('time t (s)');
ylabel('c(t)');
plot(tf-tmid,cf,'k-','LineWidth',3);
% PART 2: create "fine" time series with thus
% covariance using Fourier techniques
% amplitude spectral density of true time series
cftilde = fft(cf);
cfasd = sqrt(abs(cftilde));
% To create a time series with a specified covariannce,
% start with a times series of uniformly-distrbuted
% random noise, which has a uniform a.s.d., Fourier
% transorm it, multiply teh transform by the a.s.d
% of the specified covariance, trasnform back, and adjust
% the amplitude to the desired variance.
mf = random('Normal',0.0,1.0,Mf,1); % random numbers
mf = mf - mean(mf); % remove mean
mftilde = fft(mf); % Fourier transform of mf
mfasd = abs(mftilde); % a.s.d. of mf
mftilde = (cfasd./mfasd).*mftilde; % adjust so has a.s.d. of c
mf = real(ifft(mftilde)); % inverse Fourier transform
varmf = var(mf); % variance of mf
mf = (gamma/sqrt(varmf))*mf; % normalize mf to variance gamma^2
% PART 3: Regularly sub-sample mf to get mtrue, irregularly
% sub-sample mf to get true data dtrue, and add uncorrelated
% noise to dtrue to get observed data dobs
mtrue = mf(1:K:Mf);
Ntarget = 150; % fewer observed data than true data
j=zeros(Ntarget,1);
Nj=0;
% create list, j, of N unique, randomly chosen indices
% that are not in the sequence 1, K+1, 2K+1, ...
while (Nj<Ntarget)
i=unidrnd(Mf);
if( (mod(i-1,K)~=0) && (~ismember(i,j)) )
Nj=Nj+1;
j(Nj)=i;
end
end
j = unique(j);
j = sort(j,'ascend'); % sort for plotting purposes
N = length(j);
fprintf('target N: %d, actual N:%d\n',Ntarget,N);
td = tf(j); % times of observations
dtrue = mf(j); % irregularly sub-sample mf to get true data
% observed data is true data + uncorrelated random noise
sigmad2 = sigmad^2;
sigmadr = 1.0/sigmad;
dobs = dtrue + random('Normal', 0, sigmad, N, 1 );
% the data kernel, G, for the equation d=Gm
G = [ zeros(N,M), zeros(N,N) ];
% time series figure
figure(2);
clf;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ 0, tmax, -MPLOT, MPLOT] );
xlabel('time t (s)');
ylabel('d(t) and m(t)');
hold on;
plot(tm,mtrue,'k-','LineWidth',3);
plot(td,dtrue,'ko','LineWidth',2);
% Part 4: GPR solution for inital parameters
% (for comparison purposes)
% Since the true (w0, s) are in general unknown
% this code allows them to be assigned differently
% from the true values
w0est = wbias*w0; % w0est = 1.0*w0 is the true value
sest = sbias*s;
% build covariances via outer product
Tdd = abs( td*ones(1,N) - ones(N,1)*(td') );
Ccc = gamma2*cos(w0est*Tdd) .* exp(-sest*Tdd);
Tmd = abs( tm*ones(1,N) - ones(M,1)*(td') );
Ctc = gamma2*cos(w0est*Tmd) .* exp(-sest*Tmd);
% solve GPE problem
A = (Ccc + sigmad2*eye(N,N) );
mest = Ctc * (A\dobs);
% plot mest
plot(tm,mest,'-','LineWidth',2,'Color',[0.8,0.8,0.8]);
% plot estimated covariance
figure(1);
cest = xcorr(mest,mest)/M;
figure(1);
tc = Dtm*[(-M+1):(M-1)]';
plot(tc,cest,'k:','LineWidth',3);
% PART 5: Set up for tuning
% covariance plot
figure(3);
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ -tmid, tmid, -CPLOT, CPLOT] );
xlabel('time t (s)');
ylabel('c(t)');
plot(tf-tmid,cf,'k-','LineWidth',3);
% time series plot
figure(4);
clf;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ 0, tmax, -MPLOT, MPLOT] );
xlabel('time t (s)');
ylabel('d(t) and m(t)');
hold on;
plot(tm,mtrue,'k-','LineWidth',3);
plot(td,dtrue,'ko','LineWidth',2);
% Part 6: Tuning
% Gradient Descent Method adjustment of w0, now called q
qtrue=w0; % true parameter, q
qg = w0est; % inital guess of parameter, q
Dq = 0.025; % initial step size for gradiet descent
NITER = 30; % maximum number of iterations
qg = qg-Dq; % so first iteration starts qg
alpha = 1.0; % step size reducer
nu = 1.0; % sign of gradient
mtpri = zeros(M,1); % prior model parameters are zero
mcpri = zeros(N,1); % prior model parameters are zero
mpri = [mtpri;mcpri]; % prior model parameters are zero
tall = [tm; td]; % concatenated time vectors
qsave = zeros(NITER,1); % save estimates, for plot
for iter = [1:NITER] % loop over descent increments
qg = qg+nu*alpha*Dq; % new guess of parameter, q
% full covariance matrix and its derivative
Tall = tall*ones(1,M+N) - ones(M+N,1)*(tall');
Call = gamma2 * cos(qg*Tall) .* exp(-sest*abs(Tall));
dCalldq = -gamma2 * Tall .* sin(qg*Tall) .* exp(-sest*abs(Tall));
% break out four quandrants
Ctt = Call(1:M,1:M);
dCttdq = dCalldq(1:M,1:M);
Ctc = Call(1:M,M+1:M+N);
dCtcdq = dCalldq(1:M,M+1:M+N);
Cct = Ctc';
Ccc = Call(M+1:M+N,M+1:M+N);
dCccdq = dCalldq(M+1:M+N,M+1:M+N);
% A = Ccc + sigmad^2 I
A = Ccc + sigmad2*eye(N);
dAdq = dCccdq;
% various functions of the covariance matrix
CallI = inv(Call);
CallIsr = sqrtm(CallI);
dCallIdq = - CallI * dCalldq * CallI;
dCallIsrdq = sylvester( CallIsr, CallIsr, dCallIdq );
R = chol(Call);
logdetCall = 2.0*sum(log(diag(R)));
dlogdetCalldq = trace( CallI*dCalldq );
% GPR solution
u = A\dobs;
mtg = Ctc*u;
mcg = Ccc*u;
mg = [mtg;mcg];
dg = mg(M+1:M+N);
% derivative of predicted data
v = A\(dAdq*u);
dmestdq = [dCtcdq*u - Ctc*v; dCccdq*u - Ccc*v];
ddpredq = dmestdq(M+1:N+M);
detildedq = -sigmadr*ddpredq;
e = dobs-dg;
etilde = sigmadr*e;
% derivative of total data error
E = etilde'*etilde;
dEdp = 2*etilde'*detildedq;
% prior error and its derivative
el=(mpri-mg);
eltilde = CallIsr*el;
deltildedq = dCallIsrdq*el - CallIsr*dmestdq;
% total prior error and its derivative
L = eltilde'*eltilde;
dLdp = 2*eltilde'*deltildedq;
% function, PSI being minimized
PSI = logdetCall + E + L;
% and its deriative
dPSIdp = dlogdetCalldq + dEdp + dLdp;
nu = -sign(dPSIdp); % downhill vector
% gradient method method
if( iter==1 ) % on the first iteration do nothing but save
qglast = qg;
PSIlast = PSI;
fprintf('1: qg %e nu %e alpha %e PSI %e PSIlast %e\n', qg, nu, alpha, PSI, PSIlast);
elseif (PSI < PSIlast )
qglast = qg;
PSIlast = PSI;
fprintf('2: qg %e nu %e alpha %e PSI %e PSilast %e\n', qg, nu, alpha, PSI, PSIlast);
else % psi increased, so decrease alpha
qg = qglast;
PSI = PSIlast;
alpha = alpha/2;
fprintf('2: qg %e nu %e alpha %e PSI %e PSIlast %e\n', qg, nu, alpha, PSI, PSIlast);
end
qsave(iter) = qg;
end
pest = qg;
mest = mg;
mtest = mest(1:M);
mcest = mest(M+1:M+N);
dpre = mcest;
% plot estimated covariance
figure(3);
cest = xcorr(mtest,mtest)/M;
tc = Dtm*[(-M+1):(M-1)]';
plot(tc,cest,'k:','LineWidth',3);
% plot estimated time series
figure(4);
plot( tm, mtest, '-', 'LineWidth', 2, 'Color', [0.8,0.8,0.8]);
fprintf('qtrue %e qest %e relating error %e\n', qtrue, pest, (qtrue-pest)/qtrue );
figure(5);
clf;
hold on;
set(gca,'LineWidth',2);
set(gca,'FontSize', 14);
axis( [ 0, NITER, 0, 0.3] );
xlabel('iteration');
ylabel('w0');
plot([1:NITER]',qsave,'k-','LineWidth',2);
plot([1:NITER]',qsave,'ko','LineWidth',2);
plot([0:NITER]',w0*ones(NITER+1,1),'k:','LineWidth',2);
% edama_10_08: 2D interpolation of pressure data
% this version used griddata()
% this script
% creates simulated data (Figue 1)
% creates Delaunay triangle mesh (Figure 2)
% finds triangle enclosong an arbitrary point (Figure 3)
% linear interpolation (Figure 4)
% cubic spline interpolation
%
% Note: this script used the griddata() function, which the MatLab
% documentation says is no longer recommended. See eda08_05 for an
% alternative.
% standard set-up of x-axis, rows of grid are x
I = 41;
Dx = 1;
xmin = 0;
xmax = 40;
x = xmin + (xmax-xmin)*[0:I-1]'/(I-1);
% standard set-up of y-axis, columns of grid are y
J = 41;
Dy = 1;
ymin = 0;
ymax = 40;
y = ymin + (ymax-ymin)*[0:J-1]'/(J-1);
% grid 10% populated with data
N0 = floor(I*J/10);
rowindex0=unidrnd(I,N0,1);
colindex0=unidrnd(J,N0,1);
% remove duplicate points, because the interpolation
% function griddata() doesn't like them, by checking
% eack point against a matrix of ones and zeros
A = zeros(I,J);
rowindex=zeros(N0,1);
colindex=zeros(N0,1);
N=0;
for i = [1:N0]
if ( A(rowindex0(i),colindex0(i)) == 0 )
N=N+1;
rowindex(N)=rowindex0(i);
colindex(N)=colindex0(i);
end
A(rowindex0(i),colindex0(i)) = 1;
end
rowindex=rowindex(1:N);
colindex=colindex(1:N);
xobs = x(rowindex);
yobs = y(colindex);
N=length(xobs);
% create simulated pressure data (solution to Laplace's equation)
kappa = 0.1;
dtrue = 10*sin(kappa*xobs).*exp(-kappa*yobs);
sigmad = 0.1;
dobs = dtrue + random('Normal',0.0,sigmad,N,1);
% populate matrix A with observations
A = zeros(I,J);
for i = [1:N]
A(rowindex(i),colindex(i)) = dobs(i);
end
% plot data
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(A))-min(min(A));
if( range==0 )
range=1;
end
imagesc( [0, I], [0, J], (A-min(min(A)))/range);
hold on;
xlabel('y');
ylabel('x');
title('data');
% plot data and triangulation
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(A))-min(min(A));
if( range==0 )
range=1;
end
imagesc( [0, I], [0, J], (A-min(min(A)))/range);
hold on;
xlabel('y');
ylabel('x');
title('data + triangulation');
% Delaunay triangulation
mytri = DelaunayTri(xobs,yobs);
% XY is a NXY by 2 list of (x,y) positons of vertices
XY=mytri.X;
[NXY, i] = size(XY);
% TRI is an NTRI by 3 integer list of which three vertices
% are associated with bogen triangle
TRI=mytri.Triangulation;
[NTRI, i] = size(TRI);
for i = [1:NTRI]
v1=TRI(i,1); v2=TRI(i,2); v3=TRI(i,3);
plot( [XY(v1,2), XY(v2,2)], [XY(v1,1), XY(v2,1)], 'k-', 'LineWidth', 2);
plot( [XY(v2,2), XY(v3,2)], [XY(v2,1), XY(v3,1)], 'k-', 'LineWidth', 2);
plot( [XY(v3,2), XY(v1,2)], [XY(v3,1), XY(v1,1)], 'k-', 'LineWidth', 2);
end
% find triangle enclosing point x0
x0=[20,30];
tri0 = pointLocation(mytri, x0);
plot( x0(2), x0(1), 'ko', 'LineWidth', 2);
if( ~isnan(tri0) )
v1=TRI(tri0,1); v2=TRI(tri0,2); v3=TRI(tri0,3);
plot( [XY(v1,2), XY(v2,2)], [XY(v1,1), XY(v2,1)], 'k-', 'LineWidth', 3);
plot( [XY(v2,2), XY(v3,2)], [XY(v2,1), XY(v3,1)], 'k-', 'LineWidth', 3);
plot( [XY(v3,2), XY(v1,2)], [XY(v3,1), XY(v1,1)], 'k-', 'LineWidth', 3);
end
% grid of (x,y) values
XGRID = x * ones(1,J);
YGRID = ones(I,1) * y';
% linear interpolation
B=griddata(xobs,yobs,dobs,XGRID,YGRID,'linear');
% plot interpolated data
figure(3);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(B))-min(min(B));
if( range==0 )
range=1;
end
imagesc( [0, I], [0, J], (B-min(min(B)))/range);
hold on;
xlabel('y');
ylabel('x');
title('linear interpolation');
% cubic interpolation
C=griddata(xobs,yobs,dobs,XGRID,YGRID,'cubic');
% plot interpolated data
figure(4);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(C))-min(min(C));
if( range==0 )
range=1;
end
imagesc( [0, I], [0, J], (C-min(min(C)))/range);
hold on;
xlabel('y');
ylabel('x');
title('cubic interpolation');
% edama_10_09: 2D interpolation, alternate method
% this version uses simulated pressure data
% this script
% creates simulated data (Figue 1)
% creates Delaunay triangle mesh (Figure 2)
% finds triangle enclosing an arbitrary point (Figure 3)
% linear interpolation (Figure 4)
% natural spline interpolation
% Note: this script used the TriScatteredInterp() function, which the MatLab
% documentation says is now recommended in favor of the older, griddata()
% function implemented in eda08_04. Note, however, that TriScatteredInterp() does
% not provide cubic splines interpolation, but instead "natural neighbor
% interpolation" (which is smoother than linear interpoaltion, but which not
% as smooth as cubic interpolation). See the MatLab documentation for further details.
% rows of grid are x
I = 41;
Dx = 1;
xmin = 0;
xmax = 40;
x = xmin + (xmax-xmin)*[0:I-1]'/(I-1);
% columns of grid are y
J = 41;
Dy = 1;
ymin = 0;
ymax = 40;
y = ymin + (ymax-ymin)*[0:J-1]'/(J-1);
% grid 10% populated with data
N0 = floor(I*J/10);
rowindex0=unidrnd(I,N0,1);
colindex0=unidrnd(J,N0,1);
% remove duplicate points, because the interpolation
% function griddata() doesn't like them, by checking
% eack point against a matrix of ones and zeros
A = zeros(I,J);
rowindex=zeros(N0,1);
colindex=zeros(N0,1);
N=0;
for i = [1:N0]
if ( A(rowindex0(i),colindex0(i)) == 0 )
N=N+1;
rowindex(N)=rowindex0(i);
colindex(N)=colindex0(i);
end
A(rowindex0(i),colindex0(i)) = 1;
end
rowindex=rowindex(1:N);
colindex=colindex(1:N);
xobs = x(rowindex);
yobs = y(colindex);
N=length(xobs);
% create simulated pressure data (solution to Laplace's equation)
kappa = 0.1;
dtrue = 10*sin(kappa*xobs).*exp(-kappa*yobs);
sigmad = 0.1;
dobs = dtrue + random('Normal',0.0,sigmad,N,1);
% populate matrix A with observations
A = zeros(I,J);
for i = [1:N]
A(rowindex(i),colindex(i)) = dobs(i);
end
% plot data
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(A))-min(min(A));
if( range==0 )
range=1;
end
imagesc( [0, I], [0, J], (A-min(min(A)))/range);
hold on;
xlabel('y');
ylabel('x');
title('data');
% plot data and triangulation
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(A))-min(min(A));
if( range==0 )
range=1;
end
imagesc( [0, I], [0, J], (A-min(min(A)))/range);
hold on;
xlabel('y');
ylabel('x');
title('data + triangulation');
% Delaunay triangulation
mytri = DelaunayTri(xobs,yobs);
% XY is a NXY by 2 list of (x,y) positons of vertices
XY=mytri.X;
[NXY, i] = size(XY);
% TRI is an NTRI by 3 integer list of which three vertices
% are associated with bogen triangle
TRI=mytri.Triangulation;
[NTRI, i] = size(TRI);
for i = [1:NTRI]
v1=TRI(i,1); v2=TRI(i,2); v3=TRI(i,3);
plot( [XY(v1,2), XY(v2,2)], [XY(v1,1), XY(v2,1)], 'k-', 'LineWidth', 2);
plot( [XY(v2,2), XY(v3,2)], [XY(v2,1), XY(v3,1)], 'k-', 'LineWidth', 2);
plot( [XY(v3,2), XY(v1,2)], [XY(v3,1), XY(v1,1)], 'k-', 'LineWidth', 2);
end
% find triangle encl;osing exemplary point and plot them
x0=[20,30];
tri0 = pointLocation(mytri, x0);
plot( x0(2), x0(1), 'ko', 'LineWidth', 2);
if( ~isnan(tri0) )
v1=TRI(tri0,1); v2=TRI(tri0,2); v3=TRI(tri0,3);
plot( [XY(v1,2), XY(v2,2)], [XY(v1,1), XY(v2,1)], 'k-', 'LineWidth', 3);
plot( [XY(v2,2), XY(v3,2)], [XY(v2,1), XY(v3,1)], 'k-', 'LineWidth', 3);
plot( [XY(v3,2), XY(v1,2)], [XY(v3,1), XY(v1,1)], 'k-', 'LineWidth', 3);
end
% grid of (x,y) values
XGRID = x * ones(1,J);
YGRID = ones(I,1) * y';
% linear interpolation
F=TriScatteredInterp(xobs,yobs,dobs,'linear');
B = F(XGRID,YGRID);
% plot interpolated data
figure(3);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(B))-min(min(B));
if( range==0 )
range=1;
end
imagesc( [0, I], [0, J], (B-min(min(B)))/range);
hold on;
xlabel('y');
ylabel('x');
title('linear interpolation');
% cubic interpolation
F=TriScatteredInterp(xobs,yobs,dobs,'natural');
C = F(XGRID,YGRID);
% plot interpolated data
figure(4);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(C))-min(min(C));
if( range==0 )
range=1;
end
imagesc( [0, I], [0, J], (C-min(min(C)))/range);
hold on;
xlabel('y');
ylabel('x');
title('cubic interpolation');
% edam_10_10: two-dimensional Fourier transform
% Part 1, FFT of cosine wave, reordered so zero us in center
% Part 2, Illustration of use of symmetry
% standard set-up of wavenumber kx-axis
Nx=32;
Nkx=Nx/2+1;
Dx=1;
x=Dx*[0:Nx-1]';
xmax = Dx*(Nx-1);
kxmax=2*pi/(2*Dx);
Dkx=kxmax/(Nx/2);
kx=Dkx*[0:Nx/2,-Nx/2+1:-1]';
% standard set-up of wavenumber ky-axis
Ny=32;
Nky=Ny/2+1;
Dy=1;
y=Dy*[0:Ny-1]';
ymax = Dy*(Ny-1);
kymax=2*pi/(2*Dy);
Dky=kymax/(Ny/2);
ky=Dky*[0:Ny/2,-Ny/2+1:-1]';
% Part 1: cosine wave
theta = (pi/180)*30;
t1 = [cos(theta), sin(theta)]';;
kr1 = kxmax/4;
F = zeros(Nx,Ny);
for n = [1:Nx]
for m = [1:Ny]
F(n,m) = cos( kr1*t1(1)*x(n) + kr1*t1(2)*y(m) );
end
end
% amplitude spectral density
Ftt = fft2(F);
Fasd = ((Dx^2)/xmax) * ((Dy^2)/ymax) * sqrt((abs(Ftt).^2));
% put into more reasonable order, with zero in the center
Fasd2=zeros(Nx,Ny);
Fasd2(Nx-Nkx+1:Nx,Ny-Nky+1:Ny) = Fasd(1:Nkx,1:Nky);
Fasd2(1:Nkx-2,Ny-Nky+1:Ny) = Fasd(Nkx+1:Nx,1:Nky);
Fasd2(1:Nkx-2,1:Nky-2) = Fasd(Nkx+1:Nx,Nky+1:Ny);
Fasd2(Nkx-1:Nx,1:Nky-2) = Fasd(1:Nkx,Nky+1:Ny);
eda_draw(' ', F, 'caption F',' ',' ',' ', Fasd, 'caption Ftt',' ',' ',' ', Fasd2, 'caption Ftt' );
% Part 2: example of the use of symmetry
% new function F is noise, which has a complicated Fourier Transform
F=random('Normal',0,1,Nx,Ny);
Ftt = fft2(F);
% for real data, use symmetries to build right hand columns of
% 2D Fourier Transform from left hand columns
% First copy left Ny/2+1 columns; they are treated as known
Ftt2(:,1:Ny/2+1)=Ftt(:,1:Ny/2+1);
% then rebuild right from left columns
for m = [2:Ny/2]
mp=Ny-m+2;
Ftt2(1,mp) = conj(Ftt2(1,m));
Ftt2(2:Nx,mp) = flipud(conj(Ftt2(2:Nx,m)));
end
% invert back to F
F2 = ifft2(Ftt2);
% compute error, as a check
Ett=sum(sum(abs(Ftt2-Ftt)));
E = sum(sum(abs(F2-F)));
fprintf('Symmetry test:\n');
fprintf(' error in transform %.2f\n', Ett);
fprintf(' error in function %.2f\n', E);
eda_draw(' ', F, 'caption F', ' ', F2, 'caption F2');
% edam_10_11: Filling in gaps with fft with POCS
clear all;
% Part 1, generate true synthetic data set
% with anisotropic Gaussian autocorrelation function
% that simulates a crudely-layerd medium
% standard set-up of wavenumber kx-axis
Nx=128;
Nkx=Nx/2+1;
Dx=1;
x=Dx*[0:Nx-1]';
xmax = Dx*(Nx-1);
kxmax=2*pi/(2*Dx);
Dkx=kxmax/(Nx/2);
kx=Dkx*[0:Nx/2,-Nx/2+1:-1]';
% standard set-up of wavenumber ky-axis
Ny=128;
Nky=Ny/2+1;
Dy=1;
y=Dy*[0:Ny-1]';
ymax = Dy*(Ny-1);
kymax=2*pi/(2*Dy);
Dky=kymax/(Ny/2);
ky=Dky*[0:Ny/2,-Ny/2+1:-1]';
% wavenumber variances
% they scale a reciprocal of spatial variances
% values are chosen to make quasi-layers
sx = 5.0*Dkx;
sy = 1.0*Dky;
sx2 = sx^2;
sy2 = sy^2;
asd = zeros(Nx,Ny); % target amplitude spectral density
for ikx=[1:Nx]
for iky=[1:Ny]
psd = exp( -0.5 * (kx(ikx)^2) / sx2 - 0.5 * (ky(iky)^2) / sy2);
asd(ikx,iky) = sqrt(psd);
end
end
% start with random data and
mtrue = random('Normal',0.0,1.0,Nx,Ny);
Ftt = fft2(mtrue); % take its fft
Ftt = asd.*Ftt; % shape to target
mtrue = real(ifft2(Ftt)); % inverse fft
sigmamtrue = sqrt(var(mtrue(:))); % true variance
% normalization ratio, gives target asd the same
% variance as the fft of the data
R=sqrt(var(abs(Ftt(:)))) / sqrt(var(abs(asd(:))));
Rasd=R*asd;
% sample mtrue along vertical profiles
j = floor( [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 ]'*Ny/8.0 );
Nd = length(j);
d = zeros(Nx,Nd);
for i=[1:Nd]
d(:,i)=mtrue(:,j(i));
end
fprintf('profiles through true model\n');
eda_draw(d(:,1),d(:,2),d(:,3),d(:,4),d(:,5),d(:,6),d(:,7));
% Part 2, POCS method to fill in data gaps
% make starting model that has the right a.s.d. and variance
mest = random('Normal',0.0,1.0,Nx,Ny); % random data
Ftt = fft2(mest); % take fft
Ftt=asd.*Ftt; % shape to target
mest = real(ifft2(Ftt)); % take inverse fft
sigmamest = sqrt(var(mest(:))); % estimated variance
mest = (sigmamtrue/sigmamest)*mest; % scale to right variance
mstarting = mest; % save starting model
% POCS algorithm
Niter=50;
E=zeros(Niter,1);
for iter = [1:Niter] % loop over iterations
% enforce data are satisfied
for i=[1:Nd]
mest(:,j(i))=d(:,i);
end
% enforce a.s.d. is not greater than target
Ftt = fft2(mest);
for ikx=[1:Nx] % loop over wavenumbers
for iky=[1:Ny]
a=Rasd(ikx,iky);
b=abs(Ftt(ikx,iky));
% if asd too big, make it smaller, w/o altering phase
if( b>a )
Ftt(ikx,iky)=(a/b)*Ftt(ikx,iky);
end
end
end
mest = real(ifft2(Ftt));
% model error
e = mtrue(:)-mest(:);
E(iter) = e'*e;
end
sigmamtrue=sqrt(var(mtrue(:)));
sigmamstarting=sqrt(var(mstarting(:)));
sigmamest=sqrt(var(mest(:)));
fprintf('sqrt variances: starting %e, estimated %e, true %e\n', sigmamtrue, sigmamstarting, sigmamest);
fprintf('starting, estimated, true model\n');
eda_draw(mstarting, 'caption starting', mest, 'caption est',mtrue,'caption true');
figure(2);
clf;
set(gca,'LineWIdth',2);
set(gca,'FontSize',14);
axis([0, Niter, 0, E(1)]);
plot([1:Niter]',E,'k-', 'LineWidth',2);
xlabel('iteration, i');
ylabel('model error, E')
