In [1]:
# gdapy11_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/05/02 - Created by W. Menke, from MATLAB Version
# 2024/05/23 -  Patches to accomodate new verions of numpy, scipy, matplotlib
#               patched dot product to return scalar value
#               patched issue with gray-coding arrays

# 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, atan2   # math functions
from scipy.special import erf
import scipy.linalg as la             # linear algebra functions
import scipy.signal as sg             # signal processing
import scipy.io as scipyio            # input/output
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] = v[0:N,0];
        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 gda_FTFmul(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;


# function to perform the multiplication FT f
def gda_FTFrhs(f):
    # 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;
    return gdaFsparse.transpose()*f;

def gda_matrixmin(E):
    E1 = np.min(E,axis=0);
    k = np.argmin(E,axis=0);
    Emin = np.min(E1);
    ic = np.argmin(E1);
    ir = k[ic];
    return ir, ic, Emin;

def gda_Esurface( xc, yr, E, sigmad2 ):
# gda_Esurface: analysis of 2D error surface
# input:
#     E: matrix of errors xc increases with columns, yr with rows
#     xc, yr: corresponding vectors of (xc, yr) values; must be evenly-spaced
#     sigmad2: variance of data that went into E=e'*e with e=dobs-dpre
# output:
#    x0, y0: location of minimum in E
#     E0: value of E at minimum
#     covxy: covariance matrix
#     status: 1 on success, 0 on fail
# method: quadratic approximation

# default return values
    x0=0.0;
    y0=0.0;
    E0=0.0;
    covxy=np.zeros((2,2));
    status=0;
    Nx,i = np.shape(xc);
    Ny,i = np.shape(yr);
    Dx = xc[1,0]-xc[0,0];
    Dy = yr[1,0]-yr[0,0];

    ir,ic,Emin = gda_matrixmin(E)

    if( (ir<1) or (ir>(Ny-2)) or (ic<1) or (ic>(Nx-2)) ):
        # minimum not in interior
        return x0, y0, E0, covxy, status;

    # quadratic model with 9 data
    # E = m1 + m2*x + m3*y + m4*x*x/2 + m5*x*y + m6*y*y/2
    # dE/dx = m2 + m4*x + m5*y
    # dE/dy = m3 + m5*x + m6*y
    # at minimum
    # 0 = d1 + D2 [x,y]'  so [x,y] = -inv(D2)*d1;

    x = Dx*gda_cvec(-1, 0, 1, -1, 0, 1, -1, 0, 1);
    y = Dy*gda_cvec(-1, -1, -1, 0, 0, 0, 1, 1, 1);
    d = gda_cvec( E[ir-1,ic-1], E[ir-1,ic], E[ir-1,ic+1],
                  E[ir,ic-1],   E[ir,ic],   E[ir,ic+1],
                  E[ir+1,ic-1], E[ir+1,ic], E[ir+1,ic+1]);
    G = np.concatenate((np.ones((9,1)),x,y,
                        0.5*np.multiply(x,x),np.multiply(x,y),0.5*np.multiply(y,y)),axis=1);
    GTG = np.matmul(G.T,G);
    GTd = np.matmul(G.T,d);
    m = la.solve(GTG,GTd);
    d1 = gda_cvec(m[1,0], m[2,0] );
    D2 = np.array( [[m[3,0], m[4,0]], [m[4,0], m[5,0]]] );
    invD2 = la.inv(D2);
    v = -np.matmul(invD2,d1);
    x0 = xc[ic,0]+v[0,0];
    y0 = yr[ir,0]+v[1,0];
    #E0= m(1) +   m(2)*v(1) +     m(3)*v(2) +     m(4)*v(1)*v(1)/2 +         m(5)*v(1)*v(2) +       m(6)*v(2)*v(2)/2;
    E0 = m[0,0] + m[1,0]*v[0,0] + m[2,0]*v[1,0] + m[3,0]*v[0,0]*v[0,0]/2.0 + m[4,0]*v[0,0]*v[1,0] + m[5,0]*v[1,0]*v[1,0]/2.0;
    covxy = (sigmad2)*2.0*invD2;
    status=1;
    return x0, y0, E0, covxy, status;

# support for gray-encoded numbers.  You must execute
# the following code in the main program:
# graydict=scipyio.loadmat("..\Data\gda_gray1table.mat");
# gray1=graydict["gray1"];
# gray1index=graydict["gray1index"];
# graypow2=graydict["pow2"];

def gda_int2gray(i):
    global gray1, gray1index; 
    # look up g in a table that lists all Gray codes
    # for the integers 0-65535 in order. The gray code
    # is returned as a length-16 character string of
    # ascii zeros and ones
    return gray1[i].astype(str);

def gda_gray2int(g):
    # The gray code is length-16 character string of
    # ascii zeros and ones representing the integers 0-65535
    global gray1, gray1index, graypow2;
    # hash g to a unique integer k
    j=0;
    k=0;
    if( isinstance(g,np.ndarray) ):
        gs = '';
        for i in range(16):
            gs = gs + g[i];
    elif( isinstance(g,str) ):
        gs=g;
    else:
        raise TypeError("gda_gray2int: %s not supported" % type(g));

    for i in gs.encode():
        k=k+(i-48)*graypow2[j,0];
        j=j+1;
    # look up j in the table
    return int(gray1index[k,0]-1);

# these two functions translate between integers
# and standard (no-gray) codes.
def gda_int2bin(n):
    L = 16;
    b = '';
    m = 1;
    for j in range(L):
        k = m & n;
        if( k ):
            b = b + "1";
        else:
            b = b + "0";
        m=m<<1;
    return b[::-1];

def gda_bin2int(b):
    L = 16;
    v = 0;
    m = 1;
    br = b[::-1];
    for j in range(L):
        if( br[j]=="1"):
            v = v + m
        m=m<<1;
    return v;
In [2]:
# gdapy11_01

# Least squares fit of straight line to d(z) and z(d).
 
print("gdapy11_01");

# auxially variable z and data d
z = gda_cvec(1, 2, 3, 4);
d = gda_cvec(1, 2, 3, 5);
N,i = np.shape(z);
 
# least squares fit to d(z)
M=2;
G=np.concatenate( (np.ones((N,1)), z), axis=1 );
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,d);
mest = la.solve(GTG,GTd);
dpre = np.matmul(G,mest);

# least squares fit to  z(d)
G2= np.concatenate( (np.ones((N,1)), d), axis=1);
GTG = np.matmul(G2.T,G2);
GTd = np.matmul(G2.T,z);
mest2 = la.solve(GTG,GTd);
dpre2 = np.matmul(G2,mest2);
 
# convert model parameters for z(d) case to d(z)
mest3=np.zeros((2,1));
mest3[0,0]=-mest2[0,0]/mest2[1,0];
mest3[1,0]=1/mest2[1,0];
dpre3=np.matmul(G,mest3);

print("d(z):   intercept %f slope %f"  % (mest[0,0],  mest[1,0])  );
print("dp(zp): intercept %f slope %f" %  (mest2[0,0], mest2[1,0]) );
print("zp(dp): intercept %f slope %f" %  (mest3[0,0], mest3[1,0]) );

# plot d(z) case
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,3,1);
plt.axis( [0, 6, 0, 6] );
plt.plot( z, d, "ko", lw=3);
plt.plot( z, dpre, "r-", lw=3);
plt.xlabel("z");
plt.ylabel("d");
plt.title("d(z)");
 
# plot z(d) case
ax2=plt.subplot(1,3,2);
plt.axis( [0, 6, 0, 6] );
plt.plot( d, z, "ko", lw=3);
plt.plot( d, dpre2, "g-", lw=3);
plt.xlabel("d");
plt.ylabel("z");
plt.title("z(d)");
 
# plot z(d) and transformed z(d) cases
ax3=plt.subplot(1,3,3);
plt.axis( [0, 6, 0, 6] );
plt.plot( z, d, "ko", lw=3);
plt.plot( z, dpre3, "g-", lw=3);
plt.plot( z, dpre, "r-", lw=2);
plt.xlabel("z");
plt.ylabel("d");
plt.title("d(z) and inv z(d)");
plt.show();
print("Caption: Linear fits to z(d) to data (circles): (left) fits to d(z),");
print("(middle) fit to z(d), (right) d(z) implied by fit to z(d)");
gdapy11_01
d(z):   intercept -0.500000 slope 1.300000
dp(zp): intercept 0.457143 slope 0.742857
zp(dp): intercept -0.615385 slope 1.346154
No description has been provided for this image
Caption: Linear fits to z(d) to data (circles): (left) fits to d(z),
(middle) fit to z(d), (right) d(z) implied by fit to z(d)
In [3]:
# gdapy11_02

# example of transforming distributions
 
print("gdapy11_02");

# model parameter, m
M=2000;
mmin = 0.0;
mmax = 1.0;
Dm = (mmax-mmin)/(M-1);
m = gda_cvec( mmin + Dm*np.linspace(0,M-1,M));
 
# uniform distribution
Pm = np.ones((M,1));
Pmp = np.zeros((M,1));
 
# transformation mp = m^2
mp=m;
Pmp[1:M,0:1] = np.divide(0.5*np.ones((M-1,1)),np.sqrt(mp[1:M,0:1]));
Pmp[0,0]=0.0; # set singularity to zero
 
# check area
Am = Dm*np.sum(Pm);
Amp = Dm*np.sum(Pmp);
print("Areas: p(m): %.2f p(mp) %.2f" % (Am, Amp) );
 
# plot p(m)
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,2,1);
plt.axis( [mmin, mmax, 0, 2] );
plt.plot( m, Pm, "r-", lw=3 );
plt.xlabel("m");
plt.ylabel("p(m)");
 
# plot p(m') where m'=m^2
ax1=plt.subplot(1,2,2);
plt.axis( [mmin, mmax, 0, 2] );
plt.plot( mp[1:M,0:1], Pmp[1:M,0:1], "b-", lw=3 );
plt.xlabel("mp");
plt.ylabel("p(mp)");
plt.show();
print("Caption: (left) Pdf p(m) is uniform, Pdf m(mp) where mp=m**2 is not uniform");
         
gdapy11_02
Areas: p(m): 1.00 p(mp) 0.98
No description has been provided for this image
Caption: (left) Pdf p(m) is uniform, Pdf m(mp) where mp=m**2 is not uniform
In [4]:
# gdapy11_03
 
# An inverse problem to determine the model parameter
# in a model with an exponential function, d(z)=m1*exp(m2 z)
# using (Case 1) a linearizing transformation
# and (Case 2) Newton's Method.
 
print("gdapy11_03");

# auxiallary variable z
N=20;
zmin = 0;
zmax = 1;
Dz = (zmax-zmin)/(N-1);
z = gda_cvec( zmin + Dz*np.linspace(0,N-1,N));
 
# data, a exponential d=m1*exp(m2 z)
M=2;
mtrue = gda_cvec( 0.9, -4.0 );
dtrue = mtrue[0,0]*np.exp(mtrue[1,0]*z);
 
# additive noise
sd=0.02;
bottom= 0.02;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1));
# but don't let data become negative
kk=np.where(dobs<0.0);
i = kk[0];
if( i.any() ):
    dobs[i,0:1]=bottom;

# linearizing transformation, solve by standard least squares
# d     = m1 exp( m2 z)
# ln(d) = ln(m1) + m2 * z
G=np.zeros((N,M));
G=np.concatenate( (np.ones((N,1)), z), axis=1 );
lndobs = np.log(dobs);
GTG = np.matmul(G.T,G);
GTd = np.matmul(G.T,lndobs);
mestlin = la.solve(GTG,GTd);
mestlin[0,0] = exp(mestlin[0,0]);
dprelin = mestlin[0,0]*np.exp(mestlin[1,0]*z);
 
# Newton's Method iterative solution
# d     = m1 exp( m2 z)
# dd/dm1 = exp( m2 z)
# dd/dm2 = m1 z exp( m2 z)
 
# initial guess is previous solution
mest = mestlin;
Nit=10; # maximum number of iterations
sdmmin = 1e-5; # terminate iteration early when increment dm
               # has sqrt(variance(dm))<sdmmin
for it in range(Nit):
    # Newton's Method
    dd = dobs - (mest[0,0]*np.exp(mest[1,0]*z));
    G = np.concatenate((np.exp(mest[1,0]*z),
                    mest[0,0]*np.multiply(z,np.exp(mest[1,0]*z))),axis=1);
    GTG = np.matmul(G.T,G);
    GTd = np.matmul(G.T,dd);
    dm = la.solve(GTG,GTd)
    mest = mest+dm;
    sdmsq = np.matmul(dm.T,dm); sdmsq=sdmsq[0,0];
    sdm = sqrt(sdmsq/M);  # sqrt(variance(dm))
    if( sdm < sdmmin ):
        break; # terminate iterations early

# least squares prediction
dpre = mest[0,0]*np.exp(mest[1,0]*z);
 
fig1 = plt.figure(1,figsize=(12,5));
 
# linear plot
ax1=plt.subplot(1,2,1);
plt.axis( [zmin, zmax, 0, 1] );
plt.plot( z, dobs, "ro", lw=3 );
plt.plot( z, dtrue, "r-", lw=3 );
plt.plot( z, dprelin, "b-", lw=3 );
plt.plot( z, dpre, "g-", lw=3 );
plt.xlabel("z");
plt.ylabel("d");
 
# log plot
ax2=plt.subplot(1,2,2);
plt.axis( [zmin, zmax, -3, 0] );
plt.plot( z, np.log10(dobs), "ro", lw=3 );
plt.plot( z, np.log10(dtrue), "r-", lw=3 );
plt.plot( z, np.log10(dprelin), "b-", lw=3 );
plt.plot( z, np.log10(dpre), "g-", lw=3 );
plt.xlabel("z");
plt.ylabel("log10(d)");
plt.show();
print("Caption (left) Fit to exponentially decating data (circles), showing");
print("true data (red), linearized fit (blue), iterative fit (green). (Right)");
print("same as left, except plot is semi-logarithmic");
print("  ");
 
print("True        m1: %.3f m2: %.3f" % (mtrue[0,0], mtrue[1,0]) );
print("Newtons     m1: %.3f m2: %.3f" % (mest[0,0], mest[1,0]) );
print("Linearizing m1: %.3f m2: %.3f" % (mestlin[0,0], mestlin[1,0]) );
gdapy11_03
No description has been provided for this image
Caption (left) Fit to exponentially decating data (circles), showing
true data (red), linearized fit (blue), iterative fit (green). (Right)
same as left, except plot is semi-logarithmic
  
True        m1: 0.900 m2: -4.000
Newtons     m1: 0.910 m2: -4.072
Linearizing m1: 0.845 m2: -3.859
In [5]:
# gdapy11_04

# 1D grid search for the one-parameter linear problem d=m1*z
 
print("gdapy11_04");

# auxiliary parameter z
N = 11;
zmin = 0;
zmax = 5.0;
Dz = (zmax-zmin)/(N-1);
z = gda_cvec( zmin + Dz*np.linspace(0,N-1,N));
 
# only one model parameter, m1
M=1;
 
# linear model:  d = m1*z
mtrue=2.5; # true model
dtrue=mtrue*z; # true data
# observed data is true data plus random noise
sd=2.0;
dobs=dtrue+np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# set up grid
Mg = 101;
mmin = 0.0;
mmax = 5.0;
Dm = (mmax-mmin)/(Mg-1);
m = gda_cvec(mmin + Dm*np.linspace(0,Mg-1,Mg));
 
# tabulate error E on a grid
E=np.zeros((Mg,1));
for i in range(Mg):
    dpre = m[i,0]*z;
    e = dobs - dpre;
    myE = np.matmul(e.T,e); myE=myE[0,0];
    E[i,0] = myE;

# find point of minimum Error
Emin = np.min(E);
iEmin = np.argmin(E);
mest=m[iEmin,0];
 
# plot Error and its minimum
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [mmin, mmax, 0, 1.1*np.max(E)] );
plt.plot( m, E, "k-", lw=3 );
plt.plot( mest, Emin, "ko", lw=3);
plt.plot( gda_cvec(mest, mest), gda_cvec(0.0, np.max(E)/50.0), "k-", lw=2);
plt.xlabel("m");
plt.ylabel("E");
plt.show();
print("Caption: One dimensional grid search of error curve E(m) (black curve)");
print("for minimum value (black circle).");
gdapy11_04
No description has been provided for this image
Caption: One dimensional grid search of error curve E(m) (black curve)
for minimum value (black circle).
In [6]:
# gdapy11_05

# E(m) for several hypothetial 1D non-linear problems
 
print("gdama11_05");

# auxiliary parameter z
N = 11;
zmin = 0;
zmax = 5.0;
Dz = (zmax-zmin)/(N-1);
z = gda_cvec( zmin + Dz*np.linspace(0,N-1,N));
 
# set up 1D grid search
Mg = 501;
mmin = 0;
mmax = 5;
Dm = (mmax-mmin)/(Mg-1);
m = gda_cvec(mmin + Dm*np.linspace(0,Mg-1,Mg));
 
# only one model parameter, m1
M=1;
 
# model 1
m1true=2.5;
d1true = np.sin(5*(m1true-2.5)*z)-((m1true-2.5))*z;
sd=2.0;
d1obs=d1true+np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# tabulate error on the grid
E1=np.zeros((Mg,1));
for i in range(Mg):
    d1pre = np.sin(5*(m[i,0]-2.5)*z)-((m[i,0]-2.5))*z;
    e = d1obs - d1pre;
    myE1 = np.matmul(e.T,e); myE1=myE1[0,0];
    E1[i,0] = myE1;

# find point of minimum error
E1min = np.min(E1);
iE1min = np.argmin(E1);
m1est=m[iE1min,0];
 
fig1 = plt.figure(1,figsize=(12,5));
 
# plot E(m)
ax1=plt.subplot(2,2,1);
plt.axis( [mmin, mmax, 0.0, 1.1*np.max(E1)] );
plt.plot( m, E1, "k-", lw=3 );
plt.xlabel("m");
plt.ylabel("E");
 
# model 2
m1true=2.5;
d1true = ((np.abs((m1true-2))-0.5)*z);
sd=2.0;
d1obs=d1true+np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# tabulate error on the grid
E1=np.zeros((Mg,1));
for i in range(Mg):
    d1pre = ((abs((m[i,0]-2))-0.5)*z);
    e = d1obs - d1pre;
    myE1 = np.matmul(e.T,e); myE1=myE1[0,0]
    E1[i,0] = myE1;

# find point of minimum error
E1min = np.min(E1);
iE1min = np.argmin(E1);
m1est=m[iE1min,0];
 
# plot E(m)
ax1=plt.subplot(2,2,2);
plt.axis( [mmin, mmax, 0.0, 1.1*np.max(E1)] );
plt.plot( m, E1, "k-", lw=3 );
plt.xlabel("m");
plt.ylabel("E");
 
# model 3
m1true=2.5;
d1true = np.sin(5.0*(m1true+z));
sd=1.0;
d1obs=d1true+np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# tabulate error on the grid
E1=np.zeros((Mg,1));
for i in range(Mg):
    d1pre = np.sin(5.0*(m[i,0]+z));
    e = d1obs - d1pre;
    myE1 = np.matmul(e.T,e); myE1=myE1[0,0];
    E1[i,0] = myE1;

# find point of minumum error
E1min = np.min(E1);
iE1min = np.argmin(E1);
m1est=m[iE1min,0];
 
# plot E(m)
ax1=plt.subplot(2,2,3);
plt.axis( [mmin, mmax, 0.0, 1.1*np.max(E1)] );
plt.plot( m, E1, "k-", lw=3 );
plt.xlabel("m");
plt.ylabel("E");
 
# model 4
m1true=2.5;
t = abs(m1true-2.0) + abs(m1true-3.0);
d1true = np.exp(-np.power(t*z/10.0,2));
sd=1.0;
d1obs=d1true+np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# tabulate error on the grid
E1=np.zeros((Mg,1));
for i in range(Mg):
    t = abs(m[i,0]-2.0) + abs(m[i,0]-3.0);
    d1pre = np.exp(-np.power(t*z/10,2));
    e = d1obs - d1pre;
    myE1 = np.matmul(e.T,e); myE1=myE1[0,0];
    E1[i,0] = myE1;

# find point of minumum error
E1min = np.min(E1);
iE1min = np.argmin(E1);
m1est=m[iE1min,0];
 
# plot E(m)
ax1=plt.subplot(2,2,4);
plt.axis( [mmin, mmax, 0.0, 1.1*np.max(E1)] );
plt.plot( m, E1, "k-", lw=3 );
plt.xlabel("m");
plt.ylabel("E");
plt.show();
print("Caption: Hypothetical error surfaces");
gdama11_05
No description has been provided for this image
Caption: Hypothetical error surfaces
In [7]:
# gdapy11_06

# E(m) for several 2D non-linear problems
 
print("gdapy11_06");

# auxiliary parameter z
N = 11;
zmin = 0.0;
zmax = 1.0;
Dz = (zmax-zmin)/(N-1);
z = gda_cvec(zmin + Dz*np.linspace(0,N-1,N));
za = np.power(z,0.10);
zb = np.power(z,0.20);
 
# 2D grid of model parameters m1 and m2
# variable m1
Mg1 = 101;
m1min = 0.0;
m1max = 1.0;
Dm1 = (m1max-m1min)/(Mg1-1);
m1 = gda_cvec(m1min + Dm1*np.linspace(0,Mg1-1,Mg1));
 
# variable m2
Mg2 = 101;
m2min = 0.0;
m2max = 1.0;
Dm2 = (m2max-m2min)/(Mg2-1);
m2 = gda_cvec(m2min + Dm2*np.linspace(0,Mg2-1,Mg2));
 
# Case 1
 
# true model
m1true = 0.4;
m2true = 0.6;
 
# true data
# MATLAB(R): dtrue = ((m1true-0.2)*za-((m2true-0.3)*zb).^2).*z;
A = (m1true-0.2)*za;
B = np.power((m2true-0.3)*zb,2);
dtrue = np.multiply(A+B,z);
 
# observed data is true data plus random noise
sd=0.01;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# tabulate error on the grid, keeping track of
# minimum error along the way
E=np.zeros((Mg1,Mg2));
first=1;
for i in range(Mg1):
    for j in range(Mg2):
        A = (m1[i,0]-0.2)*za;
        B = np.power((m2[j,0]-0.3)*zb,2);
        dpre = np.multiply(A+B,z);
        e = dobs - dpre;
        myE = np.matmul(e.T,e); myE=myE[0,0];
        E[i,j] = myE;
        if( first==1 ):
            Emin=E[i,j];
            Emini=i;
            Eminj=j;
            dpresave = dpre;
            first=0;
        elif( E[i,j]<Emin ):
            Emin=E[i,j];
            Emini=i;
            Eminj=j;
            dpresave = dpre;

print("Case 1:");
print("True solution at     %f %f" % (m1true,m2true));
print("Emin=E(%d,%d): %f at %f %f" % (Emini, Eminj, Emin, m1[Emini,0], m2[Eminj,0]) );
print(" ");

# plot error surface E(m1,m2)
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m1min, m1max, m2min, m2max] );
left=m1min;
right=m1max;
bottom=m2min;
top=m2max;
dmin=np.min(E);
dmax=np.max(E);
# for checking orientation
# E=dmin*np.ones((Mg1,Mg2)); i1=4; i2=4; E[i1,i2]=dmax;
# E=dmin*np.ones((Mg1,Mg2)); i1=4; i2=Mg2-4; E[i1,i2]=dmax;
# E=dmin*np.ones((Mg1,Mg2)); i1=Mg1-4; i2=4; E[i1,i2]=dmax;
plt.imshow( np.flipud(E), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.gca().invert_yaxis();
plt.plot( m2[Eminj,0], m1[Emini,0], "wo", lw=3 );
plt.xlabel("m2");
plt.ylabel("m1");
plt.show();
print("Caption: Error surface E(m1,m2)");

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [zmin, zmax, np.min(dtrue), np.max(dtrue)] );
plt.plot( z, dtrue, "k-", lw=3 );
plt.plot( z, dobs, "ko", lw=3 );
plt.plot( z, dpresave, "r-", lw=1 );
plt.xlabel("z");
plt.ylabel("d");
plt.show();
print("Caption: True data (black curve), observed (circles), predicted (red).")
print("  ");

# Case 2
 
# true model parameters
m1true = 0.4;
m2true = 0.6;
 
# true data
dtrue = sin(12.0*pi*(m1true-0.5))*z+np.power(m2true*z,0.5);
 
# observed data is true data plus random noise
sd=0.1;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# Tabulate error E(m1,m2) on the grid, searching
# for the point of minimum error along the way
E=np.zeros((Mg1,Mg2));
first=1;
for i in range(Mg1):
    for j in range(Mg2):
        dpre = sin(12.0*pi*(m1[i,0]-0.5))*z+np.power(m2[j,0]*z,0.5);
        e = dobs - dpre;
        myE = np.matmul(e.T,e); myE=myE[0,0];
        E[i,j] = myE;
        if( first==1 ):
            Emin=E[i,j];
            Emini=i;
            Eminj=j;
            dpresave=dpre;
            first=0;
        elif( E[i,j]<Emin ):
            Emin=E[i,j];
            Emini=i;
            Eminj=j;
            dpresave=dpre;

print("Case 2:");
print("True solution at     %f %f" % (m1true,m2true));
print("Emin=E(%d,%d): %f at %f %f" % (Emini, Eminj, Emin, m1[Emini,0], m2[Eminj,0]) );

# plot error surface E(m1,m2)
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m1min, m1max, m2min, m2max] );
left=m1min;
right=m1max;
bottom=m2min;
top=m2max;
dmin=np.min(E);
dmax=np.max(E);
# for checking orientation
# E=dmin*np.ones((Mg1,Mg2)); i1=4; i2=4; E[i1,i2]=dmax;
# E=dmin*np.ones((Mg1,Mg2)); i1=4; i2=Mg2-4; E[i1,i2]=dmax;
# E=dmin*np.ones((Mg1,Mg2)); i1=Mg1-4; i2=4; E[i1,i2]=dmax;
plt.imshow( np.flipud(E), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.gca().invert_yaxis();
plt.plot( m2[Eminj,0], m1[Emini,0], "wo", lw=3 );
plt.xlabel("m2");
plt.ylabel("m1");
plt.show();
print("Caption: Error surface E(m1,m2)");

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [zmin, zmax, np.min(dtrue), np.max(dtrue)] );
plt.plot( z, dtrue, "k-", lw=3 );
plt.plot( z, dobs, "ko", lw=3 );
plt.plot( z, dpresave, "r-", lw=1 );
plt.xlabel("z");
plt.ylabel("d");
plt.show();
print("Caption: True data (black curve), observed (circles), predicted (red).")
gdapy11_06
Case 1:
True solution at     0.400000 0.600000
Emin=E(39,62): 0.000508 at 0.390000 0.620000
 
No description has been provided for this image
Caption: Error surface E(m1,m2)
No description has been provided for this image
Caption: True data (black curve), observed (circles), predicted (red).
  
Case 2:
True solution at     0.400000 0.600000
Emin=E(35,75): 0.147524 at 0.350000 0.750000
No description has been provided for this image
Caption: Error surface E(m1,m2)
No description has been provided for this image
Caption: True data (black curve), observed (circles), predicted (red).
In [8]:
# gdapy11_07
 
# Grid Search method applied to inverse problem
# d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2)
 
print("gdapy11_07");

# data are in a simple auxillary variable, x
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = gda_cvec(Dx*np.linspace(0,N-1,N));
 
# two model parameters
M=2;
 
# true model parameters
mt = gda_cvec(1.21, 1.54);
 
# d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2);
w0=20;
dtrue = np.sin(w0*mt[0,0]*x) + mt[0,0]*mt[1,0]; # true data
sd=0.4;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1)); # observed data
 
# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, xmax, 0, 4 ] );
plt.plot(x,dtrue,"k-",lw=3);
plt.plot(x,dobs,"ko",lw=3);
plt.xlabel("x");
plt.ylabel("d");
 
# Define 2D grid
L = 101;
Dm = 0.02;
m1min=0.0;
m2min=0.0;
m1a = gda_cvec(m1min+Dm*np.linspace(0,L-1,L));
m2a = gda_cvec(m2min+Dm*np.linspace(0,L-1,L));
m1max = m1a[L-1,0];
m2max = m2a[L-1,0];
 
# tabulate error E on grid.
# Note m1 varies along rows of E and m2 along columns
E = np.zeros((L,L));
for j in range(L):
    for k in range(L):
        dpre = np.sin(w0*m1a[j,0]*x) + m1a[j,0]*m2a[k,0];
        myE = np.matmul((dobs-dpre).T,(dobs-dpre)); myE=myE[0,0];
        E[j,k] = myE;

# Example of use of Esurface() function
# It refines the minimum using a quadratic approximation
# and provides an error estimate using the 2nd derivative
# note that the "column" variable is sent and returned
# before the "row" variable in this function
m2est, m1est, E0, cov21, status = gda_Esurface( m2a, m1a, E, sd**2 );
sigmam1 = sqrt( cov21[1,1] );
sigmam2 = sqrt( cov21[0,0] );

# evaluate prediction and plot it
dpre = np.sin(w0*m1est*x) + m1est*m2est;
plt.plot(x,dpre,"r-",lw=3);
plt.show();
print("Caption: Exemplary curve fit wia grid search, showing true data (black curve)");
print("observed data (circles) and predicted data (red curve).");

# plot error surface E(m1,m2)
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m2min, m2max, m1min, m1max] );
left=m1min;
right=m1max;
bottom=m2min;
top=m2max;
dmin=np.min(E);
dmax=np.max(E);
# for checking orientation
# E=dmin*np.ones((Mg1,Mg2)); i1=4; i2=4; E[i1,i2]=dmax;
# E=dmin*np.ones((Mg1,Mg2)); i1=4; i2=Mg2-4; E[i1,i2]=dmax;
# E=dmin*np.ones((Mg1,Mg2)); i1=Mg1-4; i2=4; E[i1,i2]=dmax;
plt.imshow( np.flipud(E), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.gca().invert_yaxis();
plt.plot( m2[Eminj,0], m1[Emini,0], "wo", lw=3 );
plt.xlabel("m2");
plt.ylabel("m1");
plt.plot( m2est, m1est, "wo", lw=3);
plt.plot( mt[1,0], mt[0,0], "go", lw=3);
# plot error bars
plt.plot( gda_cvec(m2est, m2est), gda_cvec(m1est-2.0*sigmam1, m1est+2.0*sigmam1), "r-", lw=2);
plt.plot( gda_cvec(m2est-2*sigmam2, m2est+2*sigmam2), gda_cvec(m1est, m1est), "r-", lw=2);
plt.show();
print("Caption: Error surface, with true solution (green circle), estimate solution");
print("(white circles) and 95 percent confidence intervals (red bars), superimposed");
gdapy11_07
No description has been provided for this image
Caption: Exemplary curve fit wia grid search, showing true data (black curve)
observed data (circles) and predicted data (red curve).
No description has been provided for this image
Caption: Error surface, with true solution (green circle), estimate solution
(white circles) and 95 percent confidence intervals (red bars), superimposed
In [9]:
# gdapy11_08

# Newton's Method solution to exemplary 1D non-linear
# inverse problem, d(z) = sin(k*m*z)*z;
# This version uses a starting guess for the model
# parameter m1 that leads to convergence to the global
# minimum of the error surface (and hence to the correct
# value of the model parameter).
 
print("gdapy11_08");

# z-axis
N = 21;
zmin = 0;
zmax = 5.0;
Dz = (zmax-zmin)/(N-1);
z = gda_cvec(zmin + Dz*np.linspace(0,N-1,N));
zp = np.power(z,0.5);
 
# define 1D grid for model parameter m1
Mg = 501;
mmin = 0.0;
mmax = 5.0;
Dm = (mmax-mmin)/(Mg-1);
m = gda_cvec(mmin + Dm*np.linspace(0,Mg-1,Mg));
 
# only one model parameter, m1
M=1;
 
# true model parameter
mtrue=3.0;
 
# true data
k=1.0;
dtrue = np.multiply(np.sin(k*mtrue*zp),z);
 
# observed data is true data plus random noise
sd=2.0;
dobs=dtrue+np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# tabulate error E(m) on the grid (for plotting purposes only)
E1=np.zeros((Mg,1));
for i in range(Mg):
    dpre = np.multiply(np.sin(k*m[i,0]*zp),z);
    e = dobs - dpre;
    myE1 = np.matmul(e.T,e); myE1=myE1[0,0];
    E1[i,0] = myE1;

# find global minimum
E1min = np.min(E1);
iE1min = np.argmin(E1);
m1est=m[iE1min,0];
 
# Newton's method
# derivative
# d = sin(k*m(i).*zp).*z;
# dd/dm = k.*zp.*z.*cos(k*m(i).*zp);
# d = d(m0) + (dd/dm)|m0 (m-m0)
 
# tabulate a quadratic version of the error
# correspinding to a linearized version of
# the inverse problem (linearized around a point m0)
E2=np.zeros((Mg,1));
m0 = 2.5;
im0 = int(floor((m0-mmin)/Dm));
d0 = np.multiply(np.sin(k*m0*zp),z);
dddm = np.multiply(k*np.multiply(zp,z),np.cos(k*m0*zp));
for i in range(Mg):
    dpre = d0 + dddm * (m[i,0]-m0); # linearized approximation
    e = dobs - dpre;
    myE2 = np.matmul(e.T,e); myE2=myE2[0,0];
    E2[i,0] = myE2;
    
# find minimum of error surface
E2min = np.min(E2);
iE2min = np.argmin(E2);
m2est=m[iE2min,0];
 
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mminp=1; mmaxp=4;
plt.axis( [mmin, mmax, 0, 1.1*np.max(E1)] );
 # true error surface, with circle at minimum, and line dropped to axis
plt.plot( m, E1, "k-", lw=3 );
plt.plot( m1est, E1min, "ko", lw=3 );
plt.plot( gda_cvec(m1est, m1est), gda_cvec(0.0, E1min), "k:", lw=2);
 
# parabiolic approximation, with circle at minimum, and line dropped to axis
plt.plot( m, E2, "r:", lw=3 );
plt.plot( m2est, E2min, "ro", lw=3 );
plt.plot( gda_cvec(m2est, m2est), gda_cvec(0, E2min), "r:", lw=2);
 
plt.plot( m0, E2[im0,0], "ko", lw=3 );
plt.plot( gda_cvec(m0, m0), gda_cvec(0.0, E2[im0,0]), "k:", lw=2);
 
plt.xlabel("m");
plt.ylabel("E");
plt.show();

print("Caption: Exemplary error surface E(m) (back curve) with minimum value at m=3");
print("(black circle) showing parabola tangegent to curve at m=2.5 (red curve), with");
print("minimum at m=3.3 (red circle).")
gdapy11_08
No description has been provided for this image
Caption: Exemplary error surface E(m) (back curve) with minimum value at m=3
(black circle) showing parabola tangegent to curve at m=2.5 (red curve), with
minimum at m=3.3 (red circle).
In [10]:
# gdapy11_09

# Newton's Method solution to exemplary 1D non-linear
# inverse problem, d(z) = sin(k*m*z)*z
# This version uses a starting guess for the model
# parameter m1 that leads to convergence to a local
# minimum of the error surface (and hence to an incorrect
# value of the model parameter)
 
print("gdapy11_09");

# z-axis
N = 21;
zmin = 0;
zmax = 5.0;
Dz = (zmax-zmin)/(N-1);
z = gda_cvec(zmin + Dz*np.linspace(0,N-1,N));
zp = np.power(z,0.5);
 
# define 1D grid for model parameter m1
Mg = 501;
mmin = 0.0;
mmax = 5.0;
Dm = (mmax-mmin)/(Mg-1);
m = gda_cvec(mmin + Dm*np.linspace(0,Mg-1,Mg));
 
# only one model parameter, m1
M=1;
 
# true model parameter
mtrue=3.0;
 
# true data
k=1.0;
dtrue = np.multiply(np.sin(k*mtrue*zp),z);
 
# observed data is true data plus random noise
sd=2.0;
dobs=dtrue+np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# tabulate error E(m) on the grid (for plotting purposes only)
E1=np.zeros((Mg,1));
for i in range(Mg):
    dpre = np.multiply(np.sin(k*m[i,0]*zp),z);
    e = dobs - dpre;
    myE1 = np.matmul(e.T,e); myE1=myE1[0,0];
    E1[i,0] = myE1;

# find global minimum
E1min = np.min(E1);
iE1min = np.argmin(E1);
m1est=m[iE1min,0];
 
# Newton's method
# derivative
# d = sin(k*m(i).*zp).*z;
# dd/dm = k.*zp.*z.*cos(k*m(i).*zp);
# d = d(m0) + (dd/dm)|m0 (m-m0)
 
# tabulate a quadratic version of the error
# correspinding to a linearized version of
# the inverse problem (linearized around a point m0)
E2=np.zeros((Mg,1));
m0 = 1.0;
im0 = int(floor((m0-mmin)/Dm));
d0 = np.multiply(np.sin(k*m0*zp),z);
dddm = np.multiply(k*np.multiply(zp,z),np.cos(k*m0*zp));
for i in range(Mg):
    dpre = d0 + dddm * (m[i,0]-m0); # linearized approximation
    e = dobs - dpre;
    myE2 = np.matmul(e.T,e); myE2=myE2[0,0];
    E2[i,0] = myE2;
    
# find minimum of error surface
E2min = np.min(E2);
iE2min = np.argmin(E2);
m2est=m[iE2min,0];
 
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mminp=1; mmaxp=4;
plt.axis( [mmin, mmax, 0, 1.1*np.max(E1)] );
 # true error surface, with circle at minimum, and line dropped to axis
plt.plot( m, E1, "k-", lw=3 );
plt.plot( m1est, E1min, "ko", lw=3 );
plt.plot( gda_cvec(m1est, m1est), gda_cvec(0.0, E1min), "k:", lw=2);
 
# parabiolic approximation, with circle at minimum, and line dropped to axis
plt.plot( m, E2, "r:", lw=3 );
plt.plot( m2est, E2min, "ro", lw=3 );
plt.plot( gda_cvec(m2est, m2est), gda_cvec(0, E2min), "r:", lw=2);
 
plt.plot( m0, E2[im0,0], "ko", lw=3 );
plt.plot( gda_cvec(m0, m0), gda_cvec(0.0, E2[im0,0]), "k:", lw=2);
 
plt.xlabel("m");
plt.ylabel("E");
plt.show();

print("Caption: Exemplary error surface E(m) (back curve) with minimum value");
print("at m=3 (black circle) showing parabola tangegent to curve at m=1");
print("(red curve), with minimum at m=3.3 (red circle).");
gdapy11_09
No description has been provided for this image
Caption: Exemplary error surface E(m) (back curve) with minimum value
at m=3 (black circle) showing parabola tangegent to curve at m=1
(red curve), with minimum at m=3.3 (red circle).
In [11]:
# gdapy11_10
# Newton's method applied to inverse problem
# d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2)

print("gdapy11_10");

# x-axis
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = gda_cvec(Dx*np.linspace(0,N-1,N));
 
# true model parameters
M=2;
mt = gda_cvec(1.21, 1.54);
 
# d(x)=g(x, m1, m2)   with   g=sin(w0*m(1)*x) + m(1)*m(2)
w0=20.0;
dtrue = np.sin(w0*mt[0,0]*x) + mt[0,0]*mt[1,0]; # true data
sd=0.4;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1)); # observed data
 
# 2D grid, for plotting purposes only
L = 101;
Dm = 0.02;
m1min=0.0;
m2min=0.0;
m1a = gda_cvec(m1min+Dm*np.linspace(0,L-1,L));
m2a = gda_cvec(m2min+Dm*np.linspace(0,L-1,L));
m1max = m1a[L-1,0];
m2max = m2a[L-1,0];
 
# tabulate error, E, on grid for plotting purposed only
E = np.zeros((L,L));
for j in range(L):
    for k in range(L):
        dpre = np.sin(w0*m1a[j,0]*x) + m1a[j,0]*m2a[k,0];
        myE = np.matmul((dobs-dpre).T,(dobs-dpre)); myE=myE[0,0];
        E[j,k] = myE;

# Newton's method, calculate derivatives
# y = sin( w0 m1  x) + m1 m2;
# dy/dm1 = w0 x cos( w0 m1 x) + m2
# dy/dm2 = m2
 
# initial guess and corresponding error
mg=gda_cvec(1.1,0.95);
dg = np.sin(w0*mg[0,0]*x) + mg[0,0]*mg[1,0];
Eg = np.matmul((dobs-dg).T,(dobs-dg)); Eg=Eg[0,0];

# save solution and minimum error as a function of iteration number
Niter=10;
Ehis=np.zeros((Niter,1));
m1his=np.zeros((Niter,1));
m2his=np.zeros((Niter,1));
Ehis[0,0]=Eg;
m1his[0,0]=mg[0,0];
m2his[0,0]=mg[1,0];
 
# iterate to improve initial guess
for k in range(1,Niter):
 
    # predict data at current guess for solution
    dg = np.sin( w0*mg[0,0]*x) + mg[0,0]*mg[1,0];
    dd = dobs-dg;
    Eg=np.matmul(dd.T,dd); Eg=Eg[0,0]; # current error
    Ehis[k,0]=Eg; # save error history
    
    # build linearized data kernel
    G = np.zeros((N,2));
    G[0:N,0:1] = w0 * np.multiply(x,np.cos(w0*mg[0,0]*x)) + mg[1,0];
    G[0:N,1:2] = mg[1,0]*np.ones((N,1));
    
    # least squares solution for improvement to m
    GTG = np.matmul(G.T,G);
    GTd = np.matmul(G.T,dd);
    dm = la.solve(GTG,GTd);
    
    # update solution
    mg = mg+dm;
    
    # save solution history
    m1his[k,0]=mg[0,0];
    m2his[k,0]=mg[1,0];

# estimated solution is from final iteration
m1est = m1his[Niter-1,0];
m2est = m2his[Niter-1,0];
 
# evaluate prediction and plot it
dpre = np.sin(w0*m1est*x) + m1est*m2est;

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0.0, xmax, 0.0, 4.0 ] );
plt.plot(x,dtrue,"k-",lw=3);
plt.plot(x,dobs,"ko",lw=3);
plt.xlabel("x");
plt.ylabel("d");
plt.plot(x,dpre,"r-",lw=3);
plt.show();
print("Caption: Exemplary curve fit with Newtons methiod, showing observed");
print("data (circles), true data (black curve) and predicted data (red curve)");

# plot error surface E(m1,m2)
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m2min, m2max, m1min, m1max] );
left=m1min;
right=m1max;
bottom=m2min;
top=m2max;
dmin=np.min(E);
dmax=np.max(E);
plt.imshow( np.flipud(E), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.gca().invert_yaxis();
plt.plot( m2his[0,0], m1his[0,0], "wo", lw=3 );
plt.plot( m2his, m1his, "w-", lw=1 );
plt.plot( mt[1,0], mt[0,0], "go", lw=2 );
plt.plot( m2est, m1est, "ro", lw=2 );
plt.xlabel("m2");
plt.ylabel("m1");
plt.show();
print("Caption: Error surface (colors), with initial guess of solution")
print("(white circle), trajectory (white line), estimated solution (red");
print("circle) and true solution (green).")

# plot history
w = gda_cvec(np.linspace(0,Niter-1,Niter));
fig3 = plt.figure(3,figsize=(12,5));

# plot m1 history
ax1=plt.subplot(3,1,1);
plt.axis( [0, Niter, 0.0, 2.0 ] );
plt.plot(w,m1his,"k-",lw=3);
plt.plot(gda_cvec(0.0,Niter),gda_cvec(mt[0,0],mt[0,0]),'r:',lw=2);
plt.xlabel("iteration i");
plt.ylabel("m1");

# plot m2 history
ax1=plt.subplot(3,1,2);
plt.axis( [0, Niter, 0.0, 2.0 ] );
plt.plot(w,m2his,"k-",lw=3);
plt.plot(gda_cvec(0.0,Niter),gda_cvec(mt[1,0],mt[1,0]),'r:',lw=2);
plt.xlabel("iteration i");
plt.ylabel("m2");

# plot E history
ax1=plt.subplot(3,1,3);
plt.axis( [0, Niter, 0.0, 1.1*np.max(Ehis) ] );
plt.plot(w,Ehis,"k-",lw=3);
plt.xlabel("iteration i");
plt.ylabel("E");
plt.show();
print("Caption: Newtons method trajectory of (top) m1, (middle) m2");
print("and (bottom) E (black curves), with true value (red).");
gdapy11_10
No description has been provided for this image
Caption: Exemplary curve fit with Newtons methiod, showing observed
data (circles), true data (black curve) and predicted data (red curve)
No description has been provided for this image
Caption: Error surface (colors), with initial guess of solution
(white circle), trajectory (white line), estimated solution (red
circle) and true solution (green).
No description has been provided for this image
Caption: Newtons method trajectory of (top) m1, (middle) m2
and (bottom) E (black curves), with true value (red).
In [12]:
# gdapy11_11

# draws a simple Normal pdf p(x1,x2)

print("gdapy11_11");

# x1-axis 
Nx1 = 101;
x1min = 0.0;
x1max = 5.0;
Dx1 = (x1max-x1min)/(Nx1-1);
x1 = gda_cvec(x1min + Dx1*np.linspace(0,Nx1-1,Nx1));
 
# x2-axis 
Nx2 = 101;
x2min = 0.0;
x2max = 5.0;
Dx2 = (x2max-x2min)/(Nx2-1);
x2 = gda_cvec(x2min + Dx2*np.linspace(0,Nx2-1,Nx2));
 
# pdf
P1=np.zeros((Nx1,Nx2));
x1bar = 2.25;
x2bar = 2.08;
bar = gda_cvec(x1bar, x2bar);
sx1 = 0.5;
sx2 = 1.0;
C1 = np.diag( gda_cvec(sx1**2, sx2**2).ravel() );
CI1 = la.inv(C1);
DC1 = abs(la.det(C1));
norm1 = (1.0/(2.0*pi)) * (1.0/sqrt(DC1));
 
for i in range(Nx1):
    for j in range(Nx2):
        x = gda_cvec(x1[i,0], x2[j,0] ) - bar;
        R2 = np.matmul(x.T,np.matmul(CI1,x)); R2=R2[0,0];
        P1[i,j] = norm1*np.exp( -0.5*R2 );
        
# plot pdf p(x1,x2)
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [x1min, x1max, x2min, x2max] );
left=x1min;
right=x1max;
bottom=x2min;
top=x2max;
dmin=0.0;
dmax=np.max(P1);
plt.imshow( np.flipud(P1), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.gca().invert_yaxis();
plt.xlabel("x2");
plt.ylabel("x1");
plt.plot( x2bar, x1bar, "wo", lw=3);
plt.plot( gda_cvec(x2bar, x2bar), gda_cvec(x1min, x1min+0.15), "w-", lw=3);
plt.plot( gda_cvec(x2min, x2min+0.15), gda_cvec(x1bar, x1bar), "w-", lw=3);
plt.show();
print("Caption: Exemplary Normal pdf p(x1,x2) with mean (white circle) superimposed");
gdapy11_11
No description has been provided for this image
Caption: Exemplary Normal pdf p(x1,x2) with mean (white circle) superimposed
In [13]:
# gdapy11_12

# Plots a simple Normal p.d.f. p(x1,x1) with a parameteric
# curve superimposed and then samples the p.d.f. under the
# curve and makes a second plot of p.d.f. as a function of
# arc-length s along the curve. In this version, the curve
# is very smooth and p(s) is bell-shaped.
 
# 2D grid in (x1,x2)

print("gdapy11_12");

# x1-axis
Nx1 = 101;
x1min = 0.0;
x1max = 5.0;
Dx1 = (x1max-x1min)/(Nx1-1);
x1 = gda_cvec(x1min + Dx1*np.linspace(0,Nx1-1,Nx1));
 
# x2-axis
Nx2 = 101;
x2min = 0.0;
x2max = 5.0;
Dx2 = (x2max-x2min)/(Nx2-1);
x2 = gda_cvec(x2min + Dx2*np.linspace(0,Nx2-1,Nx2));
 
# parameters for a Normal pdf
P1=np.zeros((Nx1,Nx2));
x1bar = 2.25;
x2bar = 2.08;
bar = gda_cvec(x1bar, x2bar);
sx1 = 0.5;
sx2 = 1.0;
C1 = np.diag( gda_cvec(sx1**2, sx2**2).ravel() );
CI1 = la.inv(C1);
DC1 = abs(la.det(C1));
norm1 = (1.0/(2.0*pi)) * (1.0/sqrt(DC1));
 
# tabulate the Normal pdf on the grid
for i in range(Nx1):
    for j in range(Nx2):
        x = gda_cvec(x1[i,0], x2[j,0]) - bar;
        R2 = np.matmul(x.T,np.matmul(CI1,x)); R2=R2[0,0];
        P1[i,j] = norm1*exp( -0.5*R2 );

# axis for parametric curve
Ns = 51;
smin = 0.0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = gda_cvec(smin + Ds*np.linspace(0,Ns-1,Ns));

# parameric curve
x2p = 1.0+s+np.sin(2.0*pi*s/smax)-2.0*np.power(s/smax,2);
x1p = s;
 
# pdf along parameteric curve
ipx1 = (np.floor((x1p-x1min)/Dx1)).astype(int);
ipx2 = (np.floor((x2p-x2min)/Dx2)).astype(int);
# insure indices in range
kk = np.where(ipx1<0);
i=kk[0];
if( i.any() ):
    ipx1[i,0:1]=0;
kk = np.where(ipx1>=Nx1);
i=kk[0];
if( i.any() ):
    ipx1[i,0:1]=Nx1-1;
         
kk = np.where(ipx2<0);
i = kk[0];
if( i.any() ):
    ipx2[i,0:1]=0;
kk=np.where(ipx2>=Nx2);
i=kk[0];
if( i.any() ):
    ipx2[i,0]=Nx2-1;
         
Ps=np.zeros((Ns,1));
# evaluate P at indices
for i in range(Ns):
    Ps[i,0] = P1[ipx1[i,0], ipx2[i,0] ];
    
    
# calculate mean of pdf p(s)
norm = Ds*np.sum(Ps);
Ps = Ps/norm;
smean = Ds*np.sum(np.multiply(s,Ps));

# maximum along curve
Pmax = np.max(Ps);
ismax =np.argmax(Ps);
x1smax = x1p[ismax,0];
x2smax = x2p[ismax,0];
smaxP = x1smax;

# plot pdf p(x1,x2)
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [x1min, x1max, x2min, x2max] );
left=x1min;
right=x1max;
bottom=x2min;
top=x2max;
dmin=0.0;
dmax=np.max(P1);
plt.imshow( np.flipud(P1), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.gca().invert_yaxis();
plt.xlabel("x2");
plt.ylabel("x1");
plt.plot( x2bar, x1bar, "wo", lw=3 );
plt.plot( gda_cvec(x2bar, x2bar), gda_cvec(x1min, x1min+0.1), "w-", lw=3 );
plt.plot( gda_cvec(x2min, x2min+0.1), gda_cvec(x1bar, x1bar), "w-", lw=3 );
plt.plot( x2p, x1p, "w-", lw=3 );
plt.plot( x2smax, x1smax, "ko", lw=4);
plt.plot( gda_cvec(x2min, x2smax), gda_cvec(x1smax, x1smax), "w:", lw=2);
plt.plot( gda_cvec(x2smax, x2smax), gda_cvec(x1min, x1smax), "w:", lw=2);
plt.show();
print("Caption: Exemplary normal pdf (colors) with mean (white circle),");
print("curve f(x)=0 (white curve), and maximum likelihood point (black");
print("circle), superimposed.");

# plot of 1D pdf (along curve)
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [smin, smax, 0, 1.1*np.max(Ps)] );
plt.plot( s, Ps, "b-", lw=3 );
plt.xlabel("s");
plt.ylabel("p(s)");
plt.plot( gda_cvec(smaxP, smaxP), gda_cvec(0.0, 0.03), "k-", lw=3);
plt.plot( smaxP, Pmax, "bo", lw=3);
plt.plot( gda_cvec(smean, smean), gda_cvec(0, 0.02), "k-", lw=3);
plt.show();
print("Caption: Probability along the curve p(x)=0, with the maximum");
print("likelihood point shown (circle and long bar) and mean (short bar).");
gdapy11_12
No description has been provided for this image
Caption: Exemplary normal pdf (colors) with mean (white circle),
curve f(x)=0 (white curve), and maximum likelihood point (black
circle), superimposed.
No description has been provided for this image
Caption: Probability along the curve p(x)=0, with the maximum
likelihood point shown (circle and long bar) and mean (short bar).
In [14]:
# gdapy11_13

# Plots a simple Normal p.d.f. p(x1,x1) with a parameteric
# curve superimposed and then samples the p.d.f. under the
# curve and makes a second plot of p.d.f. as a function of
# arc-length s along the curve. In this version, the curve
# is complicated and p(s) is bimodal.
 
# 2D grid in (x1,x2)

print("gdapy11_13");

# x1-axis
Nx1 = 101;
x1min = 0.0;
x1max = 5.0;
Dx1 = (x1max-x1min)/(Nx1-1);
x1 = gda_cvec(x1min + Dx1*np.linspace(0,Nx1-1,Nx1));
 
# x2-axis
Nx2 = 101;
x2min = 0.0;
x2max = 5.0;
Dx2 = (x2max-x2min)/(Nx2-1);
x2 = gda_cvec(x2min + Dx2*np.linspace(0,Nx2-1,Nx2));
 
# parameters for a Normal pdf
P1=np.zeros((Nx1,Nx2));
x1bar = 2.25;
x2bar = 2.08;
bar = gda_cvec(x1bar, x2bar);
sx1 = 0.5;
sx2 = 1.0;
C1 = np.diag( gda_cvec(sx1**2, sx2**2).ravel() );
CI1 = la.inv(C1);
DC1 = abs(la.det(C1));
norm1 = (1.0/(2.0*pi)) * (1.0/sqrt(DC1));
 
# tabulate the Normal pdf on the grid
for i in range(Nx1):
    for j in range(Nx2):
        x = gda_cvec(x1[i,0], x2[j,0]) - bar;
        R2 = np.matmul(x.T,np.matmul(CI1,x)); R2=R2[0,0];
        P1[i,j] = norm1*exp( -0.5*R2 );

# axis for parametric curve
Ns = 51;
smin = 0.0;
smax = 5.0;
Ds = (smax-smin)/(Ns-1);
s = gda_cvec(smin + Ds*np.linspace(0,Ns-1,Ns));

# parameric curve
a = 1.0+s+np.sin(2.0*pi*s/smax)-2.0*np.power(s/smax,2);
b = 0.3*np.exp(-0.5*(np.power(s-2.3,2)/(0.1**2)));
c = -0.3*np.exp(-0.5*(np.power(s-2.8,2)/(0.1**2)));          
x2p = a+b+c;
x1p = s;
 
# pdf along parameteric curve
ipx1 = (np.floor((x1p-x1min)/Dx1)).astype(int);
ipx2 = (np.floor((x2p-x2min)/Dx2)).astype(int);
# insure indices in range
kk = np.where(ipx1<0);
i=kk[0];
if( i.any() ):
    ipx1[i,0:1]=0;
kk = np.where(ipx1>=Nx1);
i=kk[0];
if( i.any() ):
    ipx1[i,0:1]=Nx1-1;
         
kk = np.where(ipx2<0);
i = kk[0];
if( i.any() ):
    ipx2[i,0:1]=0;
kk=np.where(ipx2>=Nx2);
i=kk[0];
if( i.any() ):
    ipx2[i,0]=Nx2-1;
         
Ps=np.zeros((Ns,1));
# evaluate P at indices
for i in range(Ns):
    Ps[i,0] = P1[ipx1[i,0], ipx2[i,0] ];
    
    
# calculate mean of pdf p(s)
norm = Ds*np.sum(Ps);
Ps = Ps/norm;
smean = Ds*np.sum(np.multiply(s,Ps));

# maximum along curve
Pmax = np.max(Ps);
ismax =np.argmax(Ps);
x1smax = x1p[ismax,0];
x2smax = x2p[ismax,0];
smaxP = x1smax;

# plot pdf p(x1,x2)
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [x1min, x1max, x2min, x2max] );
left=x1min;
right=x1max;
bottom=x2min;
top=x2max;
dmin=0.0;
dmax=np.max(P1);
plt.imshow( np.flipud(P1), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.gca().invert_yaxis();
plt.xlabel("x2");
plt.ylabel("x1");
plt.plot( x2bar, x1bar, "wo", lw=3 );
plt.plot( gda_cvec(x2bar, x2bar), gda_cvec(x1min, x1min+0.1), "w-", lw=3 );
plt.plot( gda_cvec(x2min, x2min+0.1), gda_cvec(x1bar, x1bar), "w-", lw=3 );
plt.plot( x2p, x1p, "w-", lw=3 );
plt.plot( x2smax, x1smax, "ko", lw=4);
plt.plot( gda_cvec(x2min, x2smax), gda_cvec(x1smax, x1smax), "w:", lw=2);
plt.plot( gda_cvec(x2smax, x2smax), gda_cvec(x1min, x1smax), "w:", lw=2);
plt.show();
print("Caption: Exemplary normal pdf (colors) with mean (white circle),");
print("curve f(x)=0 (white curve), and maximum likelihood point (black");
print("circle), superimposed.");

# plot of 1D pdf (along curve)
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [smin, smax, 0, 1.1*np.max(Ps)] );
plt.plot( s, Ps, "b-", lw=3 );
plt.xlabel("s");
plt.ylabel("p(s)");
plt.plot( gda_cvec(smaxP, smaxP), gda_cvec(0.0, 0.03), "k-", lw=3);
plt.plot( smaxP, Pmax, "bo", lw=3);
plt.plot( gda_cvec(smean, smean), gda_cvec(0, 0.02), "k-", lw=3);
plt.show();
print("Caption: probability alonr the curve p(x)=0, with the maximum");
print("likelihood point shown (circle and long bar) and mean (short bar).");
gdapy11_13
No description has been provided for this image
Caption: Exemplary normal pdf (colors) with mean (white circle),
curve f(x)=0 (white curve), and maximum likelihood point (black
circle), superimposed.
No description has been provided for this image
Caption: probability alonr the curve p(x)=0, with the maximum
likelihood point shown (circle and long bar) and mean (short bar).
In [15]:
# gdapy11_14

# Gradient method applied to inverse problem
# d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2)

print("gdapy11_14");

# Auxially variable x
N=40;
x = gda_cvec(np.linspace(0,N-1,N))/N;
 
# true model parameters
mt = gda_cvec(1.21, 1.54);
 
# true data
w0=20;
dtrue = np.sin(w0*mt[0,0]*x) + mt[0,0]*mt[1,0];
 
# observed data
sd=0.4;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# 2D (m1,m2) grid
M = 101;
dm = 0.02;
m1min=0.0;
m2min=0.0;
m1a = gda_cvec(m1min+dm*np.linspace(0,M-1,M));
m2a = gda_cvec(m2min+dm*np.linspace(0,M-1,M));
m1max = m1a[M-1,0];
m2max = m2a[M-1,0];
 
# tabulate E(m1,m2) on grid (for plotting purposes only)
E = np.zeros((M,M));
for j in range(M):
    for k in range(M):
        yest = np.sin(w0*m1a[j,0]*x) + m1a[j,0]*m2a[k,0];
        myE = np.matmul((dobs-yest).T,(dobs-yest)); myE=myE[0,0];
        E[j,k] = myE;

# parameters for gradient method
alpha = 0.05;
c1 = 0.0001;
c2 = 0.9;
tau = 0.5;
 
# trial solution
mgo=gda_cvec(1,1);

# save history of Error and model parameters
Niter=500;
Ehis=np.zeros((Niter,1));
m1his=np.zeros((Niter,1));
m2his=np.zeros((Niter,1));
 
# error and its gradient at the trial solution
ygo = np.sin( w0*mgo[0,0]*x) + mgo[0,0]*mgo[1,0];
Ego = np.matmul((ygo-dobs).T,(ygo-dobs)); Ego=Ego[0,0];
dydmo = np.zeros((N,2));
dydmo[0:N,0:1] = w0*np.multiply(x,np.cos(w0*mgo[0,0]*x)) + mgo[1,0];
dydmo[0:N,1:2] = mgo[1,0]*np.ones((N,1));
dEdmo = 2.0*np.matmul(dydmo.T,ygo-dobs);

Nhis=0;  # itration counter
Ehis[Nhis,0]=Ego;
m1his[Nhis,0]=mgo[0,0];
m2his[Nhis,0]=mgo[1,0];
Nhis=Nhis+1;
w=np.ones((N,1));
 
for k in range(1,Niter):

    myv = np.matmul(dEdmo.T,dEdmo); myv=myv[0,0]
    v = -dEdmo / sqrt(myv); # downhill direction
    
    # backstep
    for kk in range(10):
        mg = mgo+alpha*v;
        yg = np.sin( w0*mg[0,0]*x) + mg[0,0]*mg[1,0];
        Eg = np.matmul((yg-dobs).T,(yg-dobs)); Eg=Eg[0,0];
        dydm = np.zeros((N,2));
        dydm[0:N,0:1] = w0*np.multiply(x,np.cos(w0*mg[0,0]*x)) + mg[1,0];
        dydm[0:N,1:2] = mg[1,0]*w;
        dEdm = 2.0*np.matmul(dydm.T,yg-dobs);
        if( Eg <= (Ego + c1*alpha*np.matmul(v.T,dEdmo)) ):
            break;
        alpha = tau*alpha;
 
    # change in solution
    Dgmsq=np.matmul((mg-mgo).T,(mg-mgo)); Dgmsq=Dgmsq[0,0];
    Dmg = sqrt(Dgmsq);
    
    # update
    mgo=mg;
    ygo = yg;
    Ego = Eg;
    dydmo = dydm;
    dEdmo = dEdm;
    
    # save history
    Ehis[Nhis,0]=Eg;
    m1his[Nhis,0]=mg[0,0];
    m2his[Nhis,0]=mg[1,0];
    Nhis=Nhis+1;
    
    if( Dmg < 1.0e-6 ):
        break; # terminate iterations when change
               # in solution is suffiently small

# estimated solution is last iteration
m1est = m1his[Nhis-1,0];
m2est = m2his[Nhis-1,0];
 
# truncate to actual lenth
Ehis = Ehis[0:Nhis,0:1];
m1his=m1his[0:Nhis,0:1];
m2his=m2his[0:Nhis,0:1];
 
# evaluate prediction and plot it
dpre = np.sin(w0*m1est*x) + m1est*m2est;

# evaluate prediction and plot it
dpre = np.sin(w0*m1est*x) + m1est*m2est;

# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0.0, xmax, 0.0, 4.0 ] );
plt.plot(x,dtrue,"k-",lw=3);
plt.plot(x,dobs,"ko",lw=3);
plt.xlabel("x");
plt.ylabel("d");
plt.plot(x,dpre,"r-",lw=3);
plt.show();
print("Caption: Exemplary curve fit with Newtons methiod, showing observed");
print("data (circles), true data (black curve) and predicted data (red curve)");

# plot error surface E(m1,m2)
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m2min, m2max, m1min, m1max] );
left=m1min;
right=m1max;
bottom=m2min;
top=m2max;
dmin=np.min(E);
dmax=np.max(E);
plt.imshow( np.flipud(E), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.gca().invert_yaxis();
plt.plot( m2his[0,0], m1his[0,0], "wo", lw=3 );
plt.plot( m2his, m1his, "w-", lw=1 );
plt.plot( mt[1,0], mt[0,0], "go", lw=2 );
plt.plot( m2est, m1est, "ro", lw=2 );
plt.xlabel("m2");
plt.ylabel("m1");
plt.show();
print("Caption: Error surface (colors), with initial guess of solution")
print("(white circle), trajectory (white line), estimated solution (red");
print("circle) and true solution (green).")

# plot history
w = gda_cvec(np.linspace(0,Nhis-1,Nhis));
fig3 = plt.figure(3,figsize=(12,5));

# plot m1 history
ax1=plt.subplot(3,1,1);
plt.axis( [0, Nhis, 0.0, 2.0 ] );
plt.plot(w,m1his,"k-",lw=3);
plt.plot(gda_cvec(0.0,Niter),gda_cvec(mt[0,0],mt[0,0]),'r:',lw=2);
plt.xlabel("iteration i");
plt.ylabel("m1");

# plot m2 history
ax1=plt.subplot(3,1,2);
plt.axis( [0, Nhis, 0.0, 2.0 ] );
plt.plot(w,m2his,"k-",lw=3);
plt.plot(gda_cvec(0.0,Niter),gda_cvec(mt[1,0],mt[1,0]),'r:',lw=2);
plt.xlabel("iteration i");
plt.ylabel("m2");

# plot E history
ax1=plt.subplot(3,1,3);
plt.axis( [0, Nhis, 0.0, 1.1*np.max(Ehis) ] );
plt.plot(w,Ehis,"k-",lw=3);
plt.xlabel("iteration i");
plt.ylabel("E");
plt.show();
print("Caption: (top) Error E(i) as a function of iteration number i.");
print("(middle) Model parameter m1(i) as a function of iteration number i,");
print("with true solution (red line) for comparison. (bottom) Same as");
print("middel, except for model parameter m2");
gdapy11_14
No description has been provided for this image
Caption: Exemplary curve fit with Newtons methiod, showing observed
data (circles), true data (black curve) and predicted data (red curve)
No description has been provided for this image
Caption: Error surface (colors), with initial guess of solution
(white circle), trajectory (white line), estimated solution (red
circle) and true solution (green).
No description has been provided for this image
Caption: (top) Error E(i) as a function of iteration number i.
(middle) Model parameter m1(i) as a function of iteration number i,
with true solution (red line) for comparison. (bottom) Same as
middel, except for model parameter m2
In [16]:
# gdapy11_15
# plots the distributions of one bit-mutations
# for both normal binary and Gray code representations
# of three exemplary 16-bit integers
 
print("gdapy11_15");
 
# setup for uing gray numbers
graydict=scipyio.loadmat("..\Data\gda_gray1table.mat");
gray1=graydict["gray1"];
gray1index=graydict["gray1index"];
graypow2=graydict["pow2"];

fig3 = plt.figure(3,figsize=(12,5));

i1_list = [11111, 33333, 55555];

mysubplot=1;
for i1 in i1_list:
    ax1=plt.subplot(3,2,mysubplot);
    mysubplot=mysubplot+2;
    plt.axis( [0, 65535, 0, 1] );
    plt.plot( gda_cvec(i1,i1), gda_cvec(0.8, 1), "r-", lw=2 );
    gg1 = gda_int2bin(i1);
    g1 = np.empty((16,),dtype=np.str_);
    for i in range(16):
        g1[i]=gg1[i];
    for i in range(16):
        g2 = np.copy(g1);
        k=i;
        if( g1[k] == "0" ):
            g2[k] = "1";
        else:
            g2[k] = "0";
        i2=gda_bin2int(g2);
        plt.plot( gda_cvec(i2,i2), gda_cvec(0.2, 0.8), "k-", lw=2 );

mysubplot=2;
for i1 in i1_list:
    ax1=plt.subplot(3,2,mysubplot);
    mysubplot=mysubplot+2;
    plt.axis( [0, 65535, 0, 1] );
    plt.plot( gda_cvec(i1,i1), gda_cvec(0.8, 1), "r-", lw=2 );
    gg1 = gda_int2gray(i1);
    g1 = np.empty((16,),dtype=np.str_);
    for i in range(16):
        g1[i]=gg1[i];
    for i in range(16):
        g2 = np.copy(g1);
        k=i;
        if( g1[k] == "0" ):
            g2[k] = "1";
        else:
            g2[k] = "0";
        i2=gda_gray2int(g2);
        plt.plot( gda_cvec(i2,i2), gda_cvec(0.2, 0.8), "k-", lw=2 );

plt.show();
print("Caption: Distribution of integers (black bars) resulting from a");
print("single bit mutation of a target inter (red bar). (Left column)");
print("Standard coding, (right colummn) Gray coding.");
gdapy11_15
No description has been provided for this image
Caption: Distribution of integers (black bars) resulting from a
single bit mutation of a target inter (red bar). (Left column)
Standard coding, (right colummn) Gray coding.
In [17]:
# gdapy11_16
# Genetic Algorithm applied to inverse problem
# d(x)=g(x, m1, m2) with g = sin(w0*m(1)*x) + m(1)*m(2)
 
print("gdapy11_16");

c=5; # fitness parameter

# setup for uing gray numbers
graydict=scipyio.loadmat("..\Data\gda_gray1table.mat");
gray1=graydict["gray1"];
gray1index=graydict["gray1index"];
graypow2=graydict["pow2"];

# two model parameters
M=2;
 
# Number of generations
L=100;
Esave = np.zeros((L,1));
Msave = np.zeros((L, M));
 
# x-axis
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = gda_cvec(Dx*np.linspace(0,N-1,N));
 
# true model parameters
mt = gda_cvec(1.21, 1.54);
 
# y=f(x, m1, m2);
w0=20.0;
dtrue = np.sin(w0*mt[0,0]*x) + mt[0,0]*mt[1,0];
sd=0.4;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1));
 
# plot data
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [0, xmax, 0, 4 ] );
plt.plot(x,dtrue,"k-",lw=3);
plt.plot(x,dobs,"ko",lw=3);
plt.xlabel("x");
plt.ylabel("d");
 
# define 2D grid
# this is just for plotting purposes
# the genetic algorith does not require a grid
Lm = 101;
Dm = 0.02;
m1min=0.0;
m2min=0.0;
m1a = gda_cvec(m1min+Dm*np.linspace(0,Lm-1,Lm));
m2a = gda_cvec(m2min+Dm*np.linspace(0,Lm-1,Lm));
m1max = m1a[Lm-1,0];
m2max = m2a[Lm-1,0];
 
# compute error E on the grid
# this is just for plotting purposes
# the genetic algorith does not require a grid
E = np.zeros((Lm,Lm));
for j in range(Lm):
    for k in range(Lm):
        dpre = np.sin(w0*m1a[j,0]*x) + m1a[j,0]*m2a[k,0];
        myE = np.matmul((dobs-dpre).T,(dobs-dpre)); myE=myE[0,0];
        E[j,k] = myE;

K = 100;  # number of individuals in the population
Kold = K;
m0 = gda_cvec(0.25, 0.30); # initial guess for m1 
mL = gda_cvec(m1min, m2min); # solution > mL
mR = gda_cvec(m1max, m2max); # solution < mR
 
dpre = np.sin(w0*m0[0,0]*x) + m0[0,0]*m0[1,0]; # initial prediction
e = dobs-dpre; # individual errors
E0 = np.matmul(e.T,e); # total error at start
 
mind = np.zeros((K,M)); # solution for each individual
Eind = np.zeros((K,1));     # error for each individual
F = np.zeros((K,1)); # fitness of each individual
genes = np.empty(shape=(K,16*M),dtype=np.str_); # genes for each individual 
for i in range(M):  # loop over model parameters
    r = np.random.normal(loc=0.0,scale=0.05,size=(K,1));
    mind[0:K,i:i+1] = m0[i,0]*np.ones((K,1))+r; # all individuals initialized to initial guess
for i in range(K): # loop over individuals
    dpre = np.sin(w0*mind[i,0]*x) + mind[i,0]*mind[i,1];
    e = dobs-dpre; # individual errors
    myE = np.matmul(e.T,e); myE=myE[0,0];
    Eind[i,0]=myE; # error
Emax = np.max(Eind);
# fitness is a function of error
F=np.multiply(np.exp(-c*Eind/Emax),np.random.uniform(low=0,high=1,size=(K,1)));
  
# map solution into genes
for i in range(K): # loop over individuals
    for j in range(M): # loop over model parameters
        # scale model parameter to integer k
        k = int(floor(65535.999*(mind[i,j]-mL[j,0])/(mR[j,0]-mL[j,0])));
        g = gda_int2gray(k);
        for l in range(16):
            genes[i,j*16+l] = g[l];

# save initial solution
mindsave = mind;
 
# evolve with time
for gen in range(L):

    # each individual replicates
    mind = np.concatenate((mind,mind),axis=0);
    genes = np.concatenate((genes,genes),axis=0);
    Eind = np.concatenate((Eind,Eind),axis=0);
    K = 2*K;
 
    # mutate genes
    for i in range(K): # loop over individuals
        if( np.random.randint(low=0,high=2) == 0 ): # mutate only half of individuals
            continue;
        k = np.random.randint(low=0,high=16*M); # one random mutation
        if( genes[i,k]=="0" ):
            genes[i,k]="1";
        else:
            genes[i,k]="0";
 
    # conjugate genes; only one conjugation per generation
    j1 = np.random.randint(low=0,high=K); # random individual
    j2 = np.random.randint(low=0,high=K); # another random individual
    k = np.random.randint(low=2,high=M*16-2); # split location 2<=k<=16M-1
    # be sure to  copy!
    g1c = np.copy(genes[j1,0:M*16]);
    g2c = np.copy(genes[j2,0:M*16]);
    g1c[k:M*16] = genes[j2,k:M*16];
    g2c[k:M*16] = genes[j1,k:M*16];
    genes[j1,0:16*M]=g1c;
    genes[j2,0:16*M]=g2c;
    # update solution and error from genes
    for i in range(K): # loop over individuals
        for j in range(M): # loop over mode paramters
            g = genes[i,j*16:(j+1)*16]
            k = gda_gray2int(g);
            mind[i,j] = (mR[j,0]-mL[j,0])*(k/65535.999)+mL[j,0];  
    for i in range(K): # loop over individuals
        dpre = np.sin(w0*mind[i,0]*x) + mind[i,0]*mind[i,1];
        e = dobs-dpre;
        myE=np.matmul(e.T,e); myE=myE[0,0];
        Eind[i,0] = myE;
    Emax = np.max(Eind);
    # fitness           
    F=np.multiply(np.exp(-c*Eind/Emax),np.random.uniform(low=0,high=1,size=(K,1)));
 
    # survival of the fitness
    isort = np.flipud(np.argsort(F,axis=0)).ravel();
    isort = isort[0:Kold];
    mind=mind[isort,0:M];        
    Eind = Eind[isort,0:1];
    genes = genes[isort,0:16*M];
    K=Kold;
 
    Emin = np.min(Eind);
    k = np.argmin(Eind);
    Esave[gen,0] = Emin;
    Msave[gen,0:M]= mind[k,0:M];

# solution is lowest error individual in last generation
m1est = Msave[L-1,0];
m2est = Msave[L-1,1];

print("True solution:      %.2f %.2f" % (mt[0,0],mt[1,0]) );
print("Initial solution:   %.2f %.2f" % (m0[0,0],m0[1,0]) );
print("Estimated solution: %.2f %.2f" % (m1est,m2est) );
print("  ");

# evaluate prediction and plot it
dpre = np.sin(w0*m1est*x) + m1est*m2est;
plt.plot(x,dpre,"r-",lw=3);
plt.show();
print("Caption: True data (black curve), observed data (black circles)");
print("predicted data (red curve).");

# plot error surface E(m1,m2)
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m2min, m2max, m1min, m1max] );
left=m1min;
right=m1max;
bottom=m2min;
top=m2max;
dmin=np.min(E);
dmax=np.max(E);
plt.imshow( np.flipud(E), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.gca().invert_yaxis();
plt.plot( mindsave[0:K,1], mindsave[0:K,0], "k.", lw=1 );
plt.plot( mind[0:K,1], mind[0:K,0], "w.", lw=1 );
plt.plot( mt[1,0], mt[0,0], "go", lw=2 );
plt.plot( m2est, m1est, "ro", lw=2 );
plt.xlabel("m2");
plt.ylabel("m1");
plt.show();
print("Caption: Error surface (colors), with initial population (black dots)");
print("final population (white dots), estimated solution (red circle) and");
print("true solution (green circle)");

# plot history
w = gda_cvec(np.linspace(0,L-1,L));
fig3 = plt.figure(3,figsize=(12,5));

# plot m1 history
ax1=plt.subplot(3,1,1);
plt.axis( [0, L, 0.0, 2.0 ] );
plt.plot(w,Msave[0:L,0],"k-",lw=3);
plt.plot(gda_cvec(0.0,L),gda_cvec(mt[0,0],mt[0,0]),'r:',lw=2);
plt.xlabel("generation i");
plt.ylabel("m1");

# plot m2 history
ax1=plt.subplot(3,1,2);
plt.axis( [0, L, 0.0, 2.0 ] );
plt.plot(w,Msave[0:L,1],"k-",lw=3);
plt.plot(gda_cvec(0.0,L),gda_cvec(mt[1,0],mt[1,0]),'r:',lw=2);
plt.xlabel("generation i");
plt.ylabel("m2");

# plot E history
ax1=plt.subplot(3,1,3);
plt.axis( [0, L, 0.0, 1.1*np.max(Esave) ] );
plt.plot(w,Esave,"k-",lw=3);
plt.xlabel("generation i");
plt.ylabel("E");
plt.show();
print("Caption: (top) Error as a function of generation number (black curve), with");
print("inital error (red line). (middle) Model parametre m1 as a function of");
print("generation number (black curve), with true value shown form comparison");
print("(red line). (bottom) Same as middle except for model parameter m2.");
gdapy11_16
True solution:      1.21 1.54
Initial solution:   0.25 0.30
Estimated solution: 1.20 1.61
  
No description has been provided for this image
Caption: True data (black curve), observed data (black circles)
predicted data (red curve).
No description has been provided for this image
Caption: Error surface (colors), with initial population (black dots)
final population (white dots), estimated solution (red circle) and
true solution (green circle)
No description has been provided for this image
Caption: (top) Error as a function of generation number (black curve), with
inital error (red line). (middle) Model parametre m1 as a function of
generation number (black curve), with true value shown form comparison
(red line). (bottom) Same as middle except for model parameter m2.
In [18]:
# gdapy11_17
# Genetic algorithm applied to an "approximately
# separable inverse problem" where the M=20 model parameters
# have a natural ordering and the data are localized
# averages of the model parameters.
 
# This code compares error of two cases, A and B, 
# where the probability of mutation and recombination
# differ between the cases.  Each case has 5 trials.

print("gdapy11_17");

# setup for uing gray numbers
graydict=scipyio.loadmat("..\Data\gda_gray1table.mat");
gray1=graydict["gray1"];
gray1index=graydict["gray1index"];
graypow2=graydict["pow2"];

print("This calculation is pretty slow, so progress shown (ends at B4)");

# model parameters
M=20;

# Number of generations
L=50;
Esave = np.zeros((L,1));
Msave = np.zeros((L,M));
 
# plot error
fig1 = plt.figure(1,figsize=(12,5));
ax1=plt.subplot(1,1,1);
plt.axis( [1.0, L, 0.0, 1.5 ] );
plt.xlabel("generation");
plt.ylabel("E");
 
Ntrials=5;
 
# part 1 --------------------------------------------
Nrecom = 0;
imutate = 8;
 
for itrial in range(Ntrials):

    # fitness coefficient
    c=5;
 
    # true model parameters are linear ramp
    mt = gda_cvec(0.5+np.linspace(0,M-1,M)/(M-1));
 
    # data average of 3 neighboring model parameters
    N=M;
    x = gda_cvec(np.linspace(1,M,M));
    dtrue = (mt
             +gda_cvec(0.0,mt[0:M-1,0:1])
             +gda_cvec(0.0,0.0,mt[0:M-2,0:1]))/3.0;
    sd=0.025;
    dobs = dtrue + np.random.normal(loc=0,scale=sd,size=(N,1));
 
    K = 100;  # number of individuals in the population
    Kold = K;
    m0 = np.ones((M,1));     # initial guess for  solution
    mL = np.zeros((M,1));    # solution > mL
    mR = 2.0*np.ones((M,1)); # solution < mR
    d0 = (m0
             +gda_cvec(0.0,m0[0:M-1,0:1])
             +gda_cvec(0.0,0.0,m0[0:M-2,0:1]))/3.0;

    mind = np.zeros((K,M)); # solution for each individual
    Eind = np.zeros((K,1));    # error for each individual
    Ftind = np.zeros((K,1));   # fitness of each individual
    genes = np.empty(shape=(K,16*M),dtype=np.str_); # genes for each individual 

    # all individuals initialized to initial guess + random number
    for i in range(M):
        r = np.random.normal(loc=0.0,scale=0.025,size=(K,1));
        mind[0:K,i:i+1] = m0[i,0]*np.ones((K,1))+r; 

    # compute initial error
    for i in range(K):
        dpre = ( gda_cvec(mind[i:i+1,0:M])
            +gda_cvec(0.0,mind[i:i+1,0:M-1])
            +gda_cvec(0.0,0.0,mind[i:i+1,0:M-2])  )/3.0;
        e = dobs-dpre;
        myE = np.matmul(e.T,e); myE=myE[0,0]
        Eind[i,0] = myE;
    Emax = np.max(Eind);
    Emaxstart = Emax;
    # initial fitness
    Ftind = np.multiply( np.exp(-c*Eind/Emax),
            np.random.uniform(low=0.0,high=1.0,size=(K,1)) ); 

    # map solution into genes
    for i in range(K): # loop over individuals
        for j in range(M): # loop over model parameters
            # scale model parameter to integer k
            k = int(floor(65535.999*(mind[i,j]-mL[j,0])/(mR[j,0]-mL[j,0])));
            g = gda_int2gray(k);
            for l in range(16):
                genes[i,j*16+l] = g[l];

    for gen in range(L): # loop over generations
 
        # each individual replicates
        mind = np.concatenate((mind,mind),axis=0);
        genes = np.concatenate((genes,genes),axis=0);
        Eind = np.concatenate((Eind,Eind),axis=0);
        K = 2*K;
        
        # mutate genes
        for i in range(K): # loop over individuals
            if( np.random.randint(low=0,high=imutate) == 0 ): # mutate only some
                continue;
            k = np.random.randint(low=0,high=16*M); # one random mutation
            if( genes[i,k]=="0" ):
                genes[i,k]="1";
            else:
                genes[i,k]="0";
             
        # conjugate genes
        for rc in range(Nrecom):
            j1 = np.random.randint(low=0,high=K); # random individual
            j2 = np.random.randint(low=0,high=K); # another random individual
            k = np.random.randint(low=2,high=M*16-2); # split location 2<=k<=16M-1
            # be sure to  copy!
            g1c = np.copy(genes[j1,0:M*16]);
            g2c = np.copy(genes[j2,0:M*16]);
            g1c[k:M*16] = genes[j2,k:M*16];
            g2c[k:M*16] = genes[j1,k:M*16];
            genes[j1,0:16*M]=g1c;
            genes[j2,0:16*M]=g2c;
            
        # update solution and error from genes
        for i in range(K): # loop over individuals
            for j in range(M): # loop over mode paramters
                g = genes[i,j*16:(j+1)*16]
                k = gda_gray2int(g);
                mind[i,j] = (mR[j,0]-mL[j,0])*(k/65535.999)+mL[j,0];  
        for i in range(K): # loop over individuals
            dpre = ( gda_cvec(mind[i:i+1,0:M])
                +gda_cvec(0.0,mind[i:i+1,0:M-1])
                +gda_cvec(0.0,0.0,mind[i:i+1,0:M-2])  )/3.0;
            e = dobs-dpre;
            myE = np.matmul(e.T,e); myE=myE[0,0];
            Eind[i,0] = myE;
        Emax = np.max(Eind);
        # fitness           
        F=np.multiply(np.exp(-c*Eind/Emax),np.random.uniform(low=0,high=1,size=(K,1)));

        # survival of the fitness
        isort = np.flipud(np.argsort(F,axis=0)).ravel();
        isort = isort[0:Kold];
        mind=mind[isort,0:M];        
        Eind = Eind[isort,0:1];
        genes = genes[isort,0:16*M];
        K=Kold;
        
        Emin = np.min(Eind);
        k = np.argmin(Eind);
        Esave[gen,0] = Emin;
        Msave[gen,0:M]= mind[k,0:M];

    # solution is lowest error individual in last generation
    mest = gda_cvec(Msave[L-1,0:M]);
    dpre = ( gda_cvec(mest)
            +gda_cvec(0.0,mest[0:M-1,0])
            +gda_cvec(0.0,0.0,mest[0:M-2,0])  )/3.0;
    w = gda_cvec(np.linspace(1,L,L));
    plt.plot(w,Esave,"k-",lw=2);
    print("Done with A %d" % (itrial) );
    
# part 2 --------------------------------------------
Nrecom = 8;
imutate = 8;
 
for itrial in range(Ntrials):

    # fitness coefficient
    c=5.0;
 
    # true model parameters are linear ramp
    mt = gda_cvec(0.5+np.linspace(0,M-1,M)/(M-1));
 
    # data average of 3 neighboring model parameters
    N=M;
    x = gda_cvec(np.linspace(1,M,M));
    dtrue = (mt
             +gda_cvec(0.0,mt[0:M-1,0:1])
             +gda_cvec(0.0,0.0,mt[0:M-2,0:1]))/3.0;
    sd=0.025;
    dobs = dtrue + np.random.normal(loc=0,scale=sd,size=(N,1));
 
    K = 100;  # number of individuals in the population
    Kold = K;
    m0 = np.ones((M,1));     # initial guess for  solution
    mL = np.zeros((M,1));    # solution > mL
    mR = 2.0*np.ones((M,1)); # solution < mR
    d0 = (m0
             +gda_cvec(0.0,m0[0:M-1,0:1])
             +gda_cvec(0.0,0.0,m0[0:M-2,0:1]))/3.0;

    mind = np.zeros((K,M)); # solution for each individual
    Eind = np.zeros((K,1));    # error for each individual
    Ftind = np.zeros((K,1));   # fitness of each individual
    genes = np.empty(shape=(K,16*M),dtype=np.str_); # genes for each individual 

    # all individuals initialized to initial guess + random number
    for i in range(M):
        r = np.random.normal(loc=0.0,scale=0.025,size=(K,1));
        mind[0:K,i:i+1] = m0[i,0]*np.ones((K,1))+r; 

    # compute initial error
    for i in range(K):
        dpre = ( gda_cvec(mind[i:i+1,0:M])
            +gda_cvec(0.0,mind[i:i+1,0:M-1])
            +gda_cvec(0.0,0.0,mind[i:i+1,0:M-2])  )/3.0;
        e = dobs-dpre;
        myE =np.matmul(e.T,e); myE=myE[0,0];
        Eind[i,0] = myE;
    Emax = np.max(Eind);
    Emaxstart = Emax;
    # initial fitness
    Ftind = np.multiply( np.exp(-c*Eind/Emax),
            np.random.uniform(low=0.0,high=1.0,size=(K,1)) ); 

    # map solution into genes
    for i in range(K): # loop over individuals
        for j in range(M): # loop over model parameters
            # scale model parameter to integer k
            k = int(floor(65535.999*(mind[i,j]-mL[j,0])/(mR[j,0]-mL[j,0])));
            g = gda_int2gray(k);
            for l in range(16):
                genes[i,j*16+l] = g[l];

    for gen in range(L): # loop over generations
 
        # each individual replicates
        mind = np.concatenate((mind,mind),axis=0);
        genes = np.concatenate((genes,genes),axis=0);
        Eind = np.concatenate((Eind,Eind),axis=0);
        K = 2*K;
        
        # mutate genes
        for i in range(K): # loop over individuals
            if( np.random.randint(low=0,high=imutate) == 0 ): # mutate only some
                continue;
            k = np.random.randint(low=0,high=16*M); # one random mutation
            if( genes[i,k]=="0" ):
                genes[i,k]="1";
            else:
                genes[i,k]="0";
             
        # conjugate genes
        for rc in range(Nrecom):
            j1 = np.random.randint(low=0,high=K); # random individual
            j2 = np.random.randint(low=0,high=K); # another random individual
            k = np.random.randint(low=2,high=M*16-2); # split location 2<=k<=16M-1
            # be sure to  copy!
            g1c = np.copy(genes[j1,0:M*16]);
            g2c = np.copy(genes[j2,0:M*16]);
            g1c[k:M*16] = genes[j2,k:M*16];
            g2c[k:M*16] = genes[j1,k:M*16];
            genes[j1,0:16*M]=g1c;
            genes[j2,0:16*M]=g2c;
            
        # update solution and error from genes
        for i in range(K): # loop over individuals
            for j in range(M): # loop over mode paramters
                g = genes[i,j*16:(j+1)*16]
                k = gda_gray2int(g);
                mind[i,j] = (mR[j,0]-mL[j,0])*(k/65535.999)+mL[j,0];  
        for i in range(K): # loop over individuals
            dpre = ( gda_cvec(mind[i:i+1,0:M])
                +gda_cvec(0.0,mind[i:i+1,0:M-1])
                +gda_cvec(0.0,0.0,mind[i:i+1,0:M-2])  )/3.0;
            e = dobs-dpre;
            myE = np.matmul(e.T,e); myE=myE[0,0];
            Eind[i,0] = myE;
        Emax = np.max(Eind);
        # fitness           
        F=np.multiply(np.exp(-c*Eind/Emax),np.random.uniform(low=0,high=1,size=(K,1)));

        # survival of the fitness
        isort = np.flipud(np.argsort(F,axis=0)).ravel();
        isort = isort[0:Kold];
        mind=mind[isort,0:M];        
        Eind = Eind[isort,0:1];
        genes = genes[isort,0:16*M];
        K=Kold;
        
        Emin = np.min(Eind);
        k = np.argmin(Eind);
        Esave[gen,0] = Emin;
        Msave[gen,0:M]= mind[k,0:M];

    # solution is lowest error individual in last generation
    mest = gda_cvec(Msave[L-1,0:M]);
    dpre = ( gda_cvec(mest)
            +gda_cvec(0.0,mest[0:M-1,0])
            +gda_cvec(0.0,0.0,mest[0:M-2,0])  )/3.0;
    w = gda_cvec(np.linspace(1,L,L));
    plt.plot(w,Esave,"r-",lw=2);
    print("Done with B %d" % (itrial) );

plt.show();
print("Caption: Error as a function of generation number for five trials");
print("that omit conjugation (black curves) and five that include it (red curves).");
gdapy11_17
This calculation is pretty slow, so progress shown (ends at B4)
Done with A 0
Done with A 1
Done with A 2
Done with A 3
Done with A 4
Done with B 0
Done with B 1
Done with B 2
Done with B 3
Done with B 4
No description has been provided for this image
Caption: Error as a function of generation number for five trials
that omit conjugation (black curves) and five that include it (red curves).
In [19]:
# gdapy11_18
# Bootstrap estimate of confidence interval
# for solution computed via Newton's method
# This function with these derivatives is being solved
# d(x) = sin( w0 m1  x) + m1 m2;
# dy/dm1 = w0 x cos( w0 m1 x) + m2
# dy/dm2 = m2
 
print("gdama11_18");

# x-axis
N=40;
xmin=0;
xmax=1.0;
Dx=(xmax-xmin)/(N-1);
x = gda_cvec(Dx*np.linspace(0,N-1,N));
 
# true model parameters
mt = gda_cvec(1.21, 1.54);
 
# y=f(x, m1, m2);
w0=20.0;
dtrue = np.sin(w0*mt[0,0]*x) + mt[0,0]*mt[1,0];
sd=0.4;
dobs = dtrue + np.random.normal(loc=0.0,scale=sd,size=(N,1));

# 2D (m1,m2) grid
M = 101;
dm = 0.02;
m1min=0.0;
m2min=0.0;
m1a = gda_cvec(m1min+dm*np.linspace(0,M-1,M));
m2a = gda_cvec(m2min+dm*np.linspace(0,M-1,M));
m1max = m1a[M-1,0];
m2max = m2a[M-1,0];

# compute error E on the grid
# this is just for plotting purposes
E = np.zeros((Lm,Lm));
for j in range(Lm):
    for k in range(Lm):
        dpre = np.sin(w0*m1a[j,0]*x) + m1a[j,0]*m2a[k,0];
        myE = np.matmul((dobs-dpre).T,(dobs-dpre));  myE=myE[0,0];
        E[j,k] = myE;

# compute error E on the grid
# this is just for plotting purposes
# the genetic algorith does not require a grid
E = np.zeros((Lm,Lm));
for j in range(Lm):
    for k in range(Lm):
        dpre = np.sin(w0*m1a[j,0]*x) + m1a[j,0]*m2a[k,0];
        myE=np.matmul((dobs-dpre).T,(dobs-dpre)); myE=myE[0,0];
        E[j,k] = myE;
        
# plot error surface E(m1,m2)
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(1,1,1);
mycmap=matplotlib.colormaps['jet'];
plt.axis( [m2min, m2max, m1min, m1max] );
left=m1min;
right=m1max;
bottom=m2min;
top=m2max;
dmin=np.min(E);
dmax=np.max(E);
plt.imshow( np.flipud(E), cmap=mycmap, vmin=dmin, vmax=dmax, extent=(left,right,bottom,top) );
plt.colorbar();
plt.gca().invert_yaxis();
plt.xlabel("m2");
plt.ylabel("m1");
        
Nresamples = 1000;  # number of resampled problems to solve
m1save = np.zeros((Nresamples,1));
m2save = np.zeros((Nresamples,1));
for ir in range(Nresamples):  # loop over resampled problems
 
    # resampling with duplications of data
    # (first estimate is without resampling)
    if( ir==0 ):
        dresampled = dobs;
        xresampled = x;
    else:
        # random resampling with duplications
        rowindex = np.random.randint(low=0,high=N,size=(N,));
        xresampled = x[rowindex,0:1];
        dresampled = dobs[rowindex,0:1];
        
    # Newton's method
    # initial guess and corresponding error
    mg = gda_cvec(1.3, 1.5);
    dg = np.sin(w0*mg[0,0]*xresampled) + mg[0,0]*mg[1,0];
    Eg = np.matmul((dobs-dg).T,(dobs-dg)); Eg=Eg[0,0];
 
    # iterate to improve initial guess
    Niter=100;
    for k in range(Niter):

        # predicted data
        dg = np.sin( w0*mg[0,0]*xresampled) + mg[0,0]*mg[1,0];
        dd = dresampled-dg;
        Eg=np.matmul(dd.T,dd); Eg=Eg[0,0];
        Ehis[k,0]=Eg;
    
        # linearized data kernel
        G = np.zeros((N,2));
        G[0:N,0:1] = w0*np.multiply(xresampled,np.cos(w0*mg[0,0]*xresampled)) + mg[1,0];
        G[0:N,1:2] = mg[1,0]*np.ones((N,1));
    
        # least squares solution
        GTG = np.matmul(G.T,G);
        GTd = np.matmul(G.T,dd);
        dm = la.solve(GTG,GTd);
    
        # update
        mg = mg+dm;
    
    # save resampled solutions
    m1save[ir,0] = mg[0,0];
    m2save[ir,0] = mg[1,0];
    
    if( ir==0 ):
        # plot estimated solution
        m1est = mg[0,0];
        m2est = mg[1,0];
        plt.plot( m2est, m1est, "ro", lw=2 );
    else:
        # plot resampled solutions
        plt.plot( mg[1,0], mg[0,0], "w.", lw=1 );
plt.plot( mt[1,0], mt[0,0], "go", lw=3 );

# mean solution
m1mean = np.mean(m1save);
m2mean = np.mean(m2save);

plt.plot( m2est, m1est, "ro", lw=2 );
plt.show();
print("Caption: Error surface (colors), with bootstrap solutions (white dots),");
print("estimated solution (red circle) and true solution (green circle), superimposed.");

# histogram of m1
Lb1=50;
m1hmin=np.min(m1save);
m1hmax=np.max(m1save);
h1, e1 = np.histogram(m1save,Lb1,(m1hmin,m1hmax));  # create histogram
Nb1 = len(h1);  # lengths of histogram
Ne1 = len(e1);  # length of edges
edges1 =  gda_cvec( e1[0:Ne1-1] );   # vector of edges
bins1 = edges1 + 0.5*(edges1[1,0]-edges1[0,0]); # centers
Db1 = bins1[1,0]-bins1[0,0];
pm1 = gda_cvec(h1)/(Db1*np.sum(h1)); # empirical pdf p(m1)
Pm1 = Db1*np.cumsum(pm1,axis=0); # cumulative distribution
pm1max=np.max(pm1);
# 95 percent confidence limits
z = np.where( Pm1>0.025 );
k = z[0];
il = k[0];
qL1 = bins1[il,0];
z = np.where( Pm1>0.975);
k = z[0];
ir = k[0];
qR1 = bins1[ir,0];

# histogram of m2
Lb2=50;
m2hmin=np.min(m2save);
m2hmax=np.max(m2save);
h2, e2 = np.histogram(m2save,Lb2,(m2hmin,m2hmax));  # create histogram
Nb2 = len(h2);  # lengths of histogram
Ne2 = len(e2);  # length of edges
edges2 =  gda_cvec( e2[0:Ne2-1] );   # vector of edges
bins2 = edges2 + 0.5*(edges2[1,0]-edges2[0,0]); # centers
Db2 = bins2[1,0]-bins2[0,0];
pm2 = gda_cvec(h2)/(Db2*np.sum(h2)); # empirical pdf p(m1)
Pm2 = Db2*np.cumsum(pm2,axis=0); # cumulative distribution
pm2max=np.max(pm2);
# 95 percent confidence limits
z = np.where( Pm2>0.025 );
k = z[0];
il = k[0];
qL2 = bins2[il,0];
z = np.where( Pm2>0.975);
k = z[0];
ir = k[0];
qR2 = bins2[ir,0];

# plot error surface E(m1,m2)
fig2 = plt.figure(2,figsize=(12,5));
ax1=plt.subplot(2,1,1);
plt.axis( [m1hmin, m1hmax, 0, pm1max] );
plt.xlabel('m1');
plt.ylabel('p(m1)');
plt.plot( bins1, pm1, 'k-', 'LineWidth', 3 );
plt.plot( gda_cvec(m1mean,m1mean), gda_cvec(0.0,pm1max/10.0), "r-", lw=3 );
plt.plot( gda_cvec(qL1,qL1), gda_cvec(0.0,pm1max/20.0), "r-", lw=3 );
plt.plot( gda_cvec(qR1,qR1), gda_cvec(0.0,pm1max/20.0), "r-", lw=3 );

ax2=plt.subplot(2,1,2);
plt.axis( [m2hmin, m2hmax, 0, pm2max] );
plt.xlabel('m2');
plt.ylabel('p(m2)');
plt.plot( bins2, pm2, 'k-', 'LineWidth', 3 );
plt.plot( gda_cvec(m2mean,m2mean), gda_cvec(0.0,pm2max/10.0), "r-", lw=3 );
plt.plot( gda_cvec(qL2,qL2), gda_cvec(0.0,pm2max/20.0), "r-", lw=3 );
plt.plot( gda_cvec(qR2,qR2), gda_cvec(0.0,pm2max/20.0), "r-", lw=3 );
plt.show();
print("Caption: (Top) Empirical pdfp(m1) for model parameter m1. (black curv) with");
print("mean (large red bar) and 95 percent confidence intervals (small red bars).");
print(" ");

print("estimated m1: %f < %f < %f (95 percent confidence)" % (qL1, m1est, qR1) );
print("estimated m2: %f < %f < %f (95 percent confidence)" % (qL2, m2est, qR2) );
gdama11_18
No description has been provided for this image
Caption: Error surface (colors), with bootstrap solutions (white dots),
estimated solution (red circle) and true solution (green circle), superimposed.
No description has been provided for this image
Caption: (Top) Empirical pdfp(m1) for model parameter m1. (black curv) with
mean (large red bar) and 95 percent confidence intervals (small red bars).
 
estimated m1: 1.208547 < 1.225435 < 1.239858 (95 percent confidence)
estimated m2: 1.476239 < 1.569842 < 1.675842 (95 percent confidence)
In [ ]: