In [51]:
# gdapy09_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:
# 2023/02/07 - Created by W. Menke
# 2023/04/27 - More work ...
# 2024/05/23 - Patches to accomodate new verions of numpy, scipy, matplotlib
# patched 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
from scipy import optimize as opt
import matplotlib
# function to make a numpy N-by-1 column vector
# c=gda_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;
def gdaarrow(ax,rbar,mys,myl):
# pretty crazy way to draw an arrowhead!
rbarsq = np.matmul(rbar.T,rbar); rbarsq=rbarsq[0,0];
tangent = rbar/sqrt(rbarsq);
per1 = np.cross( tangent, gda_cvec(0.0, 0.0, 1.0), axis=0 );
per1lensq = np.matmul(per1.T,per1); per1lensq=per1lensq[0,0];
per1 = per1/sqrt(per1lensq);
per2 = np.cross( tangent, per1, axis=0 );
per2lensq = np.matmul(per2.T,per2); per2lensq=per2lensq[0,0];
per2 = per2/sqrt(per2lensq);
L = 0.05;
px=gda_cvec(0.0,rbar[0,0]);
py=gda_cvec(0.0,rbar[1,0]);
pz=gda_cvec(0.0,rbar[2,0]);
ax.plot3D( px.ravel(), py.ravel(), pz.ravel(), mys, lw=myl );
v1 = rbar - L*tangent + 0.25*L*per1;
v2 = rbar - L*tangent - 0.25*L*per1;
v3 = rbar - L*tangent + 0.25*L*per2;
v4 = rbar - L*tangent - 0.25*L*per2;
x0=gda_cvec(rbar[0,0], v1[0,0]);
y0=gda_cvec(rbar[1,0], v1[1,0]);
z0=gda_cvec(rbar[2,0], v1[2,0]);
ax.plot3D( x0.ravel(), y0.ravel(), z0.ravel(), mys, lw=myl );
x0=gda_cvec(rbar[0,0], v2[0,0]);
y0=gda_cvec(rbar[1,0], v2[1,0]);
z0=gda_cvec(rbar[2,0], v2[2,0]);
ax.plot3D( x0.ravel(), y0.ravel(), z0.ravel(), mys, lw=myl );
x0=gda_cvec(rbar[0,0], v3[0,0]);
y0=gda_cvec(rbar[1,0], v3[1,0]);
z0=gda_cvec(rbar[2,0], v3[2,0]);
ax.plot3D( x0.ravel(), y0.ravel(), z0.ravel(), mys, lw=myl );
x0=gda_cvec(rbar[0,0], v4[0,0]);
y0=gda_cvec(rbar[1,0], v4[1,0]);
z0=gda_cvec(rbar[2,0], v4[2,0]);
ax.plot3D( x0.ravel(), y0.ravel(), z0.ravel(), mys, lw=myl );
return 1;
def gdabox(ax,xmin,xmax,ymin,ymax,zmin,zmax):
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmin,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymin,ymax).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmax).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmax,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymax,ymin).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmin).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmin).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymin,ymin).ravel(), gda_cvec(zmax,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymax,ymin).ravel(), gda_cvec(zmax,zmax).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmin,xmin).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmax,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmax).ravel(), gda_cvec(ymax,ymin).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
ax.plot3D( gda_cvec(xmax,xmin).ravel(), gda_cvec(ymax,ymax).ravel(), gda_cvec(zmin,zmin).ravel(), "k-", lw=2 );
return 1;
In [52]:
# gdapy09_01
# depicts a 3D box containing a vector
print("gdapy09_01");
# draw 3D data model space vector
# independent variable x
Nx = 51;
xmin = 0.0;
xmax = 1.0;
Dx = (xmax-xmin)/(Nx-1);
x = gda_cvec( xmin + Dx*np.linspace(0,Nx-1,Nx));
# independent variable y
Ny = 51;
ymin = 0.0;
ymax = 1.0;
Dy = (ymax-ymin)/(Ny-1);
y = gda_cvec(ymin + Dy*np.linspace(0,Ny-1,Ny));
# independent variable z
Nz = 51;
zmin = 0.0;
zmax = 1.0;
Dz = (zmax-zmin)/(Nz-1);
z = gda_cvec(zmin + Dz*np.linspace(0,Nz-1,Nz));
# setup for 3D figure
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax = fig1.add_subplot(111, projection='3d');
plt.axis('off');
plt.grid(b=None);
ax.axes.set_xlim3d(left=xmin, right=xmax);
ax.axes.set_ylim3d(bottom=ymin, top=ymax);
ax.axes.set_zlim3d(bottom=zmin, top=zmax);
# 3D box
gdabox(ax,xmin,xmax,ymin,ymax,zmin,zmax);
# point in 3-space
rbar = gda_cvec(0.7, 0.8, 0.9);
# draw sphere using wireframe technique
Nu = 21;
u0 = 2.0*pi*np.linspace(0,Nu-1,Nu)/(Nu-1);
Nv = 11;
v0 = 2.0*pi*np.linspace(0,Nv-1,Nv)/(Nv-1);
u, v = np.meshgrid(u0,v0);
diameter=0.05;
x0 = diameter*np.cos(u)*np.sin(v)+rbar[0,0];
y0 = diameter*np.sin(u)*np.sin(v)+rbar[1,0];
z0 = diameter*np.cos(v)+rbar[2,0];
ax.plot_wireframe(x0, y0, z0, color="r", facecolor="r");
rbarsq = np.matmul(rbar.T,rbar); rbarsq=rbarsq[0,0];
tangent = rbar/sqrt(rbarsq);
rbar = rbar-diameter*tangent;
gdaarrow(ax,rbar,"k-",3);
plt.show();
print("Caption: Model parameters in the space of possible models, depicted as a");
print("point (red) and vector (black).");
# draw 3D data space vector
# independent variable x
Nx = 51;
xmin = 0.0;
xmax = 1.0;
Dx = (xmax-xmin)/(Nx-1);
x = gda_cvec( xmin + Dx*np.linspace(0,Nx-1,Nx));
# independent variable y
Ny = 51;
ymin = 0.0;
ymax = 1.0;
Dy = (ymax-ymin)/(Ny-1);
y = gda_cvec(ymin + Dy*np.linspace(0,Ny-1,Ny));
# independent variable z
Nz = 51;
zmin = 0.0;
zmax = 1.0;
Dz = (zmax-zmin)/(Nz-1);
z = gda_cvec(zmin + Dz*np.linspace(0,Nz-1,Nz));
# setup for 3D figure
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax = fig1.add_subplot(111, projection='3d');
plt.axis('off');
plt.grid(b=None);
ax.axes.set_xlim3d(left=xmin, right=xmax);
ax.axes.set_ylim3d(bottom=ymin, top=ymax);
ax.axes.set_zlim3d(bottom=zmin, top=zmax);
# 3D box
gdabox(ax,xmin,xmax,ymin,ymax,zmin,zmax);
# point in 3-space
rbar = gda_cvec(0.2, 0.6, 0.6);
# draw sphere using wireframe technique
Nu = 21;
u0 = 2.0*pi*np.linspace(0,Nu-1,Nu)/(Nu-1);
Nv = 11;
v0 = 2.0*pi*np.linspace(0,Nv-1,Nv)/(Nv-1);
u, v = np.meshgrid(u0,v0);
diameter=0.05;
x0 = diameter*np.cos(u)*np.sin(v)+rbar[0,0];
y0 = diameter*np.sin(u)*np.sin(v)+rbar[1,0];
z0 = diameter*np.cos(v)+rbar[2,0];
ax.plot_wireframe(x0, y0, z0, color="b", facecolor="b");
rbarsq = np.matmul(rbar.T,rbar); rbarsq = rbarsq[0,0]
tangent = rbar/sqrt(rbarsq);
rbar = rbar-diameter*tangent;
gdaarrow(ax,rbar,"k-",3);
plt.show();
print("Caption: Datum in the space of possible measurements, depicted as a");
print("point (blue) and vector (black).");
gdapy09_01
Caption: Model parameters in the space of possible models, depicted as a point (red) and vector (black).
Caption: Datum in the space of possible measurements, depicted as a point (blue) and vector (black).
In [53]:
# gdapy09_02
# draw vectors in 3D space
print("gdapy09_02");
# independent variable, x
xmin = 0.0;
xmax = 1.0;
# independent variable, y
ymin = 0.0;
ymax = 1.0;
# independent variable, z
zmin = 0.0;
zmax = 1.0;
# graph 1
# these three vectors span 3D space
m1 = gda_cvec(0.7, 0.8, 0.8);
m2 = gda_cvec(0.5, 0.7, 0.8);
m3 = gda_cvec(0.8, 0.6, 0.7);
# setup for 3D figure
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax = fig1.add_subplot(111, projection='3d');
plt.axis('off');
plt.grid(b=None);
ax.axes.set_xlim3d(left=xmin, right=xmax);
ax.axes.set_ylim3d(bottom=ymin, top=ymax);
ax.axes.set_zlim3d(bottom=zmin, top=zmax);
# 3D box
gdabox(ax,xmin,xmax,ymin,ymax,zmin,zmax);
# arrows
gdaarrow(ax,m1,"k-",3);
gdaarrow(ax,m2,"k-",3);
gdaarrow(ax,m3,"k-",3);
a1 = gda_cvec(m1[0,0], m2[0,0], m3[0,0], m1[0,0]);
a2 = gda_cvec(m1[1,0], m2[1,0], m3[1,0], m1[1,0]);
a3 = gda_cvec(m1[2,0], m2[2,0], m3[2,0], m1[2,0]);
ax.plot3D( a1.ravel(), a2.ravel(), a3.ravel(), "r:", lw=2 );
plt.show();
print("Caption: The three vectors (black) span the 3D space.");
# graph 2
# three vectors lying in a plane
m1 = gda_cvec(0.7, 0.8, 0.9);
m2 = gda_cvec(0.5, 0.7, 0.8);
m3 = gda_cvec(0.5*(m1+m2));
# setup for 3D figure
fig2 = plt.figure(2,figsize=(12,5)); # smallish figure
ax = fig2.add_subplot(111, projection='3d');
plt.axis('off');
plt.grid(b=None);
ax.axes.set_xlim3d(left=xmin, right=xmax);
ax.axes.set_ylim3d(bottom=ymin, top=ymax);
ax.axes.set_zlim3d(bottom=zmin, top=zmax);
# 3D box
gdabox(ax,xmin,xmax,ymin,ymax,zmin,zmax);
# arrows
gdaarrow(ax,m1,"k-",3);
gdaarrow(ax,m2,"k-",3);
gdaarrow(ax,m3,"k-",3);
a1 = gda_cvec(m1[0,0], m2[0,0], m3[0,0], m1[0,0]);
a2 = gda_cvec(m1[1,0], m2[1,0], m3[1,0], m1[1,0]);
a3 = gda_cvec(m1[2,0], m2[2,0], m3[2,0], m1[2,0]);
ax.plot3D( a1.ravel(), a2.ravel(), a3.ravel(), "r:", lw=2 );
plt.show();
print("Caption: The three vectors (black) lie in a plane (red),");
print("and so do not span the 3D space.");
gdapy09_02
Caption: The three vectors (black) span the 3D space.
Caption: The three vectors (black) lie in a plane (red), and so do not span the 3D space.
In [54]:
# gdapy09_03
# Householder rotation (over-determined problem)
# Note: No special effort is made to optimize the computations
print("gdapy09_03");
# random matrix and data
N=10;
M=6;
G = np.random.uniform(low=-1.0,high=1.0,size=(N,M));
mtrue = np.random.uniform(low=-1.0,high=1.0,size=(M,1));
dobs = np.matmul(G,mtrue) + np.random.normal(loc=0.0,scale=0.1,size=(N,1));
gda_draw('title G',G,'title d',dobs);
# least squares olution for reference
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mls = la.solve(GTG,GTd);
els = dobs - np.matmul(G,mls);
Els = np.matmul(els.T,els); Els=Els[0,0];
# setup for recursion
d = np.copy(dobs);
T = np.zeros((N,N,M));
for i in range(M):
a = sqrt(np.sum(np.power(G[i:N,i:i+1],2)));
if( G[i,i]<0.0 ):
a=-a;
v = np.zeros((N,1));
v[i:N,0:1] = G[i:N,i:i+1];
v[i,0] = G[i,i]-a;
T[0:N,0:N,i] = np.identity(N)-2.0*np.matmul(v,v.T)/np.matmul(v.T,v);
Ti=np.squeeze(T[0:N,0:N,i]);
G = np.matmul(Ti,G);
d = np.matmul(Ti,d);
gda_draw("title G%d" % i,G,"title d%d" % i,d);
E = np.sum(np.power(d[M:N,0:1],2));
print("Total Error: least squares %f Householder %f" % (Els, E));
print(" ");
# backsolve
m = np.zeros((M,1));
for i in reversed(range(M)):
m[i,0] = d[i,0];
for j in range(i+1,M):
m[i,0] = m[i,0] - G[i,j]*m[j,0];
m[i,0] = m[i,0]/G[i,i];
print("m-ls m-House error");
D = np.concatenate( (mls,m,np.abs(mls-m)), axis=1 )
print(D);
gdapy09_03
Total Error: least squares 0.008752 Householder 0.008752 m-ls m-House error [[ 5.65428838e-01 5.65428838e-01 2.22044605e-16] [-1.45656258e-01 -1.45656258e-01 1.94289029e-16] [-1.03632770e+00 -1.03632770e+00 4.44089210e-16] [-5.56821448e-01 -5.56821448e-01 2.22044605e-16] [ 3.62848630e-02 3.62848630e-02 5.34294831e-16] [ 6.10152859e-01 6.10152859e-01 2.22044605e-16]]
In [55]:
# gdapy09_04
# transformation that simplifies L=m'*wm*m to L=mp'*mp
print("gdapy09_04");
# simple inverse problem
N=50;
M=N;
m = np.random.normal(loc=0.0,scale=1.0,size=(M,1)); # a model chosen randomly
# prior information of flatness
epsi = 0.1;
tr = gda_cvec(1.0,-1.0,np.zeros((N-2,1)));
tc = gda_cvec(1.0,np.zeros((N-1,1)));
D = epsi*la.toeplitz(tr.ravel(),tc.ravel());
Wm = np.matmul(D.T,D);
# weighted length of m
L = np.matmul(m.T,np.matmul(Wm,m)); L=L[0,0];
# transformtion
Wmsqrt = la.sqrtm(Wm);
mp = np.matmul(Wmsqrt,m);
# unweigthed length of mp
Lp = np.matmul(mp.T,mp); Lp=Lp[0,0];
print("Weighted length of original %.4f" % (L));
print("Unweighted length of transformed %.4f" % (Lp));
gda_draw("title D",D," ", "title Wm^1/2", Wmsqrt);
print("Caption: The first difference operator D (left) is not");
print("equal to sqrt(Wm) (right), as the former is anti-symmetric");
print("and the latter symmetric.");
gdapy09_04 Weighted length of original 1.0147 Unweighted length of transformed 1.0147
Caption: The first difference operator D (left) is not equal to sqrt(Wm) (right), as the former is anti-symmetric and the latter symmetric.
In [56]:
# gdapy09_05
# solve Gm=d using Singular Value Decomposition
# with a G that averages adjacent model pararmaters
print("gdapy09_05");
# model parameters m(z) where z is an auxially variable
M=20;
zmin = 0.0;
zmax = 5.0;
Lz = zmax-zmin;
Dz = (zmax-zmin)/(M-1);
z = gda_cvec(zmin + Dz*np.linspace(0,M-1,M));
# tre model parameter is a sine wave
mtrue = np.sin( 2.0*pi*z/Lz );
# data kernel avarages L adjacent model parameters
# but wraps around from the end to the beginning of
# the model
g1 = gda_cvec(1.0, 1.0, 1.1, 1.0, 1.0);
L,i=np.shape(g1);
g=gda_cvec(g1, np.zeros((M-L,1)));
N=M;
G=np.zeros((N,M));
for i in range(N):
v = gda_cvec(np.roll(g.ravel(),i)).T;
G[i:i+1,0:M] = v;
G=np.fliplr(G);
# observed data is true data plus random noise
sigmad=0.01;
dobs = np.matmul(G,mtrue) + np.random.normal(loc=0.0,scale=sigmad,size=(N,1));
# singular value decomposition of data kernel
U, lam, VT = la.svd(G,full_matrices=False);
V = VT.T;
# selection of non-zero singular valies
p = 16;
# throw out part of decomposition associated with zero singular values
lambdap = gda_cvec(lam[0:p]);
Up=U[0:N,0:p];
Vp=V[0:M,0:p];
# natural solution
# Note: this works, but only whenlambdap is a px1 vector, o watche it!
mest = np.matmul(Vp,np.divide(np.matmul(Up.T,dobs),lambdap));
# damped minimum length solution, for comparison
e2=1e-2;
mestDML =np.matmul(G.T,la.solve(np.matmul(G,G.T)+e2*np.identity(N),dobs));
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
# plot singular values
plt.axis( [1.0, M, 0.0, np.max(lam) ] );
plt.plot( np.linspace(1,M,M), lam, "k-", lw=2 );
plt.plot( np.linspace(1,M,M), lam, 'ko', 'LineWidth', 2 );
plt.xlabel("i");
plt.ylabel("L(i)");
plt.show();
print("Caption: Singular values of the data kernel G.");
# draw dobs = G m for various m's
gda_draw("title G", G," ","title mtrue", mtrue,"="," ","title dobs",dobs," ","title mtrue",mtrue," ","title mest",mest," ","title mDML",mestDML);
print("Caption: (Left) Graphical depiction of the equation Gm=d. (right)");
print("True solution, natural solution and damped minimum length solution.");
gdapy09_05
Caption: Singular values of the data kernel G.
Caption: (Left) Graphical depiction of the equation Gm=d. (right) True solution, natural solution and damped minimum length solution.
In [57]:
# gdapy09_06
# Singular Value Decomposition of two G's
# Both G's averages adjacent model pararmaters
# but the second is a smoother average
print("gdapy09_06");
# model parameters m(z) where z is an auxially variable
M=20;
zmin = 0.0;
zmax = 5.0;
Lz = zmax-zmin;
Dz = (zmax-zmin)/(M-1);
z = gda_cvec(zmin + Dz*np.linspace(0,M-1,M));
# data kernel G1 avarages L adjacent model parameters
# but wraps around from the end to the beginning of
# the model
g1 = gda_cvec(1.0, 1.0, 1.0, 1.0, 1.0);
L,i=np.shape(g1);
g=gda_cvec(g1, np.zeros((M-L,1)));
N=M;
G1=np.zeros((N,M));
for i in range(N):
v = gda_cvec(np.roll(g.ravel(),i)).T;
G1[i:i+1,0:M] = v;
G1=np.fliplr(G1);
# the data krnel G2 is a smoother average, and does not wrap
cmin=0;
cmax=10.0*(1.0/M);
c = gda_cvec(cmin + (cmax-cmin)*np.linspace(0,M-1,M)/(M-1));
l = gda_cvec(np.linspace(1,M,M)); # vector of sequence 1, 2, 3 .. M
C = np.matmul(c,l.T); # tensor product
G2 = np.exp( -C ) - 0.9*np.exp(-np.power(2.0*C,2));
gda_draw(G1,G2);
# Singular Value Decomposition of G1
s1 = la.svd(G1,compute_uv=False);
pp = np.shape(s1);
p = pp[0];
if( (M-p)>0 ): # pad to length M (may not be necessary)
s1 = gda_cvec(S1, np.zeros((M-p),1));
# Singular Value Decomposition of G2
s2 = la.svd(G2,compute_uv=False);
pp = np.shape(s2);
p = pp[0];
if( (M-p)>0 ): # pad to length M (may not be necessary)
s2 = gda_cvec(s2, np.zeros((M-p),1));
# plot singular values, s1
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
plt.axis( [1, M, 0, np.max(s1) ] );
plt.plot( l, s1, "k-", lw=2 );
plt.plot( l, s1, "ko", lw=2 );
plt.xlabel("i");
plt.ylabel("S_i");
plt.show();
# plot singular values, s1
fig2 = plt.figure(2,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
plt.axis( [1, M, 0, np.max(s2) ] );
plt.plot( l, s2, "k-", lw=2 );
plt.plot( l, s2, "ko", lw=2 );
plt.xlabel("i");
plt.ylabel("S_i");
plt.show();
print("Caption: (Top) Singular-values of the first data kernel.");
print("Caption: (Bottom) Singular-values of the second data kernel.");
gdapy09_06
Caption: (Top) Singular-values of the first data kernel. Caption: (Bottom) Singular-values of the second data kernel.
In [58]:
# gdapy09_07
# exercising formula associated the the derivation of the SVD
print("gdapy09_07");
Etotal=0.0; # accumulated error
# exemplary sizes; requires N>M, M>p,
N=19;
M=7;
p=4;
print("N: %d M: %d p: %d" % (N,M,p) );
# random matrix with p non-zero singular values
G = np.concatenate( (np.random.normal(loc=0.0,scale=1.0,size=(p,M)), np.zeros((N-p,M))), axis=0);
# S and its eigenvalues and eigenvectors. Flip the
# order to put the singular values in descending order
S = np.concatenate( (np.concatenate((np.zeros((N,N)),G),axis=1),np.concatenate((G.T,np.zeros((M,M))),axis=1)), axis=0);
lStrue, VS = la.eigh(S);
lStrue = gda_cvec(np.real(lStrue));
VS=np.real(VS);
# check that S=VS*L*VST
e = S-np.matmul(VS,np.matmul(np.diag(lStrue.ravel()),VS.T));
E0 = np.max(np.abs(e));
Etotal = Etotal + E0;
print("Maximum deviation of S from VS*L*VST: %.6e" % (E0) );
# flip to descending order
lStrue = np.flipud(lStrue);
VS = np.fliplr(VS);
# eigenvalues and eigenvectors of GTG
# eigenvalues in ascending order
GTG = np.matmul(G.T,G);
l2V, V = la.eigh(GTG);
l2V = np.flipud(gda_cvec(np.real(l2V)));
V = np.fliplr(np.real(V));
lV = np.sqrt(np.abs(l2V));
Vp = V[0:M,0:p];
V0 = V[0:M,p:M];
Vz = np.zeros((M,N-M)); # fills out so [Vp, Vp, Vz] is MxN
lVp = lV[0:p,0:1];
# eigenvalues and eigenvectors of GGT
# eigenvalues in ascending order
GGT = np.matmul(G,G.T);
l2U, U = la.eigh(GGT);
l2U = np.flipud(gda_cvec(np.real(l2U)));
U=np.fliplr(U);
lU = np.sqrt(np.abs(l2U));
Up = U[0:N,0:p];
U0 = U[0:N,p:M];
Uz = U[0:N,M:N];
lUp = lU[0:p,0:1];
# check that GTG and GGT have the same non-zero eigenvalues
e = lVp-lUp;
El = np.max(np.abs(e));
Etotal = Etotal + El;
print("Maximum deviation of lVp from lUp: %.6e" % (El) );
# build eigenvectors of S from U and V
# MATLAB(R): W = [Up, U0, sqrt(2)*Uz, -U0, -fliplr(Up); Vp, V0, Vz, V0, fliplr(Vp)]/sqrt(2);
Wtop = np.concatenate( (Up, U0, sqrt(2.0)*Uz, -U0, -np.fliplr(Up)),axis=1);
Wbot = np.concatenate( (Vp, V0, Vz, V0, np.fliplr(Vp)), axis=1);
W = np.concatenate( (Wtop,Wbot), axis=0)/sqrt(2);
# check that W is an orthogonal matrix saitifying WWT = WTW = 2I
# the reason for the sqrt(2) is that if the u's and v's are of
# unit lenth, then you need to rescale them when you build the w's.
WWT=np.matmul(W,W.T);
WTW=np.matmul(W.T,W);
e = WWT-np.identity(M+N);
E1 = np.max(np.abs(e));
e = WTW-np.identity(M+N);
E2 = np.max(np.abs(e));
print("Maximum deviation of WWT from I: %.6e" % (E1) );
print("Maximum deviation of WTW from I: %.6e" % (E2) );
Etotal = Etotal + E1;
Etotal = Etotal + E2;
# WT*S*W should be the eigenvalues, so check that they are,
# but resort then to ensure consistent order. There is
# a subtle point here: there is no guarantee that a (u,v)
# pair generated from (GTG and GGT) corresponds to a positive
# singular value - it could be negative, in which case
# (-u,v) generates the positive singular value. So one needs
# to diddle with signs to get the singular values to get the
# list into decending order.
lSest =gda_cvec(np.diag( np.matmul(W.T,np.matmul(S,W))));
lSpre = np.concatenate( (lUp,np.zeros((N+M-2*p,1)),-lUp),axis=0);
for i in range(p):
i2 = N+M-i-1;
if( lSest[i,0] < 0.0 ):
lSest[i,0] = -lSest[i,0];
W[0:N,i:i+1] = -W[0:N,i:i+1];
lSest[i2,0] = -lSest[i2,0];
W[0:N,i2:i2+1] = -W[0:N,i2:i2+1];
# recompute eigenvalues and check in the right order
count = 0;
lSest = gda_cvec(np.diag( np.matmul(W.T,np.matmul(S,W))));
tol=1e-6;
for i in range(1,N+M):
e=lSest[i,0]-lSest[i-1,0];
if(e>tol):
count=count+1;
print("%d %e" % (i,e) );
print("%d eigenvalues are out of order" % (count) );
Etotal = Etotal + count;
# recheck orthonormality
WWT=np.matmul(W,W.T);
WTW=np.matmul(W.T,W);
e = WWT-np.identity(M+N);
E1 = np.max(np.abs(e));
e = WTW-np.identity(M+N);
E2 = np.max(np.abs(e));
print("Maximum deviation of WWT from I: %.6e (after resorting)" % (E1) );
print("Maximum deviation of WTW from I: %.6e (after resorting)" % (E2) );
Etotal = Etotal + E1;
Etotal = Etotal + E2;
# recheck eigenvalue
e = lStrue - lSest;
ElS1 = np.max(np.abs(e));
print("Maximum deviation of WT*S*W from eigenvalues: %.6e (after resorting)" % (ElS1) );
Etotal = Etotal + ElS1;
# [lUP, 0, -LUp] should be the eigenvalues, too,
# so check that they are, but resort then to ensure consistent order
# MATLAB(R): lSpre = sort(lSpre,'descend');
lSpre = np.flipud(np.sort(lSpre,axis=0));
e = lStrue - lSpre;
ElS2 = np.max(np.abs(e));
print("Maximum deviation of (U,V) from eigenvalues: %.6e" % (ElS2) );
Etotal = Etotal + ElS2;
# singular value decomposition
U, l, VT = la.svd(G,full_matrices=False);
l = gda_cvec(l);
V = VT.T;
Vp = V[0:M,0:p];
Up = U[0:N,0:p];
lp = l[0:p,0:1];
# check singular vaues
e = lp - lSpre[0:p,0:1];
ElS3 = np.max(np.abs(e));
print("Maximum deviation of SVD from eigenvalues: %.6e" % (ElS3) );
Etotal = Etotal + ElS3;
# check reconstruction
Gest = np.matmul(U,np.matmul(np.diag(l.ravel()),V.T));
e = G-Gest;
EG = np.max(np.abs(e));
print("Maximum deviation of ULVT from G: %.6e" % (EG) );
# check alternate reconstruction
Gest = np.matmul(Up,np.multiply(Vp.T,lp));
e = G-Gest;
EG = np.max(np.abs(e))
print("Maximum deviation of ULVT from G: %.6e (method 2)" % (EG));
Etotal = Etotal + EG;
# check U eigenvectors
# warning: check will not work right if any non-zero singular values
# are degenerate
Eu = 0.0;
for i in range(p):
uw = sqrt(2)*W[0:N,i:i+1];
usvd = U[0:N,i:i+1];
e = np.abs(1.0-np.abs(np.matmul(uw.T,usvd)));
Eu = Eu + e;
Eu=Eu[0,0];
print("Maximum deviation of U from Usvd: %.6e" % (Eu) );
Etotal = Etotal + Eu;
# check V eigenvectors
# warning: check will not work right if any non-zero singular values
# are degenerate
Ev = 0.0;
for i in range(p):
vw = sqrt(2)*W[N:N+M,i:i+1];
vsvd = V[0:M,i:i+1];
e = np.abs(1.0-np.abs(np.matmul(vw.T,vsvd)));
Ev = Ev + e;
Ev=Ev[0,0];
print("Maximum deviation of V from Vsvd: %.6e" % (Ev) );
Etotal = Etotal + Eu;
# if the total error is close to zero, then all tests
# succeeded and the derivation method really seems to
# reproduce the result of the SVD
print("Overall Error %e" % (Etotal) );
gdapy09_07 N: 19 M: 7 p: 4 Maximum deviation of S from VS*L*VST: 8.538944e-15 Maximum deviation of lVp from lUp: 1.554312e-15 Maximum deviation of WWT from I: 9.037840e-16 Maximum deviation of WTW from I: 4.586419e-16 0 eigenvalues are out of order Maximum deviation of WWT from I: 9.037840e-16 (after resorting) Maximum deviation of WTW from I: 4.586419e-16 (after resorting) Maximum deviation of WT*S*W from eigenvalues: 5.773160e-15 (after resorting) Maximum deviation of (U,V) from eigenvalues: 4.440892e-15 Maximum deviation of SVD from eigenvalues: 1.776357e-15 Maximum deviation of ULVT from G: 2.442491e-15 Maximum deviation of ULVT from G: 2.442491e-15 (method 2) Maximum deviation of U from Usvd: 1.110223e-15 Maximum deviation of V from Vsvd: 7.771561e-16 Overall Error 2.947145e-14
In [59]:
# gdapy09_08
# graphical interpretation of Kuhn-Tucker theorem
print("gdapy09_08");
# (x,y) grid
N=51;
M=51;
Dx=0.01;
Dy=0.01;
x=gda_cvec(Dx*np.linspace(0,N-1,N));
y=gda_cvec(Dy*np.linspace(0,M-1,M));
# a long-wavelngth, 2D sinusoid
E=-np.matmul(np.sin(x),np.sin(y).T);
# minus gradient
dEdx = np.matmul(np.cos(x),np.sin(y).T);
dEdy = np.matmul(np.sin(x),np.cos(y).T);
# gradient list
gs = 0.1; # scale factor
gi = (gda_cvec(40, 30, 40, 25, 40, 10, 10, 20, 10)-1).astype(int);
gj = (gda_cvec(40, 30, 25, 40, 10, 40, 10, 10, 20)-1).astype(int);
Ng,i = np.shape(gi);
# plot
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [ x[0,0], x[N-1,0], y[0,0], y[M-1,0]] );
left=x[0,0]; right=x[N-1,0];
top=y[M-1,0]; bottom=y[0,0];
dmax=np.max(E);
ax1.invert_yaxis();
plt.imshow(np.flipud(E), cmap=mycmap, vmin=-0.25, vmax=0.0, extent=(left,right,bottom,top) );
plt.colorbar();
for k in range(Ng):
i = gi[k,0];
j = gj[k,0];
plt.plot( x[i,0], y[j,0], "wo", lw=3 );
plt.plot( gda_cvec(x[i,0],x[i,0]+gs*dEdx[i,j]), gda_cvec(y[j,0],y[j,0]+gs*dEdy[i,j]), "w-", lw=3 );
# plot inequality constraint line
k1 = 7-1;
k2 = N-7;
plt.plot( gda_cvec(x[k1,0], x[k2,0]), gda_cvec(y[M-1,0], y[0,0]), "k-", lw=3);
# plot arrows on feasible side of equality
p1 = gda_cvec( x[k1,0], y[M-1,0] );
p2 = gda_cvec( x[k2,0], y[0,0] );
t = p1-p2;
tsq = np.matmul(t.T,t); tsq=tsq[0,0];
t = t/sqrt(tsq);
n = gda_cvec(-t[1,0], t[0,0]);
nsq = np.matmul(n.T,n); nsq=nsq[0,0];
n = n/sqrt(nsq);
sn = 0.03;
Nn = 10;
for k in range(2,Nn-1):
a = (k/(Nn-1));
p = a*p1 + (1.0-a)*p2;
plt.plot( gda_cvec(p[0,0],p[0,0]+sn*n[0,0]), gda_cvec(p[1,0],p[1,0]+sn*n[1,0]), "-", lw=2, color=[0.5,0.5,0.5] );
# tabulate the deviation between n and -gradE/|gradE| along
# the constraint line
Nn = 101;
e = np.zeros((Nn,1))
ie = np.zeros((Nn,1),dtype=int);
je = np.zeros((Nn,1),dtype=int);
for k in range(Nn):
a = (k/(Nn-1));
p = a*p1 + (1-a)*p2;
i = int(floor(p[0,0]/Dx));
j = int(floor(p[1,0]/Dy));
g = -gda_cvec(dEdx[i,j], dEdy[i,j]);
gsq = np.matmul(g.T,g); gsq=gsq[0,0];
g = g/sqrt(gsq);
e[k,0] = (n[0,0]-g[0,0])**2 + (n[1,0]-g[1,0])**2;
ie[k,0] = i;
je[k,0] = j;
# find and plot point where gradE and n are anti-parallel
emin = np.min(e);
k = np.argmin(e);
i = ie[k,0];
j = je[k,0];
plt.plot( x[i,0], y[j,0], "ko", lw=4 );
plt.xlabel("m1");
plt.ylabel("m2");
plt.show();
print("Caption: Graphical interpretation of the Kuhn-Tucker theorem, showing");
print("error surface (colors) and its gradient (while lines with dot on uphill end),");
print("constraint surface (black line) and its feasible-pointing normals (grey bars), and fessible point");
print("with minimum error.");
gdapy09_08
Caption: Graphical interpretation of the Kuhn-Tucker theorem, showing error surface (colors) and its gradient (while lines with dot on uphill end), constraint surface (black line) and its feasible-pointing normals (grey bars), and fessible point with minimum error.
In [60]:
# gdapy09_09
# non-negative least squares applied to a simple straight line problem
print("gdapy09_09");
# z-axis
N=11;
zmin=0.0;
zmax=10.0;
Dz = (zmax-zmin)/(N-1);
z = gda_cvec(np.linspace(zmin,zmax,N));
# 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;
# create synthetic data
# d = m1 + m2*z + nois0
mtrue=np.zeros((M,1));
mtrue[0,0]=-2.0;
mtrue[1,0]=1.0;
sd=2.0;
dtrue = np.matmul(G,mtrue);
dobs = dtrue+np.random.normal(loc=0.0,scale=sd,size=(N,1));
# plot data
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, -10, 20 ] );
plt.plot( z, dtrue, "r-", lw=3);
plt.plot( z, dobs, "b.", lw=3);
plt.xlabel("z");
plt.ylabel("d");
# standard least squares solution for oomparison
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mestls = la.solve(GTG,GTd);
dprels = np.matmul(G,mestls);
Els = np.matmul((dobs-dprels).T,(dobs-dprels)); Els=Els[0,0];
plt.plot( z, dprels, "g-", lw=2);
# non-negative least squares solution
vest, E = opt.nnls(G,dobs.ravel());
mest = gda_cvec(vest);
# prediced data
dpre = np.matmul(G,mest);
E = np.matmul((dobs-dpre).T,(dobs-dpre)); E=E[0,0]
plt.plot( z, dpre, "b-", lw=2);
plt.show();
print("Caption: Data (circles) fit by strsight lines, true (black), least squares");
print("(red), non-negative least squares (green).");
print(" ");
print("true m: %.2f %.2f Error %.2f" % (mtrue[0,0], mtrue[1,0], 0.0));
print("ls m: %.2f %.2f Error %.2f" % (mestls[0,0], mestls[1,0], Els));
print("nnls m: %.2f %.2f Error %.2f" % (mest[0,0], mest[1,0], E));
gdapy09_09
Caption: Data (circles) fit by strsight lines, true (black), least squares (red), non-negative least squares (green). true m: -2.00 1.00 Error 0.00 ls m: -1.89 0.89 Error 24.56 nnls m: 0.00 0.62 Error 35.79
In [61]:
# gdapy09_10
# illustration of natural and non-negative least squares
# for gravity problem
print('gdapy09_10\n');
# model is two-dimensional grid of pixels, each of which represents mass
# coorinate system, x is down, y is right, origin in upper-left
Nx = 20;
Ny = 20;
M=Nx*Ny;
# model x coordinate
xmin = 0.0;
xmax = 1.0;
Lx = xmax-xmin;
Dx = (xmax-xmin)/(Nx-1);
x = gda_cvec( xmin + Dx*np.linspace(0,Nx-1,Nx));
# model y coordinate
ymin = 0.0;
ymax = 1.0;
Ly = ymax-ymin;
Dy = (ymax-ymin)/(Ny-1);
y = gda_cvec(ymin + Dy*np.linspace(0,Ny-1,Ny));
# gravity data are on top surface, extending past edges of model by several
# model widths
N=30*Ny;
xo = -0.5; # observations well above surface of model
yomin=ymin-2*Ly;
yomax=ymax+2*Ly;
Lyo = (yomax-yomin);
Dyo = (yomax-yomin)/(N-1);
yo = gda_cvec(yomin + Dyo*np.linspace(0,N-1,N));
# data kernel for vertical component of gravity
g = 1.0; # gravitational constant (not a relistic value)
G=np.zeros((N,M));
jofq=np.zeros((M,1),dtype=int); # backpointer to x(j) value of model parameters
kofq=np.zeros((M,1),dtype=int); # backpointer to y(k) value of model parameters
for i in range(N): # loop over observation points
# loop over masses
q = 0; # counts model parameters
for j in range(Nx): # x
for k in range(Ny): # y
jofq[q,0]=j;
kofq[q,0]=k;
D2 = (x[j,0]-xo)**2 + (y[k,0]-yo[i,0])**2; # distance squared
ct = (x[j,0]-xo) / sqrt(D2); # cosine of vertical angle
G[i,q] = ct*g/D2;
q = q+1;
# true model parameters; just something cooked up to look interesting
wx = gda_cvec(np.linspace(1,Nx,Nx)-Nx/2);
wy = gda_cvec(np.linspace(1,Ny,Ny)-Ny/2);
A = np.matmul(wx,np.ones((1,Ny)))/Nx;
B = np.matmul(np.ones((Nx,1)),wy.T)/Ny;
mxytrue = np.matmul( A , B );
mxytrue = mxytrue - np.min(mxytrue);
mxytrue = mxytrue / np.max(mxytrue) + np.random.uniform(low=0.0,high=0.1,size=(Nx,Ny));
# convert to vector
mtrue=np.zeros((M,1));
for q in range(M):
mtrue[q,0] = mxytrue[ jofq[q,0], kofq[q,0] ];
# convert back to image as a check
#mxytrue2=np.zeros((Nx,Ny));
#for q in range(M):
# mxytrue2[ jofq[q,0], kofq[q,0] ] = mtrue[q,0];
#gda_draw(mxytrue," ",mxytrue2);
# create and plot synthetic data
sigmad=5;
dobs = np.matmul(G,mtrue) + np.random.normal(loc=0.0,scale=sigmad,size=(N,1));
# calculate and plot singular values
U, lam, VT = la.svd(G,full_matrices=False);
lam = gda_cvec(lam);
Nlam, i = np.shape(lam);
V = VT.T;
p=15;
lambdap=lam[0:p,0:1];
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
plt.axis( [1, Nlam, 0.0, np.max(lam) ] );
plt.plot( gda_cvec(np.linspace(1,Nlam,Nlam)), lam, "k-", lw=2 );
plt.plot( gda_cvec(np.linspace(1,Nlam,Nlam)), lam, "ko", lw=2 );
plt.plot( gda_cvec(np.linspace(1,p,p)), lambdap, "r-", lw=2 );
plt.plot( gda_cvec(np.linspace(1,p,p)), lambdap, "ro", lw=2 );
plt.xlabel("i");
plt.ylabel("L(i)");
plt.show();
print("Caption: Singular values (black), with p retained values highlighted (red)");
fig2 = plt.figure(2,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
plt.axis( [yomin, yomax, 0.0, np.max(dobs) ] );
plt.plot( yo, dobs, "k-", lw=3 );
plt.xlabel("y");
plt.ylabel("data d(y)");
# natural solution
p=15;
Up=U[0:N,0:p];
Vp=V[0:M,0:p];
lambdap=lam[0:p,0:1];
mestN = np.matmul(Vp,np.divide(np.matmul(Up.T,dobs),lambdap));
# convert to image
mxyestN=np.zeros((Nx,Ny));
for q in range(M):
mxyestN[ jofq[q,0], kofq[q,0] ] = mestN[q,0];
# natural solution, smaller p
p=4;
Up=U[0:N,0:p];
Vp=V[0:M,0:p];
lambdap=lam[0:p,0:1];
mestN2 = np.matmul(Vp,np.divide(np.matmul(Up.T,dobs),lambdap));
# convert to image
mxyestN2=np.zeros((Nx,Ny));
for q in range(M):
mxyestN2[ jofq[q,0], kofq[q,0] ] = mestN2[q,0];
# predicted data
dpreN = np.matmul(G,mestN);
plt.plot( yo, dpreN, "ro", lw=2 );
eN = dobs-dpreN;
EN=np.matmul(eN.T,eN); EN=EN[0,0];
# non-negative splution
t1, t2 = opt.nnls(G,dobs.ravel());
mestNN = gda_cvec(t1);
# convert to image
mxyestNN=np.zeros((Nx,Ny));
for q in range(M):
mxyestNN[ jofq[q,0], kofq[q,0] ] = mestNN[q,0];
# predicted data
dpreNN = np.matmul(G,mestNN);
plt.plot( yo, dpreNN, "g.", lw=2 );
plt.show();
eNN = dobs-dpreNN;
ENN=np.matmul(eNN.T,eNN); ENN=ENN[0,0];
print("Caption: Observed data (black), data predicted from natural solution (red),");
print("data predicted from non-negative solution (green)");
gda_draw(" ",mxytrue," ",mxyestN2," ",mxyestN," ",mxyestNN);
print("Caption: (Left) True solution, (2nd) Natural solution with small p,");
print("(3rd) Natural solution with large p, (right) Non-negative solution.");
print(" ");
print("Total Error: Natural %f Non-Negative %f" % (EN,ENN) );
gdapy09_10
Caption: Singular values (black), with p retained values highlighted (red)
Caption: Observed data (black), data predicted from natural solution (red), data predicted from non-negative solution (green)
Caption: (Left) True solution, (2nd) Natural solution with small p, (3rd) Natural solution with large p, (right) Non-negative solution. Total Error: Natural 13874.321107 Non-Negative 13948.785672
In [62]:
# gdapy09_11
# depiction of feasible areas, for 4 cases
print("gdapy09_11");
# m1-axis
m1min=0.0;
m1max=2.0;
Nm1 = 41;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec( m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
# m2-axis
m2min=0.0;
m2max=2.0;
Nm2 = 41;
Dm2 =(m2max-m2min)/(Nm2-1);
m2 = gda_cvec( m2min + Dm2*np.linspace(0,Nm2-1,Nm2));
X = np.matmul( m1, np.ones((1,Nm2)) );
Y = np.matmul( np.ones((Nm1,1)), m2.T );
# three inequality constraints
# m1-m2>=0 [1 -1] [0]
# 0.5*m1+m2>=1 [0.5 1] [1]
# m1>=0.2 [1 0] [0.2]
# in the last constraint, the 0.2 can be diddled
# right now this constraint does not affect the
# feasible region, but changing 0.2 to 1.2 changes
# the region considerably
H = np.array( [[1.0, -1.0], [0.5, 1.0], [1.0, 0.0] ] );
NN, MH = np.shape(H);
h = gda_cvec(0.0, 1.0, 0.2);
p=0;
C1 = np.zeros((Nm1,Nm2));
for i in range(Nm1):
for j in range(Nm2):
m = gda_cvec( m1[i,0], m2[j,0] );
myv = np.matmul(H[p:p+1,0:MH],m); myv=myv[0,0]
C1[i,j]=myv-h[p,0];
C1a = (C1>=0)+0.001;
p=1;
C2 = np.zeros((Nm1,Nm2));
for i in range(Nm1):
for j in range(Nm2):
m = gda_cvec( m1[i,0], m2[j,0] );
myv = np.matmul(H[p:p+1,0:MH],m); myv=myv[0,0];
C2[i,j]=myv-h[p,0]
C2a = (C2>=0)+0.001;
p=2;
C3 = np.zeros((Nm1,Nm2));
for i in range(Nm1):
for j in range(Nm2):
m = gda_cvec( m1[i,0], m2[j,0] );
myv = np.matmul(H[p:p+1,0:MH],m); myv=myv[0,0];
C3[i,j]=myv-h[p,0];
C3a = (C3>=0)+0.001;
CT = ((C1>=0)&(C2>=0)&(C3>=0))+0.001;
Gp = np.concatenate((H,h),axis=1).T;
dp = gda_cvec( np.zeros((MH,1)), 1.0);
t1, t2 = opt.nnls(Gp,dp.ravel());
mp = gda_cvec(t1);
ep = dp - np.matmul(Gp,mp);
Ne,i = np.shape(ep);
m = -ep[0:Ne-1,0:1]/ep[Ne-1,0];
im1est = int(floor((m[0,0]-m1min)/Dm1));
jm2est = int(floor((m[1,0]-m2min)/Dm2));
CT[im1est,jm2est]=0.5;
gda_draw("title (A)",C1a," ","title (B)", " ", C2a, "title (C)", C3a, "title (D)",CT);
print("Caption: (Left three frames) Three inequality constraints, esch dividing");
print("the (m2,m2) plane into feasible (brown) and infeasible (blue( region. (Right)");
print("Union of the three inequality constraints, with green dot the feasible point");
print("that minimizes |m|**2");
print(" ");
e=np.matmul(H,m)-h;
for i in range(3):
print("Constraint %d: Hm-h = %f" % (i, e[i,0]) );
gdapy09_11
Caption: (Left three frames) Three inequality constraints, esch dividing the (m2,m2) plane into feasible (brown) and infeasible (blue( region. (Right) Union of the three inequality constraints, with green dot the feasible point that minimizes |m|**2 Constraint 0: Hm-h = 0.000000 Constraint 1: Hm-h = 0.000000 Constraint 2: Hm-h = 0.466667
In [64]:
# gdapy09_12
# least squares with inequality constraints
# exemplary stright line problem
print("gdapy09_12");
# three inequality constraints
# m1-m2>=0 [1 -1] [0]
# 0.5*m1+m2>=1 [0.5 1] [1]
# m1>=0.2 [1 0] [0.2]
# in the last constraint, the 0.2 can be diddled
# right now this constraint does not affect the
# feasible region, but changing 0.2 to 1.2 cnages
# the region considerably
H = np.array( [[1.0, -1.0], [0.5, 1.0], [1.0, 0.0] ] );
NN, MH = np.shape(H);
h = gda_cvec(0.0, 1.0, 0.2);
# straight line problem
# in feasible region: mtrue = [1, 0.8]';
# out of feasbile region mtrue = [1, 1.3]';
mtrue = gda_cvec(1.0, 0.2);
N=11;
z = gda_cvec(np.linspace(0,N-1,N)/(N-1));
G = np.concatenate( (np.ones((N,1)), z), axis=1);
dtrue = np.matmul(G,mtrue);
sd=0.05;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1));
# solve by ordinary least squares, just for check
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mestls = la.solve(GTG,GTd);
dprels = np.matmul(G,mestls);
print("Ordinary least squares:");
print(" mtrue %f %f" % (mtrue[0,0], mtrue[1,0]) );
print(" mest %f %f" % (mestls[0,0], mestls[1,0]) );
print(" ");
# m1-axis
m1min=0.0;
m1max=2.0;
Nm1 = 41;
Dm1 = (m1max-m1min)/(Nm1-1);
m1 = gda_cvec( m1min + Dm1*np.linspace(0,Nm1-1,Nm1));
# m2-axis
m2min=0.0;
m2max=2.0;
Nm2 = 41;
Dm2 =(m2max-m2min)/(Nm2-1);
m2 = gda_cvec( m2min + Dm2*np.linspace(0,Nm2-1,Nm2));
X = np.matmul( m1, np.ones((1,Nm2)) );
Y = np.matmul( np.ones((Nm1,1)), m2.T );
# Error surface, for plotting purposes
E = np.zeros((Nm1,Nm2));
for i in range(Nm1):
for j in range(Nm2):
m = gda_cvec( m1[i,0], m2[j,0] );
e = dobs - np.matmul(G,m);
myE = np.matmul(e.T,e); myE=myE[0,0];
E[i,j]=myE;
Emin = np.min(E);
Emax = np.max(E);
p=0;
C1 = np.zeros((Nm1,Nm2));
for i in range(Nm1):
for j in range(Nm2):
m = gda_cvec( m1[i,0], m2[j,0] );
myv = np.matmul(H[p:p+1,0:MH],m); myv=myv[0,0];
C1[i,j]=myv-h[p,0];
C1a = (C1>=0)+0.001;
p=1;
C2 = np.zeros((Nm1,Nm2));
for i in range(Nm1):
for j in range(Nm2):
m = gda_cvec( m1[i,0], m2[j,0] );
myv = np.matmul(H[p:p+1,0:MH],m); myv=myv[0,0];
C2[i,j]=myv-h[p,0]
C2a = (C2>=0)+0.001;
p=2;
C3 = np.zeros((Nm1,Nm2));
for i in range(Nm1):
for j in range(Nm2):
m = gda_cvec( m1[i,0], m2[j,0] );
myv = np.matmul(H[p:p+1,0:MH],m); myv=myv[0,0];
C3[i,j]=myv-h[p,0]
C3a = (C3>=0)+0.001;
CT = ((C1>=0)&(C2>=0)&(C3>=0))+0.001;
# plot true solution by making dot in CT
im1est = int(floor((mtrue[0,0]-m1min)/Dm1));
jm2est = int(floor((mtrue[1,0]-m2min)/Dm2));
CT[im1est,jm2est]=0.7;
E[im1est,jm2est]=0.7*(Emax-Emin);
# SVD
# note that G must be overdetermined so V=Vp
# economy size (flag 0) SVD so U is Up
Up, lam, VpT = la.svd(G,full_matrices=False);
lam = gda_cvec(lam);
Vp = VpT.T;
rlambda = np.reciprocal(lam);
Lpi = np.diag(rlambda.ravel());
# transformation 1
Hp = -np.matmul(H,np.matmul(Vp,Lpi));
NHp, MHp = np.shape(Hp);
hp = h + np.matmul(Hp,np.matmul(Up.T,dobs));
# transformation 2
# MATLAB(R): Gp = [Hp, hp]';
Gp = np.concatenate( (Hp,hp), axis=1).T
# MATLAB(R): dp = [zeros(1,length(Hp(1,:))), 1]';
dp = gda_cvec( np.zeros((MHp,1)), 1.0 );
t1, t2 = opt.nnls(Gp,dp.ravel());
mpp = gda_cvec(t1);
ep = dp - np.matmul(Gp,mpp);
Ne,i = np.shape(ep);
mp = -ep[0:Ne-1,0:1]/ep[Ne-1,0];
# take mp back to m
# MATLAB(R): mest = Vp*Lpi*(Up'*dobs-mp);
mest = np.matmul( Vp, np.matmul(Lpi, (np.matmul(Up.T,dobs)-mp)));
dpre = np.matmul(G,mest);
print("Constrained least squares:");
print(" mtrue %f %f" % (mtrue[0,0], mtrue[1,0]) );
print(" mest %f %f" % (mest[0,0], mest[1,0]) );
print(" ");
# plot estimated solution by making dot in CT
im1est = int(floor((mest[0,0]-m1min)/Dm1));
jm2est = int(floor((mest[1,0]-m2min)/Dm2));
CT[im1est,jm2est]=0.5;
E[im1est,jm2est]=0.5*(Emax-Emin);
gda_draw("title (A)",CT," ","title (B)",E);
print("Caption: (Left) Feasible region (brown) with error surface (colors), unconstrained");
print("solution (orange dot) and constrained solution (green dot). (Right) Same as");
print("left, but with feasible region not shown.");
fig1 = plt.figure(1,figsize=(12,5)); # smallish figure
ax1 = plt.subplot(1, 1, 1);
plt.axis( [0.0, 1.0, 0.0, 2.0 ] );
plt.plot( z, dobs, "ko", lw=3 );
plt.plot( z, dtrue, "k-", lw=3);
plt.plot( z, dprels, "r-", lw=3);
plt.plot( z, dpre, "g-", lw=3);
plt.xlabel("z");
plt.ylabel("d(z)");
plt.show();
print("Caption: Straight line fit, showing true data (black line), observed data");
print("(circles), unconstrained least squares fit (red), and constrained fit (green),");
print("using least squares with inequality constraint.");
gdapy09_12 Ordinary least squares: mtrue 1.000000 0.200000 mest 1.016964 0.146775 Constrained least squares: mtrue 1.000000 0.200000 mest 0.826255 0.586872
Caption: (Left) Feasible region (brown) with error surface (colors), unconstrained solution (orange dot) and constrained solution (green dot). (Right) Same as left, but with feasible region not shown.
Caption: Straight line fit, showing true data (black line), observed data (circles), unconstrained least squares fit (red), and constrained fit (green), using least squares with inequality constraint.
In [ ]: