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
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
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
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
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
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
Caption: Error surface E(m1,m2)
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
Caption: Error surface E(m1,m2)
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
Caption: Exemplary curve fit wia grid search, showing true data (black curve) observed data (circles) and predicted data (red curve).
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
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
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
Caption: Exemplary curve fit with Newtons methiod, showing observed data (circles), true data (black curve) and predicted data (red curve)
Caption: Error surface (colors), with initial guess of solution (white circle), trajectory (white line), estimated solution (red circle) and true solution (green).
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
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
Caption: Exemplary normal pdf (colors) with mean (white circle), curve f(x)=0 (white curve), and maximum likelihood point (black circle), superimposed.
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
Caption: Exemplary normal pdf (colors) with mean (white circle), curve f(x)=0 (white curve), and maximum likelihood point (black circle), superimposed.
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
Caption: Exemplary curve fit with Newtons methiod, showing observed data (circles), true data (black curve) and predicted data (red curve)
Caption: Error surface (colors), with initial guess of solution (white circle), trajectory (white line), estimated solution (red circle) and true solution (green).
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
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
Caption: True data (black curve), observed data (black circles) predicted data (red curve).
Caption: Error surface (colors), with initial population (black dots) final population (white dots), estimated solution (red circle) and true solution (green circle)
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
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
Caption: Error surface (colors), with bootstrap solutions (white dots), estimated solution (red circle) and true solution (green circle), superimposed.
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 [ ]: