In [1]:
# edapy_07_00 clear all variables and import vatious modules
# History
# 2022/06/26 - checked with most recent Python & etc. software
# 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)
%reset -f
import os
from datetime import date
from math import exp, pi, sin, sqrt, floor, ceil
import numpy as np
import scipy.sparse.linalg as las
from scipy import sparse
import scipy.linalg as la
import scipy.signal as sg
from matplotlib import pyplot as plt
import matplotlib.cm as cm
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;
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;
def GLSFilterMul(v):
# this function is used by the bicongugate gradient solver to solve the
# geneneralized least squares problem Fm=f with F=[G;eH] and f=[d;e*h],
# where G is a topplitz matrix built from the filter g
# Note that "H" must be sparse and called "edaHsparse" and the
# thefilter "g" bust be called edafilterg and must be a column vector
# and that the parameter "e" must be called eda_e
global eda_e;
global edaHsparse;
global edafilterg;
# rearrange as column-vector
s = np.shape(v);
M = s[0];
if(len(s)==1): # convert to column-vector
vv = np.zeros((s[0],1));
vv[:,0] = v;
else:
vv=v;
N, i = np.shape(edafilterg);
temp1 = np.zeros((N+M-1,1));
temp1[:,0] = np.convolve(edafilterg.ravel(),vv.ravel());
a = np.zeros((N,1));
a[:,0] = temp1[0:N,0];
b = eda_e*edaHsparse * v;
temp2 = np.zeros((N+N-1,1));
temp2[:,0] = np.convolve(
(np.flipud(edafilterg)).ravel(),a.ravel());
a2 = temp2[N-1:N+M-1,0];
b2 = eda_e*edaHsparse.transpose()*b;
# a2+b2 = FT F v = GT G v + (e^**2) HT H v
return (a2+b2);
def GLSFilterMulRHS(g,dobs,e,H,h,M):
# companion function to GLSFilterMul()
# creates RHS of generalized least squares equattion
# set up right hand side, F'f = GT qobs + e HT h
# note that H must be sparse
N,i = np.shape(dobs);
temp = eda_cvec( np.convolve(np.flipud(g).ravel(),dobs.ravel()) );
FTfa = eda_cvec( temp[N-1:N+M-1,0] );
FTfb = (e**2)*H.T*h;
FTf=FTfa+FTfb;
return (FTf);
# 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));
# 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;
In [2]:
# edapy_07_01: least squares estimation of filter input
# example using impulse response of a heat-producing layer
# with solution by damped least squares
# Note: symbol q for theta, temperture
# This is the non-sparse implementation.
# generic time setup, except make two time arrays.
# one in seconds for calulations, one in days for plotting
M=1024; # length of time series
N=M; # number of model parameters equal number of data
Dtdays = 0.5; # sampling interval in days
tdays = eda_cvec(np.linspace(0,Dtdays*(N-1),N)); # time vector (days)
Dtseconds = Dtdays*3600*24; # samping interval in seconds
tseconds = Dtseconds*tdays; # time vector (seconds)
# materal constants
rho = 1500; # density, kg/m3
cp = 800; # heat capacity, J / kg-K
kappa = 1.25e-6; # thermal diffusivity, m2/s
z = 1.0; # depth of observation, m
# impulse response, which is a complicated analytic formula
# time vector without zero
tsecondswo0 = eda_cvec( tseconds[1:M,0:1] );
P1 = (1.0/(rho*cp));
P2 = (Dtseconds/sqrt(2*pi));
P3 = np.reciprocal( np.sqrt(2*kappa*tsecondswo0));
P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tsecondswo0) );
g = eda_cvec( 0.0, P1 * P2 * np.multiply(P3 ,P4) );
# true model parameters, heat production
t1 = Dtdays*N/4; # time of heat pulse 2
t2 = Dtdays*7*N/16; # time of heat pulse 2
sd = Dtdays*N/256; # width of pulses
Ah=10; # amplitude of pulses
P5 = np.exp(-np.power(tdays-t1,2) / (2*sd**2) ); # gaussian pulse 1
P6 = np.exp(-np.power(tdays-t2,2) / (2*sd**2) ); # gaussian pulse 2
htrue = Ah*P5 + 0.5*Ah*P6; # sum pulses
# predict true temperature by convolving g and htrue
# note that np.convolve() expects Nx0 arrays, so the
# Nx1 columnn-vectors must be "raveled" before call
# to np.convolve, and then turned back into a column-vector
# by a call to eda_cvec()
qtrue = eda_cvec(np.convolve(g.ravel(),htrue.ravel()));
qtrue = np.copy(qtrue[0:N,0:1]);# keep only first N values
# plot impulse response
fig1 = plt.figure(1,figsize=(8,4));
ax1 = plt.subplot(3,1,1);
plt.axis([tdays[0,0], tdays[M-1,0], 0, 1.1*np.max(g)] );
plt.plot(tdays[0:M],g,'k-');
plt.xlabel('time t, days');
plt.ylabel('g(t)');
# plot true heat production
ax1 = plt.subplot(3,1,2);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:M,0:1],htrue,'k-');
plt.xlabel('time t, days');
plt.ylabel('htrue(t)');
# plot true temperature
ax1 = plt.subplot(3,1,3);
plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qtrue)] );
plt.plot(tdays[0:N,0:1],qtrue,'k-');
plt.xlabel('time t, days');
plt.ylabel('qtrue(t)');
plt.show();
# Toeplitz matrix version of impulse response
Gcol = eda_cvec( g, np.zeros((N-M,1)));
Grow = eda_cvec( g[0,0], np.zeros((M-1,1)));
G = la.toeplitz( Gcol, Grow.T );
# temperature = G times heat_production
qtrue2 = np.matmul(G,htrue);
# create simulated observations, qobs = qtrue+noise
sigmad=0.001*np.amax(qtrue); # noise is 0.1% amplitude of signal
qobs=qtrue+np.random.normal(0,sigmad,(N,1));
# plot simulated observations
fig2 = plt.figure(2,figsize=(8,4));
ax1 = plt.subplot(3,1,1);
plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );
plt.plot(tdays[0:N,0:1],qobs,'k-');
plt.xlabel('time t, days');
plt.ylabel('qobs(t)');
# solve for h given qobs and G using damped least squares
GTG = np.matmul(G.T,G); # G-transpose times G
GTGmax = np.max(np.abs(GTG)); # maximum value of G
e2=1e-4*GTGmax; # damping is a fraction of GTG
hest1 = la.solve(GTG+e2*np.identity(M), np.matmul(G.T,qobs) );
# plot true heat production
ax1 = plt.subplot(3,1,2);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:N,0:1],htrue,'k-');
plt.xlabel('time t, days');
plt.ylabel('htrue(t)');
# plot estimated heat production
ax1 = plt.subplot(3,1,3);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:N,0:1],hest1,'k-');
plt.xlabel('time t, days');
plt.ylabel('hest1(t)');
plt.show();
In [3]:
# edapy_07_02: estimation of filter input, with smoothing
# using impulse repsonse of a heat-producing layer
# with solution via generalized least squares
# with the prior information of smoothness
# Note: symbol q for theta, temperture.
# This is the non-sparse implementation.
# generic time setup, except make two time arrays.
# one in seconds for calulations, one in days for plotting
M=1024; # length of time series
N=M; # number of model parameters equal number of data
Dtdays = 0.5; # sampling interval in days
tdays = eda_cvec(np.linspace(0,Dtdays*(N-1),N)); # time vector (days)
Dtseconds = Dtdays*3600*24; # samping interval in seconds
tseconds = Dtseconds*tdays; # time vector (seconds)
# materal constants
rho = 1500.0; # density, kg/m3
cp = 800.0; # heat capacity, J / kg-K
kappa = 1.25e-6; # thermal diffusivity, m2/s
z = 1.0; # depth of observation, m
# impulse response, which is a analysic formula that's
# fairly complicated, so break it into three pieces,
# P1, P2 and P3
# impulse respone is g, but time=0 will be handled separately
tsecondswo0 = eda_cvec( tseconds[1:M,0] ); # time vector without zero
P1 = (1.0/(rho*cp));
P2 = (Dtseconds/sqrt(2*pi));
P3 = np.reciprocal( np.sqrt(2*kappa*tsecondswo0));
P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tsecondswo0) );
g = eda_cvec( 0.0, P1 * P2 * np.multiply(P3 ,P4) );
#true model parameters, heat production
t1 = Dtdays*N/4; # time of heat pulse 1
t2 = Dtdays*7*N/16; # time of heat pulse 2
sd = Dtdays*N/256; # width of pulses
Ah = 10.0; # amplitude of pulses
P5 = np.exp(-np.power(tdays-t1,2) / (2*sd**2) ); # gaussian pulse 1
P6 = np.exp(-np.power(tdays-t2,2) / (2*sd**2) ); # gaussian pulse 2
htrue = Ah*P5 + 0.5*Ah*P6; # sum pulses
# Toeplitz matrix version of impulse response
Gcol = eda_cvec( g, np.zeros((N-M,1)));
Grow = eda_cvec( g[0,0], np.zeros((M-1,1)));
G = la.toeplitz( Gcol, Grow.T );
# temperature = G times heat_production
qtrue2 = np.matmul(G,htrue);
# create simulated observations, qobs = qtrue+noise
# noise leve is 1% of maximum temperature
sigmad=0.01*np.max(np.abs(qtrue2));
qobs=qtrue2+np.random.normal(0,sigmad,(N,1));
# solve for h given qobs and g using damped least squared
GTG = np.matmul(G.T,G); # G_transpose time G
GTGmax = np.amax(GTG); # maximum value of GTG
e2=1e-4*GTGmax; # damping is fraction of GTG
# damped least squares solution
hest2 = la.solve(GTG+e2*np.identity(M), np.matmul(G.T,qobs) );
# generalized least squares with F = [G; H] and f = [qobs, 0]
# second derivative smoothing in interior of interval; nothing on ends
# H is a Toeplitz matrix with diaginals 1, -2 and 1. Note
# that he factor of 1/(Dt^2) is omitted because the r.h.s. is zero.
# The damping cnsnstant e = sigma_h/sigma_d determines the
# relative weighting of Gm=d and Hm=h. It is being tuned
# here by trail and error to achieve the desired degree of smoothness
e=10*np.max(np.abs(G));# e is fraction of maximum value of G
Hrow = eda_cvec( 1.0, -2.0, 1.0, np.zeros((M-3,1)));
Hcol = eda_cvec( 1.0, np.zeros((M-3,1)));
H = la.toeplitz( Hcol, Hrow.T ); # construct H
h = np.zeros((M-2,1)); # since prior info is that 2nd derivative is zero, h=0
F = np.concatenate( (G,e*H), axis=0 ); # combine G and H into F
f = np.concatenate( (qobs,e*h), axis=0 ); # combine qobs and h into f
# deneralized least squares solution
FTF = np.matmul(F.T,F); # F_transpose times F
FTf = np.matmul(F.T,f); # F_transpose times f
hest3 = la.solve(FTF, FTf); # generalized least squares solution
# figure 1
fig1 = plt.figure(1,figsize=(8,4));
# plot observations
ax1 = plt.subplot(3,1,1);
plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );
plt.plot(tdays[0:N,0:1],qobs,'k-');
plt.xlabel('time t, days');
plt.ylabel('htrue(t)');
plt.title('Damped Least Squares');
# plot true heat production
ax1 = plt.subplot(3,1,2);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:N,0:1],htrue,'k-');
plt.xlabel('time t, days');
plt.ylabel('htrue(t)');
# plot estimated heat production
ax1 = plt.subplot(3,1,3);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:N,0:1],hest2,'k-');
plt.xlabel('time t, days');
plt.ylabel('hest(t)');
plt.show();
# figure 2
fig2 = plt.figure(2,figsize=(8,4));
# plot observations
ax1 = plt.subplot(3,1,1);
plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );
plt.plot(tdays[0:N,0:1],qobs,'k-');
plt.xlabel('time t, days');
plt.ylabel('htrue(t)');
plt.title('Generalized Least Squares with 2nd Derivative Smoothing');
# plot true heat production
ax1 = plt.subplot(3,1,2);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:N,0:1],htrue,'k-');
plt.xlabel('time t, days');
plt.ylabel('htrue(t)');
# plot estimated heat production
ax1 = plt.subplot(3,1,3);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:N,0:1],hest3,'k-');
plt.xlabel('time t, days');
plt.ylabel('hest1(t)');
plt.show();
In [4]:
# edapy_07_03: sparse, least sq estimation of filter input
# using impulse repsonse of a heat-producing layer
# with solution via generalized least squares
# with the prior information of smoothness
# Note: symbol q for theta, temperture
# This is the sparse implementation.
# generic time setup, except make two time arrays.
# one in seconds for calulations, one in days for plotting
M=1024; # length of time series
N=M; # number of model parameters equal number of data
Dtdays = 0.5; # sampling interval in days
tdays = eda_cvec(np.linspace(0,Dtdays*(N-1),N)); # time vector (days)
Dtseconds = Dtdays*3600*24; # samping interval in seconds
tseconds = Dtseconds*tdays; # time vector (seconds)
# materal constants
rho = 1500.0; # density, kg/m3
cp = 800.0; # heat capacity, J / kg-K
kappa = 1.25e-6; # thermal diffusivity, m2/s
z = 1.0; # depth of observation, m
# impulse response, which is a analysic formula that's
# fairly complicated, so break it into three pieces,
# P1, P2 and P3
# impulse respone is g, but time=0 will be handled separately
tsecondswo0 = eda_cvec( tseconds[1:M,0] ); # time vector without zero
P1 = (1.0/(rho*cp));
P2 = (Dtseconds/sqrt(2*pi));
P3 = np.reciprocal( np.sqrt(2*kappa*tsecondswo0));
P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tsecondswo0) );
g = eda_cvec( 0.0, P1 * P2 * np.multiply(P3 ,P4) );
edafilterg = g; # edafilterg is a global variable defined in GLSFilterMul()
# true heat production, sum of two pulses
t1 = Dtdays*N/4; # time of first pulse
t2 = Dtdays*7*N/16; # time of second pulse
sd = Dtdays*N/256; # variance of pulse
Ah=10; # amplitude of pulse
# sum of two gaussian pulses
P5 = np.exp(-np.power(tdays-t1,2) / (2*sd**2) );
P6 = np.exp(-np.power(tdays-t2,2) / (2*sd**2) );
htrue = Ah*P5 + 0.5*Ah*P6;
# predict true temperature by convolving g and htrue
# note that np.convolve() expects Nx0 arrays, so the
# Nx1 columnn-vectors must be "raveled" before call
# to np.convolve, and then turned back into a column-vector
# by a call to eda_cvec()
qtrue = eda_cvec(np.convolve(g.ravel(), htrue.ravel()));
qtrue = eda_cvec(qtrue[0:N,0:1]); # keep only first N values
# create noisy simulated observations, qobs = qtrue+noise
sigmad=0.01*np.max(np.abs(qtrue)); # fraction of signal
qobs=qtrue+np.random.normal(0,sigmad,(N,1))
# for comparison purposes, solve for h given qobs and G using damped least squares
# Toeplitz matrix version of impulse response
Gcol = eda_cvec( g, np.zeros((N-M,1)));
Grow = eda_cvec( g[0,0], np.zeros((M-1,1)));
G = la.toeplitz( Gcol, Grow.T );
# predict true temperature from true heat prodcution
qtrue2 = np.matmul(G,htrue);
# solve with damped least squares
GTG = np.matmul(G.T,G); # G_transpose times G
GTGmax = np.max(np.abs(GTG)); # maximum value of GTG
e2=1e-4*GTGmax; # damping is fraction of maximum value of GTG
# damped least squares ssolution
hest2 = la.solve(GTG+e2*np.identity(M), np.matmul(G.T,qobs) );
# generalized leasr squares with F = [G; H] and f = [qobs, 0]
# second derivative smoothing in interior of interval; nothing on ends
eda_e=10.0*np.amax(G); # note damoing is a global variable
# Implement prior information of smoothness (2nd derivative = 0).
# Note that a factor of 1/(Dt**2) is omitted because the r.h.s. is zero.
# H is Toeplitz, so use sparse.spdiags() to create a sparse version of it.
# The [0,1,2] means the non-zero diagonals are the main diagonal
# (named 0) and two immediateley the right of it (the 1 and 2), The
# (D,-2.0*D,D) are the corresponding values of the diagonals.
# The sparse matrix is K by M.
K=M-2;
D = np.ones((1,M));
edaHsparse = sparse.spdiags(
np.concatenate((D,-2.0*D,D),axis=0), [0, 1, 2], K, M);
h = np.zeros((M-2,1));
# set up right hand side, F'f = GT qobs + e HT h
FTf = GLSFilterMulRHS(edafilterg,qobs,eda_e,edaHsparse,h,M)
# define linear operator needed for conjugate gradienet solver
LO=las.LinearOperator(shape=(M,M),matvec=GLSFilterMul,rmatvec=GLSFilterMul);
# solve for estimated model parameters using congugate gradient algrorithm
hest3 = np.zeros((M,1));
q=las.cg(LO,FTf,rtol=1e-12, maxiter=(3*(N+M)+100));
# note that q is a tupple, the first element of which is the solution
hest3 = eda_cvec(q[0].astype(float));
# figure 1
fig1 = plt.figure(1,figsize=(8,4));
# plot observations
ax1 = plt.subplot(3,1,1);
plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );
plt.plot(tdays[0:N],qobs,'k-');
plt.xlabel('time t, days');
plt.ylabel('htrue(t)');
plt.title('Damped Least Squares');
# plot true heat production
ax1 = plt.subplot(3,1,2);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:N],htrue,'k-');
plt.xlabel('time t, days');
plt.ylabel('htrue(t)');
# plot estimated heat production
ax1 = plt.subplot(3,1,3);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:N],hest2,'k-');
plt.xlabel('time t, days');
plt.ylabel('hest(t)');
plt.show();
# figure 2
fig2 = plt.figure(2,figsize=(8,4));
# plot observations
ax1 = plt.subplot(3,1,1);
plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );
plt.plot(tdays[0:N],qobs,'k-');
plt.xlabel('time t, days');
plt.ylabel('htrue(t)');
plt.title('Generalized Least Squares with 2nd Derivative Smoothing');
# plot true heat production
ax1 = plt.subplot(3,1,2);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:N],htrue,'k-');
plt.xlabel('time t, days');
plt.ylabel('htrue(t)');
# plot estimated heat production
ax1 = plt.subplot(3,1,3);
plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
plt.plot(tdays[0:N],hest3,'k-');
plt.xlabel('time t, days');
plt.ylabel('hest1(t)');
plt.show();
In [5]:
# edapy_07_04: filter estimation problem
# example of estimating a filter r that predicts river
# discharge d from precipitation data g; that is, g*r=d
# in this example, the filter length is chosen to be
# M<<N, to embody the idea that the response of the river
# (meaning the length M of the filter r) is short
M=25; # length of filter r
# load merged discharge and precipitation data
# column 1: time in inter count of days starting 1/1/2002
# column 2: Hudson River discharge in cfs
# column 3: Albany NY precipitation in inches
D = np.genfromtxt('../Data/precip_discharge_data.txt', delimiter='\t')
[Nlong, K]=D.shape;
tlong = eda_cvec(D[:,0]); # time in days
dlong = eda_cvec(D[:,1]/35.3146); # discharge in cubic meters per second
glong = eda_cvec(D[:,2]*25.4) # precipitation, im millimeters pper day;
Dt = tlong[1,0]-tlong[0,0]; # sampling interval
# in this example, we use only a 101 day portion of the
# dataset, mainly so that when we plot it, storm events
# are clearly visible
# select small segment of time, discharge and precip data
N=101;
istart = 900-1;
t = eda_cvec( tlong[istart:istart+N,0] );
d = eda_cvec( dlong[istart:istart+N,0] );
g = eda_cvec( glong[istart:istart+N,0] );
tmin = t[0,0]; # start time
tmax = t[N-1,0]; # end time
edafilterg = g; # edafilterg is a global variable defined in GLSFilterMul()
# figure 1: discharge and precipitation data
fig1 = plt.figure(1,figsize=(8,4));
# plot observations
ax1 = plt.subplot(2,1,1);
plt.axis([tmin, tmax, -20, 60 ] );
plt.plot(t,g,'k-');
plt.xlabel('time t, days');
plt.ylabel('precipitation (mm/day)');
ax1 = plt.subplot(2,1,2);
plt.axis([tmin, tmax, -1000, 2000 ] );
plt.plot(t,d,'k-');
plt.xlabel('time t, days');
plt.ylabel('discharge (m3/s)');
plt.show();
# Given the filter equation g*r=d, solve for r given d and g
# using genralized least squares
# Note that the convolution relationship g*r=d is equivalent
# to the matrix equation G r = d where G is the (g*) operator,
# Since we are using filterfun(), we do not need to construct
# G. We do, however, need to implement the prior information
# equation Hr=h.
# This code can toggle between the prior information of r
# being small, its first derivative being small and its second
# derivative being small, depending on whether DERIVATIVE is 0,
# 1 or 2. (However, the results turn out not to be especially
# sensitive to the choice).
DERIVATIVE=1;
if( DERIVATIVE==2 ):
# The prior information equation H*r = h is implemented
# for the smoothness condition of the second derivative
# being close to zero
eda_e=2*np.max(np.abs(g)); # the coefficient was manually tuned
# to produce a reasonably smooth filter
K=M-2;
L=N+K;
# second derivative operator; no end conditions
D = np.ones((1,M));
edaHsparse = sparse.spdiags( np.concatenate((D,-2.0*D,D),axis=0), [0, 1, 2], K, M);
h = np.zeros((K,1));
elif (DERIVATIVE==1):
# The prior information equation H*r = h is implemented
# for the smoothness condition of the first derivative
# being close to zero
eda_e=4.0*np.max(np.abs(g)); # the coefficient was manually tuned
# % to produce a reasonably smooth filter
K=M-1;
L=N+K;
D = np.ones((1,M));
edaHsparse = sparse.spdiags( np.concatenate((D,-D),axis=0), [0, 1], K, M);
h = np.zeros((K,1));
else: # DERIVATIVE == 0
# The prior information equation H*r = h is implemented
# for the smallness condition of the filter is close to zero
eda_e=0.1*np.max(np.abs(g)); # the coefficient was manually tuned
# to produce a reasonably small filter
K=M;
L=N+K;
D = np.ones((1,M));
edaHsparse = sparse.spdiags( D, [0], K, M);
h = np.zeros((K,1));
# set up right hand side, F'f = GT qobs + e HT h
FTf = GLSFilterMulRHS(edafilterg,d,eda_e,edaHsparse,h,M);
# define linear operator needed for conjugate gradienet solver
LO=las.LinearOperator(shape=(M,M),matvec=GLSFilterMul,rmatvec=GLSFilterMul);
# solve for estimated model parameters using congugate gradient algrorithm
# note las.cg() returns a tupple, the first element of which is the solution
q=las.cg(LO,FTf,rtol=1e-12, maxiter=(3*(N+M)+100));
r = eda_cvec( q[0].astype(float) );
# filter time axis
tr = eda_cvec( Dt*np.linspace(0,M-1,M) );
# predicted discharge, convolve and truncate to length N
dpre = eda_cvec( np.convolve(g.ravel(),r.ravel()) );
dpre = dpre[0:N,0:1];
e = d-dpre;
E = np.matmul(e.T,e); E=E[0,0];
print("total error E: %e" % (E));
# figure 2: filter and predicted discharge
fig2 = plt.figure(2,figsize=(8,4));
# plot observations
ax1 = plt.subplot(2,1,1);
plt.axis([0, M, -5, 10 ] );
plt.plot(tr,r,'k-');
plt.plot([tr[0,0], tr[M-1,0]],[0,0],'k:');
plt.xlabel('time t, days');
plt.ylabel('filter r(t)');
ax1 = plt.subplot(2,1,2);
plt.axis([tmin, tmax, -1000, 2000 ] );
plt.plot(t,d,'k:');
plt.plot(t,dpre,'k-');
plt.xlabel('time t, days');
plt.ylabel('discharge (m3/s)');
plt.show();
total error E: 2.316860e+06
In [6]:
# edapy_07_05: prediction error filter (pef)
# this examples uses the for Neuse Riverr Hydrograph
# this code uses generalized least squares to implement
# the prior information that pef(1)=(-1). Two solutions
# are calculated, one via standard matrix algebra,
# (FTF)\(FTf) and the other via the biconjugate gradient method
# read in data
D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
[N, K]=D.shape;
t = eda_cvec( D[:,0] ); # time in days
d = eda_cvec( D[:,1] ); # discharge in cubic meters per second
Dt = t[1,0]-t[0,0]; # sampling interval
tmin = t[0,0]; # start time
tmax = t[N-1,0]; # end time
# sizes of various arrays
M=100; # length of pef is user tunable
K=1; # one piece of prior information
L=N+K; # total number of equations
# for a prediction error filter, pef, the data is the filter, g
g = d;
edafilterg = d; # edafilterg is a global variable defined in GLSFilterMul()
# prior information that first element of the pef is (-1)
eda_e=(10^1)*np.amax(g); # large weight
H=np.zeros((K,M)); # the pror information equation, Hm=h, such that m(1)=(-1);
H[0,0]=1.0;
h=eda_cvec( (-1.0) );
# Data Kernel (needed only for standard matrix algebra calculation
# data kernel is Toeoplitz, so use standard setup
# Toeplitz matrix version of impulse response
Gcol = g;
Grow = eda_cvec( g[0,0], np.zeros((M-1,1)));
G = la.toeplitz( Gcol, Grow.T );
# for pef, Gm=rhs=0
# F needed only for standard matrix algebra calculation
F = np.concatenate( (G, eda_e*H), axis=0 );
rhs = np.zeros((N,1));
f = np.concatenate( (rhs, eda_e*h), axis=0 );
# standard GLS solution
FTF = np.matmul( F.T, F ); # F_transpose times F
FTf = np.matmul( F.T, f ); # F_transpose times f
pef1 = la.solve(FTF, FTf ); # least squares solution
# now set up for biconjugate gradient method
# Hm=h
eda_e=(10^1)*np.amax(g); # large weight
# edaHsparse is global variable defined in GLSFilterMul()
edaHsparse = sparse.diags([1.0], [0], shape=(1, M));
h=np.zeros((1,1));
h[0,0]=(-1.0);
# FT f can be done analytucally:
# FTf = GT 0 + e HT e h = 0 + [e 0 0 ... 0]T [-e] = [-e^2 0 0 ... 0]T
FTf = np.zeros((M,1));
FTf[0,0] = -(eda_e**2);
# define linear operator needed for conjugate gradienet solver
LO=las.LinearOperator(shape=(M,M),matvec=GLSFilterMul,rmatvec=GLSFilterMul);
# solve for estimated model parameters using congugate gradient algrorithm
r = np.zeros((M,1));
# conjugate gradient solver returns soln vector as first element of a tupple
q=las.cg(LO,FTf,rtol=1e-12, maxiter=(3*(N+M)+100));
pef2 = eda_cvec(q[0].astype(float));
# compute prediction error by colvolving g with pef, then truncate to length N
temp = np.convolve(pef2.ravel(),g.ravel());
perror = eda_cvec( temp[0:N] );
# short time axis
tpef = t[0:M,0:1];
# plot pef
fig1 = plt.figure(1,figsize=(8,4));
# plot prediction error filters
ax1 = plt.subplot(2,1,1);
plt.axis([ 0, Dt*(M-1), np.min(pef1), np.max(pef1)] );
plt.plot(tpef,pef1,'k-');
plt.xlabel('time t, days');
plt.ylabel('pef(t) (std)');
ax1 = plt.subplot(2,1,2);
plt.axis([ 0, Dt*(M-1), np.min(pef1), np.max(pef1)] );
plt.plot(tpef,pef2,'k-');
plt.xlabel('time t, days');
plt.ylabel('pef(t) (std)');
plt.show();
# plot hydrograph data and prediction error for first year
fig2 = plt.figure(2,figsize=(8,4));
Np = 365;
ax1 = plt.subplot(2,1,1);
plt.axis([tmin, tmin+Dt*Np, -5000, 15000]);
plt.plot(t[0:Np,0],d[0:Np,0],'k-');
plt.xlabel('time t, days');
plt.ylabel('d(t)');
ax1 = plt.subplot(2,1,2);
plt.axis([tmin, tmin+Dt*Np, -5000, 15000]);
plt.plot(t[0:Np,0],perror[0:Np,0],'k-');
plt.xlabel('time t, days');
plt.ylabel('e(t)');
plt.show();
In [9]:
# edapy_07_06 Script supporting CribSheet 7.2
# Least squares filter estimation of filter m
# obeying g*m=d
# standard set-up of time axis, t
N=101; # number of data
Dt = 1.0; # sampling interval
t = eda_cvec(np.linspace(0,Dt*(N-1),N));
tmin = t[0,0]; # start time
tmax = t[N-1,0]; # end time
# make some synthetic timeseries g
# draw from an exponential distribution so its all positive
g = np.random.exponential(1,(N,1));
# edafilterg is a global variable defined in GLSFilterMul()
edafilterg = g;
# exemplary true filter is decaying exponential
M=10; # length of filter
c=5.0*Dt; # decay constant
mtrue = eda_cvec(0.3*np.exp(-t[0:M,0]/c));
# the true data is g covolved with m; then truncate to length N
dtrue = eda_cvec(np.convolve(g.ravel(),mtrue.ravel()));
dtrue = dtrue[0:N,0:1];
# the observed data is the true data plus noise
# the noise is Normally distributed with zero mean and variance sigmad^2
sigmad = 0.1; # standard deviation of noise
noise = np.random.normal(0,sigmad,(N,1)); # random noise
d = dtrue + noise; # data = signal + noise
# now invert for the filter using the observed data and the time series g
# prior infoirmation of smoothness, with no end conditions
eda_e = 0.001; # 1 damping factor = sigma_h / sigma_d
# make Toeplitz matrix of second derivatives
K = M-2;
D = np.ones((1,M))/(Dt**2);
DIAGS = np.concatenate((D,-2.0*D,D),axis=0);
edaHsparse = sparse.spdiags(DIAGS, [0, 1, 2], K, M);
h = np.zeros((M-2,1)); # prior info that second derivative is zero,
# set up right hand side, F'f = GT qobs + e HT h
FTf = GLSFilterMulRHS(edafilterg,d,eda_e,edaHsparse,h,M);
# define linear operator needed for conjugate gradienet solver
LO=las.LinearOperator(shape=(M,M),matvec=GLSFilterMul,rmatvec=GLSFilterMul);
# solve for estimated model parameters using congugate gradient algrorithm
# returned value 1 is a tuple whose first element is the solution vector
q=las.cg(LO,FTf,rtol=1e-12, maxiter=(3*(N+M)+100));
mest = eda_cvec(q[0].astype(float));
# predicted data by colvolving mest and g, then truncat to length N
dpre = eda_cvec(np.convolve(g.ravel(),mest.ravel()));
dpre = dpre[0:N,0:1]; # trucate
# compute the error
e = d-dpre; # individual errors
E = np.matmul(e.T, e); E=E[0,0]; # total error
# Figure 1
fig1 = plt.figure(1,figsize=(8,4));
ax1 = plt.subplot(2,1,1);
plt.axis([tmin, tmax, -5, 5] );
plt.plot(t,g,'k-');
plt.xlabel('t');
plt.ylabel('g(t)');
ax1 = plt.subplot(2,1,2);
plt.axis([tmin, tmax, -5, 5] );
plt.plot(t,d,'k-',lw=1);
plt.plot(t,dpre,'k:',lw=3);
plt.xlabel('t');
plt.ylabel('d(t)');
plt.show();
# Figure 2
fig1 = plt.figure(1,figsize=(8,4));
ax1 = plt.subplot(1,1,1);
plt.axis([t[0,0], t[M-1,0], 0, 0.5] );
plt.plot( t[0:M], mtrue.ravel(),'k-');
plt.plot( t[0:M], mest.ravel(), '-', color = (0.8,0.8,0.8));
plt.xlabel('t');
plt.ylabel('m(t)');
plt.show();
In [10]:
# edapy_07_07: inverse filter via z-transforms
# use z-transforms to:
# check if a filter is minumum phase
# represent it as a cascade of length-2 filters
# compute the inverse filter
# define a filter
g = eda_cvec( [6, 5.5, 5.3, 5, 4, 3, 2] );
N,K = np.shape(g);
gN = g[N-1,0];
# find roots of g
pol = (np.flipud(g)).ravel(); # pol is a list of coefficients
r = np.roots(pol); # r is a list of roots
# check size of roots
badroots=0;
for j in range(N-1):
if( abs(r[j]) <= 1 ):
badroots=badroots+1;
printstr = "%d roots not outside the unit circle" % (badroots);
print(printstr);
# rebuild filter from its roots
gp = np.zeros((1,1), dtype=complex);
gp[0,0]=gN;
for j in range(N-1):
temp = np.convolve(gp.ravel(),[-r[j],1]);
K = len(temp);
gp = np.zeros((K,1), dtype=complex);
gp[:,0] = temp;
# there should be no imiginary part
gp = np.real(gp);
# check that it worked by computing error
E1 = np.matmul( (g-gp).T , (g-gp) ); E1=E1[0,0];
printstr = "error of filter rebuilt from its roots: %f" % (E1);
print(printstr);
# construct inverse filter, gi, of length Ni
Ni = 50;
gi = np.zeros((Ni,1),dtype=complex);
gi[0,0]=1/gN;
# filter cascade, one filter per loop
for j in range(N-1):
# construct inverse filter of a length-2 filter
temp = np.zeros((Ni,1),dtype=complex);
temp[0,0] = 1;
for k in range(1,Ni):
temp[k,0]=temp[k-1,0]/r[j];
temp = -temp/r[j];
temp2 = eda_cvec(np.convolve(gi.ravel(),temp.ravel()));
gi[0:Ni,0:1]=temp2[0:Ni,0:1];
# delete imaginary part (which should be zero)
gi = np.real(gi);
# test quality of inverse filter; truncate to length Ni
gig=eda_cvec(np.convolve(g.ravel(),gi.ravel()));
gig=eda_cvec( gig[0:Ni,0:1] );
# impulse function
d = np.zeros((Ni,1));
d[0,0]=1.0;
E2 = np.matmul( (gig-d).T, (gig-d) ); E2=E2[0,0]; # error
printstr = "error of inverse filter: %f" % (E2);
print(printstr);
# plot results
fig1 = plt.figure(1,figsize=(8,4));
# plot filter
ax1 = plt.subplot(3,1,1);
plt.axis([0, Ni, np.min(g), np.max(g)] );
gx = np.concatenate( (g, np.zeros((Ni-N,1))), axis=0);
plt.plot( np.linspace(0,Ni-1,Ni), gx.ravel(),'k-');
plt.xlabel('element j');
plt.ylabel('g(j)');
ax1 = plt.subplot(3,1,2);
plt.axis([0, Ni, np.min(gi), np.max(gi)] );
plt.plot( np.linspace(0,Ni-1,Ni), gi.ravel(),'k-');
plt.xlabel('element j');
plt.ylabel('ginv(j)');
ax1 = plt.subplot(3,1,3);
plt.axis([0, Ni, -1.2, 1.2] );
plt.plot( np.linspace(0,Ni-1,Ni), gig.ravel(),'k-');
plt.xlabel('element j');
plt.ylabel('[g*ginv](j)');
0 roots not outside the unit circle error of filter rebuilt from its roots: 0.000000 error of inverse filter: 0.000000
In [11]:
# edapy_07_08, recursive smoothing filter, by for-loops
# example of recursive smoothing filter
# this version used for-loops, for demonstration purposes
# the method in edapy_07_09 is prefered
# standard set-up of time axis, t
N=101; # number of data
Dt=1.0; # sampling interval
t = eda_cvec( Dt*np.linspace(0,N-1,N) ); # time vector
# create two synthetic datasets
# h1, a unit spike
# h2, random noise with zero mean and unit variance
h1=np.zeros((N,1)); # h1 is a time series with a spike in its middle
No2 = floor(N/2); # middle at N/2
h1[No2,0]=1.0; # spike has unit amplitude
h2=np.random.normal(0.0,1.0,(N,1)); # h2 is a time series of random noise
# apply recursive filter to h1 to yield q1
q1=np.zeros((N,1));
q1[0,0]=0.5*h1[0,0];
for j in range(2,N):
q1[j,0]=0.5*(h1[j,0]+q1[j-1,0]);
# apply recursive filter to h2 to yield q2
q2=np.zeros((N,1));
q2[0,0]=0.5*h2[0,0];
for j in range(2,N):
q2[j,0]=0.5*(h2[j,0]+q2[j-1,0]);
# plot results
fig1 = plt.figure(1,figsize=(8,4));
ax1 = plt.subplot(2,1,1);
plt.axis( [t[0,0], t[N-1,0], -1.1, 1.1] );
plt.plot( t, h1, 'k.');
plt.plot( t, q1, '-', color=(0.7,0.7,0.7));
plt.xlabel('time t');
plt.ylabel('h1(t) and q1(t)');
ax1 = plt.subplot(2,1,2);
plt.axis([t[0,0], t[N-1,0], -5, 5]);
plt.plot( t, h2, 'k.');
plt.plot( t, q2, '-', color=(0.7,0.7,0.7));
plt.xlabel('time t');
plt.ylabel('h2(t) and q2(t)');
plt.show();
displaystr = "sum of elements of the unit spike: %f" % (np.sum(h1));
print(displaystr);
displaystr = "sum of elements of the filtered spike: %f" % (np.sum(q1));
print(displaystr);
sum of elements of the unit spike: 1.000000 sum of elements of the filtered spike: 1.000000
In [12]:
# edapy_07_09: recursive smoothing filter, by sg.lfilter()
# example of recursive smoothing filter
# this version used scipy.signal.lfilter()
# standard set-up of time axis, t
N=101; # number of data
Dt=1.0; # sampling interval
t = eda_cvec(np.linspace(0,Dt*(N-1),N)); # time vector
# create two synthetic datasets
# h1, a unit spike
# h2, random noise with zero mean and unit variance
h1=np.zeros((N,1)); # h1 is a time series with a spike in its middle
No2 = floor(N/2); # middle position
h1[No2,0]=1.0; # unit spike
h2=np.random.normal(0.0,1.0,(N,1)); # h2 is a random time series
Nu = 2; # length of filter u
Nv = 2; # length of filter v
u = eda_cvec( [0.5, 0.0] ); # filter u
v = eda_cvec([1.0, -0.5] ); # filter v
# filter to h1 to q1
q1 = eda_cvec(sg.lfilter(u.ravel(),v.ravel(),h1.ravel()));
q1 = eda_cvec( q1[0:N,0:1] ); # truncate to length N
# apply filter to h2 to q2
q2 = eda_cvec(sg.lfilter(u.ravel(),v.ravel(),h2.ravel()));
q2 = eda_cvec( q2[0:N,0:1] ); # truncate to length N
# plot results
fig1 = plt.figure(1,figsize=(8,4));
ax1 = plt.subplot(2,1,1);
plt.axis( [t[0,0], t[N-1,0], -1.1, 1.1] );
plt.plot( t, h1, 'k.');
plt.plot( t, q1, '-', color=(0.7,0.7,0.7));
plt.xlabel('time t');
plt.ylabel('h1(t) and q1(t)');
ax1 = plt.subplot(2,1,2);
plt.axis([t[0,0], t[N-1,0], -5, 5]);
plt.plot( t, h2, 'k.');
plt.plot( t, q2, '-', color=(0.7,0.7,0.7));
plt.xlabel('time t');
plt.ylabel('h2(t) and q2(t)');
plt.show();
displaystr = "sum of elements of the unit spike: %f" % (np.sum(h1));
print(displaystr);
displaystr = "sum of elements of the filtered spike: %f" % (np.sum(q1));
print(displaystr);
sum of elements of the unit spike: 1.000000 sum of elements of the filtered spike: 1.000000
In [ ]: