In [15]:
# gdapy10_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/15 - Translated from the MATLAB version by W. Menke
# 2024/05/23 - Patches to accomodate new verions of numpy, scipy, matplotlib
# patched dot product to return scalar value
# 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, log, log1p # math functions
import scipy.linalg as la # linear algebra functions
import scipy.optimize as opt # optimization functions such as linear programming
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 [16]:
# gdapy10_01
# comparizon of Normal and Exponential pdfs
print("gdapy10_01");
# d-axis
N=101;
dmin = -5.0;
dmax = 5.0;
Dd = (dmax-dmin)/(N-1);
d = gda_cvec(np.linspace(dmin,dmax,N));
# exponential pdf
dbar=0.0;
sd=1.0;
pE = (1/sqrt(2)) * (1/sd)* np.exp(-sqrt(2)*np.abs(d-dbar)/sd );
AEcheck = Dd*np.sum(pE);
dbarEcheck = Dd*np.sum(np.multiply(d,pE));
sdEcheck = sqrt(Dd*np.sum(np.multiply(np.power(d-dbarEcheck,2),pE)));
# Normal pdf
pN = (1/sqrt(2*pi)) * (1/sd)* np.exp( -0.5*(np.power(d-dbar,2)/(sd**2)));
ANcheck = Dd*np.sum(pN);
dbarNcheck = Dd*np.sum(np.multiply(d,pN));
sdNcheck = sqrt(Dd*np.sum(np.multiply(np.power(d-dbarNcheck,2),pN)));
# area, stddev to make sure it's what I expect
print("exact A %f dbar %f sd %f" % (1, dbar, sd));
print("PE A %f dbar %f sd %f" % (AEcheck, dbarEcheck, sdEcheck));
print("PN A %f dbar %f sd %f" % (ANcheck, dbarNcheck, sdNcheck));
# plot pdf
fig1 = plt.figure(1,figsize=(12,5));
ax1 = plt.subplot(1,1,1);
plt.axis( [dmin, dmax, 0, 1] );
plt.plot(d, pE, "r-", lw=3 );
plt.plot(d, pN, "b-", lw=3 );
plt.xlabel("d");
plt.ylabel("p(d)");
plt.show();
# realizations of two-sided exponential distribution, using one-sided
# MatLab function
Nr = 10000;
mu = sd/sqrt(2);
dra = 2*np.random.randint(low=0,high=2,size=(Nr,1)) - 1.0; # random sign
drb = np.random.exponential(scale=mu,size=(Nr,1)); # exp pdf
dr = dbar + np.multiply(dra,drb);
Lb=100;
h, e = np.histogram(dr,Lb,(dmin,dmax)); # create histogram
Nbin = len(h); # lengths of histogram
Ne = len(e); # length of edges
h = gda_cvec(h); # convert h1 to column-vector
edges = gda_cvec(e[0:Ne-1]); # convert e to column-vector
bins = edges + 0.5*(edges[1,0]-edges[0,0]); # centers of bins
Db = bins[1,0]-bins[0,0]; # size of bins
hmax=np.max(h); # maximum counts
h = h / (Db*np.sum(h));
# histogram
Nb=100;
h, e = np.histogram(dr,Nb,(dmin,dmax)); # create histogram
Nbin = len(h); # lengths of histogram
Ne = len(e); # length of edges
h = gda_cvec(h); # convert h1 to column-vector
edges = gda_cvec(e[0:Ne-1]); # convert e to column-vector
bins = edges + 0.5*(edges[1,0]-edges[0,0]); # centers of bins
Db = bins[1,0]-bins[0,0]; # size of bins
pdest=h/(np.sum(h)*Db);
# max likelihood estimate of mean and stddev
dbarest = np.median(dr);
sdest = (sqrt(2)/Nr)*np.sum(np.abs(dr-dbarest));
# plot histogram
fig1 = plt.figure(1,figsize=(12,5));
ax1 = plt.subplot(1,1,1);
plt.axis( [-5, 5, 0, max(pE)] );
plt.plot(bins,pdest,"k-",lw=3);
plt.plot(d, pE, "r-",lw=2);
plt.xlabel("d");
plt.ylabel("p(d)");
plt.show();
# some further info for verification purposes
print("Maximum liklihood estimates for exponential distribution");
print("mean %f sd %f" % (dbarest, sdest) );
gdapy10_01 exact A 1.000000 dbar 0.000000 sd 1.000000 PE A 1.000875 dbar 0.000000 sd 0.986603 PN A 1.000000 dbar 0.000000 sd 0.999994
Maximum liklihood estimates for exponential distribution mean 0.001554 sd 1.002096
In [17]:
# gdapy10_02
# comparizon of L1 and L2 error in estimating the mean of data
print("gdapy10_02");
fig1 = plt.figure(1,figsize=(12,5));
# this is a one-dimensional grid search
# define limits of m, the mean of the data
M=1001;
mmin = 0;
mmax = 2;
Dm = (mmax-mmin)/(M-1);
m = gda_cvec(np.linspace(mmin,mmax,M));
# randomly generate the data
N=12;
d = np.random.uniform( low=mmin, high=mmax, size=(N,1) );
# L1 norm with even number of data
E1 = np.zeros((M,1));
for i in range(M):
e=d-(m[i,0]*np.ones((N,1)));
E1[i,0]=np.sum(np.abs(e));
# location of each of the data on the mmin-to-mmax interval
# allows to interogate what E[idata] is at the values of
# m where m equals one of the data
idata = (np.floor((d-mmin)/Dm)).astype(int);
# plot E(m) for this case
ax1 = plt.subplot(1,3,1);
mminp=0.0;
mmaxp=2.0;
plt.axis( [mminp, mmaxp, 0, 10] );
plt.plot(m, E1, "r-", lw=3 );
plt.plot( m[idata,0], E1[idata,0], "ko", lw=3 );
plt.xlabel("m");
plt.ylabel("E(m)");
# L1 norm with odd number of data, using the
# same data as previously, except omitting last point
N=N-1;
d = np.copy(d[0:N,0:1]);
E2 = np.zeros((M,1));
for i in range(M):
e=d-(m[i,0]*np.ones((N,1)));
E2[i,0]=np.sum(np.abs(e));
# location of each of the data on the mmin-to-mmax interval
# allows to interogate what E[idata] is at the values of
# m where m equals one of the data
idata = (np.floor((d-mmin)/Dm)).astype(int);
# plot E(m) for this case
plt.subplot(1,3,2);
mminp=0.0;
mmaxp=2.0;
plt.axis( [mminp, mmaxp, 0, 10] );
plt.plot(m, E2, "r-", lw=3 );
plt.plot(m[idata,0], E2[idata,0], "ko", lw=3 );
plt.xlabel("m");
plt.ylabel("E(m)");
# L2 norm
E3 = np.zeros((M,1));
for i in range(M):
e=d-(m[i,0]*np.ones((N,1)));
E3[i,0]=sqrt(np.sum(np.power(e,2)));
# location of each of the data on the mmin-to-mmax interval
# allows to interogate what E[idata] is at the values of
# m where m equals one of the data
idata = (np.floor((d-mmin)/Dm)).astype(int);
# plot E(m) for this case
plt.subplot(1,3,3)
mminp=0;
mmaxp=2;
plt.axis( [mminp, mmaxp, 0, 10] );
plt.plot(m, E3, "r-", lw=3 );
plt.plot(m[idata,0], E3[idata,0], "ko", lw=3 );
plt.xlabel("m");
plt.ylabel("sqrt E(m)");
plt.show();
gdapy10_02
In [18]:
# gdapy10_03
#
# L1 norm inverse problem, overdetermined case
# solved by transformation to a linear programming problem
print("gdapy10_03");
# z-axis
N=20;
zmin = 0.0;
zmax = 1.0;
Dz = (zmax-zmin)/(N-1);
z = gda_cvec(np.linspace(zmin,zmax,N));
# set up for cubic fit
M=4;
mtrue = gda_cvec( 4.0, 3.0, 2.0, 1.0 );
G = np.concatenate( ( np.ones((N,1)), z, np.power(z,2), np.power(z,3)), axis=1);
dtrue = np.matmul(G,mtrue);
sd = 1.0 * np.ones((N,1));
dobs=np.zeros((N,1));
# add exponentially-distributed random noise
for i in range(N):
r=np.random.exponential(scale=sd[i,0]/sqrt(2)); # random on interval [0,inf]
s=2.0*np.random.randint(low=0,high=2)-1.0; # random sign
dobs[i,0] = dtrue[i,0] + r*s; # noise now on interval [-inf,inf]
# outlier
dobs[N-1,0]=5.0;
# least squares solution (sure easier!)
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mls = la.solve(GTG,GTd);
dls = np.matmul(G,mls);
# L1 solution
# set up for linear programming problem
# min f*x subject to A x <= b, Aeq x = beq
# variables
# m = mp - mpp
# x = [mp', mpp', alpha', x', xp']
# with mp, mpp length M and alpha, x, xp, length N
L = 2*M+3*N;
x = np.zeros((L,1));
# f is length L
# minimize sum aplha(i)/sd(i)
f = np.zeros((L,1));
f[2*M:2*M+N,0:1]=np.reciprocal(sd);
# equality constraints
Aeq = np.zeros((2*N,L));
beq = np.zeros((2*N,1));
# first equation G(mp-mpp)+x-alpha=d
Aeq[0:N,0:M] = G;
Aeq[0:N,M:2*M] = -G;
Aeq[0:N,2*M:2*M+N] = -np.identity(N);
Aeq[0:N,2*M+N:2*M+2*N] = np.identity(N);
beq[0:N,0:1] = dobs;
# second equation G(mp-mpp)-xp+alpha=d
Aeq[N:2*N,0:M] = G;
Aeq[N:2*N,M:2*M] = -G;
Aeq[N:2*N,2*M:2*M+N] = np.identity(N);
Aeq[N:2*N,2*M+2*N:2*M+3*N] = -np.identity(N);
beq[N:2*N,0:1] = dobs;
# inequality constraints A x <= b
# part 1: everything positive
A = np.zeros((L+2*M,L));
b = np.zeros((L+2*M,1));
A[0:L,0:L] = -np.identity(L);
b[0:L,0:1] = np.zeros((L,1));
# part 2; mp and mpp have an upper bound. Note:
# Without this constraint, a potential problem is that
# mp and mpp are individually unbounded, though their
# difference, m=mp-mpp, is not. Thus there might be cases
# where the algroithm strays to very large mp and mpp.
A[L:L+2*M,0:L] = np.eye(N=2*M,M=L);
mupperbound=10.0*max(abs(mls)); # might need to be adjusted for problem at hand
b[L:L+2*M,0:1] = mupperbound;
# solve linear programming problem
res = opt.linprog(c=f,A_ub=A,b_ub=b,A_eq=Aeq,b_eq=beq);
x = gda_cvec(res.x);
fmin = -res.fun;
status = res.success;
if( status ):
print("linear programming succeeded" );
else:
print("linear programming failed" );
mest = x[0:M,0:1] - x[M:2*M,0:1];
dpre = np.matmul(G,mest);
# plot true, observed & predcted data
fig1 = plt.figure(1,figsize=(12,5));
ax1 = plt.subplot(1,1,1);
plt. axis( [zmin, zmax, 0, max(dtrue)+5 ] );
plt.plot( z, dtrue, "r-", lw=2);
plt.plot( z, dpre, "g-", lw=3);
plt.plot( z, dls, "b-", lw=2);
plt.plot( z, dobs, "ro", lw=2);
plt.xlabel("z");
plt.ylabel("d");
plt.show();
gdapy10_03 linear programming succeeded
In [19]:
# gdapy10_04
# example of solving an L1 problem by iterating
# a weighted solution to an L2 problem. The example chosen here is
# [d1 d2 ... dN]' = [1 1 ... 1]' m
# since we know that the solution is the mean when the dobs is
# Normally-distributed, that is, mest=mean(dobs) and that the
# solution is the median when dobs is exponentially-distributed,
# that is mest=median(dobs). The quantity f = (dmean-mest)/(dmean-dmedian)
# is used as a quality of solution indicator. It measures how
# much closer the solution is to the median than the mean.
# The scripts creates the empirical probability density function p(f).
# Note: The script is a bit sluggish because the L1 problem
# is being solved many time to generate a histogram
print("gdapy10_04");
# Allows you to toggle between different ways of initializing
# the solution.
USEL2 = 1;
USETRUE = 0;
USEZERO = 0;
USEONE = 0;
# Allows you to toggle between simpe formula for weigths and
# Frommlet & Nuel's (2016) more complicated (but equivalent and
# more numerically stable version)
USEFN = 1;
# true solution
mtrue = gda_cvec(0.0);
Nreal = 100; # number of realizations
f = np.zeros((Nreal,1)); # (mest-mmedian)/(mmean-mmedian)
for j in range(Nreal): # loop over realizations
# exponentially-distributed noise
N = 101; # must be odd so that median is unique
# generate exponentially-distributed random noise
smp = 0.1;
mu = smp/sqrt(2);
rsign = 2*np.random.randint(low=0,high=2,size=(N,1))-1; # +1 or -1
rval = np.random.exponential(scale=mu,size=(N,1));
dn = np.multiply(rsign,rval);
# set up problem
M=1;
G = np.ones((N,1)); # data kernel says d(i) = 1 m(1)
dtrue =np.matmul(G,mtrue); # true data
dobs = dtrue + dn; # observed data is true data plus noise
dmean = np.mean(dobs); # sample mean is poor estimate of m(1)
dmedian = np.median(dobs); # median is good estimate of m(1)
H = abs( dmean - dmedian ); # for plotting purposes
Ni=50; # number of iterations of solution
if( j==1 ): # plot solution as a function of iteration,
# but only for first realization
fig1 = plt.figure(1,figsize=(12,5));
ax1 = plt.subplot(1,1,1);
plt.axis( [0, Ni+1, min(dmean,dmedian)-H/4, max(dmean,dmedian)+H/4] );
plt.plot( gda_cvec(0.0,Ni), gda_cvec(dmean,dmean), "k:", lw=2 );
plt.plot( gda_cvec(0.0,Ni), gda_cvec(dmedian,dmedian), "k-", "LineWidth", 2 );
plt.xlabel("iteration k");
plt.ylabel("m(k)");
plt.title("median (solid), mean (dashed)");
# setup for iterative solution
delta = 1e-5;
gamma = 2.0;
n = 1;
# initial solution
if( USEL2 == 1 ):
W = np.identity(N);
GTG = np.matmul( G.T, np.matmul(W,G) );
GTd = np.matmul( G.T, np.matmul(W,dobs) );
mest = la.solve(GTG,GTd); # L2 solution
elif( USETRUE == 1 ):
mest = mtrue; # true solution
elif( USEZERO == 1 ):
mest = 0; # zero solution
elif( USEONE == 1 ):
mest = np.ones((M,1)); # unity solution
for i in range(Ni): # iterative solution
e = dobs-np.matmul(G,mest); # prediction error
ae = np.abs(e);
if( USEFN ): # Frommlet and Nuels (2016) weight formula
w = np.zeros((N,1));
for k in range(N):
if( ae[k,0]>=delta ):
w[k,0]=(delta**(n-2))*exp(((n-2)/gamma)*log1p((ae[k,0]/delta)**gamma));
else:
w[k,0]=(ae[k,0]**(n-2))*exp(((n-2)/gamma)*log1p((delta/ae[k,0])**gamma));
else: # simple weight formula
# MATLABw = ((ae.^gamma) + (delta^gamma) ).^((n-2)/gamma); % weight matrix
wx=np.power(ae,gamma)+np.ones((N,1))*(delta**gamma)
w=np.power(wx,(n-2)/gamma); # weight matrix
W = np.diag(w.ravel());
GTG = np.matmul(G.T, np.matmul(W,G) );
GTd = np.matmul(G.T, np.matmul(W,dobs) );
mest = la.solve(GTG,GTd); # L2 solution
mestscalar=mest[0,0];
if( j==1 ): # plot, but only first realization
plt.plot( i, mest[0,0], 'ko', 'LineWidth', 2 );
# f is a quality of solution indicator
# the solution should be the median, not the mean
# f is normalized closeness of solution from medium
f[j,0] = (dmean-mestscalar)/(dmean-dmedian);
# end loop over realizations
# histogram
binleft = (-2.0);
binright = 2.0;
Nb=100;
h, e = np.histogram(f,Nb,(binleft,binright)); # create histogram
Nbin = len(h); # lengths of histogram
Ne = len(e); # length of edges
h = gda_cvec(h); # convert h1 to column-vector
edges = gda_cvec(e[0:Ne-1]); # convert e to column-vector
bins = edges + 0.5*(edges[1,0]-edges[0,0]); # centers of bins
Db = bins[1,0]-bins[0,0]; # size of bins
p=h/(np.sum(h)*Db); # probability density for f
fig2 = plt.figure(2,figsize=(12,5));
ax1 = plt.subplot(1,1,1);
pmax=np.max(p);
plt. axis( [binleft, binright, 0.0, pmax ] );
plt.plot( gda_cvec(0.0,0.0), gda_cvec(0.0,pmax/10), "k:", lw=2 );
plt.plot( gda_cvec(1.0,1.0), gda_cvec(0.0,pmax/10), "k:", lw=2 );
plt.plot( bins, p, "k-", lw=3 );
plt.xlabel("f");
plt.ylabel("p(f)");
plt.title("probability density function p(f)");
plt.show();
gdapy10_04
In [20]:
# gdama10_05
# Linf norm inverse problem, overdetermined case
# solved by transformation to a linear programming problem
print("gdapy10_05");
# z-axis
N=20;
zmin = 0.0;
zmax = 1.0;
Dz = (zmax-zmin)/(N-1);
z = gda_cvec(np.linspace(zmin,zmax,N));
# set up for cubic fit
M=4;
mtrue = gda_cvec( 4.0, 3.0, 2.0, 1.0 );
G = np.concatenate( (np.ones((N,1)), z, np.power(z,2), np.power(z,3)), axis=1 );
dtrue = np.matmul(G,mtrue);
sd = 1.0 * np.ones((N,1));
dobs=np.zeros((N,1));
for i in range(N):
rv = np.random.exponential(scale=sd[i,0]/sqrt(2))
rs = 2*np.random.randint(low=0,high=2)-1;
dobs[i,0] = dtrue[i,0] + rs*rv;
# outlier
dobs[ int(0.37*N), 0 ] = 2.5;
# least squares solution (sure eassier!)
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,dobs);
mls = la.solve(GTG,GTd);
dls = np.matmul(G,mls);
# Linf solution
# set up for linear programming problem
# min f*x subject to A x <= b, Aeq x = beq
# variables
# m = mp - mpp
# x = [mp', mpp', alpha', x', xp']
# with mp, mpp length M; alpha length 1, x, xp, length N
L = 2*M+1+2*N;
x = np.zeros((L,1));
# f is length L
# minimize alpha
f = np.zeros((L,1));
f[2*M,0]=1.0;
# equality constraints
Aeq = np.zeros((2*N,L));
beq = np.zeros((2*N,1));
# first equation G(mp-mpp)+x-alpha=d
# MATLAB Aeq(1:N,1:M)=G;
# MATLAB Aeq(1:N,M+1:2*M) = -G;
# MATLAB Aeq(1:N,2*M+1:2*M+1) = -1./sd;
# MATLAB Aeq(1:N,2*M+1+1:2*M+1+N) = eye(N,N);
# MATLAB beq(1:N) = dobs;
Aeq[0:N,0:M] = G;
Aeq[0:N,M:2*M] = -G;
Aeq[0:N,2*M:2*M+1] = -1./sd;
Aeq[0:N,2*M+1:2*M+1+N] = np.identity(N);
beq[0:N,0:1] = dobs;
# second equation G(mp-mpp)-xp+alpha=d
# MATLAB Aeq(N+1:2*N,1:M) = G;
# MATLAB Aeq(N+1:2*N,M+1:2*M) = -G;
# MATLAB Aeq(N+1:2*N,2*M+1:2*M+1) = 1./sd;
# MATLAB Aeq(N+1:2*N,2*M+1+N+1:2*M+1+2*N) = -eye(N,N);
# MATLAB beq(N+1:2*N) = dobs;
Aeq[N:2*N,0:M] = G;
Aeq[N:2*N,M:2*M] = -G;
Aeq[N:2*N,2*M:2*M+1] = 1./sd;
Aeq[N:2*N,2*M+1+N:2*M+1+2*N] = -np.identity(N);
beq[N:2*N,0:1] = dobs;
# inequality constraints A x <= b
# part 1: everything positive
A = np.zeros((L+2*M,L));
b = np.zeros((L+2*M,1));
# MATLAB A(1:L,:) = -eye(L,L);
# MATLAB b(1:L) = zeros(L,1);
A[0:L,0:L] = -np.identity(L);
b[0:L,0:1] = np.zeros((L,1));
# part 2; mp and mpp have an upper bound. Note:
# Without this constraint, a potential problem is that
# mp and mpp are individually unbounded, though their
# difference, m=mp-mpp, is not. Thus there might be cases
# where the algroithm strays to very large mp and mpp.
# MATLAB A(L+1:L+2*M,:) = eye(2*M,L);
# MATLAB b(L+1:L+2*M) = mupperbound;
A[L:L+2*M,0:L] = np.eye(N=2*M,M=L);
mupperbound=10.0*np.max(np.abs(mls)); # might need to be adjusted for problem at hand
b[L:L+2*M,0:L] = mupperbound*np.ones((2*M,1));
# solve linear programming problem
res = opt.linprog(c=f,A_ub=A,b_ub=b,A_eq=Aeq,b_eq=beq);
x = gda_cvec(res.x);
fmin = -res.fun;
status = res.success;
if( status ):
print("linear programming succeeded" );
else:
print("linear programming failed" );
mest = x[0:M,0:1] - x[M:2*M,0:1];
dpre = np.matmul(G,mest);
# plot
fig1 = plt.figure(1,figsize=(12,5));
ax1 = plt.subplot(1,1,1);
plt.axis( [zmin, zmax, 0, np.max(dtrue) ] );
plt.plot( z, dtrue, "r-", lw=2);
plt.plot( z, dpre, "g-", lw=3);
plt.plot( z, dls, "b-", lw=2);
plt.plot( z, dobs, "ro", lw=2);
plt.xlabel("z");
plt.ylabel("d");
plt.show();
print("soluton");
print("mtrue %.2f %.2f %.2f %.2f" % (mtrue[0,0],mtrue[1,0],mtrue[2,0],mtrue[3,0]));
print("mls %.2f %.2f %.2f %.2f" % (mls[0,0],mls[1,0],mls[2,0],mls[3,0]));
print("mLinf %.2f %.2f %.2f %.2f" % (mest[0,0],mest[1,0],mest[2,0],mest[3,0]));
gdapy10_05 linear programming succeeded
soluton mtrue 4.00 3.00 2.00 1.00 mls 4.11 0.53 11.51 -6.38 mLinf 5.99 -19.83 52.94 -28.07
In [21]:
# gdapy10_06
# An example of solving a minimum length problem under
# the L0 and L1 norms, by the itertive reweighting algorithm
# The example chosen here is
# d(t) = g(t) * m(t) where * is convolution
# The convolution is represented by a Toeplitz matrix G
print("gdama10_06");
# t-axis
N=100;
M=N;
t = gda_cvec(np.linspace(0,N-1,N));
# a gaussian pulse g(t)
t0 = 10.0;
s = 2.5;
s2 = s**2;
g = np.exp(-np.power(t-t0,2)/(2*s2)); # gaussian g(t)
# matrix equivalent to convolution by g
G = la.toeplitz(g,gda_cvec(g[0,0],np.zeros((N-1,1))).T);
# true solution is a collection of spikes
mtrue = np.zeros((M,1));
mtrue[4,0]=1.0; mtrue[19,0]=0.5; mtrue[39,0]=0.25;
dtrue = np.matmul(G,mtrue); # true data
gda_draw("title G", G, "title m", mtrue, "=", "title d", dtrue);
# observed data is true data plus random noise
sd = 0.05;
dobs = dtrue+np.random.normal(loc=0,scale=sd,size=(N,1));
# L2 solution
mu = 0.01;
w = np.ones((M,1));
GTG = np.matmul(G.T,G)+mu*np.diag(w.ravel());
GTd = np.matmul(G.T,dobs);
mL2 = la.solve(GTG,GTd);
dL2 = np.matmul(G,mL2);
e = dobs-dL2;
EL2 = np.matmul(e.T,e); EL2=EL2[0,0];
# setup for L1 norm
n = 1;
mu = 0.01;
delta = 1e-5;
gamma = 2;
w = np.ones((M,1));
for j in range(10): # reweighting algorithm
GTG = np.matmul(G.T,G)+mu*np.diag(w.ravel());
GTd = np.matmul(G.T,dobs);
mest = la.solve(GTG,GTd);
dpre = np.matmul(G,mest);
e = dobs-dpre;
E = np.matmul(e.T,e); E=E[0,0];
am = np.abs(mest); # Frommlet & Nuel's (2016) weight formula
for k in range(M):
if( am[k,0]>=delta ):
w[k,0]=(delta**(n-2))*exp(((n-2)/gamma)*log1p((am[k,0]/delta)**gamma));
else:
w[k,0]=(am[k,0]**(n-2))*exp(((n-2)/gamma)*log1p((delta/am[k,0])**gamma));
# save solution
mL1 = np.copy(mest);
dL1 = np.copy(dpre);
EL1 = E;
# setup for L0 norm (actually 0.1 norm)
n=0.1;
mu = 0.01;
epsi = 1.0e-6;
w = np.ones((M,1));
for j in range(10): # reweighting algorithm
GTG = np.matmul(G.T,G)+mu*np.diag(w.ravel());
GTd = np.matmul(G.T,dobs);
mest = la.solve(GTG,GTd);
dpre = np.matmul(G,mest);
e = dobs-dpre;
E = np.matmul(e.T,e); E=E[0,0];
am = np.abs(mest); # Frommlet & Nuel's (2016) weight formula
for k in range(M):
if( am[k,0]>=delta ):
w[k,0]=(delta**(n-2))*exp(((n-2)/gamma)*log1p((am[k,0]/delta)**gamma));
else:
w[k,0]=(am[k,0]**(n-2))*exp(((n-2)/gamma)*log1p((delta/am[k,0])**gamma));
# save solution
mL0 = np.copy(mest);
dL0 = np.copy(dpre);
EL0 = E;
print("EL2/N %f EL1/N %f EL0/N %f sd %f" % (EL2, EL1, EL0, sd) );
fig1 = plt.figure(1,figsize=(12,5));
ax1 = plt.subplot(3,1,1); # plot g(t)
plt.axis( [0.0, N-1, -0.5, 1.5] );
plt.xlabel("t (s)");
plt.ylabel("g(t)");
plt.title("gtrue (black)");
plt.plot( t, g, "k-", lw=4);
ax1 = plt.subplot(3,1,3); # plot d(t)
plt.axis( [0.0, N-1, -0.5, 1.5] );
plt.xlabel("t (s)");
plt.ylabel("d(t)");
plt.title("dobs (black), dL0 (red), dL1 (blue) dL2 (green)");
plt.plot( t, dobs, "ko", lw=1);
plt.plot( t, dL2, "g-", lw=5);
plt.plot( t, dL1, "b-", lw=3);
plt.plot( t, dL0, "r-", lw=1);
ax1 = plt.subplot(3,1,2); # plot m(t)
plt.axis( [0.0, N-1, -0.5, 1.0] );
plt.xlabel("t (s)");
plt.ylabel("m(t)");
plt.title("mtrue (black), L2 (green) L1 (blue), L0 (red)");
plt.plot( t, mtrue, "k-", lw=4);
plt.plot( t, mL2, "g-", lw=4);
plt.plot( t, mL1, "b-", lw=4);
plt.plot( t, mL0, "r-", lw=2);
plt.show();
fig2 = plt.figure(2,figsize=(12,5));
ax1 = plt.subplot(1,1,1);
plt.axis( [0.0, 10, -0.5, 1.0] );
plt.xlabel("t (s)");
plt.ylabel("m(t)");
plt.title("mtrue (black), L2 (green) L1 (blue), L0 (red)");
plt.plot( t, mtrue, "k-", lw=7);
plt.plot( t, mL2, "g-", lw=5);
plt.plot( t, mL1, "b-", lw=3);
plt.plot( t, mL0, "w-", lw=2);
plt.plot( t, mL0, "r-", lw=1);
plt.show();
gdama10_06
EL2/N 0.189479 EL1/N 0.213675 EL0/N 0.275237 sd 0.050000
In [22]:
# gdapy10_07
# An example of solving a minimum gradient problem under
# the L0 and L1 norms using the iterative reweighting algorithm
# (The L1 solution is the TV solution).
# The example chosen here is
# d(t) = g(t) * m(t) where * is convolution
# The convolution is represented by a Toeplitz matrix G
print("gdama10_07");
# t-axis
N=100;
M=N;
t = gda_cvec(np.linspace(0.0,N-1,N));
# gaussian pulse g(t)
t0 = 10.0;
s = 2.5;
s2 = s**2;
g = np.exp(-np.power(t-t0,2)/(2*s2) ); # gaussian g(t)
# matrix equivalent to convolution by g
G = la.toeplitz(g,gda_cvec(g[0,0],np.zeros((N-1,1))).T);
# true solution is a boxcar function
mtrue = np.zeros((N,1));
mtrue[29:60,0:1]=np.ones((31,1));
# true data
dtrue = np.matmul(G,mtrue);
gda_draw("title G", G, "title m", mtrue, "=", "title d", dtrue);
# observed data is true data plus random noise
sd = 0.05;
dobs = dtrue+np.random.normal(loc=0.0,scale=sd,size=(N,1));
# first difference operator
D = la.toeplitz( gda_cvec(-1.0, np.zeros((M-2,1))), gda_cvec(-1.0,1.0,np.zeros((M-2,1))).T );
# L2 solution
mu = 0.01;
w = np.ones((M-1,1));
GTG = np.matmul(G.T,G) + mu * np.matmul( D.T, np.matmul(np.diag(w.ravel()), D));
GTd = np.matmul(G.T,dobs);
mL2= la.solve(GTG,GTd);
dL2 = np.matmul(G,mL2);
e = dobs-dpre;
EL2 = np.matmul(e.T,e); EL2=EL2[0,0];
# setup for L1 norm
n = 1;
delta = 1e-5;
gamma = 2;
mu = 0.01;
w = np.ones((M-1,1));
for j in range(10): # reweighting algorithm
GTG = np.matmul(G.T,G) + mu * np.matmul( D.T, np.matmul(np.diag(w.ravel()), D));
GTd = np.matmul(G.T,dobs);
mest = la.solve(GTG,GTd);
dpre = np.matmul(G,mest);
e = dobs-dpre;
E = np.matmul(e.T,e); E=E[0,0];
v = np.matmul(D,mest);
av = np.abs(v); # Frommlet & Nuel's (2016) weight formula
for k in range(M-1):
if( av[k,0]>=delta ):
w[k,0]=(delta**(n-2))*exp(((n-2)/gamma)*log1p((av[k,0]/delta)**gamma));
else:
w[k,0]=(av[k,0]**(n-2))*exp(((n-2)/gamma)*log1p((delta/av[k,0])**gamma));
mL1 = np.copy(mest);
dL1 = np.copy(dpre);
EL1 = E;
# setup for L0 norm (actually 0.1 norm)
n=0.1;
mu = 0.01;
epsi = 1.0e-6;
w = np.ones((M-1,1));
for j in range(10): # reweighting algorithm
GTG = np.matmul(G.T,G) + mu * np.matmul( D.T, np.matmul(np.diag(w.ravel()), D));
GTd = np.matmul(G.T,dobs);
mest = la.solve(GTG,GTd);
dpre = np.matmul(G,mest);
e = dobs-dpre;
E = np.matmul(e.T,e); E=E[0,0];
v = np.matmul(D,mest);
av = np.abs(v); # Frommlet & Nuel's (2016) weight formula
for k in range(M-1):
if( av[k,0]>=delta ):
w[k,0]=(delta**(n-2))*exp(((n-2)/gamma)*log1p((av[k,0]/delta)**gamma));
else:
w[k,0]=(av[k,0]**(n-2))*exp(((n-2)/gamma)*log1p((delta/av[k,0])**gamma));
mL0 = np.copy(mest);
dL0 = np.copy(dpre);
EL0 = E;
# print out sumamry statistics
print("EL2/N %f EL1/N %f EL0/N %f sd %f" % (EL2, EL1, EL0, sd) );
fig1 = plt.figure(1,figsize=(12,5));
ax1 = plt.subplot(3,1,1); # plot g(t);
plt.axis( [0.0, N-1, -0.5, 1.5] );
plt.xlabel("t (s)");
plt.ylabel("g(t)");
plt.title("gtrue (black)");
plt.plot( t, g, "k-", lw=4);
ax1 = plt.subplot(3,1,3); # plot d(t);
plt.axis( [0, N-1, -5, 15] );
plt.xlabel("t (s)");
plt.ylabel("d(t)");
plt.title("dtrue (black), dL2 (green), dL1 (blue), DL0 (red)");
plt.plot( t, dobs, "ko", lw=1);
plt.plot( t, dL2, "g-", lw=7);
plt.plot( t, dL1, "b-", lw=5);
plt.plot( t, dL0, "r-", lw=3);
plt.subplot(3,1,2); # plot m(t)
plt.axis( [0.0, N-1, -0.5, 1.5] );
plt.xlabel("t (s)");
plt.ylabel("m(t)");
plt.title("mtrue (black), L2 (green) L1 (blue), L0 (red)");
plt.plot( t, mtrue, "k-", lw=4);
plt.plot( t, mL2, "g-", lw=4);
plt.plot( t, mL1, "b-", lw=4);
plt.plot( t, mL0, "r-", lw=2);
plt.show();
# close up of edge of boxcar
fig2 = plt.figure(2,figsize=(12,5));
plt.axis( [25.0, 35.0, -0.5, 1.5] );
plt.xlabel("t (s)");
plt.ylabel("m(t)");
plt.title("mtrue (black), L2 (green) L1 (blue), L0 (red)");
plt.plot( t, mtrue, "k-", lw=7);
plt.plot( t, mL2, "g-", lw=5);
plt.plot( t, mL1, "b-", lw=3);
plt.plot( t, mL0, "y-", lw=2);
plt.plot( t, mL0, "r-", lw=1);
plt.show();
gdama10_07
EL2/N 1093.051567 EL1/N 0.148608 EL0/N 0.188344 sd 0.050000
In [ ]: