In [2]:
# edapy_10_00 clear all variables and import vatious modules
# History
# 2022/06/26 - checked with most recent Python & etc. software
# 2023/07/13 - fixed some issues, incl colormap, loglof plt, eda_cvec()
# 2023/10/04 - fixed issue with GPR covariance
# 2024/05/26 - Patches to accomodate new verions of numpy, scipy, matplotlib
# patched dot product to return scalar value, and similar issues
# associated with depreciation of equivalence of 1x1 matrix and scalar
# las.bicg keyword tol changes to rtol
# las.bicg output cast to float (not sure this is realy necessary)
# corrected error concerning GPE covariance in edapy_10_05, see errata
%reset -f
import os
from datetime import date
from math import exp, pi, sin, cos, tan, sqrt, floor, ceil, log
import numpy as np
import scipy.sparse.linalg as las
import scipy.interpolate as ip
import scipy.spatial as sp
from scipy import sparse
import scipy.linalg as la
import scipy.signal as sg
import matplotlib
from matplotlib import pyplot as plt
from matplotlib.colors import ListedColormap
# eda_draw function makes a "pictorial matrix equation"
# arguments are vectors, matrices and strings
# which are plotted in the order that the appear
# except that strings starting with 'title ' are plotted
# under the subseqeunt matrix or vector
# always returns a status of 1
def eda_draw(*argv):
bw = np.zeros((256,4));
v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;
bw[:,0] = v;
bw[:,1] = v;
bw[:,2] = v;
bw[:,3] = np.ones(256);
bwcmap = ListedColormap(bw);
# size of plot
W = 16;
H = 4;
fig1 = plt.figure(1);
# figsize width and height in inches
fig1.set_size_inches(W,H);
ax1 = plt.subplot(1,1,1);
plt.axis([0, W, -H/2, H/2]);
plt.axis('off');
LM = W/6; # matrix width and heoght
LV = W/40; # vector width
FS = 0.12; # character width
TO = 0.4; # title vertical offset
SP = 0.2; # space between objects
LS = 0.2; # leading space
p = LS; # starting x-position
istitle=0; # flags presence of a title
for a in argv:
if isinstance(a,np.ndarray):
sh = np.shape(a);
if len(sh) == 1: # conversion to nx1 array
n = sh[0];
m = 1;
ap = a;
a = np.zeros((n,1));
a[:,0] = ap;
else:
n = sh[0];
m = sh[1];
if m==1:
pold=p;
left=p;
right=p+LV;
bottom=-LM/2;
top=LM/2;
plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );
p = p+LV;
pm = (p+pold)/2;
if istitle:
plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');
istitle=0;
p = p+SP;
else:
pold=p;
left=p;
right=p+LM;
bottom=-LM/2;
top=LM/2;
plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );
p = p+LM;
pm = (p+pold)/2;
if istitle:
plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');
istitle=0;
p = p+SP;
elif isinstance(a,str):
ns = len(a);
istitle=0;
if( ns>=6 ):
if 'title ' in a[0:6]:
istitle=1;
titlestr=a[6:];
if( istitle != 1):
plt.text(p,0,a);
p = p + ns*FS + SP;
plt.show();
return 1;
def FTFmul(v):
# this function is used by the bicongugate gradient solver to solve the geneneralized least
# squares problem Fm=f. Note that "F" must be called "edaFsparse".
global edaFsparse;
s = np.shape(v);
if(len(s)==1):
vv = np.zeros((s[0],1));
vv[:,0] = v;
else:
vv=v;
temp = edaFsparse*vv;
return edaFsparse.transpose()*temp;
def GLSFilterMul(v):
# this function is used by the bicongugate gradient solver to solve the
# geneneralized least squares problem Fm=f with F=[G;H] and f=[d,h],
# where G is a topplitz matrix built from the filter g
# Note that "H" must be sparese and called "edaHsparse" and that
# the filter "g" bust be called edafilterg and must be a column vector
global edaHsparse;
global edafilterg;
# rearrange as column-vector
s = np.shape(v);
M = s[0];
if(len(s)==1):
vv = np.zeros((s[0],1));
vv[:,0] = v;
else:
vv=v;
N, i = np.shape(edafilterg);
# implements (GT G + HT H) v
temp1 = np.zeros((N+M-1,1));
temp1[:,0] = np.convolve(edafilterg.ravel(),vv.ravel()); # G v is of length N
a = np.zeros((N,1));
a[:,0] = temp1[0:N,0];
b = edaHsparse * v; # H v is of length K
temp2 = np.zeros((N+N-1,1));
temp2[:,0] = np.convolve((np.flipud(edafilterg)).ravel(),a.ravel()); # GT (G v) is of length M
a2 = temp2[N-1:N+M-1,0];
b2 = edaHsparse.transpose()*b; # HT (H v) is of length M
# FT F v = GT G v + HT H v
return (a2+b2);
def GLSautoMul(v):
# this function is used by the bicongugate gradient solver to solve the
# geneneralized least squares problem Fm=f with F=[G;H] and f=[d,h],
# where GTG is built from the autocorrelation of the filter g
global edaHsparse;
global edaautoc;
# rearrange as column-vector
s = np.shape(v);
M = s[0];
if(len(s)==1):
vv = np.zeros((s[0],1));
vv[:,0] = v;
else:
vv=v;
# FT F v = GT G v + HT H v
GTGv=np.zeros((M,1));
for i in range(M):
u = np.concatenate( (np.flipud(edaautoc[0:i+1,0]), edaautoc[1:M-i,0]), axis=0 );
GTGv[i,0] = np.matmul( u.T, vv );
Hv = edaHsparse*vv;
HTHv = edaHsparse.transpose()*Hv;
return (GTGv + HTHv);
def chebyshevfilt(d, Dt, flow, fhigh):
#chebyshevfilt
# chebyshev IIR bandpass filter
# d - input array of data
# Dt - sampling interval
# flow - low pass frequency, Hz
# fhigh - high pass frequency, Hz
# dout - output array of data
# u - the numerator filter
# v - the denominator filter
# these filters can be used again using dout=filter(u,v,din);
# make sure input timeseries is a column vector
s = np.shape(d);
N = s[0];
if(N==1):
dd = np.zeros((N,1));
dd[:,0] = d;
else:
dd=d;
# sampling rate
rate=1/Dt;
# ripple parameter, set to ten percent
ripple=0.1;
# normalise frequency
fl=2.0*flow/rate;
fh=2.0*fhigh/rate;
# center frequency
cf = 4 * tan( (fl*pi/2) ) * tan( (fh*pi/2) );
# bandwidth
bw = 2 * ( tan( (fh*pi/2) ) - tan( (fl*pi/2) ) );
# ripple parameter factor
rpf = sqrt((sqrt((1.0+1.0/(ripple*ripple))) + 1.0/ripple));
a = 0.5*(rpf-1.0/rpf);
b = 0.5*(rpf+1.0/rpf);
u=np.zeros((5,1));
v=np.zeros((5,1));
theta = 3*pi/4;
sr = a * cos(theta);
si = b * sin(theta);
es = sqrt(sr*sr+si*si);
tmp= 16.0 - 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw - 4.0*bw*cf*sr + cf*cf;
v[0,0] = 1.0;
v[1,0] = 4.0*(-16.0 + 8.0*bw*sr - 2.0*bw*cf*sr + cf*cf)/tmp;
v[2,0] = (96.0 - 16.0*cf - 8.0*es*es*bw*bw + 6.0*cf*cf)/tmp;
v[3,0] = (-64.0 - 32.0*bw*sr + 8.0*bw*cf*sr + 4.0*cf*cf)/tmp;
v[4,0] = (16.0 + 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw + 4.0*bw*cf*sr + cf*cf)/tmp;
tmp = 4.0*es*es*bw*bw/tmp;
u[0,0] = tmp;
u[1,0] = 0.0;
u[2,0] = -2.0*tmp;
u[3,0] = 0.0;
u[4,0] = tmp;
dout = sg.lfilter(u.ravel(),v.ravel(),dd.ravel());
return (dout,u,v);
# function to make a numpy N-by-1 column vector
# c=eda_cvec(v1, v2, ...) from a list of several
# array-like entities v1, v2, including a number
# a list of numbers, a tuple of numbers, an N-by-0 np array
# and a N-by-1 np array. The function
# also insures that, if v is an np array, that
# c is a copy, as contrasted to a view, of it
# It promotes integers to floats, and integers
# and floats to complex, by context.
# This version concatenates many argments,
# whereas c=eda_cvec1(v1) takes just one argiment.
# I recommend always using eda_cvec(v1, v2, ...)
def eda_cvec(*argv):
t = int;
Nt = 0;
for a in argv:
v = eda_cvec1(a);
N,M = np.shape(v);
Nt = Nt + N;
if( N==0 ):
continue; # skip vector of zero length
if (t==int) and isinstance(v[0,0],float):
t=float;
elif isinstance(v[0,0],complex):
t=complex;
w = np.zeros((Nt,1),dtype=t);
Nt = 0;
for a in argv:
v = eda_cvec1(a);
N,M = np.shape(v);
w[Nt:Nt+N,0] = v[0:N,0];
Nt = Nt + N;
return w;
# function to make a numpy N-by-1 column vector
# c=eda_cvec1(v) from entity v that is array-like,
# including a number, a list of numbers, a tuple
# of numbers, an N-by-0 np array and a N-by1 np array.
# It promotes integers to floats, and integers
# and floats to complex, by context. The function
# also insures that, if v is an np array, that
# c is a copy, as contrasted to a view, of it.
# This version takes just one input argment.
# whereas c=eda_cvec(v1,v2,...) concatenates
# many argiments.
def eda_cvec1(v):
if isinstance(v, int) or isinstance(v, np.int32):
w = np.zeros((1,1),dtype=int);
w[0,0] = v;
return w;
elif isinstance(v, float):
w = np.zeros((1,1),dtype=float);
w[0,0] = v;
return w;
elif isinstance(v, complex):
w = np.zeros((1,1),dtype=complex);
w[0,0] = v;
return w;
elif isinstance(v, np.ndarray):
s = np.shape(v);
if len(s) == 1:
return np.copy(np.reshape(v,(s[0],1)));
else:
[r,c]=s;
if( c==1 ):
return(np.copy(v));
elif(r==1):
return(np.copy(v.T));
else:
raise TypeError("eda_cvec: %d by %d ndarray not allowed" % (r, c));
elif isinstance(v, list):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(v, np.int32):
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("eda_cvec: list contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = v;
return w;
elif isinstance(v, tuple):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(v, np.int32):
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("eda_cvec: tuple contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = v;
return w;
else:
raise TypeError("eda_cvec: %s not supported" % type(v));
def eda_sign(x):
# sign function; apparently Python doesn't have one
if x>0.0:
s=1.0;
elif x<0.0:
s=(-1.0);
else:
s=0.0;
return(s);
# function to perform the multiplication y = (Ccc + sigmad2 I) v
def eda_CpImul(v):
global eda_N;
global eda_td;
global eda_sigmad2;
global eda_gamma2;
global eda_w0est;
global eda_sest;
# rearrange as column-vector
s = np.shape(v);
M = s[0];
if(len(s)==1):
vv = np.zeros((s[0],1));
vv[:,0] = v;
else:
vv=v;
# Note (C + sigmad2 I) v implies yi = SUMj Cij vj + sigmad2 vi
y = np.zeros((eda_N,1));
for i in range(eda_N):
y[i,0]=eda_sigmad2*vv[i,0];
for j in range(eda_N):
Tabs = abs(eda_td[i,0]-eda_td[j,0]);
Cij = eda_gamma2*cos(eda_w0est*Tabs)*exp(-eda_sest*Tabs);
y[i,0]=y[i,0]+Cij*vv[j,0];
return y;
In [2]:
# edapy_10_01, fitting a high-order polynomial
# illustrate problem associated with interpolation by
# fitting a high-order polynomial
# times
N=11;
t = eda_cvec( [0.01, 0.11, 0.24, 0.38, 0.42, 0.57, 0.62, 0.77, 0.81, 0.91, 0.99] );
# simple function
d = t/2.0 + eda_cvec([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
fig1 = plt.figure(1,figsize=(10,6));
plt.axis([0, 1, -5, 5]);
plt.plot( t, d, 'ko');
plt.xlabel('time t (s)');
plt.ylabel('d(t)');
# linear representation of polynomial relationship
G=np.zeros((N,N));
for i in range(N):
G[:,i]=np.power(t,i).ravel();
# solve for coefficients
m = la.solve(G,d);
# evaluate at equally-spaced values of t
Np=101;
tp = eda_cvec(np.linspace(0.0,1.0,Np));
Gp = np.zeros((Np,N));
for i in range(N):
Gp[:,i] = np.power( tp, i ).ravel();
dp = np.matmul(Gp,m);
# plot simple function
plt.plot(tp,dp,'k-');
plt.show();
In [3]:
# edapy_10_02, linear interpolation
# times
N=11;
t = eda_cvec( [0.01, 0.11, 0.24, 0.38, 0.42, 0.57, 0.62, 0.77, 0.81, 0.91, 0.99] );
# simple function
d = t/2 + eda_cvec( [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
fig1 = plt.figure(1,figsize=(10,6));
plt.axis([0, 1, -5, 5]);
plt.plot( t, d, 'ko');
plt.xlabel('time t (s)');
plt.ylabel('d(t)');
# qually-spaced values of t
Np=101;
tp = eda_cvec( np.linspace(0.0,1.0,Np) );
# linear interpolation
f = ip.interp1d(t.ravel(),d.ravel(),axis=0,
kind='linear',fill_value='extrapolate');
dp = eda_cvec(f(tp));
# plot simple function
plt.plot(tp,dp,'k-');
plt.show();
In [4]:
# edapy_10_03, cubic interpolation
# times
N=11;
t = eda_cvec( [0.01, 0.11, 0.24, 0.38, 0.42, 0.57, 0.62, 0.77, 0.81, 0.91, 0.99] );
# simple function
d = t/2.0 + eda_cvec( [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
fig1 = plt.figure(1,figsize=(10,6));
plt.axis([0, 1, -5, 5]);
plt.plot( t, d, 'ko');
plt.xlabel('time t (s)');
plt.ylabel('d(t)');
# qually-spaced values of t
Np=101;
tp = eda_cvec( np.linspace(0,1,Np) );
# cubic interpolation
f = ip.interp1d(t.ravel(),d.ravel(),axis=0,
kind='cubic',fill_value='extrapolate');
dp=eda_cvec(f(tp));
# plot simple function
plt.plot(tp,dp,'k-');
plt.show();
In [5]:
# edapy_10_04: make correlated random time series
# 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=eda_cvec(np.linspace(0, tmax, N));
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==1:
c = gamma2*np.exp(-0.5*s*np.power((t-tmid),2));
elif CTYPE==0:
c = gamma2*np.multiply( np.cos(w0*(t-tmid)), np.exp(-s*abs(t-tmid)) );
else:
c = gamma2*np.exp(-s*np.abs(t-tmid));
print("gamma: %.4f s: %.4f" % (gamma, s) );
# realize a time series with a specified covariance
ctilde = np.fft.fft(c,axis=0); # Fourier transform of c
casd = np.sqrt(np.abs(ctilde)); # amp. spec. density of time series
# make random numbers
d = np.random.normal(loc=0.0,scale=1.0,size=(N,1));
d = d - np.mean(d); # remove mean
dtilde = np.fft.fft(d,axis=0); # Fourier transform of d
dasd = np.abs(dtilde); # a.s.d. of d
# adjust so has a.s.d. of c
dtilde = np.multiply(np.divide(casd,dasd),dtilde);
d = np.real(np.fft.ifft(dtilde,axis=0)); # inverse Fourier transform
vard = np.var(d); # variance of d
d = (gamma/sqrt(vard))*d; # normalize to variance gamma^2
# estimated covariance
cest = eda_cvec(np.correlate(d.ravel(),d.ravel(),mode='full') )/N;
# plot this covariance
fig1 = plt.figure(1,figsize=(10,6));
plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);
plt.xlabel('time lag t (s)');
plt.ylabel('c(t)');
plt.plot( t-tmid,c,'k-',lw=3);
tc = Dt*eda_cvec(np.linspace(-N+1,N-1,2*N-1));
plt.plot( tc,cest,'k:','LineWidth',3);
# plot time series
fig2 = plt.figure(2,figsize=(10,6));
plt.axis([ 0, tmax, -MPLOT, MPLOT]);
plt.xlabel('time t (s)');
plt.ylabel('d(t)');
plt.plot( t,d,'k-',lw=3);
plt.show();
gamma: 1.0000 s: 0.0784
In [3]:
# edapy_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)|
# 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*eda_cvec(np.linspace(0,tmax,M));
# 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*eda_cvec(np.linspace(0,Mf-1,Mf));
# 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; # variance
cf = gamma2*np.multiply(np.cos(w0*(tf-tmid)),np.exp(-s*np.abs(tf-tmid)));
print("gamma: %.4f w0: %.4f s: %.4f" % (gamma, w0, s) );
# PART 2: create "fine" time series with thus
# covariance using Fourier techniques
# amplitude spectral density of true covariance
cftilde = np.fft.fft(cf,axis=0);
cfasd = np.sqrt(np.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.
# random numbers
mf = np.random.normal(loc=0.0,scale=1.0,size=(Mf,1));
mf = mf - np.mean(mf); # remove mean
mftilde = np.fft.fft(mf,axis=0); # Fourier transform of mf
mfasd = np.abs(mftilde); # a.s.d. of mf
# adjust so has a.s.d. of c
mftilde = np.multiply(np.divide(cfasd,mfasd),mftilde);
# inverse Fourier transform
mf = np.real(np.fft.ifft(mftilde,axis=0));
varmf = np.var(mf); # variance of mf
# normalize mf to variance gamma^2
mf = (gamma/sqrt(varmf))*mf;
# 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 = np.copy(mf[0:Mf:K]);
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=np.zeros((Ntarget,1),dtype=int);
Nj=0;
while (Nj<Ntarget):
i=np.random.randint(low=0,high=Mf-1);
if( (((i-1)%K) != 0) and (not np.isin(i,j)) ):
j[Nj,0]=i;
Nj=Nj+1;
# sort for plotting purposes
j = eda_cvec(np.sort(j.ravel()));
j = eda_cvec(np.unique(j.ravel()));
N,i = np.shape(j);
print('target N: %d, actual %d' % (Ntarget,N) );
# times of observations
td = np.copy(tf[j.ravel()]);
# irregularly sub-sampled data
dtrue = np.copy(mf[j.ravel()]);
# observed data is true data + uncorrelated random noise
sigmad2 = sigmad**2;
dobs = dtrue+np.random.normal(loc=0,scale=sigmad,size=(N,1));
# the data kernel, G, for the equation d=Gm
G = np.concatenate( (np.zeros((N,M)), np.zeros((N,N))), axis=1 );
# 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 = np.zeros((N,N)); # covariance matrix, C(control,control)
for i in range(N):
for j in range(N):
T = abs(td[i,0]-td[j,0]);
Ccc[i,j] = gamma2*cos(w0est*T)*exp(-sest*T);
Ctc = np.zeros((M,N)); # covariance matrix, C(target,control)
for i in range(M):
for j in range(N):
T = abs(tm[i,0]-td[j,0]);
Ctc[i,j] = gamma2*cos(w0est*T)*exp(-sest*T);
if( PLOTCONF ):
Ctt = np.zeros((M,M)); # covariance matrix, C(target,target)
for i in range(M):
for j in range(M):
T = abs(tm[i,0]-tM[j,0]);
Ctt[i,j] = gamma2*cos(w0est*T)*exp(-sest*T);
else: # use outer product
# build covariances via outer product
Tdd = np.abs(np.matmul(td,np.ones((1,N)))-np.matmul(np.ones((N,1)),(td.T)));
Ccc = gamma2*np.multiply(np.cos(w0est*Tdd),np.exp(-sest*Tdd));
Tmd = np.abs(np.matmul(tm,np.ones((1,N)))-np.matmul(np.ones((M,1)),(td.T)));
Ctc = gamma2*np.multiply(np.cos(w0est*Tmd),np.exp(-sest*Tmd));
if( PLOTCONF ):
Tmm = np.abs(np.matmul(tm,np.ones((1,M)))-np.matmul(np.ones((M,1)),(tm.T)));
Ctt = gamma2*np.multiply(np.cos(w0est*Tmm),np.exp(-sest*Tmm));
# solve GPE problem
A = (Ccc + sigmad2*np.identity(N) );
mest = np.matmul(Ctc,la.solve(A,dobs));
# estimated covariance
cest = eda_cvec(
np.correlate(mest.ravel(),mest.ravel(),mode='full'))/M;
# plot this covariance
fig1 = plt.figure(1,figsize=(10,6));
plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);
plt.xlabel('time lag t (s)');
plt.ylabel('c(t)');
plt.plot( tf-tmid,cf,'k-',lw=3);
tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));
plt.plot(tc,cest,'k:',lw=3);
plt.show();
# time series figure
fig2 = plt.figure(2,figsize=(10,6));
plt.axis( [ 0, tmax, -MPLOT, MPLOT] );
plt.xlabel('time t (s)');
plt.ylabel('d(t) and m(t)');
plt.plot(tm,mtrue,'k-',lw=3);
plt.plot(td,dtrue,'ko',lw=2);
plt.plot(tm,mest,'-',lw=2,color=[0.8,0.8,0.8]);
if( PLOTCONF ): # plot 95% confidence intervals
# 2024/05/26 corrected for error, see errata for page 328
# correct formula is
# Cttest = Ctt - Ctc A-1 Cct with A = Ccc + sigmad2 I
Cttest = Ctt - np.matmul( Ctc, la.solve(A,Ctc.T) );
sigmamest = np.sqrt(eda_cvec(np.diag(Cttest)));
plt.plot(tm,mest-2*sigmamest,'k-',lw=1);
plt.plot(tm,mest+2*sigmamest,'k-',lw=1);
plt.show();
print("Note: This verion corrects an error on page 328 of the book; See the errata.");
gamma: 1.0000 w0: 0.1478 s: 0.0784 target N: 50, actual 50
Note: This verion corrects an error on page 328 of the book; See the errata.
In [4]:
# edapy_10_06, GPE using bicg()
# this example uses synthetic data
# and a gamma2*cos(w0T)*exp(-sT) covariance function,
# with T = |t(i)-t(j)|
# 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
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*eda_cvec(np.linspace(0,tmax,M));
# 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*eda_cvec(np.linspace(0,Mf-1,Mf));
# 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; # variance
cf = gamma2*np.multiply(np.cos(w0*(tf-tmid)),np.exp(-s*np.abs(tf-tmid)));
print("gamma: %.4f w0: %.4f s: %.4f" % (gamma, w0, s) );
# PART 2: create "fine" time series with thus
# covariance using Fourier techniques
# amplitude spectral density of true time series
cftilde = np.fft.fft(cf,axis=0);
cfasd = np.sqrt(np.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.
# random numbers
mf = np.random.normal(loc=0.0,scale=1.0,size=(Mf,1));
mf = mf - np.mean(mf); # remove mean
mftilde = np.fft.fft(mf,axis=0); # Fourier transform of mf
mfasd = np.abs(mftilde); # a.s.d. of mf
# adjust so has a.s.d. of c
mftilde = np.multiply(np.divide(cfasd,mfasd),mftilde);
# inverse Fourier transform
mf = np.real(np.fft.ifft(mftilde,axis=0));
varmf = np.var(mf); # variance of mf
# normalize mf to variance gamma^2
mf = (gamma/sqrt(varmf))*mf;
# 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 = np.copy(mf[0:Mf:K]);
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=np.zeros((Ntarget,1),dtype=int);
Nj=0;
while (Nj<Ntarget):
i=np.random.randint(low=0,high=Mf-1);
if( (((i-1)%K) != 0) and (not np.isin(i,j)) ):
j[Nj,0]=i;
Nj=Nj+1;
# sort for plotting purposes
j = eda_cvec(np.sort(j.ravel()));
j = eda_cvec(np.unique(j.ravel()));
N, i =np.shape(j);
print('target N: %d, actual %d' % (Ntarget, N) );
# times of observations
td = np.copy(tf[j.ravel()]);
# irregularly sub-sampled data
dtrue = np.copy(mf[j.ravel()]);
# observed data is true data + uncorrelated random noise
sigmad2 = sigmad**2;
dobs = dtrue+np.random.normal(loc=0,scale=sigmad,size=(N,1));
# 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;
# set globals
eda_N=N;
eda_td=np.copy(td);
eda_sigmad2=sigmad2;
eda_gamma2=gamma2;
eda_w0est=w0est;
eda_sest=sest;
# define linear operator needed for conjugate gradienet solver
LO=las.LinearOperator(shape=(N,N),
matvec=eda_CpImul,rmatvec=eda_CpImul);
# solve using congugate gradient algorithm
u = np.zeros((N,1));
q=las.bicg(LO,dobs,rtol=1e-12,maxiter=(3*(N+M)+100));
u = eda_cvec( q[0].astype(float) );
mest = np.zeros((M,1));
for i in range(M):
mest[i,0]=0.0;
for j in range(N):
Tabs = abs(tm[i,0]-td[j,0]);
Cij = gamma2*cos(w0est*Tabs)*exp(-sest*Tabs);
mest[i,0]=mest[i,0]+Cij*u[j,0];
dpre = np.zeros((N,1));
for i in range(N):
dpre[i,0]=0.0;
for j in range(N):
Tabs = abs(td[i,0]-td[j,0]);
Cij = gamma2*cos(w0est*Tabs)*exp(-sest*Tabs);
dpre[i,0]=dpre[i,0]+Cij*u[j,0];
# estimated covariance
cest = eda_cvec(
np.correlate(mest.ravel(),mest.ravel(),mode='full'))/M;
# plot this covariance
fig1 = plt.figure(1,figsize=(10,6));
plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);
plt.xlabel('time lag t (s)');
plt.ylabel('c(t)');
plt.plot( tf-tmid,cf,'k-',lw=3);
tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));
plt.plot(tc,cest,'k:',lw=3);
plt.show();
# time series figure
fig2 = plt.figure(2,figsize=(10,6));
plt.axis( [ 0, tmax, -MPLOT, MPLOT] );
plt.xlabel('time t (s)');
plt.ylabel('d(t) and m(t)');
plt.plot(tm,mtrue,'k-',lw=3);
plt.plot(td,dtrue,'ko',lw=2);
plt.plot(tm,mest,'-',lw=2,color=[0.8,0.8,0.8]);
plt.show();
gamma: 1.0000 w0: 0.1478 s: 0.0784 target N: 50, actual 50
In [16]:
# edapy_10_07, GPR Tuning Example
# this example uses synthetic data
# and a gamma2*cos(w0T)*exp(-sT) covariance function,
# with T = |t(i)-t(j)|
# 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=eda_cvec(np.linspace(0,tmax,M));
# 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=eda_cvec(np.linspace(0,tmax,Mf));
# 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; # midpoint
s = Dtm/60.0; # decay rate
gamma2 = gamma**2; # variance
cf = gamma2*np.multiply(np.cos(w0*(tf-tmid)),np.exp(-s*np.abs(tf-tmid)));
print("gamma: %.4f w0: %.4f s: %.4f" % (gamma, w0, s) );
# PART 2: create "fine" time series with thus
# covariance using Fourier techniques
# amplitude spectral density of true time series
cftilde = np.fft.fft(cf,axis=0);
cfasd = np.sqrt(np.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.
# random numbers
mf = np.random.normal(loc=0.0,scale=1.0,size=(Mf,1));
mf = mf - np.mean(mf); # remove mean
mftilde = np.fft.fft(mf,axis=0); # Fourier transform of mf
mfasd = np.abs(mftilde); # a.s.d. of mf
# adjust so has a.s.d. of c
mftilde = np.multiply(np.divide(cfasd,mfasd),mftilde);
# inverse Fourier transform
mf = np.real(np.fft.ifft(mftilde,axis=0));
varmf = np.var(mf); # variance of mf
# normalize mf to variance gamma^2
mf = (gamma/sqrt(varmf))*mf;
# 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 = np.copy(mf[0:Mf:K]);
Ntarget = 150; # 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=np.zeros((Ntarget,1),dtype=int);
Nj=0;
while (Nj<Ntarget):
i=np.random.randint(low=0,high=Mf-1);
if( ((i%K) != 0) and (not np.isin(i,j)) ):
j[Nj,0]=i;
Nj=Nj+1;
# sort for plotting purposes
j = eda_cvec(np.sort(j.ravel()));
j = eda_cvec(np.unique(j));
N, i =np.shape(j);
print('target N: %d, actual %d' % (Ntarget, N) );
# times of observations
td = np.copy(tf[j.ravel()]);
# irregularly sub-sampled data
dtrue = np.copy(mf[j.ravel()]);
# observed data is true data + uncorrelated random noise
sigmad2 = sigmad**2;
sigmadr = 1.0/sigmad;
dobs = dtrue+np.random.normal(loc=0,scale=sigmad,size=(N,1));
# the data kernel, G, for the equation d=Gm
G = np.concatenate( (np.zeros((N,M)), np.zeros((N,N))), axis=1 );
# PART 4: Solve GPR problem for initial guess of parameter
w0est = wbias*w0; # w0est = 1.0*w0 is the true value
sest = sbias*s;
# build covariances via outer product
Tdd = np.abs(np.matmul(td,np.ones((1,N)))-np.matmul(np.ones((N,1)),(td.T)));
Ccc = gamma2*np.multiply(np.cos(w0est*Tdd),np.exp(-sest*Tdd));
Tmd = np.abs(np.matmul(tm,np.ones((1,N)))-np.matmul(np.ones((M,1)),(td.T)));
Ctc = gamma2*np.multiply(np.cos(w0est*Tmd),np.exp(-sest*Tmd));
# solve GPE problem
A = (Ccc + sigmad2*np.identity(N) );
mest = np.matmul(Ctc,la.solve(A,dobs));
# estimated covariance
cest = eda_cvec(np.correlate(mest.ravel(),mest.ravel(),mode='full'))/M;
# plot startng covariance
fig1 = plt.figure(1,figsize=(10,6));
plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);
plt.xlabel('time lag t (s)');
plt.ylabel('c(t)');
plt.plot( tf-tmid,cf,'k-',lw=3);
tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));
plt.plot(tc,cest,'k:',lw=3);
plt.show();
# time series figure
fig2 = plt.figure(2,figsize=(10,6));
plt.axis( [ 0, tmax, -MPLOT, MPLOT] );
plt.xlabel('time t (s)');
plt.ylabel('d(t) and m(t)');
plt.plot(tm,mtrue,'k-',lw=3);
plt.plot(td,dtrue,'ko',lw=2);
plt.plot(tm,mest,'-',lw=2,color=[0.8,0.8,0.8]);
plt.show();
# Part 5: Set-up for 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 = np.zeros((M,1)); # prior model parameters are zero
mcpri = np.zeros((N,1)); # prior model parameters are zero
# concatenate into overall prior vector mpri and times tall
mpri = np.concatenate((mtpri,mcpri),axis=0); # prior model pars
tall = np.concatenate((tm,td),axis=0); # corresponding times
qsave = np.zeros((NITER,1)); # save estimates, for plot
# Part 6: Tuning
for iter in range(NITER): # loop over descent increments
qg = qg+nu*alpha*Dq; # new guess of parameter, q
# full covariance matrix and its derivative
Tall = np.matmul(tall,np.ones((1,M+N))) - np.matmul(np.ones((M+N,1)),tall.T);
Call = gamma2*np.multiply(np.cos(qg*Tall),np.exp(-sest*np.abs(Tall)));
dCalldq = -gamma2*np.multiply(
np.multiply(Tall,np.sin(qg*Tall)),np.exp(-sest*abs(Tall)));
# break out four quandrants; note these are views, not arrays,
# so we need to be careful not to overwrite Ctt
Ctt = Call[0:M,0:M];
dCttdq = dCalldq[0:M,0:M];
Ctc = Call[0:M,M:M+N];
dCtcdq = dCalldq[0:M,M:M+N];
Cct = Ctc.T;
Ccc = Call[M:M+N,M:M+N];
dCccdq = dCalldq[M:M+N,M:M+N];
# A = Ccc + sigmad^2 I
A = Ccc + sigmad2*np.identity(N);
dAdq = dCccdq;
# various functions of the covariance matrix
CallI = la.inv(Call);
CallIsr = la.sqrtm(CallI);
dCallIdq = -np.matmul(np.matmul(CallI,dCalldq), CallI);
dCallIsrdq = la.solve_sylvester( CallIsr, CallIsr, dCallIdq );
R = la.cholesky(Call);
logdetCall = 2.0*np.sum(np.log(np.diag(R)));
dlogdetCalldq = np.trace(np.matmul(CallI,dCalldq));
# GPR solution
u=la.solve(A,dobs);
mtg = np.matmul(Ctc,u);
mcg = np.matmul(Ccc,u);
mg = np.concatenate( (mtg,mcg), axis=0 );
dg = eda_cvec(mg[M:M+N]);
# derivative of predicted data
v = la.solve( A, np.matmul(dAdq,u) );
dmestdq = np.concatenate((
np.matmul(dCtcdq,u) - np.matmul(Ctc,v),
np.matmul(dCccdq,u) - np.matmul(Ccc,v)),
axis=0 );
ddpredq = np.copy(dmestdq[M:N+M]);
detildedq = -sigmadr*ddpredq;
e = dobs-dg;
etilde = sigmadr*e;
# derivative of total data error
E = np.matmul(etilde.T,etilde); E=E[0,0];
dEdp = 2.0*np.matmul(etilde.T,detildedq);
# prior error and its derivative
el=(mpri-mg);
eltilde = np.matmul(CallIsr,el);
deltildedq = np.matmul(dCallIsrdq,el)-np.matmul(CallIsr,dmestdq);
# total prior error and its derivative
L = np.matmul(eltilde.T,eltilde); L=L[0,0];
dLdp = 2.0*np.matmul(eltilde.T,deltildedq);
# function, PSI being minimized
PSI = logdetCall + E + L;
# and its deriative
dPSIdp = dlogdetCalldq + dEdp + dLdp;
nu = -eda_sign(dPSIdp); # downhill vector
# gradient method method
if( iter==0 ): # on the first iteration do nothing but save
qglast = qg;
PSIlast = PSI;
print("1: qg %.4f nu %.2f alpha %e PSI %e PSIlast %e" % (qg,nu,alpha,PSI,PSIlast));
elif (PSI < PSIlast ):
qglast = qg;
PSIlast = PSI;
print("2: qg %.4f nu %.2f alpha %e PSI %e PSilast %e" % (qg,nu,alpha,PSI,PSIlast));
else: # psi increased, so decrease alpha
qg = qglast;
PSI = PSIlast;
alpha = alpha/2.0;
print("2: qg %.4f nu %.2f alpha %e PSI %e PSIlast %e" % (qg,nu,alpha,PSI,PSIlast));
qsave[iter,0] = qg;
pest = qg;
mest = mg;
mtest = np.copy(mest[0:M,0]);
mcest = np.copy(mest[M:M+N,0]);
dpre = mcest;
# estimated covariance
cest = eda_cvec(
np.correlate(mtest.ravel(),mtest.ravel(),mode='full'))/M;
# plot final estimated covariance
fig3 = plt.figure(3,figsize=(10,6));
plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);
plt.xlabel('time lag t (s)');
plt.ylabel('c(t)');
plt.plot( tf-tmid,cf,'k-',lw=3);
tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));
plt.plot(tc,cest,'k:',lw=3);
plt.show();
# plot final time series
fig4 = plt.figure(4,figsize=(10,6));
plt.axis( [ 0, tmax, -MPLOT, MPLOT] );
plt.xlabel('time t (s)');
plt.ylabel('d(t) and m(t)');
plt.plot(tm,mtrue,'k-',lw=3);
plt.plot(td,dtrue,'ko',lw=2);
plt.plot(tm,mtest,'-',lw=2,color=[0.8,0.8,0.8]);
plt.show();
print("qtrue %e qest %e relating error %e" % (qtrue, pest, (qtrue-pest)/qtrue) );
# plot p as a function of iteration number
fig5 = plt.figure(5,figsize=(10,6));
plt.axis( [ 0, NITER, 0, 0.3] );
plt.xlabel('iteration');
plt.ylabel('w0');
itvec = eda_cvec(np.linspace(1,NITER,NITER));
plt.plot(itvec,qsave,'k-',lw=2);
plt.plot(itvec,qsave,'ko',lw=2);
plt.plot(itvec,w0*np.ones((NITER,1)),'k:',lw=2);
gamma: 1.0000 w0: 0.0738 s: 0.0167 target N: 150, actual 150
1: qg 0.2213 nu -1.00 alpha 1.000000e+00 PSI -1.997520e+03 PSIlast -1.997520e+03 2: qg 0.1963 nu -1.00 alpha 1.000000e+00 PSI -2.081769e+03 PSilast -2.081769e+03 2: qg 0.1713 nu -1.00 alpha 1.000000e+00 PSI -2.165762e+03 PSilast -2.165762e+03 2: qg 0.1463 nu -1.00 alpha 1.000000e+00 PSI -2.243025e+03 PSilast -2.243025e+03 2: qg 0.1213 nu -1.00 alpha 1.000000e+00 PSI -2.306142e+03 PSilast -2.306142e+03 2: qg 0.0963 nu -1.00 alpha 1.000000e+00 PSI -2.347793e+03 PSilast -2.347793e+03 2: qg 0.0713 nu -1.00 alpha 1.000000e+00 PSI -2.362663e+03 PSilast -2.362663e+03 2: qg 0.0713 nu 1.00 alpha 5.000000e-01 PSI -2.362663e+03 PSIlast -2.362663e+03 2: qg 0.0713 nu -1.00 alpha 2.500000e-01 PSI -2.362663e+03 PSIlast -2.362663e+03 2: qg 0.0713 nu 1.00 alpha 1.250000e-01 PSI -2.362663e+03 PSIlast -2.362663e+03 2: qg 0.0713 nu -1.00 alpha 6.250000e-02 PSI -2.362663e+03 PSIlast -2.362663e+03 2: qg 0.0698 nu 1.00 alpha 6.250000e-02 PSI -2.362669e+03 PSilast -2.362669e+03 2: qg 0.0698 nu -1.00 alpha 3.125000e-02 PSI -2.362669e+03 PSIlast -2.362669e+03 2: qg 0.0698 nu 1.00 alpha 1.562500e-02 PSI -2.362669e+03 PSIlast -2.362669e+03 2: qg 0.0702 nu 1.00 alpha 1.562500e-02 PSI -2.362677e+03 PSilast -2.362677e+03 2: qg 0.0705 nu -1.00 alpha 1.562500e-02 PSI -2.362679e+03 PSilast -2.362679e+03 2: qg 0.0705 nu 1.00 alpha 7.812500e-03 PSI -2.362679e+03 PSIlast -2.362679e+03 2: qg 0.0705 nu -1.00 alpha 3.906250e-03 PSI -2.362679e+03 PSIlast -2.362679e+03 2: qg 0.0704 nu 1.00 alpha 3.906250e-03 PSI -2.362679e+03 PSilast -2.362679e+03 2: qg 0.0704 nu -1.00 alpha 1.953125e-03 PSI -2.362679e+03 PSIlast -2.362679e+03 2: qg 0.0704 nu 1.00 alpha 9.765625e-04 PSI -2.362679e+03 PSIlast -2.362679e+03 2: qg 0.0704 nu -1.00 alpha 4.882812e-04 PSI -2.362679e+03 PSIlast -2.362679e+03 2: qg 0.0704 nu 1.00 alpha 2.441406e-04 PSI -2.362679e+03 PSIlast -2.362679e+03 2: qg 0.0705 nu 1.00 alpha 2.441406e-04 PSI -2.362679e+03 PSilast -2.362679e+03 2: qg 0.0705 nu -1.00 alpha 2.441406e-04 PSI -2.362679e+03 PSilast -2.362679e+03 2: qg 0.0705 nu 1.00 alpha 1.220703e-04 PSI -2.362679e+03 PSIlast -2.362679e+03 2: qg 0.0705 nu -1.00 alpha 6.103516e-05 PSI -2.362679e+03 PSIlast -2.362679e+03 2: qg 0.0705 nu -1.00 alpha 6.103516e-05 PSI -2.362679e+03 PSilast -2.362679e+03 2: qg 0.0705 nu 1.00 alpha 6.103516e-05 PSI -2.362679e+03 PSilast -2.362679e+03 2: qg 0.0705 nu -1.00 alpha 3.051758e-05 PSI -2.362679e+03 PSIlast -2.362679e+03
qtrue 7.377517e-02 qest 7.045576e-02 relating error 4.499361e-02
In [17]:
# edapy_10_08, 2D interpolation of pressure data
# 2D interpolation of simulated pressure data
# this script
# creates simulated data (Figue 1a)
# creates Delaunay triangle mesh (Figure 1b)
# finds triangle enclosing an arbitrary point (Figure 2)
# linear interpolation (Figure 3)
DOCUBIC = 0; # linear interpolation on 0, cubic on 1
# rows of grid are x
I = 41; # number of x's in grid
Dx = 1.0; # sampling interval
xmin = 0.0; # start position
xmax = 40.0; # end position
# vector of x's
x = xmin + (xmax-xmin)*eda_cvec(np.linspace(0,I-1,I)/(I-1));
# columns of grid are y
J = 41; # number of y's in grid
Dy = 1.0; # sampling interval
ymin = 0.0; # start position
ymax = 40.0; # end position
# vector of y's
y = ymin + (ymax-ymin)*eda_cvec(np.linspace(0,J-1,J)/(J-1));
# grid 10% populated with data
N0 = floor(I*J/10); # number of data
rowindex0=np.random.randint(I,size=(N0,1)); # row index of data
colindex0=np.random.randint(J,size=(N0,1)); # col index of data
# remove duplicate points, because the interpolation
# function doesn't seem to like them, by checking
# eack point against a matrix F of ones and zeros
F = np.zeros((I,J),dtype=int);
rowindex1=np.zeros((N0,1),dtype=int);
colindex1=np.zeros((N0,1),dtype=int);
N=0;
for i in range(N0):
irow = rowindex0[i,0];
icol = colindex0[i,0];
if ( F[irow,icol] == 0 ):
rowindex1[N,0]=irow;
colindex1[N,0]=icol;
N=N+1;
F[irow,icol] = 1;
rowindex=eda_cvec(rowindex1[0:N,0]);
colindex=eda_cvec(colindex1[0:N,0]);
# corresponding x and y values of data
xobs = x[rowindex,0];
yobs = y[colindex,0];
# create simulated pressure data (solution to Laplace's equation)
kappa = 0.1;
dtrue = 10.0*np.multiply( np.sin(kappa*xobs), np.exp(-kappa*yobs) );
sigmad = 0.1;
dobs = dtrue + np.random.normal(0.0,sigmad,(N,1) );
# populate matrix A with observations
A = np.zeros((I,J));
for i in range(N):
A[rowindex[i,0],colindex[i,0]] = dobs[i,0];
# B&W colormap
nc = 256;
bw = np.zeros((nc,4));
v = eda_cvec(0.9*(nc - np.linspace( 0, nc-1, nc ))/(nc-1));
bw[0:nc,0:1] = v;
bw[0:nc,1:2] = v;
bw[0:nc,2:3] = v;
bw[0:nc,3:4] = np.ones((nc,1));
bwcmap = ListedColormap(bw);
fig1 = plt.figure(1,figsize=(16,6));
ax1 = plt.subplot(1,2,1);
plt.axis([ymin, ymax, xmin, xmax]);
left=ymin;
right=ymax;
bottom=xmin;
top=xmax;
plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );
plt.xlabel('y');
plt.ylabel('x');
plt.title('data');
# Delaunay triangulation
XY = np.concatenate( (xobs,yobs), axis=1 );
# XY is a N by 2 array of (x,y) positions of vertices
mytri = sp.Delaunay( XY );
# TRI is a NTRI by 3 array of the indices of
# the vertices associated with each triangle
TRI = mytri.simplices;
[NTRI, i] = np.shape(TRI);
# plot
ax1 = plt.subplot(1,2,2);
plt.axis([ymin, ymax, xmin, xmax]);
left=ymin;
right=ymax;
bottom=xmin;
top=xmax;
plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );
for i in range(NTRI):
v1=TRI[i,0];
v2=TRI[i,1];
v3=TRI[i,2];
plt.plot( [XY[v1,1], XY[v2,1]], [XY[v1,0], XY[v2,0]], 'k-');
plt.plot( [XY[v2,1], XY[v3,1]], [XY[v2,0], XY[v3,0]], 'k-');
plt.plot( [XY[v3,1], XY[v1,1]], [XY[v3,0], XY[v1,0]], 'k-');
plt.xlabel('y');
plt.ylabel('x');
plt.title('triangles');
plt.show();
fig2 = plt.figure(2,figsize=(16,6));
ax1 = plt.subplot(1,1,1);
plt.axis([ymin, ymax, xmin, xmax]);
left=ymin;
right=ymax;
bottom=xmin;
top=xmax;
plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );
plt.xlabel('y');
plt.ylabel('x');
plt.title('one triangle enclosing a specified point');
# find triangle enclosing point x0
x0=[20.5,30.5];
tri0 = mytri.find_simplex( x0 );
k = tri0;
if( k != -1 ):
plt.plot( x0[1], x0[0], 'ko');
v1=TRI[k,0];
v2=TRI[k,1];
v3=TRI[k,2];
plt.plot( [XY[v1,1], XY[v2,1]], [XY[v1,0], XY[v2,0]], 'k-');
plt.plot( [XY[v2,1], XY[v3,1]], [XY[v2,0], XY[v3,0]], 'k-');
plt.plot( [XY[v3,1], XY[v1,1]], [XY[v3,0], XY[v1,0]], 'k-');
plt.show();
# grid of (x,y) values
XGRID = np.matmul( x, np.ones((1,J)) );
YGRID = np.matmul( np.ones((I,1)), y.T );
# interpolation
if( DOCUBIC):
# linear interpolation
BB=ip.griddata( np.concatenate((xobs,yobs),axis=1),
dobs,(XGRID,YGRID), method='cubic', fill_value=0.0);
B = BB[:,:,0];
else:
# cubic intepolation
CC=ip.griddata( np.concatenate((xobs,yobs),axis=1),
dobs,(XGRID,YGRID), method='linear', fill_value=0.0);
C = CC[:,:,0];
B = C;
# plot
fig3 = plt.figure(3,figsize=(16,6));
ax1 = plt.subplot(1,2,1);
plt.axis([ymin, ymax, xmin, xmax]);
left=ymin;
right=ymax;
bottom=xmin;
top=xmax;
plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );
plt.xlabel('y');
plt.ylabel('x');
plt.title('data');
ax1 = plt.subplot(1,2,2);
plt.axis([ymin, ymax, xmin, xmax]);
left=ymin;
right=ymax;
bottom=xmin;
top=xmax;
plt.imshow( np.flipud(B), cmap=bwcmap, vmin=np.min(B), vmax=np.max(B), extent=(left,right,bottom,top) );
plt.xlabel('y');
plt.ylabel('x');
plt.title('interpolated data');
plt.show();
In [12]:
# edapy_10_09, is not provided, since the MATHLAB
# version edama_10_09 addresses an issue that does
# not arise in Python
In [18]:
# edapy_10_10, two-dimensional Fourier transform
# illustrate two-dimensional fourier transform
# Part 1, FFT of cosine wave, reordered so zero us in center
# Part 2, Illustration of use of symmetry
# standard distance/wavenumber set up
Nx=32; # number of x's, must be even
Nxo2 = floor(Nx/2); # integer version of Nx/2
Nkx=Nxo2+1; # number of non-negative x-wavenumbers kx
Dx=1.0; # sampling interval in x
x = eda_cvec(Dx*np.linspace(0,Nx-1,Nx)); # vector of x's
xmax = Dx*(Nx-1); # ending x
kxmax=2*pi/(2*Dx); # nyquist kx
Dkx=kxmax/(Nx/2); # sampling in kx
# full (pos and neg) x-wavenumber vector
kx=eda_cvec( Dkx * np.concatenate( ( np.linspace(0,Nxo2,Nkx), np.linspace(-Nxo2+1,-1,Nxo2-1) ), axis=0 ) );
Ny=32; # number of y's, must be even
Nyo2 = floor(Ny/2); # integer version of Ny/2
Nky=Nyo2+1; # number of non-negative y-wavenumbers ky
Dy=1.0; # sampling interval in y
y = eda_cvec(Dy*np.linspace(0,Ny-1,Ny)); # vector of y's
ymax = Dy*(Ny-1); # ending y
kymax=2*pi/(2*Dy); # nyquist ky
Dky=kymax/(Ny/2); # sampling in ky
# full (pos and neg) y-wavenumber vector
ky=eda_cvec(Dky * np.concatenate( ( np.linspace(0,Nyo2,Nky), np.linspace(-Nyo2+1,-1,Nyo2-1) ), axis=0 ));
# Part 1: cosine wave with radial wavenumber kr1 in the theta diection
theta = (pi/180.0)*30.0;
t1 = eda_cvec( [cos(theta), sin(theta)] );
kr1 = kxmax/4.0;
F = np.zeros((Nx,Ny));
for n in range(Nx):
for m in range(Ny):
F[n,m] = cos( kr1*t1[0,0]*x[n,0] + kr1*t1[1,0]*y[m,0] );
eda_draw('title f(x,y)', F);
# amplitude spectral density
Ftt = np.fft.fft2(F);
norm = ((Dx**2)/xmax) * ((Dy**2)/ymax);
Fasd = norm * np.sqrt( np.power( np.abs(Ftt), 2));
# put into more reasonable order, with zero in the center
Fasd2=np.zeros((Nx,Ny));
Fasd2[Nx-Nkx:Nx,Ny-Nky:Ny] = Fasd[0:Nkx,0:Nky];
Fasd2[0:Nkx-2,Ny-Nky:Ny] = Fasd[Nkx:Nx,0:Nky];
Fasd2[0:Nkx-2,0:Nky-2] = Fasd[Nkx:Nx,Nky:Ny];
Fasd2[Nkx-2:Nx,0:Nky-2] = Fasd[0:Nkx,Nky:Ny];
# Part 2: example of the use of symmetry
# new function F is noise, which has a complicated Fourier Transform
F=np.random.normal(0.0,1.0,(Nx,Ny) );
Ftt = np.fft.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=np.zeros((Nx,Ny),dtype=complex);
Ftt2[:,0:Nky] = Ftt[:,0:Nky];
# then rebuild right from left columns
for m in range(1,Nyo2):
mp=Ny-m;
Ftt2[0,mp] = np.conj(Ftt2[0,m]);
Ftt2[1:Nx,mp] = np.flipud(np.conj(Ftt2[1:Nx,m]));
# invert back to F
F2 = np.real(np.fft.ifft2(Ftt2));
eda_draw('Part 1', 'title ASD', Fasd, 'title ASD rearranged', Fasd2, 'Part 2', 'title F2', F, 'title reconstructed F2', F2);
# compute error, as a check
Ett=np.sum(np.abs(Ftt2-Ftt));
E = np.sum(np.abs(F2-F));
print('Symmetry test:');
print(" error in transform %.4f" % (Ett) );
print(" error in function %.4f" % (E) );
Symmetry test:
error in transform 0.0000
error in function 0.0000
In [19]:
# edapy_10_11: Filling in gaps with fft using POCS
# Part 1, generate true synthetic data set
# with Gaussian autocorrelation function
# that simulates a crudely-layerd medium
# standard distance/wavenumber set up
Nx=128; # number of x's, must be even
Nxo2 = floor(Nx/2); # integer version of Nx/2
Nkx=Nxo2+1; # number of non-negative x-wavenumbers kx
Dx=1.0; # sampling interval in x
x = eda_cvec(Dx*np.linspace(0,Nx-1,Nx)); # vector of x's
xmax = Dx*(Nx-1); # ending x
kxmax=2*pi/(2*Dx); # nyquist kx
Dkx=kxmax/(Nx/2); # sampling in kx
# full (pos and neg) x-wavenumber vector
kx=eda_cvec( Dkx * np.concatenate( ( np.linspace(0,Nxo2,Nkx), np.linspace(-Nxo2+1,-1,Nxo2-1) ), axis=0 ) );
Ny=128; # numer of y's, must be even
Nyo2 = floor(Ny/2); # integer version of Ny/2
Nky=Nyo2+1; # number of non-negative y-wavenumbers ky
Dy=1.0; # sampling interval in y
y = eda_cvec(Dy*np.linspace(0,Ny-1,Ny)); # vector of y's
ymax = Dy*(Ny-1); # ending y
kymax=2*pi/(2*Dy); # nyquist ky
Dky=kymax/(Ny/2); # sampling in ky
# full (pos and neg) y-wavenumber vector
ky=eda_cvec(Dky * np.concatenate( ( np.linspace(0,Nyo2,Nky), np.linspace(-Nyo2+1,-1,Nyo2-1) ), axis=0 ));
# 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 = np.zeros((Nx,Ny)); # target amplitude spectral density
for ikx in range(Nx):
for iky in range(Ny):
psd=exp(-0.5*(kx[ikx,0]**2)/sx2-0.5*(ky[iky,0]**2)/sy2);
asd[ikx,iky] = sqrt(psd);
# start with random data and
mtrue = np.random.normal(loc=0.0,scale=1.0,size=(Nx,Ny));
Ftt = np.fft.fft2(mtrue); # take its fft
Ftt = np.multiply(asd,Ftt); # shape to target
mtrue = np.real(np.fft.ifft2(Ftt)); # inverse fft
sigmamtrue = sqrt(np.var(mtrue)); # true variance
# normalization ratio, gives target asd the same
# variance as the fft of the data
R=sqrt(np.var(np.abs(Ftt))) / sqrt(np.var(np.abs(asd)));
Rasd=R*asd;
# sample mtrue along vertical profiles
j = np.floor( eda_cvec( 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 )*Ny/8.0).astype(int);
Nd,i = np.shape(j);
d = np.zeros((Nx,Nd));
for i in range(Nd):
k = j[i,0];
d[0:Nx,i:i+1]=mtrue[0:Nx,k:k+1];
print('profiles through true model');
eda_draw(d[0:Nx,0:1],d[0:Nx,1:2],d[0:Nx,2:3],d[0:Nx,3:4],d[0:Nx,4:5],d[0:Nx,5:6],d[0:Nx,6:7]);
# Part 2, POCS method to fill in data gaps
# make starting model that has the right a.s.d. and variance
mest = np.random.normal(loc=0.0,scale=1.0,size=(Nx,Ny)); # random data
Ftt = np.fft.fft2(mest); # take fft
Ftt=np.multiply(asd,Ftt); # shape to target
mest = np.real(np.fft.ifft2(Ftt)); # take inverse fft
sigmamest = sqrt(np.var(mest)); # estimated variance
mest = (sigmamtrue/sigmamest)*mest; # scale to right variance
mstarting = np.copy(mest); # save starting model
# POCS algorithm
Niter=50;
E=np.zeros((Niter,1));
for iter in range(Niter): # loop over iterations
# enforce data are satisfied
for i in range(Nd):
k = j[i,0];
mest[0:Nx,k:k+1]=d[0:Nx,i:i+1];
# enforce a.s.d. is not greater than target
Ftt = np.fft.fft2(mest);
for ikx in range(Nx): # loop over wavenumbers
for iky in range(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];
mest = np.real(np.fft.ifft2(Ftt));
# model error
e = eda_cvec( mtrue.ravel()-mest.ravel() );
myE = np.matmul(e.T,e); myE=myE[0,0];
E[iter,0] = myE;
sigmamtrue=sqrt(np.var(mtrue));
sigmamstarting=sqrt(np.var(mstarting));
sigmamest=sqrt(np.var(mest));
print('sqrt variances: starting %e, estimated %e, true %e' % (sigmamtrue, sigmamstarting, sigmamest));
print('starting, estimated, true model');
eda_draw('title starting',mstarting, 'title est', mest,'title true', mtrue);
fig2 = plt.figure(2,figsize=(16,6));
ax1 = plt.subplot(1,1,1);
plt.axis([0, Niter, 0, E[0,0] ]);
plt.plot( eda_cvec(np.linspace(1,Niter,Niter)), E, 'k-', lw=2);
plt.xlabel('iteration, i');
plt.ylabel('model error, E')
plt.show();
profiles through true model
sqrt variances: starting 4.378116e-02, estimated 4.378116e-02, true 3.372636e-02 starting, estimated, true model
In [ ]: