In [1]:
# gdapy05_00 clears all variables and import various modules
# you must run this cell first, before running any of others
# note that these examples use inverse theory techniques that
# are only covered in later chapters of the book
# History:
# 2022/11/25 - Created by W. Menke
# 2024/05/23 - Patches to accomodate new verions of numpy, scipy, matplotlib
# patches dot product to return scalar value
# las.bicg keyword tol changes to rtol
# A reset deletes any previously-created variables from memory
%reset -f
# import various packages
import numpy as np # matrices & etc
from matplotlib import pyplot as plt # general plotting
from matplotlib import patches # plotting shapes
from datetime import date # date and time arithmetic
from math import exp, pi, sin, cos, tan, sqrt, floor, log10, nan # math functions
import scipy.linalg as la # linear algebra functions
import os # operating system
import matplotlib.cm as cm
from matplotlib.colors import ListedColormap
import scipy.signal as sg
import scipy.stats as st
import scipy.sparse.linalg as las
from scipy import sparse
import matplotlib
# 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 gda_cvec(*argv):
t = int;
Nt = 0;
for a in argv:
v = gda_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 = gda_cvec1(a);
N,M = np.shape(v);
w[Nt:Nt+N,0:1] = v; # patch 20230418 was #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=gda_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=gda_cvec(v1,v2,...) concatenates
# many argiments.
def gda_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("gda_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(vi, np.int32): #patch v->vi 20230418
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("gda_cvec: list contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = np.array(v); # patch v -> np.array(v)
return w;
elif isinstance(v, tuple):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("gda_cvec: tuple contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = np.array(list(v)); # patch v -> np.array(list(v));
return w;
else:
raise TypeError("gda_cvec: %s not supported" % type(v));
# gda_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 gda_draw(*argv):
DOCOLOR=True;
if( DOCOLOR ):
bwcmap = matplotlib.colormaps['jet'];
else:
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;
# bandpass filter, used in seismological example, but hand
# in a variety of settings involving time series
def gda_chebyshevfilt(d, Dt, flow, fhigh):
# (dout,u,v)=gda_chebyshevfilt(d, Dt, flow, fhigh);
# 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 (gda_cvec(dout),gda_cvec(u),gda_cvec(v));
# function to perform the multiplication FT (F v)
def gdaFTFmul(v):
# this function is used by bicongugate gradient to
# solve the generalized least squares problem Fm=F.
# Note that "F" must be called "gdaFsparse".
global gdaFsparse;
s = np.shape(v);
if(len(s)==1):
vv = np.zeros((s[0],1));
vv[:,0] = v;
else:
vv=v;
# weirdly, "*" multiplies sparse matrices just fine
temp = gdaFsparse*vv;
return gdaFsparse.transpose()*temp;
In [2]:
# gdapy05_01
# data Resolution Matrix example
print("gdapy05_01");
# z-axis
N=101;
zmin=0.0;
zmax=10.0;
Dz = (zmax-zmin)/(N-1);
z = gda_cvec(np.linspace(zmin,zmax,N));
# create synthetic data
# d = a + b*z + noise
a=2.0;
b=1.0;
sd=0.5;
dtrue = a+b*z;
dobs = dtrue+np.random.normal(loc=0.0,scale=sd,size=(N,1));
# set up data kernel
M=2;
G=np.zeros((N,M));
G[0:N,0:1] = np.ones((N,1));
G[0:N,1:2] = z;
# generalized inverse
GTG = np.matmul(G.T,G);
GMG = la.solve(GTG,G.T);
# least squars solution
mest = np.matmul(GMG,dobs);
# predicted data
dpre = np.matmul(G,mest);
# compute and plot data resolution matrix
Nres = np.matmul(G,GMG);
gda_draw("title dpre",dpre," = ","title N",Nres,"title dobs",dobs);
# plot data
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
plt.axis( [zmin, zmax, a+b*zmin, a+b*zmax ] );
plt.plot(z,dobs,"ko");
plt.plot(z,dpre,"r-");
plt.xlabel('z');
plt.ylabel('d');
plt.show();
gdapy05_01
In [3]:
# gdapy05_02
# Model Resolution Matrix example
print("gdapy05_02");
# z-axis
M=101;
zmin=0.0;
zmax=10.0;
Dz=(zmax-zmin)/(M-1);
z=gda_cvec(np.linspace(zmin,zmax,M));
# model, m(z), moztly zero but a few spikes
mtrue = np.zeros((M,1));
mtrue[ 5-1,0]=1.0;
mtrue[10-1,0]=1.0;
mtrue[20-1,0]=1.0;
mtrue[50-1,0]=1.0;
mtrue[90-1,0]=1.0;
# experiment: exponential smoothing of model
N=80;
cmin=0.00;
cmax=7.90;
Dc=(cmax-cmin)/(N-1);
c = gda_cvec(np.linspace(cmin,cmax,N));
# data kernel, note use of outer product
G = np.exp(-np.matmul(c,z.T));
# true data and synthetic observed data
sd=1.0e-10;
dtrue = np.matmul(G,mtrue);
dobs = dtrue + np.random.normal(loc=0,scale=sd,size=(N,1));
# damped minimum length solution
epsilon=1e-6;
GGT = np.matmul(G,G.T) + epsilon*np.identity(N);
GGTINV = la.inv(GGT);
GMG = np.matmul( G.T, GGTINV ); # generalized inverse
mest = np.matmul(GMG, dobs); # estimated solution
dpre = np.matmul(G,mest);
Rres = np.matmul(GMG, G); # model resolution matrix
Nres = np.matmul(G, GMG); # data resolution matrix
# plot model resolution matrix
gda_draw('title mest',mest,' = ','title R',Rres,'title mtrue',mtrue);
# plot data kernel
gda_draw('title dpre',dpre,' = ','title G',G,'title mest',mest);
# plot generalizd inverse
gda_draw('title mest',mest,' = ','title GMG',GMG,'title dobs',dobs);
# plot model resolution matrix
gda_draw('title mest',mest,' = ','title R',Rres,'title mtrue',mtrue);
# plot data resolution matrix
gda_draw('title dpre',dpre,' = ','title N',Nres,'title dobs',dobs);
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
# plot scale
pmmin=-0.2;
pmmax=1.0;
plt.axis( [zmin, zmax, pmmin, pmmax ] );
plt.plot( z, mtrue, "r-", 'LineWidth', 3);
plt.plot( z, mest, "b-", 'LineWidth', 3);
plt.xlabel("z");
plt.ylabel("m");
gdapy05_02
In [4]:
# gdapy05_03
# least squares fit to two cases of synthetic data
# distinguished by the spacing of the z's
print("gdapy05_03");
# CASE 1: date clumped in middle of interval
# z-axis
N=10;
zmin=0;
zmax=10;
z = np.sort(np.random.uniform(low=zmin+4*(zmax-zmin)/10.0,high=zmin+6*(zmax-zmin)/10.0,size=(N,1)));
# create synthetic observed data
# d = a + b*z + noise
a=5.0;
b=0.5;
sd=1.0;
sd2 = sd**2;
dobs = a+b*z+np.random.normal(loc=0.0,scale=sd,size=(N,1));
#least squares fit
M=2;
G = np.zeros((N,M));
G[0:N,0:1]=np.ones((N,1));
G[0:N,1:2] = z;
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mest = la.solve(GTG,GTd);
Cm = sd2*la.inv(GTG); # covariance of model parameters
sm = gda_cvec( sqrt(Cm[0,0]), sqrt(Cm[1,1]) ); # sqrt(variance)
# predicted data at a new set of z's
No=10;
zeval = gda_cvec(np.linspace(zmin,zmax,No));
# data kernel
Go = np.zeros((No,M));
Go[0:No,0:1] = np.ones((No,1));
Go[0:No,1:2] = zeval;
deval = np.matmul(Go,mest); # predicted data
Cdeval =np.matmul(Go,np.matmul(Cm,Go.T)); # its covariance
sdeval = gda_cvec(np.sqrt(np.diag(Cdeval))); # its sqrt(variance)
# plot
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 2, 1);
pdmin=0;
pdmax=15;
# plot observed data, predicted data and +/- one-sigma error bounds
plt.axis( [zmin, zmax, pdmin, pdmax ] );
plt.plot( zeval, deval+sdeval, "b-", lw=3);
plt.plot( zeval, deval-sdeval, "b-", lw=3);
plt.plot( zeval, deval, "r-", lw=3);
plt.plot( z, dobs, "ko", lw=3);
for i in range(N):
plt.plot( gda_cvec(z[i,0], z[i,0]), gda_cvec(dobs[i,0]-sd, dobs[i,0]+sd), "k-", lw=2);
plt.xlabel('z');
plt.ylabel('d');
# CASE 2: date clumped in middle of interval
# z-axis
N=10;
zmin=0;
zmax=10;
z = np.sort(np.random.uniform(low=zmin,high=zmax,size=(N,1)));
# create synthetic observed data
# d = a + b*z + noise
a=5.0;
b=0.5;
sd=1.0;
sd2 = sd**2;
dobs = a+b*z+np.random.normal(loc=0.0,scale=sd,size=(N,1));
#least squares fit
M=2;
G = np.zeros((N,M));
G[0:N,0:1]=np.ones((N,1));
G[0:N,1:2] = z;
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mest = la.solve(GTG,GTd);
Cm = sd2*la.inv(GTG); # covariance of model parameters
sm = gda_cvec( sqrt(Cm[0,0]), sqrt(Cm[1,1]) ); # sqrt(variance)
# predicted data at a new set of z's
No=10;
zeval = gda_cvec(np.linspace(zmin,zmax,No));
# data kernel
Go = np.zeros((No,M));
Go[0:No,0:1] = np.ones((No,1));
Go[0:No,1:2] = zeval;
deval = np.matmul(Go,mest); # predicted data
Cdeval =np.matmul(Go,np.matmul(Cm,Go.T)); # its covariance
sdeval = gda_cvec(np.sqrt(np.diag(Cdeval))); # its sqrt(variance)
# plot
ax1 = plt.subplot(1, 2, 2);
pdmin=0;
pdmax=15;
# plot observed data, predicted data and +/- one-sigma error bounds
plt.axis( [zmin, zmax, pdmin, pdmax ] );
plt.plot( zeval, deval+sdeval, "b-", lw=3);
plt.plot( zeval, deval-sdeval, "b-", lw=3);
plt.plot( zeval, deval, "r-", lw=3);
plt.plot( z, dobs, "ko", lw=3);
for i in range(N):
plt.plot( gda_cvec(z[i,0], z[i,0]), gda_cvec(dobs[i,0]-sd, dobs[i,0]+sd), "k-", lw=2);
plt.xlabel('z');
plt.ylabel('d');
gdapy05_03
In [5]:
# gdapy05_04
# verify bordering method solution
print("gdapy05_04");
# random (M-1)x(M-1) symmetric matrix
M=6;
S = np.random.normal(loc=0.0,scale=1.0,size=(M-1,M-1));
S = 0.5*(S+S.T);
u = np.random.normal(loc=0.0,scale=1.0,size=(M-1,1));
z = np.array([[0.0]]);
# concateate and calculate inverse
np.concatenate((S,u),axis=1)
np.concatenate((u.T,z),axis=1)
Z =np.concatenate( (np.concatenate((S,u),axis=1), np.concatenate((u.T,z),axis=1)), axis=0 );
ZINV = la.inv(Z);
# Bordering method
SI = la.inv(S);
d = np.matmul( u.T, np.matmul(SI,u) );
b = np.matmul(SI,u)/d;
c = -1.0/d; # pain! make a 1x1 array
# computationally, if a matrix A is known to be symmetric
# then its best to compute it with a manifestly symmetric formula
# A = (eye(M-1)-b*u')*SI; # form in book, is not manifestly symmetric
# A = SI-b*u'*SI # multiply out SI
# A = SI-(SI*u)*(u'*SI)/d # insert u
v = np.matmul(SI,u); # shorthand for SI*u
A = SI-np.matmul(v,v.T)/d; # this form is manifestly symmetric
ZZINV =np.concatenate( (np.concatenate((A,b),axis=1), np.concatenate((b.T,c),axis=1)), axis=0 );
print("full inverse");
print(ZINV);
print(" ");
print("bordering inverse");
print(ZZINV);
print(" ");
E = np.max(np.abs(ZINV-ZZINV)); # error
print("maximum deviation in inverse: %e" % (E) );
gdapy05_04 full inverse [[-0.19124657 -0.43051428 -0.67231826 -0.54769213 0.14647913 0.13439947] [-0.43051428 -0.61635259 -0.51826484 -0.80063137 0.20362706 0.6599572 ] [-0.67231826 -0.51826484 0.65080454 -0.92965892 -0.70159972 0.33712145] [-0.54769213 -0.80063137 -0.92965892 -1.59628528 -0.34097289 1.3495111 ] [ 0.14647913 0.20362706 -0.70159972 -0.34097289 0.56014299 0.17850021] [ 0.13439947 0.6599572 0.33712145 1.3495111 0.17850021 -0.91126664]] bordering inverse [[-0.19124657 -0.43051428 -0.67231826 -0.54769213 0.14647913 0.13439947] [-0.43051428 -0.61635259 -0.51826484 -0.80063137 0.20362706 0.6599572 ] [-0.67231826 -0.51826484 0.65080454 -0.92965892 -0.70159972 0.33712145] [-0.54769213 -0.80063137 -0.92965892 -1.59628528 -0.34097289 1.3495111 ] [ 0.14647913 0.20362706 -0.70159972 -0.34097289 0.56014299 0.17850021] [ 0.13439947 0.6599572 0.33712145 1.3495111 0.17850021 -0.91126664]] maximum deviation in inverse: 6.661338e-16
In [6]:
# gdapy05_05
# Model Resolution Matrix example, Backus-Gilbert solution
print("gdapy05_05");
# z-axis
M=101;
zmin=0.0;
zmax=10.0;
Dz=(zmax-zmin)/(M-1);
z=gda_cvec(np.linspace(zmin,zmax,M));
# model, m(z), moztly zero but a few spikes
mtrue = np.zeros((M,1));
mtrue[ 5-1,0]=1.0;
mtrue[10-1,0]=1.0;
mtrue[20-1,0]=1.0;
mtrue[50-1,0]=1.0;
mtrue[90-1,0]=1.0;
# experiment: exponential smoothing of model
N=80;
cmin=0.00;
cmax=7.90;
Dc=(cmax-cmin)/(N-1);
c = gda_cvec(np.linspace(cmin,cmax,N));
# build data kernel; note use of outer product
G = np.exp(-np.matmul(c,z.T)); # data kernel
# create synthetic observed data
sd=0.0;
dtrue = np.matmul(G,mtrue);
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1));
# construct Backus-Gilbert solution row-wise
GMG = np.zeros((M,N));
u = np.matmul( G, np.ones((M,1)) );
CHECKS = 0; # one 1, checks calculation of S
for k in range(M):
# note that S is a symmetric NxN matrix
S = np.matmul(np.matmul(G,np.diag(np.power(np.arange(M)-k,2))),G.T);
# thi code checks S against a brute-force calculation
if(CHECKS):
St= np.zeros((N,N));
for i in range(N):
for j in range(N):
tmp=0;
for el in range(M):
tmp=tmp+((el-k)**2)*G[i,el]*G[j,el];
St[i,j]=tmp;
E = np.max(np.abs(S-St));
print(k,E);
epsilon=1e-6; # add a tiny bit of damping
Sp = S+epsilon*np.identity(N);
# uSpinv = np.matmul(u.T,la.inv(Sp));
uSpinv = la.solve(Sp,u).T; # this way faster than commented out line
GMG[k:k+1,0:N] = uSpinv / np.matmul(uSpinv,u); # generalized inverse
mest = np.matmul(GMG, dobs ); # estimated model parameters
dpre = np.matmul(G, mest); # predicted data
Rres = np.matmul(GMG, G); # model resolution matrix
Nres = np.matmul(G, GMG); # data resolution matrix
# plot data kernel
gda_draw('title dpre',dpre,' = ','title G',G,'title mest',mest);
# plot generalizd inverse
gda_draw('title mest',mest,' = ','title GMG',GMG,'title dobs',dobs);
# plot model resolution matrix
gda_draw('title mest',mest,' = ','title R',Rres,'title mtrue',mtrue);
# plot data resolution matrix
gda_draw('title dpre',dpre,' = ','title N',Nres,'title dobs',dobs);
# plot
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
pmmin=-0.2;
pmmax=1.0;
# plot true and estimated model
plt.axis( [zmin, zmax, pmmin, pmmax ] );
plt.plot( z, mtrue, "r-", lw=3);
plt.plot( z, mest, "b-", lw=3);
plt.xlabel("z");
plt.ylabel("m");
plt.show();
gdapy05_05
In [7]:
# gdmpy05_06
# trade off of resolution and variance
print("gdapy05_06");
# z-axis
M=101;
zmin=0.0;
zmax=10.0;
Dz=(zmax-zmin)/(M-1);
z=gda_cvec(np.linspace(zmin,zmax,M));
# experiment: exponential smoothing of model
N=80;
cmin=0.00;
cmax=7.90;
Dc=(cmax-cmin)/(N-1);
c = gda_cvec(np.linspace(cmin,cmax,N));
# build data kernel; note use of outer product
G = np.exp(-np.matmul(c,z.T)); # data kernel
# Part 1: Backus-Gilbert
# I hand-tuned the points on trade-off curve to make the plot look nice
alphavec = gda_cvec(0.999, 0.995, 0.99, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001);
iahalf=3;
Na,t=np.shape(alphavec); # number of points on the trade-off curve
spreadR1=np.zeros((Na,1)); # tabulate spread of resolution
sizeC1=np.zeros((Na,1)); # tabulate size of covariance
CHECKR=0; # on true, checks calculation of resolution matrix
for a in range(Na):
alpha = alphavec[a,0];
# construct Backus-Gilbert solution row-wise
GMG = np.zeros((M,N));
u = np.matmul(G,np.ones((M,1)));
for k in range(M):
S = np.matmul(np.matmul(G,np.diag(np.power(np.arange(M)-k,2))),G.T);
Sp = alpha*S + (1.0-alpha)*np.identity(N);
uSpinv = la.solve(Sp,u).T;
GMG[k:k+1,0:N] = uSpinv / np.matmul(uSpinv,u); # generalized inverse
Cm1 = np.matmul(GMG,GMG.T); # unit covariance matrix
R1 = np.matmul(GMG,G); # model resolution matrix
sizeC1[a,0] = np.sum(np.diag(Cm1)); # size of covariacne
# spread of resolution
v = gda_cvec(np.arange(M));
o = np.ones((M,1));
spreadR1[a,0] = np.sum(np.multiply(np.power(np.matmul(v,o.T)-np.matmul(o,v.T),2),np.power(R1,2)));
# the hard way, to check result
if( CHECKR ):
spreadR=0;
for i in range(M):
for j in range(M):
spreadR = spreadR + ((i-j)**2)*(R1[i,j]**2);
print(log10(spreadR),log10(spreadR1[a,0]),abs(spreadR-spreadR1[a,0]));
# plot
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 2, 1);
# plot trade-off curve
plt.axis( [2.0, 3.5, -3.0, 5.0 ] );
plt.plot( np.log10(spreadR1), np.log10(sizeC1), "k-", lw=4);
x = np.log10( gda_cvec( spreadR1[0,0], spreadR1[iahalf,0], spreadR1[Na-1,0]) );
y = np.log10( gda_cvec( sizeC1[0,0], sizeC1[iahalf,0], sizeC1[Na-1,0]) );
plt.plot(x, y, "ko", lw=2);
plt.xlabel("log10 spread(R)");
plt.ylabel("log10 size(Cm)");
# Part 2: Damped minimum length
# I hand-tuned the points on trade-off curve to make the plot look nice
alphavec = gda_cvec(0.9999, 0.999, 0.995, 0.99, 0.5, 0.4, 0.3, 0.2, 0.1, 0.05, 0.01);
iahalf=4;
Na,tmp=np.shape(alphavec); # number of points on the trade-of curve
spreadR2=np.zeros((Na,1)); # tabulate spread of reolution
sizeC2=np.zeros((Na,1)); # tabulate size of covariance
for a in range(Na):
alpha = alphavec[a,0];
GMG = np.matmul( (alpha*G.T), la.inv(alpha*np.matmul(G,G.T) + (1.0-alpha)*np.identity(N))); # generalized inverse
Cm2 = np.matmul(GMG,GMG.T); # model covariance
R2 = np.matmul(GMG,G); # model resolution
sizeC2[a,0] = np.sum(np.diag(Cm2)); # size of covariance
spreadR2[a,0] = np.sum(np.power(R2-np.identity(M),2)); # spread of resulution
# plot trade-off curve
plt.subplot(1,2,2);
plt.axis( [90, 100, -3, 4 ] );
plt.plot( spreadR2, np.log10(sizeC2), "k-", lw=4);
x = gda_cvec( spreadR2[0,0], spreadR2[iahalf,0], spreadR2[Na-1,0] );
y = np.log10( gda_cvec( sizeC2[0,0], sizeC2[iahalf,0], sizeC2[Na-1,0] ));
plt.plot( x , y, "ko", lw=2);
plt.xlabel("spread(R)");
plt.ylabel("log10 size(Cm)");
plt.show();
gdapy05_06
In [8]:
# gdapy05_07
# simple figures of rays crossing a box
print("gdapy05_07");
# Part 1: coarse grid
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1=plt.subplot(1,1,1);
ax1.axis('equal')
ax1.axis('off')
# box has (x,y)-axes
xmin=0.0;
xmax=2.0;
ymin=0.0;
ymax=2.0;
plt.plot( gda_cvec(xmin, xmax, xmax, xmin, xmin), gda_cvec(ymin, ymin, ymax, ymax, ymin), "k-", lw=2);
# box contains a circular object
Nth=101;
R=0.9;
x0 = 1.0;
y0 = 1.0;
th = gda_cvec(np.linspace(0.0,2.0*pi,Nth));
x = x0+R*np.sin(th);
y = y0+R*np.cos(th);
plt.axis( [xmin, xmax, ymin, ymax] );
# axis off;
# axis equal;
plt. plot( x, y, "m-", lw=4 );
# rays are generated randomly
Nr=21;
x1 = np.random.uniform(low=xmin, high=xmax, size=(Nr,1) );
x2 = np.random.uniform(low=xmin, high=xmax, size=(Nr,1) );
for i in range(Nr):
plt.plot( gda_cvec(x1[i,0], x2[i,0]), gda_cvec(ymin, ymax), "b-", lw=2);
y1 = np.random.uniform( low=ymin, high=ymax, size=(Nr,1) );
y2 =np. random.uniform( low=ymin, high=ymax, size=(Nr,1) );
for i in range(Nr):
plt.plot( gda_cvec(xmin, xmax), gda_cvec(y1[i,0], y2[i,0]), "b-", lw=2);
# (x,y) grid is superimosed on the box
Ng = 7;
xg = gda_cvec(np.linspace(xmin,xmax,Ng));
yg = gda_cvec(np.linspace(ymin,ymax,Ng));
for i in range(1,Ng-1):
plt.plot( gda_cvec(xg[0,0], xg[Ng-1,0]), gda_cvec(yg[i,0], yg[i,0]), "k-", lw=2 );
for i in range(1,Ng-1):
plt.plot( gda_cvec(xg[i,0], xg[i,0]), gda_cvec(yg[0,0], yg[Ng-1,0]), "k-", lw=2 );
plt.show();
# Part 2, finer grid
fig2 = plt.figure(2,figsize=(12,5)); # smallish figure
ax2=plt.subplot(1,1,1);
ax2.axis('equal')
ax2.axis('off')
# box has (x,y)-axes
xmin=0.0;
xmax=2.0;
ymin=0.0;
ymax=2.0;
plt.plot( gda_cvec(xmin, xmax, xmax, xmin, xmin), gda_cvec(ymin, ymin, ymax, ymax, ymin), "k-", lw=2);
# box contains a circular object
Nth=101;
R=0.9;
x0 = 1.0;
y0 = 1.0;
th = gda_cvec(np.linspace(0.0,2.0*pi,Nth));
x = x0+R*np.sin(th);
y = y0+R*np.cos(th);
plt.axis( [xmin, xmax, ymin, ymax] );
# axis off;
# axis equal;
plt. plot( x, y, "m-", lw=4 );
# rays are generated randomly
Nr=21;
x1 = np.random.uniform(low=xmin, high=xmax, size=(Nr,1) );
x2 = np.random.uniform(low=xmin, high=xmax, size=(Nr,1) );
for i in range(Nr):
plt.plot( gda_cvec(x1[i,0], x2[i,0]), gda_cvec(ymin, ymax), "b-", lw=2);
y1 = np.random.uniform( low=ymin, high=ymax, size=(Nr,1) );
y2 =np. random.uniform( low=ymin, high=ymax, size=(Nr,1) );
for i in range(Nr):
plt.plot( gda_cvec(xmin, xmax), gda_cvec(y1[i,0], y2[i,0]), "b-", lw=2);
# (x,y) grid is superimosed on the box
Ng = 21;
xg = gda_cvec(np.linspace(xmin,xmax,Ng));
yg = gda_cvec(np.linspace(ymin,ymax,Ng));
for i in range(1,Ng-1):
plt.plot( gda_cvec(xg[0,0], xg[Ng-1,0]), gda_cvec(yg[i,0], yg[i,0]), "k-", lw=2 );
for i in range(1,Ng-1):
plt.plot( gda_cvec(xg[i,0], xg[i,0]), gda_cvec(yg[0,0], yg[Ng-1,0]), "k-", lw=2 );
plt.show();
gdapy05_07
In [8]:
# gdapy05_08
# Constructing just one column of the resolution matrix
# Note:
# A row of the resolution matrix R tells you the contributions
# that all the true model parameters make to a single estimated
# model parameter; that is, the estimated model parameter is
# a weighted average of the true model parameter, where the
# elements of the row gives the weights.
#
# On the other hand, a column of the resolution matrix tells
# you all the estimated model parameters that are influenced
# by a true model parameter. The column gives the "blur"
# of estimated model parameters that would be observed if
# the true model consisted of just a single spike.
#
# The distinction between row and column is important,
# since R is, in general, not symmetric.
#
# inverse problem considered here is an acoustic tomography
# problem, where the observations are along just rows and
# columns
print("gdapy05_08");
# grid of unknowns is Lx by Ly
Lx = 20;
Ly = 20;
M = Lx*Ly;
# observations only along rows and columns
N=Lx+Ly;
# build index tables for convenience
iofk=np.zeros((M,1),dtype=int); # backward index, i(k)
jofk=np.zeros((M,1),dtype=int); # backward index, j(k)
kofij=np.zeros((Lx,Ly),dtype=int); # forward index, k(i,j)
k=0; # counter
for i in range(Lx):
for j in range(Ly):
iofk[k,0]=i;
jofk[k,0]=j;
kofij[i,j]=k;
k=k+1;
# observations across rows
G=np.zeros((N,M));
for i in range(Lx):
for j in range(Ly):
k=kofij[i,j]; # (ix,iy) maps into scalar index j
G[i,k]=1.0;
# observations across columns
for j in range(Ly):
for i in range(Lx):
k=kofij[i,j]; # (ix,iy) maps into scalar index j
G[j+Lx,k]=1.0;
# note that this problem is actually mix-determined
# since the sum of all the horizontal ray traveltimes
# equals the sum of all the vertical ray traveltimes
# so use the damped minimum-length solution when
# computing the solution
epsi = 0.0001;
epsi2=epsi**2;
GMG = np.matmul(G.T, la.inv(np.matmul(G,G.T)+epsi2*np.identity(N)));
# compute the complete resolution matrix for comparison
# note that R is symmetric in this case: R = GMG*G =
# G'/(G*G'+(epsi^2)*eye(N,N)) * G is symmetric, since
# the inverse of a symmetric matrix is symmetric
Rres = np.matmul(GMG,G);
# pull out just one column of the resolution matrix, one
# corresponding to a true model parameter near the middle of
# the model volume, at (ix, iy)
ix=floor(Lx/2)-1;
jy=floor(Lx/2)-1;
icolB=kofij[ix,jy];
RicolB=np.zeros((Lx,Ly));
for k in range(M):
RicolB[iofk[k,0],jofk[k,0]]=Rres[k,icolB];
# now construct just that column
# model parameter with unity in that col
mk = np.zeros((M,1));
mk[icolB,0] = 1.0;
# data it predicts
dk=np.matmul(G,mk);
# solve inverse problem, interpret the result as
# a column of the resolution matrix.
# rk = GT x;
rk = np.matmul(GMG,dk);
# reorganize to 2D physical model space
Rk=np.zeros((Lx,Ly));
for k in range(M):
Rk[iofk[k,0],jofk[k,0]]=rk[k,0];
gda_draw('title col from whole R',RicolB,'title col alone',Rk);
gdapy05_08
In [9]:
# gdapy05_09
# checkerboard test, same tomography problem as previous script
print("gdapy05_09");
# grid of unknowns is Lx by Ly
Lx = 20;
Ly = 20;
M = Lx*Ly;
# observations only along rows and columns
N=Lx+Ly;
# build index tables for convenience
iofk=np.zeros((M,1),dtype=int); # backward index, i(l)
jofk=np.zeros((M,1),dtype=int); # backward index, j(k)
kofij=np.zeros((Lx,Ly),dtype=int); # forward index, k(i,j)
k=0; # counter
for i in range(Lx):
for j in range(Ly):
iofk[k,0]=i;
jofk[k,0]=j;
kofij[i,j]=k;
k=k+1;
# true model parameters are checkerboard
mtrue = np.zeros((M,1));
for i in range(0,Lx,4):
for j in range(0,Ly,4):
k=kofij[i,j];
mtrue[k,0]=1.0;
# observations across rows
G=np.zeros((N,M));
for i in range(Lx):
for j in range(Ly):
k=kofij[i,j]; # (ix,iy) maps into scalar index j
G[i,k]=1.0;
# observations across columns
for j in range(Ly):
for i in range(Lx):
k=kofij[i,j]; # (ix,iy) maps into scalar index j
G[j+Lx,k]=1.0;
dtrue = np.matmul(G,mtrue);
# Generalized inverse
epsi = 0.0001;
epsi2=epsi**2;
GMG = np.matmul(G.T, la.inv(np.matmul(G,G.T)+epsi2*np.identity(N)));
# estimated data
mest=np.matmul(GMG,dtrue);
Strue=np.zeros((Lx,Ly)); # true checkerboard
Sest =np.zeros((Lx,Ly)); # recovered checkerboard
for k in range(M):
Strue[iofk[k,0],jofk[k,0]]=mtrue[k,0];
Sest[iofk[k,0],jofk[k,0]]=mest[k,0];
print("checkerboard test");
gda_draw('title Strue',Strue,'title Sest',Sest);
gdapy05_09 checkerboard test
In [10]:
# gdapy05_10
# Backus-Gilbert matrix S for 2D weight function
print("gdapy05_09");
# (x,y) grid
Lx=101;
xmin=0.0;
xmax=xmin+Lx-1;
x = gda_cvec(np.linspace(xmin,xmax,Lx));
Ly=101;
ymin=0.0;
ymax=ymin+Ly-1;
y = gda_cvec(np.linspace(ymin,ymax,Ly));
M = Lx*Ly;
# build index tables for unfolding/refolding matrix to vector
iofk=np.zeros((M,1),dtype=int); # backward index, i(k)
jofk=np.zeros((M,1),dtype=int); # backward index, j(k)
kofij=np.zeros((Lx,Ly),dtype=int); # forward index, k(i,j)
k=0; # counter
for i in range(Lx):
for j in range(Ly):
iofk[k,0]=i;
jofk[k,0]=j;
kofij[i,j]=k;
k=k+1;
# coordinates of model parameters
xm = np.zeros((M,1));
ym = np.zeros((M,1));
for i in range(Lx):
for j in range(Ly):
k = kofij[i,j];
xm[k,0] = x[i,0];
ym[k,0] = y[j,0];
# random data kernel
N=10;
G = np.random.normal(loc=0.0,scale=1.0,size=(N,M));
k=kofij[50,30]; # target row of S
wk=np.power(xm-xm[k,0],2)+np.power(ym-ym[k,0],2); # weights
S = np.matmul(G,np.matmul(np.diag(wk.ravel()),G.T));
gda_draw('title S',S);
# display weights on (x,y) grid
W = np.zeros((Lx,Ly));
for k in range(M):
W[iofk[k,0],jofk[k,0]] = wk[k,0];
gda_draw('title W',W);
gdapy05_09
In [ ]: