In [1]:
# edapy_05_00 clear all variables and import vatious modules
# this cell must be exectuted first, before any of the others
# History
# 2022/06/26 - checked with most recent Python & etc. software
# 2023/07/13 - fixed some issues, incl colormap, loglof plt, eda_cvec()
# 2024/05/26 - Patches to accomodate new verions of numpy, scipy, matplotlib
# patched dot product to return scalar value
# las.bicg keyword tol changes to rtol
# las.bicg output cast to float (not sure this is realy necessary)
# switched index arrays to type int from np.int32
%reset -f
import os
from datetime import date
from math import exp, pi, sin, sqrt, floor, ceil
import numpy as np
import scipy.linalg as la
import scipy.sparse.linalg as las
from scipy import sparse
import matplotlib
from matplotlib import pyplot as plt
from matplotlib.colors import ListedColormap
# eda_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 eda_draw(*argv):
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;
# function to perform the multiplication FT (F v)
def FTFmul(v):
# this function is used by the bicongugate gradient solver to solve the geneneralized least
# squares problem Fm=f. Note that "F" must be called "edaFsparse".
global edaFsparse;
s = np.shape(v);
if(len(s)==1):
vv = np.zeros((s[0],1));
vv[:,0] = v;
else:
vv=v;
temp = edaFsparse*vv;
return edaFsparse.transpose()*temp;
# 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 eda_cvec(*argv):
t = int;
Nt = 0;
for a in argv:
v = eda_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 = eda_cvec1(a);
N,M = np.shape(v);
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=eda_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=eda_cvec(v1,v2,...) concatenates
# many argiments.
def eda_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("eda_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(v, np.int32):
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("eda_cvec: list contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = v;
return w;
elif isinstance(v, tuple):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(v, np.int32):
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("eda_cvec: tuple contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = v;
return w;
else:
raise TypeError("eda_cvec: %s not supported" % type(v));
In [3]:
# edapy_05_01: application of Bayes theorem to least squares
# set up vectors m1 and m2
L = 40;
Dm = 1.0;
m1 = eda_cvec( np.linspace(0,Dm*(L-1),L) );
m2 = eda_cvec( np.linspace(0,Dm*(L-1),L) );
# make normal distribution pp(m)
m1bar=10;
m2bar=20;
Cpm=[[6,0],[0,10]];
CpmI=la.inv(Cpm);
norm=1/(2*pi*sqrt(la.det(Cpm)));
Pp=np.zeros((L,L));
dm = np.zeros((2,1));
for i in range(L):
for j in range(L):
dm = eda_cvec( [m1[i,0]-m1bar, m2[j,0]-m2bar] );
Pij = norm*np.exp( -0.5 * np.matmul( np.matmul( dm.T, CpmI ), dm )); Pij=Pij[0,0];
Pp[i,j]=Pij;
# search grid for maximum
k = np.argmax( Pp, axis=0 );
Pptemp = np.max( Pp, axis=0 );
j = np.argmax(Pptemp);
i=k[j];
m1bar=m1[i];
m2bar=m2[j];
# make normal distribution pdgm(m)
# with a one-row data kernel that the
# difference of two model parameters is d1
G = np.zeros((1,2));
G[0,:] = [1, -1];
d1 = 0;
Cd = 3;
CdI = 1/Cd;
norm = 1/sqrt(2*pi*Cd);
Pdgm = np.zeros((L,L));
m = np.zeros((2,1));
for i in range(L):
for j in range(L):
m[:,0] = [m1[i,0], m2[j,0]];
dd = d1 - np.matmul( G, m );
Pij = norm*np.exp( -0.5 * CdI * np.matmul( dd.T,dd ) ); Pij=Pij[0,0];
Pdgm[i,j]=Pij;
# multiply the two to get Pmgd
Pmgd = np.multiply( Pp , Pdgm);
# search grid for maximum
k = np.argmax( Pmgd, axis=0 ); # row indices of Pmgd with the largest value in each col
Pmgdtemp = np.max( Pmgd, axis=0 ); # maximums in each column
j = np.argmax(Pmgdtemp); # indice to the column with the overall maximum
i=k[j]; # row of the overall maximum
m1est=m1[i]; # m1 position on the grid of the maximum
m2est=m2[j]; # m2 position on the grid of the maximum
print('m1: prior', m1bar, 'estimate', m1est);
print('m2: prior', m2bar, 'estimate', m2est);
# plot result
eda_draw('title pp(d)',Pp,'title p(d|m)',Pdgm,'title p(m|d)',Pmgd);
m1: prior [10.] estimate [13.] m2: prior [20.] estimate [15.]
In [4]:
# edapy_05_02: illustrate the product of twos Normal p.d.f.s
# set up vectors m1 and m2
L = 40;
Dd = 1.0;
d1 = eda_cvec( np.linspace(0,Dd*(L-1),L) );
d2 = eda_cvec( np.linspace(0,Dd*(L-1),L) );
# first Normal p.d.f.
d1bar=15; # d1 mean
d2bar=15; # d2 mean
C=[[10,8],[8,10]]; # covariance
CI=la.inv(C); # its inverse
norm=1/(2*pi*sqrt(la.det(C))); # normalization factor for 2D Normal pdf
pa=np.zeros((L,L)); # matrix to hold p.d.f.
dd = np.zeros((2,1)); # vector to hold difference between value and mean
for i in range(L):
for j in range(L):
dd[:,0] = [d1[i,0]-d1bar, d2[j,0]-d2bar];
pij = norm*np.exp( -0.5 * np.matmul( np.matmul( dd.T, CI ), dd )); pij=pij[0,0];
pa[i,j]=pij;
# second Normal p.d.f.
d1bar=15;
d2bar=25;
C=[[10,-8],[-8,10]];
CI=la.inv(C);
norm=1/(2*pi*sqrt(la.det(C)));
pb=np.zeros((L,L));
dd = np.zeros((2,1));
for i in range(L):
for j in range(L):
dd[:,0] = [d1[i,0]-d1bar, d2[j,0]-d2bar];
pij = norm*np.exp( -0.5 * np.matmul( np.matmul( dd.T, CI ), dd )); pij=pij[0,0];
pb[i,j]=pij;
# product pc = pa * pb
pc=np.zeros((L,L));
pc = np.multiply(pa, pb);
# normalize to unit area
norm = (Dd**2)*np.sum(pc); # sum is approximation for integral
pc = pc/norm;
# plot results
eda_draw('title pA',pa,'title pB',pb,'title pC', pc);
In [5]:
# edapy_05_03: straight-line fit with unequal variance
# illustrate least squares fitting of a line
# when the data have unequal variance
# this version uses G and d
# set of vector of x's
N = 51;
N2 = floor((N-1)/2);
Dx = 2;
xmin = 0.0;
xmax = xmin+Dx*(N-1);
x = eda_cvec( np.linspace(xmin, xmax, N) );
# prepare some test data hat consists of
# two sections, that differ slightly in
# interrcept and slope
M=2;
a = 1.0;
b = 2.0;
dobs = a + b*x;
Da = 1;
Db = 0.3;
dobs[N2:N,0:1] = dobs[N2:N,0:1] + Da + Db*x[N2:N,0:1];
# Part 1: variance, high for thr first section low for the second section
sd = np.zeros((N,1));
sd[0:N2,0:1] = 20.0;
sd[N2:N,0:1] = 2.0;
sd2i = np.power(sd,-2);
Cdi = np.diag(sd2i.ravel());
# set up the data kernel
G=np.zeros((N,2));
G[:,0]=np.ones((N,1)).ravel();
G[:,1]=x.ravel();
# solve for estimated model parameters
FTF = np.matmul( G.T , np.matmul( Cdi, G ) );
FTf = np.matmul( G.T , np.matmul( Cdi, dobs ) );
mest = la.solve(FTF,FTf);
# predict data
dpre = mest[0] + mest[1] * x;
# plot data
fig1 = plt.figure(1,figsize=(8,6));
ax1 = plt.subplot(2,1,1);
plt.axis([x[0,0], x[N-1,0], 0, 300])
plt.plot(x,dobs,'ko');
# 95% error bars
for p in range(N):
plt.plot( [x[p], x[p]], [dobs[p]-2*sd[p], dobs[p]+2*sd[p]], 'k-', lw=2 );
plt.plot(x,dpre,'k-');
plt.xlabel('x');
plt.ylabel('d');
titlestr = "intercept %.2f slope %.2f" % (mest[0,0], mest[1,0]);
plt.title(titlestr);
# Part 2: variance, high for thr first section low for the second section
sd = np.zeros((N,1));
sd[0:N2,0:1] = 2.0;
sd[N2:N,0:1] = 20.0;
sd2i = np.power(sd,-2);
Cdi = np.diag(sd2i.ravel());
# set up the data kernel
G[0:N,0:1]=np.ones((N,1));
G[0:N,1:2]=x;
# solve for estimated model parameters
FTF = np.matmul( G.T , np.matmul( Cdi, G ) );
FTf = np.matmul( G.T , np.matmul( Cdi, dobs ) );
mest = la.solve(FTF,FTf);
# predict data
dpre = mest[0,0] + mest[1,0] * x;
# plot data
ax1 = plt.subplot(2,1,2);
plt.axis([x[0,0], x[N-1,0], 0, 300])
plt.plot(x,dobs,'ko');
for p in range(N):
plt.plot( [x[p], x[p]], [dobs[p]-2*sd[p], dobs[p]+2*sd[p]], 'k-', lw=2 );
plt.plot(x,dpre,'k-');
plt.xlabel('x');
plt.ylabel('d');
titlestr = "intercept %.2f slope %.2f" % (mest[0,0], mest[1,0]);
plt.title(titlestr);
plt.show();
In [6]:
# edapy_05_04: Straight-line fit with diagonal variance
# illustrate least squares fitting of a line
# when the data have unequal variance
# this version uses F and f
# set of vector of x's
N = 51;
N2 = floor((N-1)/2);
Dx = 2;
xmin = 0.0;
xmax = xmin+Dx*(N-1);
x = eda_cvec( np.linspace(xmin, xmax, N) );
M=2; # number of model parametres
# prepare some test data hat consists of
# two sections, that differ slightly in
# interrcept and slope
a = 1.0;
b = 2.0;
dobs = a + b*x;
Da = 1;
Db = 0.3;
dobs[N2:N,0] = np.copy( dobs[N2:N,0] + Da + Db*x[N2:N,0] );
# Part 1: variance, high for thr first section low for the second section
sd = np.zeros((N,1));
sd[0:N2,0] = 20;
sd[N2:N,0] = 2;
# set up the data kernel
G=np.zeros((N,2));
G[0:N,0:1]=np.ones((N,1));
G[0:N,1:2]=x;
# set up F and f
F = np.zeros((N,M));
for i in range(N):
F[i:i+1,0:M] = G[i:i+1,0:M] / sd[i,0];
f = np.divide( dobs, sd );
# solve for estimated model parameters
FTF = np.matmul( F.T, F );
FTf = np.matmul( F.T, f );
mest = la.solve(FTF,FTf);
# predict data
dpre = mest[0] + mest[1] * x;
# plot data
fig1 = plt.figure(1,figsize=(8,6));
ax1 = plt.subplot(2,1,1);
plt.axis([x[0,0], x[N-1,0], 0, 300])
plt.plot(x,dobs,'ko');
# 95% error barx
for p in range(N):
plt.plot( [x[p], x[p]], [dobs[p]-2*sd[p], dobs[p]+2*sd[p]], 'k-', 'LineWidth', 2 );
plt.plot(x,dpre,'k-');
plt.xlabel('x');
plt.ylabel('d');
titlestr = "intercept %.2f slope %.2f" % (mest[0,0], mest[1,0]);
plt.title(titlestr);
plt.show();
In [7]:
# edapy_05_05: using smoothness to fill in data gaps
# simple example of smoothness applied
# to the problem of filling in data gaps
# this version uses non-sparse arrays
# the model parameters are the values of a sinusoidal
# function evenly spaced along the x axis
M=101;
Dx=1.0;
Dx2 = Dx**2;
xmin=0.0;
xmax = xmin+Dx*(M-1);
x = eda_cvec( np.linspace(xmin, xmax, M) );
# the data are the values the function measured on
# only a subset of points along the x axis
N=40;
# randomly choose N of the M points to have data
rowindex = np.sort( np.random.randint(M,size=N) ); # N random integers between 1 and M
xd = eda_cvec( x[rowindex,:] ); # xd are the x's of the N data
# data is sinusoidal
freq=0.02;
d = np.sin(2*pi*freq*xd);
# prior variances
# 1% error of signal level of 1 i 0.1
sigmadi = 1.0/0.01;
# maximum slope of about (1-0)/10 = 0.1
sigmaD2i = 1.0/0.1;
# maximum 2nd der about (0-2*1-0)/(20*20) = 0.005
sigmaD1i = 1.0/0.005;
# build F and f
K=M;
F=np.zeros((N+K,M));
fobs=np.zeros((N+K,1));
nextrow = 0;
# the model parameter equals the observed data
for p in range(N):
F[nextrow,rowindex[p]] = sigmadi;
fobs[nextrow,0]=sigmadi*d[p,0];
nextrow=nextrow+1;
# smoothness (small second derivative) in the interior
for p in range(M-2):
F[nextrow,p] = sigmaD2i/Dx2;
F[nextrow,p+1] = -2*sigmaD2i/Dx2;
F[nextrow,p+2] = sigmaD2i/Dx2;
fobs[nextrow,0]=0.0;
nextrow=nextrow+1;
# flatness (small first derivative) at the left end
F[nextrow,0]=-sigmaD1i*Dx;
F[nextrow,1]=sigmaD1i*Dx;
fobs[nextrow,0]=0;
nextrow=nextrow+1;
# flatness (small first derivative) at the right end
F[nextrow,M-2]=-sigmaD1i*Dx;
F[nextrow,M-1]=sigmaD1i*Dx;
fobs[nextrow,0]=0;
nextrow=nextrow+1;
print("built %d rows out of %d" % (nextrow,N+K));
# solve generalized least squares equation
FTF = np.matmul( F.T, F );
FTf = np.matmul( F.T, fobs );
mest = np.zeros((M,1));
mest[0:M,0:1] = la.solve(FTF,FTf);
# error
e = d-mest[rowindex]; # individual error
E = np.matmul( e.T, e ); E=E[0,0]; # overall error
# plot data
fig1 = plt.figure(1,figsize=(8,6));
ax1 = plt.subplot(2,1,1);
plt.axis([x[0,0], x[M-1,0], -1.5, 1.5])
plt.plot(xd,d,'ko');
plt.plot(x,mest,'k-');
plt.xlabel('x');
plt.ylabel('d');
plt.show();
built 141 rows out of 141
In [8]:
# eda05_06: sparse matrix from list of non-zero elements
# build lists
L=7;
M=6;
Frow = eda_cvec( 0, 2, 3, 5);
Fcol = eda_cvec( 0, 2, 0, 4);
Fval = eda_cvec( 1.0, 2.0, 3.0, 4.0);
KK, tmp = np.shape(Fval);
# build sparse matrix
edaFsparse = sparse.coo_matrix((Fval[0:KK,0],
(Frow[0:KK,0],Fcol[0:KK,0])), shape=(L,M) );
M=sparse.csr_matrix.todense(edaFsparse);
print(M);
[[1. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0.] [0. 0. 2. 0. 0. 0.] [3. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 4. 0.] [0. 0. 0. 0. 0. 0.]]
In [11]:
# edapy_05_07: sparse matrix implementation of Fm=f
# simple example of smoothness applied
# to the problem of filling in data gaps
# this version uses sparse arrays
# the model parameters are the values of a sinusoidal
# function evenly spaced along the x axis
M=101;
Dx=1.0;
Dx2 = Dx**2;
xmin=0.0;
xmax = xmin+Dx*(M-1);
x = eda_cvec( np.linspace(xmin, xmax, M) );
# the data are the values the function measured on
# only a subset of points along the x axis
N=40;
# randomly choose N of the M points to have data
rowindex = np.sort( np.random.randint(M,size=N) ); # N random integers between 1 and M
xd = eda_cvec( x[rowindex,:] ); # x's of the data
# data is sinusoidal
freq=0.02;
dobs = np.sin(2*pi*freq*xd);
# prior variances
# 1% error of signal level of 1 i 0.1
sigmadi = 1.0/0.01;
# maximum slope of about (1-0)/10 = 0.1
sigmaD2i = 1.0/0.1;
# maximum 2nd der about (0-2*1-0)/(20*20) = 0.005
sigmaD1i = 1.0/0.005;
# set up f
L=N+M; # rows in F, N data, and N priors
fobs=np.zeros((L,1));
# set up F. In this impletation, the non-zero
# elements of F are put first in a vector Fval, with
# two auxilary arrays Frow and Fcol
K = int;
K = N + 3*(M-2) + 2*2 + 100; # needs to be at least as large as the number of non-zero elments in F
Frow = np.zeros((K,1),dtype=np.int32); # the row index vector is integer
Fcol = np.zeros((K,1),dtype=np.int32); # the col index vector is integer
Fval = np.zeros((K,1),dtype=np.double); # the value vector is double precision
# a trick for keeping tracks of rows is to use a row counter
# an increment it after each row is added. The counters
# are defined as integers.
nextrow = int; # next row in F; note is an integer
nextrow = 0; # starting value of zero
KK = int; # next element in Fval; note is an integer
KK = 0; # starting value of zero
# the first N rows of F are that the model parameter
# equals the observed data
for p in range(N):
# F[nextrow,rowindex[p]] = 1;
Frow[KK,0]=nextrow;
Fcol[KK,0]=rowindex[p];
Fval[KK,0]=sigmadi;
KK=KK+1;
# f[nextrow]=d[p];
fobs[nextrow,0]=sigmadi*dobs[p,0];
nextrow=nextrow+1;
# implement the smoothness constraints in the interior
# of the x interval
for p in range(M-2):
# F[nextrow,p] = shi/Dx2;
Frow[KK,0]=nextrow;
Fcol[KK,0]=p;
Fval[KK,0]=sigmaD2i/Dx2;
KK=KK+1;
# F[nextrow,p+1] = -2*shi/Dx2;
Frow[KK,0]=nextrow;
Fcol[KK,0]=p+1;
Fval[KK,0]=-2*sigmaD2i/Dx2;
KK=KK+1;
# F[nextrow,p+2] = shi/Dx2;
Frow[KK,0]=nextrow;
Fcol[KK,0]=p+2;
Fval[KK,0]=sigmaD2i/Dx2;
KK=KK+1;
# f[nextrow,0]=0.0;
fobs[nextrow,0]=0.0;
nextrow=nextrow+1;
# implement flatness constraint at left end of the x interval
# F[nextrow,0]=-shi*Dx;
Frow[KK,0]=nextrow;
Fcol[KK,0]=0;
Fval[KK,0]=-sigmaD1i*Dx;
KK=KK+1;
# F[nextrow,1]=shi*Dx;
Frow[KK,0]=nextrow;
Fcol[KK,0]=1;
Fval[KK,0]=sigmaD1i*Dx;
KK=KK+1;
# f[nextrow,0]=0;
fobs[nextrow,0]=0;
nextrow=nextrow+1;
# implement flatness constraint at right end of the x interval
# F[nextrow,M-2]=-shi*Dx;
Frow[KK,0]=nextrow;
Fcol[KK,0]=M-2;
Fval[KK,0]=-sigmaD1i*Dx;
KK=KK+1;
# F[nextrow,M-1]=shi*Dx;
Frow[KK,0]=nextrow;
Fcol[KK,0]=M-1;
Fval[KK,0]=sigmaD1i*Dx;
KK=KK+1;
# f[nextrow]=0;
fobs[nextrow]=0;
nextrow=nextrow+1;
print("%d non-zero-elements, expected %d" % (KK, N+3*(M-2)+2*2) );
print("%d rows, expected %d" % (nextrow, N+M) );
# actual rows in F
# the "F" matrix must be called "edaFsparse" in order to correctly interface
# witk the FYFmul() function called by the biconjugate gradient solver
# Note the reduced length of the arrays Fval, Frow, Fcol
edaFsparse = sparse.coo_matrix((Fval[0:KK,0], (Frow[0:KK,0],
Fcol[0:KK,0])), shape=((N+M),M) );
del Fval; # delete no loger needed arrays
del Frow;
del Fcol;
# define linear operator needed for conjugate gradienet solver
LO=las.LinearOperator(shape=(M,M),
matvec=FTFmul,rmatvec=FTFmul);
# solve for eusing congugate gradient algorithm
mest = np.zeros((M,1));
q=las.bicg(LO,edaFsparse.T*fobs,rtol=1e-12, maxiter=(3*(N+M)+100));
mest = eda_cvec( q[0].astype(float) );
# error
e = dobs-mest[rowindex]; # individual errors
E = np.matmul( e.T, e ); E=E[0,0]; # overall error
print('E:',E);
# plot data
fig1 = plt.figure(1,figsize=(8,6));
ax1 = plt.subplot(2,1,1);
plt.axis([x[0,0], x[M-1,0], -1.5, 1.5])
plt.plot(xd,dobs,'ko');
plt.plot(x,mest,'k-');
plt.xlabel('x');
plt.ylabel('d');
plt.show();
341 non-zero-elements, expected 341 141 rows, expected 141 E: 1.9744793530986667e-05
In [13]:
# edapy_05_08: folding out a matrix into a vector
I = 4; # number of rows of a matrix A
J = 3; # number of coluns of a matrix A
M = I*J; # total number of elements
# make a matrix A
A=np.zeros((I,J));
A[:,:] = [ [1, 5, 9], [2, 6, 10], [3, 7, 11], [4, 8, 12]];
print('A\n', A);
# make a Mx1 vector m
m=np.zeros((M,1));
# PART 1: using a formula
print('\nPart 1: with a formula');
# copy elements of A[i,j] into m[k]
for i in range(I):
for j in range(J):
k = i*J+j; # row-wise unfold
m[k,0] = A[i,j];
print('m\n', m);
# copy the elements of m[k] into a matrix B[i,j]
B=np.zeros((I,J));
for k in range(M):
i = floor(k/J); # row-wise fold
j = k-i*J;
B[i,j]=m[k,0];
print('B\n',B);
# PART 2: with index arrays
print('\nPart 2: with index arrays')
# the idea is to define a set of "pointer arrays"
# that translate between (i,j) and k. The arrays
# are built during one traverse of (i,j). One
# advanage is that no "formulas" are needed.
i_of_k = np.zeros((M,1), dtype=np.int32); # i given k
j_of_k = np.zeros((M,1), dtype=np.int32); # j given k
k_of_ij = np.zeros((I,J), dtype=np.int32); # k given (i,j)
k=0; # initialize counter
for i in range(I):
for j in range(J):
i_of_k[k,0]=i;
j_of_k[k,0]=j;
k_of_ij[i,j]=k;
k=k+1;
M=k;
# copy A[i,j] to m[k]
m=np.zeros((M,1));
for i in range(I):
for j in range(J):
k = k_of_ij[i,j];
m[k,0] = A[i,j];
print('m\n', m);
# copy m[k] to B[i,j]
B=np.zeros((I,J));
for k in range(M):
i = i_of_k[k,0];
j = j_of_k[k,0];
B[i,j]=m[k,0];
print('B\n',B);
A [[ 1. 5. 9.] [ 2. 6. 10.] [ 3. 7. 11.] [ 4. 8. 12.]] Part 1: with a formula m [[ 1.] [ 5.] [ 9.] [ 2.] [ 6.] [10.] [ 3.] [ 7.] [11.] [ 4.] [ 8.] [12.]] B [[ 1. 5. 9.] [ 2. 6. 10.] [ 3. 7. 11.] [ 4. 8. 12.]] Part 2: with index arrays m [[ 1.] [ 5.] [ 9.] [ 2.] [ 6.] [10.] [ 3.] [ 7.] [11.] [ 4.] [ 8.] [12.]] B [[ 1. 5. 9.] [ 2. 6. 10.] [ 3. 7. 11.] [ 4. 8. 12.]]
In [29]:
# edapy_05_09: fill in data gaps in 2D using diff eqn
# use generealized least square to fill in
# data gaps for data that satisfy Laplace's
# equation in the in the interior of a grid,
# and the normal derivative is zero on the
# edges of the grid
# load the data. These data are only 'simulated'
# (also called "synthetic") and were created by the
# script eda05_09
# load (x, y, pressure) data from file
D = np.genfromtxt('../data/pressure.txt', delimiter='\t')
[N, K]=D.shape;
xobs = eda_cvec( D[:,0] );
yobs = eda_cvec( D[:,1] );
dobs = eda_cvec( D[:,2] );
# assign variance of data
sigmad = 0.1;
rsigmad = 1/sigmad;
rsigmad2 = rsigmad**2;
# rows of grid are x
I = 41;
Dx = 1;
xmin = 0;
xmax = Dx*(I-1);
x = eda_cvec( np.linspace(xmin, xmax, I) );
rDx = 1/Dx;
rDx2 = rDx**2;
# columns of grid are y
J = 41;
Dy = 1;
ymin = 0;
ymax = Dy*(J-1);
y = eda_cvec( np.linspace(ymin, ymax, J) );
rDy = 1/Dy;
rDy2 = rDy**2;
# associate observations with points on the grid
# by scaling the range (xmin, xmax) to (1, N)
rowindex = np.zeros( (N,1), dtype=np.int32 );
rowindex[:,0] = np.floor(I*(xobs-xmin)/(xmax-xmin)+0.5).ravel();
colindex = np.zeros( (N,1), dtype=np.int32 );
colindex[:,0] = np.floor(J*(yobs-ymin)/(ymax-ymin)+0.5).ravel();
# don't let any data fall outside the grid. In this
# case, we just put it on the edge of the grid. The
# alternative is to throw it away.
for i in range(N):
if( rowindex[i,0]<0 ):
rowindex[i,0]=1;
elif( rowindex[i,0]>=I ):
rowindex[i,0]=I-1;
if( colindex[i,0]<0 ):
colindex[i,0]=0;
elif( colindex[i,0]>=J ):
colindex[i,0]=J-1;
# copy the data to the grid A for later use
Aobs = np.zeros((I,J));
for i in range(N):
Aobs[rowindex[i,0],colindex[i,0]] = dobs[i,0];
# one model parameter per grid point
M = I*J;
# Fm=f is LxM with N rows of data plus M rows of prior information
L = N+M;
# F will have no more than K elements
K = N + 5*M + 100;
# r.h.s. of Fm=F
fobs = np.zeros((L,1));
# reminders on how to move from grid to vector & vice-versa
# m(k) = A(i,j);
# k = i*J+j;
# i = floor(k/J);
# j = k-i*J;
# B(i,j)=m(k);
# in this implementation, I build vectors on zon-zero elements and
# then change them into a sparse matrix at the end
Frow = np.zeros((K,1),dtype=int);
Fcol = np.zeros((K,1),dtype=int);
Fval = np.zeros((K,1));
# build F and f
# nextrow is used to keep track of the current row of F and f
nextrow = int;
nextrow = 0;
# KK is used to keep track of the current row of Fvec, Fcol, Frow
KK = int;
KK = 0;
print(type(KK),KK);
# k=dindex(i): data i corresponds to row k of m
dindex = np.zeros((N,1),dtype=int);
# first part of Fm=f: the data
for i in range(N):
k = rowindex[i,0]*J+colindex[i,0];
# F[q,k] = rsigmad;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmad;
KK = KK+1;
fobs[nextrow,0] = rsigmad*dobs[i,0];
dindex[i:0] = k;
nextrow = nextrow+1;
# error of prior information
sigmah = 1; # was 1000
rsigmah = 1/sigmah;
# second part, laplaces equation in interior
# second derivative is a (1,-2,1)/Dx^2 in x plus a (1,-2,1)/Dy^2 in y
for i in range(1,I-1): # Matlab [2:I-1]
for j in range(1,J-1): # Matlab [2:J-1]
# k = (i-1)*J+j; and F(q,k) = rsigmah*(-2*rDx2-2*rDy2);
k = i*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*(-2*rDx2-2*rDy2);
KK=KK+1;
# k = ((i-1)-1)*J+j; and F(q,k) = rsigmah*rDx2;
k = (i-1)*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDx2
KK=KK+1;
# k = ((i+1)-1)*J+j; and F(q,k) = rsigmah*rDx2;
k = (i+1)*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDx2;
KK=KK+1;
# k = (i-1)*J+(j-1); and F(q,k) = rsigmah*rDy2;
k = i*J+(j-1);
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDy2;
KK=KK+1;
# k = (i-1)*J+(j+1); and F(q,k) = rsigmah*rDy2;
k = i*J+(j+1);
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDy2;
KK=KK+1;
fobs[nextrow,0] = 0.0;
nextrow = nextrow+1;
# third part: top edge
for j in range(1,J-1): # MATLAB [2:J-1]
i = 0;
# k = (i-1)*J+j; and F(q,k) = -rsigmah*rDx;
k = i*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = -rsigmah*rDx;
KK=KK+1;
# k = ((i+1)-1)*J+j; and F(q,k) = rsigmah*rDx;
k = (i+1)*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDx;
KK=KK+1;
fobs[nextrow,0] = 0.0;
nextrow = nextrow+1;
# fourth part: bottom edge
for j in range(1,J-1): # MATLAB [2:J-1]
i = I-1;
# k = (i-1)*J+j; and F(q,k) = rsigmah*rDx;
k = i*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDx;
KK=KK+1;
# k = ((i-1)-1)*J+j; and F(q,k) = -rsigmah*rDx;
k = (i-1)*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = -rsigmah*rDx;
KK=KK+1;
fobs[nextrow,0] = 0.0;
nextrow = nextrow+1;
# fifth part: left edge
for i in range(1,I-1): # MATLAB [2:I-1]:
j = 0;
# k = (i-1)*J+j; and F(q,k) = -rsigmah*rDy;
k = i*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = -rsigmah*rDy;
KK=KK+1;
# k = (i-1)*J+(j+1); and F(q,k) = rsigmah*rDy;
k = i*J+(j+1);
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDy;
KK=KK+1;
fobs[nextrow,0] = 0.0;
nextrow = nextrow+1;
# sixth part: right edge
for i in range(1,I-1): # MATLAB [2:I-1]
j = J-1;
# k = (i-1)*J+j; and F(q,k) = rsigmah*rDy;
k = i*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDy;
KK=KK+1;
# k = (i-1)*J+(j-1); and F(q,k) = -rsigmah*rDy;
k = i*J+(j-1);
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = -rsigmah*rDy;
KK=KK+1;
fobs[nextrow,0] = 0.0;
nextrow = nextrow+1;
# seventh part: four corners
# top left
rDz = 1/sqrt((Dx^2)+(Dy^2));
i=0;
j=0;
# k = (i-1)*J+j; and F(q,k) = -rsigmah*rDz;
k = i*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = -rsigmah*rDz;
KK=KK+1;
# k = ((i+1)-1)*J+(j+1); F(q,k) = rsigmah*rDz;
k=(i+1)*J+(j+1);
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDz;
KK=KK+1;
fobs[nextrow,0]=0.0;
nextrow=nextrow+1;
# top right
i=0;
j=J-1;
# k = (i-1)*J+j; and F(q,k) = rsigmah*rDz;
k = i*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDz;
KK=KK+1;
# k = ((i+1)-1)*J+(j-1); F(q,k) = -rsigmah*rDz;
k=(i+1)*J+(j-1);
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = -rsigmah*rDz;
KK=KK+1;
fobs[nextrow,0]=0.0;
nextrow=nextrow+1;
# bottom left
i=I-1;
j=0;
# k = (i-1)*J+j; and F(q,k) = -rsigmah*rDz;
k = i*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = -rsigmah*rDz;
KK=KK+1;
# k = ((i-1)-1)*J+(j+1); and F(q,k) = rsigmah*rDz;
k = (i-1)*J+(j+1);
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDz;
KK=KK+1;
fobs[nextrow,0]=0.0;
nextrow=nextrow+1;
# bottom right
i=I-1;
j=0;
# k = (i-1)*J+j; and F(q,k) = rsigmah*rDz;
k = i*J+j;
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = rsigmah*rDz;
KK=KK+1;
# k = ((i-1)-1)*J+(j-1); and F(q,k) = -rsigmah*rDz;
k = (i-1)*J+(j-1);
Frow[KK,0] = nextrow;
Fcol[KK,0] = k;
Fval[KK,0] = -rsigmah*rDz;
KK=KK+1;
fobs[nextrow,0]=0.0;
nextrow=nextrow+1;
# actual rows in F
LL = int;
LL = nextrow;
print("%d non-zero elements" % (KK));
print("%d rows, expected %d" % (LL, (N+M)));
# the "F" matrix must be called "edaFsparse" in order to correctly interface
# witk the FYFmul() function called by the biconjugate gradient solver
# Note the reduced length of the arrays Fval, Frow, Fcol
edaFsparse = sparse.coo_matrix((Fval[0:KK,0], (Frow[0:KK,0], Fcol[0:KK,0])), shape=(LL,M) );
del Fval; # delete the arrays used to build the sparse matrix
del Frow;
del Fcol;
# define linear operator needed for conjugate gradient solver
LO=las.LinearOperator(shape=(M,M),matvec=FTFmul,rmatvec=FTFmul);
# solve for estimated model parameters using congugate gradient algrorithm
mest = np.zeros((M,1));
q=las.bicg(LO,edaFsparse.transpose()*fobs,rtol=1e-6, maxiter=(3*(N+M)+100));
mest = eda_cvec( q[0].astype(float) );
# reorganize estimated model parameters back into grid
Apre = np.zeros((I,J));
for k in range(M):
i = floor(k/J);
j = k-i*J;
Apre[i,j] = mest[k,0];
# plot results
eda_draw('title p-obs',Aobs,'title p-pre',Apre);
# compute weighted prediction error Ed
dpre=np.zeros((N,1));
for i in range(N):
k=dindex[i,0];
dpre[i,0] = mest[k,0];
e = dobs-dpre;
Ed = np.matmul( e.T, e ); Ed=Ed[0,0]
sigmadpos = sqrt(Ed/N);
print('posterior errror in data',sigmadpos);
<class 'int'> 0 8093 non-zero elements 1849 rows, expected 1849
posterior errror in data 2.341640515787999
In [30]:
# edapy_05_10: make simulated (or "synthetic") pressure data
# note: in order not to accidentally overwrite the data file
# in the distribution, this code only writes the data to
# mypressure.txt not pressure.txt; you must rename it to
# have is used by eda05_08
# rows of grid are x
I = 41;
Dx = 1;
xmin = 0;
xmax = Dx*(I-1);
xobs = eda_cvec( np.linspace(xmin, xmax, I) );
# columns of grid are y
J = 41;
Dy = 1;
ymin = 0;
ymax = Dy*(J-1);
yobs = eda_cvec( np.linspace(ymin, ymax, J) );
# randomly choose N of the I*J points to have data
N = floor(I*J/10);
rowindex = np.random.randint(I,size=(N,1));
xd = eda_cvec( xobs[rowindex[:,0],:] ); # x coords of data
colindex = np.random.randint(J,size=(N,1));
yd = eda_cvec( yobs[colindex[:,0],:] ); # y coords of data
# synthetoc data is chosen to be a solution to Laplace's equation
kappa = 0.1;
dtrue = eda_cvec(
10*np.multiply(np.sin(kappa*xd),np.exp(-kappa*yd)) );
# synthetic observed data are true data plus random noise
sigmad = 0.1;
noise = np.random.normal(loc=0.0, scale=sigmad, size=(N,1));
dobs = dtrue + noise;
# populate a matrix A with data and display it
A = np.zeros((I,J));
for i in range(N):
A[rowindex[i],colindex[i]] = dobs[i,0];
eda_draw('title pressure',A);
# save matrix to file of tab-separated values
# note that the filename has been changed to "pressure2.txt"
# so that the original file "pressure.txr" is not overwritten
# If you want to investigate the performance of eda05_08 with
# different choices of data, change the filename in eda05_08 to
# "pressure2.txt"
DD = np.concatenate( (xd,yd,dobs),axis=1); # concatenate (x,y,data) vectors into matrix
np.savetxt("../data/mypressure.txt",DD, delimiter="\t"); # write the matrix to a file
In [ ]: