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
No description has been provided for this image
No description has been provided for this image
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
No description has been provided for this image
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
No description has been provided for this image
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
No description has been provided for this image
No description has been provided for this image
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
No description has been provided for this image
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
No description has been provided for this image
EL2/N 0.189479 EL1/N 0.213675 EL0/N 0.275237 sd 0.050000
No description has been provided for this image
No description has been provided for this image
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
No description has been provided for this image
EL2/N 1093.051567 EL1/N 0.148608 EL0/N 0.188344 sd 0.050000
No description has been provided for this image
No description has been provided for this image
In [ ]: