In [36]:
# gdapy04_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
# erroneous call to sum, hould be np.sum
# 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 [37]:
# gdapy04_01
# least squares fit to synthetic data
print("gdapy04_01");
# z-axis of randomly-chosen points
N=30;
zmin=0.0;
zmax=10.0;
# randomly chosen zs
z = np.sort(np.random.uniform(low=zmin,high=zmax,size=(N,1)));
# synthetic data
# d = a + b*z + noise
a=2.0;
b=1.0;
sd=1.0;
sd2=sd**2;
dtrue = a+b*z;
dobs = dtrue+np.random.normal(loc=0.0,scale=sd,size=(N,1));
# least squares fit
# construct data kernel
M=2;
G = np.zeros((N,M));
G[0:N,0:1] = np.ones((N,1));
G[0:N,1:2] = z;
# alternatively, this single commands can be used
# G = np.concatenate( (np.ones((N,1)), z), axis=1);
# least squares solution
GTG = np.matmul(G.T,G);
mest = la.solve(GTG,np.matmul(G.T,dobs));
covm = sd2 * la.inv(GTG);
# predicted data
dpre = np.matmul(G,mest);
e = dobs - dpre;
emax = np.max(np.abs(e));
iemax = np.argmax(np.abs(e));
# plot
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
# plot scale
pdmin=0.0;
pdmax=15.0;
# plot of data & straight line fit
ax1 = plt.subplot(1, 2, 1);
plt.axis( [zmin, zmax, pdmin, pdmax ] );
plt.plot( z, dobs, "ro", lw=2);
plt.plot( z, dpre, "b-", lw=2);
plt.xlabel("z");
plt.ylabel("d");
# plot illustrating definition of prediction error
ax2 = plt.subplot(1, 2, 2);
plt.axis( [zmin, zmax, pdmin, pdmax ] );
i=iemax;
plt.plot( z[i,0], dobs[i,0], "ro", lw=2);
plt.plot( z, dpre, "b-", lw=2);
plt.plot( gda_cvec(z[i,0],z[i,0]), gda_cvec(pdmin,dpre[i,0]), "k:", lw=2);
plt.plot( gda_cvec(zmin,z[i,0]), gda_cvec(dpre[i,0],dpre[i,0]), "k:", lw=2);
plt.plot( gda_cvec(zmin,z[i,0]), gda_cvec(dobs[i,0],dobs[i,0]), "k:", lw=2);
plt.xlabel('z');
plt.ylabel('d');
plt.show();
gdapy04_01
In [38]:
# gdapy04_02
# comparison of norms
print("gdapy04_02");
# z's
N=31;
zmin=0.0;
zmax=10.0;
Dz=(zmax-zmin)/(N-1);
z=gda_cvec(np.linspace(zmin,zmin+Dz*(N-1),N));
# random errors
e = np.random.uniform(low=-1,high=1,size=(N,1));
# variery of powers of errors
e1 = np.abs(e);
E1 = np.sum(e1);
e2 = np.power(np.abs(e),2);
E2 = sqrt(np.sum(e1));
e10 = np.power(np.abs(e),10);
E10 = np.sum(e10)**0.1;
print("Errors E1 %f E2 %f E10 %f" % (E1, E2, E10) );
# plot of errors
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
# plot scale
pemin=-1;
pemax=1;
# plot error
ax1 = plt.subplot(4, 1, 1);
plt.axis( [zmin, zmax, pemin, pemax ] );
# improvise bar chart
for i in range(N-1):
plt.plot( gda_cvec(z[i,0],z[i,0],z[i+1,0],z[i+1,0]),
gda_cvec(0.0, e[i,0], e[i,0], 0.0), "k-", lw=2 );
plt.xlabel("z");
plt.ylabel("e");
plt.plot( gda_cvec(zmin,zmax), gda_cvec(0.0,0.0), "k-", lw=2 );
# plot |error|
ax2 = plt.subplot(4, 1, 2);
plt.axis( [zmin, zmax, pemin, pemax] );
# improvise bar chart
for i in range(N-1):
plt.plot( gda_cvec(z[i,0],z[i,0],z[i+1,0],z[i+1,0]),
gda_cvec(0.0, e1[i,0], e1[i,0], 0.0), "k-", lw=2 );
plt.xlabel("z");
plt.ylabel("|e|");
plt.plot( gda_cvec(zmin,zmax), gda_cvec(0.0,0.0), "k-", lw=2 );
# plot |error|^2
ax3 = plt.subplot(4, 1, 3);
plt.axis( [zmin, zmax, pemin, pemax] );
# improvise bar chart
for i in range(N-1):
plt.plot( gda_cvec(z[i,0],z[i,0],z[i+1,0],z[i+1,0]),
gda_cvec(0.0, e2[i,0], e2[i,0], 0.0), "k-", lw=2 );
plt.xlabel("z");
plt.ylabel("|e|**2");
plt.plot( gda_cvec(zmin,zmax), gda_cvec(0.0,0.0), "k-", lw=2 );
# plot |error|**10
ax4 = plt.subplot(4, 1, 4);
plt.axis( [zmin, zmax, pemin, pemax] );
# improvise bar chart
for i in range(N-1):
plt.plot( gda_cvec(z[i,0],z[i,0],z[i+1,0],z[i+1,0]),
gda_cvec(0.0, e10[i,0], e10[i,0], 0.0), "k-", lw=2 );
plt.xlabel("z");
plt.ylabel("|e|**10");
plt.plot( gda_cvec(zmin, zmax), gda_cvec(0.0,0.0), "k-", lw=2 );
plt.show();
gdapy04_02 Errors E1 14.081021 E2 3.752469 E10 1.082987
In [39]:
# gdapy04_03
# Ln fits of straight line to to synthetic data
print("gdapy04_03");
print("Note: This script is pretty slow.");
# z-axis of randomly-chosen points
N=30;
zmin=0.0;
zmax=10.0;
# randomly chosen zs
z = np.sort(np.random.uniform(low=zmin,high=zmax,size=(N,1)));
# linear model d = a + b*z + noise
a=2.0;
b=1.0;
sd=0.5;
dobs = a+b*z+np.random.normal(loc=0.0,scale=sd,size=(N,1));
# one terrible outlier
z[N-1,0]=zmax;
dobs[N-1,0]=4.0;
# grid (a=intercept, b=slope)
Na=1001;
amin=-20;
amax=20;
Da = (amax-amin)/(Na-1);
Nb=1001;
bmin=-6.0;
bmax=6.0;
Db = (bmax-bmin)/(Nb-1);
# populate grid with errors
E1 = np.zeros((Na,Nb));
E2 = np.zeros((Na,Nb));
Einf = np.zeros((Na,Nb));
for i in range(Na):
for j in range(Nb):
ao = amin+Da*i;
bo = bmin+Db*j;
dpre = ao + bo*z;
e = dobs-dpre;
abse = np.abs(e);
E1[i,j]=np.sum(abse);
E2[i,j]=np.sum(np.power(abse,2));
Einf[i,j]=np.sum(np.power(abse,20)); # cheating; using large but finite power
# Important! Examine the error surface to make sure that
# the minimum is within the grid!
print("error surfaces");
gda_draw("title L1",np.log10(E1)," ",'title L2',np.log10(E2)," ","title Linf",np.log10(Einf));
# find minimum error
# find the (r1,c1) for which E1(r,c) is minimum
k = np.argmin( E1, axis=0 );
Etemp = np.min( E1, axis=0 );
c1 = np.argmin(Etemp);
r1=k[c1];
E1min = E1[r1,c1]; # minimum error
# predict data
a1 = amin+Da*r1;
b1 = bmin+Db*c1;
dpre1 = a1 + b1*z;
# find the (r2,c2) for which E2(r2,c2) is minimum
k = np.argmin( E2, axis=0 );
Etemp = np.min( E2, axis=0 );
c2 = np.argmin(Etemp);
r2=k[c2];
E2min = E2[r2,c2]; # minimum error
# predict data
a2 = amin+Da*r2;
b2 = bmin+Db*c2;
dpre2 = a2 + b2*z;
# find the (rinf,cinf) for which Einf(rinf,cinf) is minimum
k = np.argmin( Einf, axis=0 );
Etemp = np.min( Einf, axis=0 );
cinf = np.argmin(Etemp);
rinf=k[cinf];
Einfmin = E2[rinf,cinf]; # minimum error
# predict data
ainf = amin+Da*rinf;
binf = bmin+Db*cinf;
dpreinf = ainf + binf*z;
# write our intecept & slope
print("1: a %.3f b %.3f" % (a1, b1) );
print("2: a %.3f b %.3f" % (a2, b2) );
print("inf: a %.3f b %.3f" % (ainf, binf) );
# plot results
fig2 = plt.figure(2,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
# plot bounds
pdmin=0.0;
pdmax=15.0;
# plot data and fits
plt.axis( [zmin, zmax, pdmin, pdmax ] );
plt.plot( z, dobs, "ko", lw=2);
plt.plot( z, dpre1, "r-", lw=2);
plt.plot( z, dpre2, "g-", lw=2);
plt.plot( z, dpreinf, "b-", lw=2);
plt.xlabel("z");
plt.ylabel("d");
plt.show();
gdapy04_03 Note: This script is pretty slow. error surfaces
1: a 1.800 b 1.020 2: a 2.240 b 0.888 inf: a 4.760 b 0.336
In [40]:
# gdapy04_04
# comparison of long-tailed and short tailed pdf's
print("gdapy04_04");
# d-axis
Dd = 0.1;
N = 101;
d = gda_cvec(np.linspace(0.0,Dd*(N-1),N));
dmin=d[0,0];
dmax=d[N-1,0];
# short tailed pdf: Normal pdf
dbar=dmax/2.0;
d2=np.power(d-dbar,2);
sd = 1.0;
sd2 = sd**2;
p1 = np.exp(-0.5*d2/sd2)/(sqrt(2*pi)*sd);
A1 = Dd*np.sum(p1);
# long-tailed distribution: Cauchy–Lorentz distribution
g = 1.0;
g2 = g**2;
p2 = np.divide(1.0,(pi*g*(1.0+(d2/g2))));
A2 = Dd*np.sum(p2);
print("check on areas: true %.4f est1 %.4f est2 %.4f" % (1.0, A1, A2));
print("(but missing some tail on long-tailed pdf)");
fig2 = plt.figure(2,figsize=(12,5)); # smallish figure
# plot long-tailed pdf
# long-tailed pdf
plt.subplot(1,2,2);
plt.axis( [dmin, dmax, 0, 0.5 ] );
plt.plot(d,p1,"r-",lw=3);
plt.xlabel("d");
plt.ylabel("p(d)");
# short-tailed pdf
plt.subplot(1,2,1);
plt.axis( [dmin, dmax, 0.0, 0.5 ] );
plt.plot(d,p2,"b-",lw=3);
plt.xlabel("d");
plt.ylabel("p(d)");
plt.show();
gdapy04_04 check on areas: true 1.0000 est1 1.0000 est2 0.8756 (but missing some tail on long-tailed pdf)
In [41]:
# gdapy04_05
# Kepler's 3rd law example, period^2 = radius^3
print("gdapy04_05");
# read data
D = np.genfromtxt("../data/planetary.txt", delimiter="\t");
N, K = np.shape(D);
radius=D[0:N,0:1]/1.0e9; # work inunits of 10^9 km
period=D[0:N,1:2]/1000.0; # and units of 10^3 days
z = period; # shorthand for auxiliary variable, period
# take radius^3 to be the observation
# and period to be the auxillary variable
# (which is ok, since time can be measured so accurately)
dobs = np.power(radius,3);
M=3;
G = np.zeros((N,M));
G[0:N,0:1] = np.ones((N,1));
G[0:N,1:2] = z;
G[0:N,2:3] = np.power(z,2);
# least squares solution
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mest = la.solve(GTG,GTd);
# predicted data
dpre = np.matmul(G,mest);
# error
e = dobs-dpre;
E = np.matmul(e.T,e); E=E[0,0];
sd2 = E/(N-M); # posterior variance
sd = sqrt(sd2);
# covariance and confidence limits
covm = sd2*la.inv(GTG);
sm0 = sqrt(covm[0,0]);
sm1 = sqrt(covm[1,1]);
sm2 = sqrt(covm[2,2]);
print("mest0: %.4f +/- %.4f (95%%)" % (mest[0,0], 2.0*sm0) );
print("mest1: %.4f +/- %.4f (95%%)" % (mest[1,0], 2.0*sm1) );
print("mest2: %.4f +/- %.4f (95%%)" % (mest[2,0], 2.0*sm2) );
# plot smooth parabola, so use lots of z's
periodeval=gda_cvec(np.linspace(0.0,100.0,101));
deval=mest[0,0]+mest[1,0]*periodeval+mest[2,0]*np.power(periodeval,2);
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1=plt.subplot(2,1,1);
plt.axis( [0, 100, 0, 250 ]);
plt.plot(period,dobs,"ro",lw=3);
plt.plot(periodeval,deval,"k-",lw=3);
plt.xlabel('period, Kdays');
plt.ylabel('radius^3, Tm^3');
ax2=plt.subplot(2,1,2);
plt.axis( [0, 100, -0.5, 0.5 ]);
plt.plot(period,e,"ro",lw=3);
plt.plot( gda_cvec(0.0, 100), gda_cvec(0.0, 0.0), "k-",lw=2);
plt.xlabel('period, Kdays');
plt.ylabel('error, Tm^3');
plt.show();
gdapy04_05 mest0: -0.0253 +/- 0.0913 (95%) mest1: 0.0100 +/- 0.0090 (95%) mest2: 0.0250 +/- 0.0001 (95%)
In [44]:
# gdapy04_06
#
# planar fit to depths of earthquakes in the Kurile subduction zone
print("gdapy04_06");
# read data
D = np.genfromtxt("../data/kurile_eqs.txt", delimiter="\t");
N, K = np.shape(D);
lat=D[0:N,0:1];
lon=D[0:N,1:2];
depth=D[0:N,2:3];
# convert to km
x = 111.12*np.cos((pi/180)*np.mean(lat))*(lon-np.min(lon));
y = 111.12*(lat-np.min(lat));
dobs=(-depth);
z = depth;
# bounds
xmin = 0.0;
xmax = 1200.0;
ymin = 0.0;
ymax = 1200.0;
zmin = -700.0;
zmax = 0.0;
# three-d plot of data
fig1 = plt.figure(1,figsize=(10,10));
ax1 = fig1.add_subplot(111, projection='3d')
ax1.axes.set_xlim3d(left=xmin, right=xmax);
ax1.axes.set_ylim3d(bottom=ymin, top=ymax);
ax1.axes.set_zlim3d(bottom=zmin, top=zmax);
ax1.set_xlabel("x");
ax1.set_ylabel("y");
ax1.set_zlabel("z");
# edapy_08_21, 3D plot of clusters in Atlantic Rocks dataset
# solve inverse problem of fitting a plane to data
M=3;
G = np.zeros((N,M));
G[0:N,0:1] = np.ones((N,1));
G[0:N,1:2] = x;
G[0:N,2:3] = y;
# solve inverse problem
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mest = la.solve(GTG,GTd);
# predict data
dpre=np.matmul(G,mest);
# display fit as mesh in 3D space
K=31;
xx=gda_cvec(np.linspace(xmin,xmax,K));
yy=gda_cvec(np.linspace(ymin,ymax,K));
X,Y=np.meshgrid(xx, yy);
Z = mest[0,0]+mest[1,0]*X+mest[2,0]*Y;
# clip
CLIPMETHOD=1;
if(CLIPMETHOD):
zminv=zmin;
zmaxv=zmax;
for ix in range(K):
for iy in range(K):
if( Z[ix,iy]<zmin ):
Z[ix,iy]=zminv;
elif( Z[ix,iy]>zmax ):
Z[ix,iy]=zmaxv;
else:
Z = np.ma.masked_where(Z<=zmin, Z, copy=True);
Z = np.ma.masked_where(Z>=zmax, Z, copy=True);
surf = ax1.plot_surface(X, Y, Z, cmap=cm.jet, lw=0, antialiased=False)
ax1.scatter( xs=x, ys=y, zs=dobs, color=(0,0,0) );
ax1.view_init(1.0, 45.0);
# add colorbar
fig1.colorbar(surf, shrink=0.4, aspect=10)
plt.show();
# plot observed and predicted depths
fig2 = plt.figure(2,figsize=(7,7));
plt.axis( [-750.0, 0.0, -750.0, 0.0 ]);
plt.xlabel("observed depth (km)");
plt.ylabel("predicted depth (km)")
plt.plot(dobs,dpre,"ko",lw=2);
plt.plot( gda_cvec(-750.0,0.0), gda_cvec(-750.0,0.0), "r-", lw=1 );
plt.show();
gdapy04_06
In [45]:
# gdapy04_07
# Seismic Reflection example
print("gdapy04_07");
# note that functions are used to simplify the struture of the
# code. As an example of defining a function, consider a function
# "myfcn" that takes A, B and returns C=2A and D=2B. It is defined as:
# def myfcn(A,B):
# C = 2*A;
# D = 2*B;
# return C,D;
# funtion to compute "P Down" refl coefficients
def gda_PD( iP, r, a, b ):
# input
# iP: angle of incidence of P wave in degrees
# r: density, kg/m3
# a: compressional velocity, m/s
# b: shear velocity, m/s
# ouput
# Upgoing reflected P wave sensitivity, PDPUr, PDPUa, PDPUb
# Upgoing converted S wave sensitivity, PDSUr, PDSUa, PDSUb
# with respect to (Dr/r), (Da/a) and (Db/b)
# see AKi and Richards (2009, their equation 5.46)
# Warning: waves must not be evanescent
DTOR = pi/180.0;
# squares of velocities
b2 = b**2;
a2 = a**2;
# ray parameter, p
# p = sin(iP)/a
# sin2(iP)=a2*p2
# cos2(iP)=1-a2*p2
p = sin(DTOR*iP)/a;
p2 = p**2;
cosiP2 = 1.0-a2*p2;
cosiP = sqrt(cosiP2);
cosiPoa = cosiP/a;
# S wave has same ray parameter, p = sin(jS)/b
# sin(jS) = b*p;
cosjS2 = 1.0 - b2*p2;
cosjS = sqrt(cosjS2);
cosjSob = cosjS/b;
# upgoing P wave sensitivities
PDPUr = 0.5*(1-4*b2*p2);
PDPUa = 1.0/(2.0*cosiP2);
PDPUb = -4.0*b2*p2;
C = (-p*a)/(2*cosjS);
PDSUr = C*( 1.0 - 2.0*b2*p2 + 2.0*b2*cosiPoa*cosjSob );
PDSUa = 0.0;
PDSUb = C*(-4*b2)*(p2-cosiPoa*cosjSob);
return PDPUr, PDPUa, PDPUb, PDSUr, PDSUa, PDSUb;
# funtion to compute "P Down" refl coeffocients
def gda_SD( jS, r, a, b ):
# input
# jS: angle of incidence of S wave in degrees
# r: density, kg/m3
# a: compressional velocity, m/s
# b: shear velocity, m/s
# ouput
# Upgoing reflected P wave sensitivity, SDPUr, SDPUa, SDPUb
# Upgoing converted S wave sensitivity, SDSUr, SDSUa, SDSUb
# with respect to (Dr/r), (Da/a) and (Db/b)
# see AKi and Richards (2009, their equation 5.46)
# Warning: waves must not be evanescent
DTOR = pi/180;
# squares of velocities
b2 = b**2;
a2 = a**2;
# ray parameter, p
# p = sin(jS)/b
# sin2(jS)=a2*p2
# cos2(jS)=1-a2*p2
p = sin(DTOR*jS)/b;
p2 = p**2;
cosjS2 = 1.0-b2*p2;
cosjS = sqrt(cosjS2);
cosjSob = cosjS/b;
# P wave has same ray parameter, p = sin(iP)/a
# sin(iP) = a*p;
cosiP2 = 1.0 - a2*p2;
cosiP = sqrt(cosiP2);
cosiPoa = cosiP/a;
# upgoing P wave sensitivities
C = (-p*a)/(2*cosjS);
PDSUr = C*( 1.0 - 2.0*b2*p2 + 2.0*b2*cosiPoa*cosjSob );
PDSUa = 0.0;
PDSUb = C*(-4.0*b2)*(p2-cosiPoa*cosjSob);
D = (cosjS/a)*(b/cosiP);
SDPUr = D*PDSUr;
SDPUa = D*PDSUa;
SDPUb = D*PDSUb;
SDSUr = -0.5*(1-4*b2*p2);
SDSUa = 0.0;
SDSUb = -( (1.0/(2*cosjS2)) - 4.0*b2*p2 );
return SDPUr, SDPUa, SDPUb, SDSUr, SDSUa, SDSUb;
# setup q-axis (q = theta = ange of incidence)
Nq = 21;
qmin = 0.0;
qmax = 25.0;
Dq = (qmax-qmin)/(Nq-1);
q = gda_cvec(np.linspace(qmin, qmax, Nq));
# medium properties
a = 5000; # compressional velocity
b = 3000; # shear velocity
r = 2500; # density
# interface properties
Dror = 0.07; # normalized jump in compressional velocity
Daoa = 0.10; # normalized jump in shear velocity
Dbob = 0.11; # normalized jump in density velocity
mtrue = gda_cvec(Dror, Daoa, Dbob); # group into model parameters
N=4*Nq; # number of data, 4 reflection coefficients times Nq angles
M=3; # number of model parameters
G = np.zeros((N, M)); # data kernel
dtrue = np.zeros((N, 1)); # true data
d = np.zeros((N, 1)); # used for consistncy check
for iq in range(Nq): # loop over angles of incidence
PDPUr, PDPUa, PDPUb, PDSUr, PDSUa, PDSUb = gda_PD( q[iq,0], r, a, b );
PDPU = PDPUr*Dror + PDPUa*Daoa + PDPUb*Dbob;
PDSU = PDSUr*Dror + PDSUa*Daoa + PDSUb*Dbob;
SDPUr, SDPUa, SDPUb, SDSUr, SDSUa, SDSUb = gda_SD( q[iq,0], r, a, b );
SDPU = SDPUr*Dror + SDPUa*Daoa + SDPUb*Dbob;
SDSU = SDSUr*Dror + SDSUa*Daoa + SDSUb*Dbob;
G[iq:iq+1,0:M] = gda_cvec(PDPUr, PDPUa, PDPUb).T;
G[iq+Nq:iq+Nq+1,0:M] = gda_cvec(PDSUr, PDSUa, PDSUb).T;
G[iq+2*Nq:iq+2*Nq+1,0:M] = gda_cvec(SDPUr, SDPUa, SDPUb).T;
G[iq+3*Nq:iq+3*Nq+1,0:M] =gda_cvec(SDSUr, SDSUa, SDSUb).T;
d[iq,0] = PDPU;
d[iq+Nq,0] = PDSU;
d[iq+2*Nq,0] = SDPU;
d[iq+3*Nq,0] = SDSU;
# reflection coefficients
dtrue = np.matmul(G,mtrue);
# check consitency of d and G*mtrue true reflection coefficients
e = dtrue-d;
E = np.matmul(e.T,e); E=E[0,0];
print("consistency check on Gm=d: %e should be close to zero" % (E) );
print(" ");
# plot the reflection coefficents as a function of angle
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
plt.axis( [qmin, qmax, -0.1, 0.1] );
plt.plot( gda_cvec(qmin, qmax), gda_cvec(0,0), "k:", lw=2 );
plt.plot( q, dtrue[0:Nq,0:1], "k-", lw=3 );
plt.plot( q, dtrue[Nq:2*Nq,0:1], "r-", lw=3 );
plt.plot( q, dtrue[2*Nq:3*Nq,0:1], "g-", lw=3 );
plt.plot( q, dtrue[3*Nq:4*Nq], "b-", lw=3 );
# make synthetic noisy data and plot it
sigmad = 0.005;
dobs = dtrue + np.random.normal(loc=0.0,scale=sigmad,size=(N,1));
plt.plot( q, dobs[0:Nq,0:1], "k-", lw=2 );
plt.plot( q, dobs[Nq:2*Nq,0:1], "r-", lw=2 );
plt.plot( q, dobs[2*Nq:3*Nq,0:1], "g-", lw=2 );
plt.plot( q, dobs[3*Nq:4*Nq,0:1], "b-", lw=2 );
# least squares estimate of the layer interface properties
GTG = np.matmul(G.T,G);
GTd= np.matmul(G.T,dobs);
GTGinv = la.inv(GTG);
mest = la.solve(GTG,GTd); # least squares soln
dpre = np.matmul(G,mest); # predicted data
e = dobs-dpre; # error vector
E = np.matmul(e.T,e); E=E[0,0]; # least squares error
sigmadest2 = E/(N-M); # posterior variance
sigmadest = sqrt(sigmadest2);
# print covariance matrix
print("sigmad true %.4f est %.4f" % (sigmad, sigmadest));
print(" ");
covm = sigmadest2*GTGinv; # covariance of the estimates
print("cov m");
print("%.5f %.5f %.5f" % (covm[0,0], covm[0,1], covm[0,2]) );
print("%.5f %.5f %.5f" % (covm[1,0], covm[1,1], covm[1,2]) );
print("%.5f %.5f %.5f" % (covm[2,0], covm[2,1], covm[2,2]) );
print(" ");
sigmam1 = sqrt(covm[0,0]);
sigmam2 = sqrt(covm[1,1]);
sigmam3 = sqrt(covm[2,2]);
print("m0: true %.2f est %.2f +/- %.2f (95%%)" % (mtrue[0,0], mest[0,0], 2*sigmam1) );
print("m1: true %.2f est %.2f +/- %.2f (95%%)" % (mtrue[1,0], mest[1,0], 2*sigmam2) );
print("m2: true %.2f est %.2f +/- %.2f (95%%)" % (mtrue[2,0], mest[2,0], 2*sigmam3) );
plt.plot( q, dpre[0:Nq,0:1], "ko", lw=2 );
plt.plot( q, dpre[Nq:2*Nq,0:1], "ro", lw=2 );
plt.plot( q, dpre[2*Nq:3*Nq,0:1], "go", lw=2 );
plt.plot( q, dpre[3*Nq:4*Nq,0:1], "bo", lw=2 );
plt.xlabel("theta");
plt.ylabel("reflection coefficient data");
plt.show();
fig2 = plt.figure(2,figsize=(12,5)); # smallish figure
plt.subplot(1,3,1);
plt.axis( [0.0, 1.0, 0.0, 0.25] );
plt.plot( gda_cvec(0.0, 1.0), gda_cvec(mtrue[0,0],mtrue[0,0]), "k-", lw=5 );
plt.plot( gda_cvec(0.0, 1.0), gda_cvec(mest[0,0],mest[0,0]), "r-", lw=3 );
plt.plot( gda_cvec(0.5, 0.5), gda_cvec(mest[0,0]-2*sigmam1,mest[0,0]+2*sigmam1), "r-", lw=3 );
plt.title("Dror");
plt.subplot(1,3,2);
plt.axis( [0, 1, 0, 0.25] );
plt.plot( gda_cvec(0.0, 1.0), gda_cvec(mtrue[1,0],mtrue[1,0]), "k-", lw=5 );
plt.plot( gda_cvec(0.0, 1.0), gda_cvec(mest[1,0],mest[1,0]), "r-", lw=3 );
plt.plot( gda_cvec(0.5, 0.5), gda_cvec(mest[1,0]-2*sigmam2,mest[1,0]+2*sigmam2), "r-", lw=3 );
plt.title("Daoa");
plt.subplot(1,3,3);
plt.axis( [0, 1, 0, 0.25] );
plt.plot( gda_cvec(0.0, 1.0), gda_cvec(mtrue[2,0],mtrue[2,0]), "k-", lw=5 );
plt.plot( gda_cvec(0.0, 1.0), gda_cvec(mest[2,0],mest[2,0]), "r-", lw=3 );
plt.plot( gda_cvec(0.5, 0.5), gda_cvec(mest[2,0]-2*sigmam3,mest[2,0]+2*sigmam3), "r-", lw=3 );
plt.title("Dbob");
plt.show();
gdapy04_07 consistency check on Gm=d: 0.000000e+00 should be close to zero sigmad true 0.0050 est 0.0050 cov m 0.00007 -0.00007 -0.00007 -0.00007 0.00008 0.00008 -0.00007 0.00008 0.00009 m0: true 0.07 est 0.09 +/- 0.02 (95%) m1: true 0.10 est 0.08 +/- 0.02 (95%) m2: true 0.11 est 0.09 +/- 0.02 (95%)
In [46]:
# gdapy04_08
# simple exampe of setting up and solving a sparse problem
# Fm=h where F is an anti-identity matrix
print("gdapy04_08");
# F is L by M
L=5;
M=5;
# Define row-index, column-index, value arrays for sparse matrix
P=0; # element counter
Pmax = 100; # maximum number of elements
Fr = np.zeros((Pmax,1),dtype=int);
Fc = np.zeros((Pmax,1),dtype=int);
Fv = np.zeros((Pmax,1));
# create the element arrays corresponding to an anti-identity matrix
for i in range(L): # row i
j=M-i-1; # column j
Fr[P,0]=i; # row index
Fc[P,0]=j; # column index
Fv[P,0]=1.0; # value
P=P+1; # increment element number
print("the element arrays");
print("P: ",P);
print("rF cF vF");
print(np.concatenate((Fr[0:P,0:1],Fc[0:P,0:1],Fv[0:P,0:1]),axis=1));
print(" ");
gdaFsparse=sparse.coo_matrix((Fv[0:P,0],
(Fr[0:P,0],Fc[0:P,0])),shape=(L,M));
del Fr; # delete no longer needed arrays
del Fc;
del Fv;
print("the anti-identity matrix gdaFsparse");
print(gdaFsparse.todense());
print(" ");
# true model parameters are [0, 1, 2, 3, 4]
mtrue = gda_cvec(np.linspace(0.0,1.0*(M-1),M));
# true h from
htrue = gdaFsparse*mtrue;
hobs = htrue;
print("hobs");
print(hobs);
print(" ");
# define linear operator needed for biconjugate gradient solver
LO=las.LinearOperator(shape=(M,M),
matvec=gdaFTFmul,rmatvec=gdaFTFmul);
# solve using conjugate gradient algorithm
tol = 1e-12; # tolerance
maxit = 3*(L+M); # maximum number of iterations allowed
FTh = gdaFsparse.T*hobs;
q=las.bicg(LO,FTh,rtol=tol,maxiter= maxit);
mest = gda_cvec(q[0]);
print("mtrue dtrue mest");
print(np.concatenate((mtrue,htrue,mest),axis=1));
print(" ");
e = mtrue-mest
E = np.matmul(e.T,e);
print("error in solution: ", E);
gdapy04_08 the element arrays P: 5 rF cF vF [[0. 4. 1.] [1. 3. 1.] [2. 2. 1.] [3. 1. 1.] [4. 0. 1.]] the anti-identity matrix gdaFsparse [[0. 0. 0. 0. 1.] [0. 0. 0. 1. 0.] [0. 0. 1. 0. 0.] [0. 1. 0. 0. 0.] [1. 0. 0. 0. 0.]] hobs [[4.] [3.] [2.] [1.] [0.]] mtrue dtrue mest [[0. 4. 0.] [1. 3. 1.] [2. 2. 2.] [3. 1. 3.] [4. 0. 4.]] error in solution: [[0.]]
In [47]:
# gdapy04_09
# Example of building a sparse matrix row-wise, where each row is created
# as a full-vector and then then only the non-zero elements
# transferred to the sparse matrix. This strategy can be
# used when it's difficult to figure out in advance the number
# of non-zero elements in the row
print("gdapy04_09");
# F is LxM
L=6;
M=6;
# Define row-index, column-index, value arrays for sparse matrix
P=0; # element counter
Pmax = 100; # maximum number of elements
Fr = np.zeros((Pmax,1),dtype=int);
Fc = np.zeros((Pmax,1),dtype=int);
Fv = np.zeros((Pmax,1));
# process matric row-wise
for i in range(L): # row i
Frow = np.zeros((M,1)); # one row of matrix
# In this example, I cook up a matrix that is the
# an identify matrix modified by by setting either
# one or two randomly chosen columns to a random number
Frow[i,0]=1.0; # identity matrix
Ncols = np.random.randint(low=1,high=2);
for j in range(Ncols):
irow = np.random.randint(low=0,high=M);
vrow = np.random.normal(loc=0.0,scale=1.0);
Frow[irow,0]=vrow;
# process row, transferring non-zero values to the element arrays
for j in range(M):
if( Frow[j,0] != 0.0 ):
Fr[P,0]=i; # row index
Fc[P,0]=j; # column index
Fv[P,0]=Frow[j,0]; # value
P=P+1; # increment element number
print("P: ",P);
print("rF cF vF");
print(np.concatenate((Fr[0:P,0:1],Fc[0:P,0:1],Fv[0:P,0:1]),axis=1));
print(" ");
gdaFsparse=sparse.coo_matrix((Fv[0:P,0],
(Fr[0:P,0],Fc[0:P,0])),shape=(L,M));
del Fr; # delete no longer needed arrays
del Fc;
del Fv;
print("The F matrix gdaFsparse");
print(gdaFsparse.todense());
print(" ");
mtrue = gda_cvec(np.linspace(0.0,1.0*(L-1),L));
htrue = gdaFsparse*mtrue;
hobs = htrue;
print("hobs");
print(hobs);
print(" ");
# define linear operator needed for biconjugate gradient solver
LO=las.LinearOperator(shape=(M,M),
matvec=gdaFTFmul,rmatvec=gdaFTFmul);
# solve using conjugate gradient algorithm
tol = 1e-12; # tolerance
maxit = 3*(L+M); # maximum number of iterations allowed
FTh = gdaFsparse.T*hobs;
q=las.bicg(LO,FTh,rtol=tol,maxiter= maxit);
mest = gda_cvec(q[0]);
print("mtrue htrue mest");
print(np.concatenate((mtrue,htrue,mest),axis=1));
print(" ");
e = mtrue-mest;
E = np.matmul(e.T,e);
print("error in solution: ", E);
gdapy04_09 P: 11 rF cF vF [[ 0. 0. 1. ] [ 0. 4. -1.89709761] [ 1. 1. 0.03935712] [ 2. 2. 1. ] [ 2. 4. 0.14314966] [ 3. 3. 1. ] [ 3. 4. 2.60304265] [ 4. 1. -0.48255317] [ 4. 4. 1. ] [ 5. 2. 1.11457505] [ 5. 5. 1. ]] The F matrix gdaFsparse [[ 1. 0. 0. 0. -1.89709761 0. ] [ 0. 0.03935712 0. 0. 0. 0. ] [ 0. 0. 1. 0. 0.14314966 0. ] [ 0. 0. 0. 1. 2.60304265 0. ] [ 0. -0.48255317 0. 0. 1. 0. ] [ 0. 0. 1.11457505 0. 0. 1. ]] hobs [[-7.58839042] [ 0.03935712] [ 2.57259864] [13.41217059] [ 3.51744683] [ 7.22915011]] mtrue htrue mest [[ 0.00000000e+00 -7.58839042e+00 -9.51543693e-15] [ 1.00000000e+00 3.93571176e-02 1.00000000e+00] [ 2.00000000e+00 2.57259864e+00 2.00000000e+00] [ 3.00000000e+00 1.34121706e+01 3.00000000e+00] [ 4.00000000e+00 3.51744683e+00 4.00000000e+00] [ 5.00000000e+00 7.22915011e+00 5.00000000e+00]] error in solution: [[1.20218798e-26]]
In [48]:
# gdapy04_10
#
# straight line problem
# least squares fit to synthetic data
# using two solution methods
# standard solver
# biconjugate gradient solver
print("gdapy04_10");
# z-axis of randomly-chosen points
N=30;
zmin=0.0;
zmax=10.0;
# randomly chosen zs
z = np.sort(np.random.uniform(low=zmin,high=zmax,size=(N,1)));
# synthetic data
# d = a + b*z + noise
a=2.0;
b=1.0;
print(" true: intercept %.4f slope: %.4f" % (a, b));
sd=1.0;
sd2=sd**2;
dtrue = a+b*z;
dobs = dtrue+np.random.normal(loc=0.0,scale=sd,size=(N,1));
# least squares fit
# construct full data kernel
M=2;
G = np.zeros((N,M));
G[0:N,0:1] = np.ones((N,1));
G[0:N,1:2] = z;
# least squares solution
GTG = np.matmul(G.T,G);
mest1 = la.solve(GTG,np.matmul(G.T,dobs));
print("non-sparse: intercept %.4f slope: %.4f" % (mest1[0,0], mest1[1,0]));
# predicted data
dpre1 = np.matmul(G,mest1);
# the sparse code uses a nomenclature where
# the data kernel is called edaFsparse and is
# LxM in size with P non-zero elements. In
# this case, edaFsparse, is not really sparse
# but the method work fine, anyway
# construce spare data kernel
# Define row-index, column-index, value arrays for sparse matrix
L = N;
P=0; # element counter
Pmax = N*M+2; # maximum number of elements
Fr = np.zeros((Pmax,1),dtype=int);
Fc = np.zeros((Pmax,1),dtype=int);
Fv = np.zeros((Pmax,1));
# build element arrays for sparse matrix
for i in range(N): # row i
Fr[P,0]=i; # row index
Fc[P,0]=0; # column index
Fv[P,0]=1.0; # value
P=P+1; # increment element number
Fr[P,0]=i; # row index
Fc[P,0]=1; # column index
Fv[P,0]=z[i,0]; # value
P=P+1; # increment element number
# build sparse matrix
gdaFsparse=sparse.coo_matrix((Fv[0:P,0],
(Fr[0:P,0],Fc[0:P,0])),shape=(L,M));
del Fr; # delete no longer needed arrays
del Fc;
del Fv;
# define linear operator needed for biconjugate gradient solver
LO=las.LinearOperator(shape=(M,M),
matvec=gdaFTFmul,rmatvec=gdaFTFmul);
# solve using conjugate gradient algorithm
tol = 1e-12; # tolerance
maxit = 3*(L+M); # maximum number of iterations allowed
FTd = gdaFsparse.T*dobs;
q=las.bicg(LO,FTd,rtol=tol,maxiter= maxit);
mest2 = gda_cvec(q[0]);
print(" sparse: intercept %.4f slope: %.4f" % (mest2[0,0], mest2[1,0]));
# predicted data
dpre2 = gdaFsparse*mest2;
# plot
fig2 = plt.figure(2,figsize=(12,5)); # smallish figure
pdmin=0.0;
pdmax=15.0;
plt.axis( [zmin, zmax, pdmin, pdmax ] );
plt.plot( z, dobs, "ko", lw=2);
plt.plot( z, dpre1, "r-", lw=3);
plt.plot( z, dpre2, "g--", lw=2);
plt.xlabel("z");
plt.ylabel("d");
plt.show();
gdapy04_10
true: intercept 2.0000 slope: 1.0000
non-sparse: intercept 2.7235 slope: 0.9372
sparse: intercept 2.7235 slope: 0.9372
In [49]:
# gdapy04_11
# fill in missing observations in times series
# using prior inforamerion of smoothness
print("gdapy04_11");
# trength of prior information
epsi = 1.e-2;
# x-axis
N=20;
M=101;
Dx = 1.0;
a = epsi/Dx;
Dx2 = Dx**2;
b = epsi/Dx2;
x = gda_cvec(np.linspace(0.0,Dx*(M-1),M))
xmin=x[0,0];
xmax=x[M-1,0];
# true model is a sine wave
mtrue = np.sin(3.0*pi*x/xmax);
# select N uniue random points (xd,dtrue) of data
tmp=np.random.permutation(np.linspace(0,M-1,M).astype(int));
nd = np.sort(tmp[0:N]);
xd = x[nd,0:1];
dtrue = mtrue[nd,0:1];
dobs = dtrue;
# plot the data and the predicted line
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
plt.axis( [xmin, xmax, -1.1, 1.1 ] );
plt.plot(x,mtrue,"k-",lw=4);
plt.plot(xd,dtrue,"ro",lw=2);
plt.xlabel("x");
plt.ylabel("d and m");
# The array F is L=N+M by M
# the first N rows are the data
# the next M-2 rows are second differences of the interior
# the final two rows are first difference of the ends
L = N+M; # row in F, f
P=0;
Pmax = 3*L;
Fr = np.zeros((Pmax,1),dtype=int);
Fc = np.zeros((Pmax,1),dtype=int);
Fv = np.zeros((Pmax,1));
f = np.zeros((L,1));
nL = 0; # row counter
# N data rows; that is, matrix G
for i in range(N):
j=nd[i]; # j-th model parameter is associated with i-th wata
Fr[P,0]=nL; # row index
Fc[P,0]=j; # column index
Fv[P,0]=1.0; # element of G is identity matrix
P=P+1; # increment element number
f[nL,0]=dobs[i,0]; # data
nL=nL+1; # incrememt row-of-F counter
# M-2 second differnece rows; that is middle M-2 rows of H
for i in range(1,M-1):
j = i-1;
Fr[P,0]=nL;
Fc[P,0]=j;
Fv[P,0]=b;
P=P+1;
j = i;
Fr[P,0]=nL;
Fc[P,0]=j;
Fv[P,0]=-2.0*b;
P=P+1;
j = i+1;
Fr[P,0]=nL;
Fc[P,0]=j;
Fv[P,0]=b;
P=P+1;
f[nL,0]=0.0; # right hand side
nL=nL+1; # incrememt row-of-F counter
# 2 first derivative rows; that is top and bottom rows of H
j=0;
Fr[P,0]=nL; # row index
Fc[P,0]=j; # column index
Fv[P,0]=-a; # first difference operator
P=P+1; # increment element number
f[nL]=0.0; # right hand side
Fr[P,0]=nL; # row index
Fc[P,0]=j+1; # column index
Fv[P,0]=a; # first difference operator
P=P+1; # increment element number
f[nL,0]=0.0; # right hand side
nL=nL+1; # incrememt row-of-F counter
j=M-2;
Fr[P,0]=nL; # row index
Fc[P,0]=j; # column index
Fv[P,0]=-a; # first difference operator
P=P+1; # increment element number
f[nL]=0.0; # right hand side
Fr[P,0]=nL; # row index
Fc[P,0]=j+1; # column index
Fv[P,0]=a; # first difference operator
P=P+1; # increment element number
f[nL,0]=0.0; # right hand side
nL=nL+1; # incrememt row-of-F counter
print("After build of Fm=h");
print("Expected elements P=N+3*(M-2)+2*2 = ",N+3*(M-2)+2*2," got ",P);
print("Expected nL=N+M = ",N+M," got ", nL);
# Define row-index, column-index, value arrays for sparse matrix
gdaFsparse=sparse.coo_matrix((Fv[0:P,0],
(Fr[0:P,0],Fc[0:P,0])),shape=(L,M));
del Fr; # delete no longer needed arrays
del Fc;
del Fv;
# define linear operator needed for biconjugate gradient solver
LO=las.LinearOperator(shape=(M,M),
matvec=gdaFTFmul,rmatvec=gdaFTFmul);
# solve using conjugate gradient algorithm
tol = 1e-14; # tolerance
maxit = 3*(L+M); # maximum number of iterations allowed
FTf = gdaFsparse.T*f;
q=las.bicg(LO,FTf,rtol=tol,maxiter= maxit);
mest = gda_cvec(q[0]);
plt.plot(x,mest,"g-",lw=2);
gdapy04_11 After build of Fm=h Expected elements P=N+3*(M-2)+2*2 = 321 got 321 Expected nL=N+M = 121 got 121
In [50]:
# gdapy04_12
# Example of smoothing operators
# using the formula of Menke and Eilon (2015)
print("gdapy04_12");
# x-axis
N=501;
Dx = 0.02;
xmin = -N*Dx/2;
xmax = xmin + Dx*(N-1)
x = gda_cvec(np.linspace(xmin,xmax,N) );
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1=plt.subplot(2,1,1);
ep = 0.3;
epi = 1/ep;
m1 = (epi/2)*np.exp(-epi*np.abs(x));
amp = np.max(np.abs(m1));
plt.axis( [xmin, xmax, -1.1*amp, 1.1*amp ] );
plt.plot( gda_cvec(xmin, xmax), gda_cvec(0.0,0.0), "k:", lw=1 );
plt.plot( x, m1, "k-", lw=3 );
ep = 0.6;
epi = 1/ep;
m1 = (epi/2)*np.exp(-epi*np.abs(x));
plt.plot( x, m1, "r-", lw=3 );
plt.xlabel("x");
plt.ylabel("m");
ax2=plt.subplot(2,1,2);
ep = 0.3;
epi = 1/ep;
a = sqrt(2*ep);
V = (a**3)/(8*ep*ep);
m2 = np.multiply( V*np.exp(-np.abs(x)/a), np.cos(np.abs(x)/a)+np.sin(np.abs(x)/a) );
amp = np.max(np.abs(m2));
plt.axis( [xmin, xmax, -1.1*amp, 1.1*amp ] );
plt.plot( gda_cvec(xmin,xmax), gda_cvec(0.0,0.0), "k:", lw=1 );
plt.plot( x, m2, "k-", lw=3 );
ep = 0.6;
epi = 1/ep;
a = sqrt(2*ep);
V = (a**3)/(8*ep*ep);
m2 = np.multiply( V*np.exp(-np.abs(x)/a), np.cos(abs(x)/a)+np.sin(abs(x)/a) );
plt.plot( x, m2, "r-", lw=3 );
plt.xlabel("x");
plt.ylabel("m");
plt.show();
gdapy04_12
In [51]:
# gdapy04_13
#
# constrained least squares fit to synthetic data
print("gdapy04_13");
# z-axus
N=30;
zmin=0;
zmax=10;
z = np.sort(np.random.uniform(low=zmin,high=zmax,size=(N,1)));
# synthetic data
# d = a + b*z + noise
a=2.0;
b=1.0;
sd=0.5;
dobs = a+b*z+np.random.normal(loc=0.0,scale=sd,size=(N,1));
# constraint
zp = 8;
dp = 6;
# data kernel G is NxM, equation i Gm=dobs
M=2;
G = np.zeros((N,2));
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);
# constraint matrix H is PxM , equation is Hm=h
P=1;
H = np.array( [[1.0, zp]] );
h = np.zeros((P,1));
h[0,0] = dp;
# combined matrix is (N+M)x, A [m,lambda] = b
LL = np.concatenate((GTG,H),axis=0);
RR = np.concatenate((H.T,np.zeros((P,P))),axis=0)
A = np.concatenate((LL,RR), axis=1)
b = np.concatenate((GTd, h), axis=0);
v = la.solve(A,b);
mest = np.copy( v[0:M,0:1] );
# predicted data
dpre = np.matmul(G,mest);
e = dobs - dpre;
emax = np.max(e);
iemax = np.argmax(e);
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1=plt.subplot(1,1,1);
# plot bounds
pdmin=0;
pdmax=15;
# plot observed and presicted data, constained point
plt.axis( [zmin, zmax, pdmin, pdmax ] );
plt.plot( z, dobs, "ro", lw=2);
plt.plot( zp, dp, "bo", lw=4);
plt.plot( z, dpre, "g-", lw=3);
plt.xlabel("z");
plt.ylabel("d");
# Weighted least squares fit as a check.
# The idea is to set the weight on the
# prior information to a very large number
# so that it functions as a hard constraint.
wd = 1.0;
wh = 1.0e6;
FF = np.concatenate( (G*wd, H*wh), axis=0 );
ff = np.concatenate( (dobs*wd, h*wh), axis=0 );
mest2 = la.solve( np.matmul(FF.T,FF), np.matmul(FF.T,ff) );
dpre2 = np.matmul(G,mest2);
plt.plot( z, dpre2, "y--", lw=1);
plt.plot( zp, dp, "bo", lw=4);
plt.show();
# Explicit version of solution
mls = la.solve(GTG,GTd);
Z = np.matmul(H,la.solve(GTG,H.T));
Dh = h-np.matmul(H,mls);
q = la.solve(Z,Dh);
mex = mls + la.solve(GTG,np.matmul(H.T,q));
D = np.concatenate((mest,mex,mest-mex),axis=1);
print('comparison with explicit solution');
print('mest mex mest-mes');
print(D);
gdapy04_13
comparison with explicit solution mest mex mest-mes [[ 2.75996088e+00 2.75996088e+00 -4.44089210e-16] [ 4.05004889e-01 4.05004889e-01 -5.55111512e-17]]
In [52]:
# gdapy04_14
# two examples of variance
print("gdapy04_14");
# z-axis
N=101;
zmin=0.0;
zmax=100.0;
Dz=(zmax-zmin)/(N-1);
z=gda_cvec(np.linspace(zmin,zmax,N));
# number of model parameters = number of model data
M=N;
# true model all ones
mtrue = np.ones((M,1));
# two different problems
# 1: sum of first K blocks
# 2: sum of pairs of blocks
X = np.ones((N,1))
Y = np.zeros((1,N));
Y[0,0]=1.0;
G1=la.toeplitz( X, Y );
d1true = np.matmul(G1,mtrue);
X = np.zeros((N,1))
X[0,0]=1.0;
X[1,0]=1.0;
Y = np.zeros((1,N));
Y[0,0]=1.0;
G2=la.toeplitz(X, Y);
d2true = np.matmul(G2,mtrue);
sd=0.1; # standard deviation of data
# synthetic observed data
d1obs = d1true + np.random.normal(loc=0.0,scale=sd,size=(N,1));
d2obs = d2true + np.random.normal(loc=0.0,scale=sd,size=(N,1));
# least squares solutions
G1TG1 = np.matmul(G1.T,G1);
G2TG2 = np.matmul(G2.T,G2);
G1Td1 = np.matmul(G1.T,d1obs);
G2Td2 = np.matmul(G2.T,d2obs);
m1est = la.solve(G1TG1,G1Td1);
m2est = la.solve(G2TG2,G2Td2);
# variances
C1 = (sd**2) * la.inv(G1TG1);
sm1 = gda_cvec( np.sqrt(np.diag(C1)) );
C2 = (sd**2) * la.inv(G2TG2);
sm2 = gda_cvec(np.sqrt(np.diag(C2)) );
# plot
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
# plot bounds
pmmin=-3.0;
pmmax=3.0;
psmin=0.0;
psmax=1.0;
# plot solutions
ax1=plt.subplot(2,1,1);
plt.axis( [zmin, zmax, pmmin, pmmax ] );
plt.plot( z, mtrue, "k-", lw=2 );
plt.ylabel("m");
plt.xlabel("z");
plt.plot( gda_cvec(zmin,zmax), gda_cvec(0.0,0.0), "k-", lw=2 );
plt.plot( z, m2est, "b-", lw=3 );
plt.plot( z, m1est, "k-", lw=3 );
plt.plot( z, m1est, "r-", lw=3 );
# plot sqrt(variance)
ax2=plt.subplot(2,1,2);
plt.axis( [zmin, zmax, psmin, psmax ] );
plt.plot( z, sm1, "r-", lw=3 );
plt.plot( z, sm2, "b-", lw=3 );
plt.ylabel("sm");
plt.xlabel("z");
plt.show();
gdapy04_14
In [53]:
# gdapy04_15
# Examine error surface of intercept and slope in a straight
# line fit to synthetic data
print("gdapy04_15");
# z-axis
N=30;
zmin=-5.0;
zmax=5.0;
z = np.sort(np.random.uniform(low=zmin,high=zmax,size=(N,1)));
# PART 1: data evenly spread along interval
# d = a + b*z + noise
a=2.0;
b=1.0;
sd=0.51;
dobs = a+b*z+np.random.normal(loc=0,scale=sd,size=(N,1));
# intercept grid
Na=101;
amin=0;
amax=4;
Da = (amax-amin)/(Na-1);
a = gda_cvec(np.linspace(amin, amax, Na));
# slope grid
Nb=101;
bmin=0;
bmax=4;
Db = (bmax-bmin)/(Nb-1);
b = gda_cvec(np.linspace(bmin, bmax, Nb));
# populate grid with errors
EA = np.zeros((Na,Nb));
for i in range(Na):
for j in range(Nb):
ao = amin+Da*(i-1);
bo = bmin+Db*(j-1);
dpre = ao + bo*z;
e = dobs-dpre;
myE = np.matmul(e.T,e); myE=myE[0,0];
EA[i,j]=myE;
# minimum error, E1min = E(r1,c1)
k = np.argmin( EA, axis=0 );
Etemp = np.min( EA, axis=0 );
c1 = np.argmin(Etemp);
r1=k[c1];
EAmin = EA[r1,c1];
EAmax = np.max(EA);
# predict data
a1 = amin+Da*(r1-1);
b1 = bmin+Db*(c1-1);
dpre = a1 + b1*z;
# covariance calculation
# least squares formula
G = np.zeros((N,2));
G[0:N,0:1] = np.ones((N,1));
G[0:N,1:2] = z;
C1 = (sd**2)*la.inv(np.matmul(G.T,G));
print("case 1: data straddles origin");
print(" method 1: sm1 %.4f sm2 %.4f cov %.4f" % (sqrt(C1[0,0]), sqrt(C1[1,1]), C1[0,1] ));
# curvature of error suface formula
j=1;
d2Eda2 = (EA[r1+j,c1]-2*EA[r1,c1]+EA[r1-j,c1])/((j*Da)**2);
d2Edb2 = (EA[r1,c1+j]-2*EA[r1,c1]+EA[r1,c1-j])/((j*Db)**2);
d2Edadb = (EA[r1+j,c1+j]-EA[r1+j,c1-j]-EA[r1-j,c1+j]+EA[r1-j,c1-j])/(4.0*j*Da*j*Db);
DA=np.zeros((2,2));
DA[0,0]=d2Eda2;
DA[0,1]=d2Edadb;
DA[1,0]=d2Edadb;
DA[1,1]=d2Edb2;
C2 = (sd**2)*la.inv(DA/2.0);
print(" method 2: sm1 %.4f sm2 %.4f cov %.4f" % (sqrt(C2[0,0]), sqrt(C2[1,1]), C2[0,1] ));
print(" ");
# plot data and fit
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [-5, 5, -10, 10] );
plt.plot( z, dobs, "ro", lw=3);
plt.plot( z, dpre, "b-", lw=2);
plt.xlabel("distance z");
plt.ylabel("data d");
plt.show();
# plot error surface
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [bmin, bmax, amin, amax] );
left=bmin;
right=bmax;
bottom=amax;
top=amin;
plt.imshow( EA, cmap=mycmap, vmin=EAmin, vmax=EAmax, extent=(left,right,bottom,top) );
plt.colorbar(label='error E'); # draw color bar
plt.gca().invert_yaxis();
plt.plot( b1, a1, "wo", lw=2);
X, Y = np.meshgrid(b,a);
DeltaE=(EAmax-EAmin)/100.0;
plt.contour(X,Y,EA,levels=[EAmin+DeltaE],colors=['white'],linewidths=3);
plt.xlabel("slope m2");
plt.ylabel("intercept m1");
plt.show();
# PART 2: data bunched up at end of interval
# z-axis
N=30;
zmin=0.0;
zmax=5.0;
z = np.sort(np.random.uniform(low=zmin,high=zmax,size=(N,1)));
# d = a + b*z + noise
a=2.0;
b=1.0;
sd=0.51;
dobs = a+b*z+np.random.normal(loc=0,scale=sd,size=(N,1));
# intercept grid
Na=101;
amin=0;
amax=4;
Da = (amax-amin)/(Na-1);
a = gda_cvec(np.linspace(amin, amax, Na));
# slope grid
Nb=101;
bmin=0;
bmax=4;
Db = (bmax-bmin)/(Nb-1);
b = gda_cvec(np.linspace(bmin, bmax, Nb));
# populate grid with errors
EA = np.zeros((Na,Nb));
for i in range(Na):
for j in range(Nb):
ao = amin+Da*(i-1);
bo = bmin+Db*(j-1);
dpre = ao + bo*z;
e = dobs-dpre;
myE = np.matmul(e.T,e); myE=myE[0,0];
EA[i,j]=myE;
# minimum error, E1min = E(r1,c1)
k = np.argmin( EA, axis=0 );
Etemp = np.min( EA, axis=0 );
c1 = np.argmin(Etemp);
r1=k[c1];
EAmin = EA[r1,c1];
EAmax = np.max(EA);
# predict data
a1 = amin+Da*(r1-1);
b1 = bmin+Db*(c1-1);
dpre = a1 + b1*z;
# covariance calculation
# least squares formula
G = np.zeros((N,2));
G[0:N,0:1] = np.ones((N,1));
G[0:N,1:2] = z;
C1 = (sd**2)*la.inv(np.matmul(G.T,G));
print("case 1: data straddles origin");
print(" method 1: sm1 %.4f sm2 %.4f cov %.4f" % (sqrt(C1[0,0]), sqrt(C1[1,1]), C1[0,1] ));
# curvature of error suface formula
j=1;
d2Eda2 = (EA[r1+j,c1]-2*EA[r1,c1]+EA[r1-j,c1])/((j*Da)**2);
d2Edb2 = (EA[r1,c1+j]-2*EA[r1,c1]+EA[r1,c1-j])/((j*Db)**2);
d2Edadb = (EA[r1+j,c1+j]-EA[r1+j,c1-j]-EA[r1-j,c1+j]+EA[r1-j,c1-j])/(4.0*j*Da*j*Db);
DA=np.zeros((2,2));
DA[0,0]=d2Eda2;
DA[0,1]=d2Edadb;
DA[1,0]=d2Edadb;
DA[1,1]=d2Edb2;
C2 = (sd**2)*la.inv(DA/2.0);
print(" method 2: sm1 %.4f sm2 %.4f cov %.4f" % (sqrt(C2[0,0]), sqrt(C2[1,1]), C2[0,1] ));
print(" ");
# plot data and fit
fig3 = plt.figure(3,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [-5, 5, -10, 10] );
plt.plot( z, dobs, "ro", lw=3);
plt.plot( z, dpre, "b-", lw=2);
plt.xlabel("distance z");
plt.ylabel("data d");
plt.show();
# plot error surface
fig3 = plt.figure(3,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [bmin, bmax, amin, amax] );
left=bmin;
right=bmax;
bottom=amax;
top=amin;
plt.imshow( EA, cmap=mycmap, vmin=EAmin, vmax=EAmax, extent=(left,right,bottom,top) );
plt.colorbar(label='error E'); # draw color bar
plt.gca().invert_yaxis();
plt.plot( b1, a1, "wo", lw=2);
X, Y = np.meshgrid(b,a);
DeltaE=(EAmax-EAmin)/100.0;
plt.contour(X,Y,EA,levels=[EAmin+DeltaE],colors=['white'],linewidths=3);
plt.xlabel("slope m2");
plt.ylabel("intercept m1");
plt.show();
gdapy04_15
case 1: data straddles origin
method 1: sm1 0.0965 sm2 0.0324 cov -0.0008
method 2: sm1 0.0965 sm2 0.0324 cov -0.0008
case 1: data straddles origin
method 1: sm1 0.1904 sm2 0.0689 cov -0.0114
method 2: sm1 0.1904 sm2 0.0689 cov -0.0114
In [ ]: