In [7]:
# edapy_13_00 clear all variables and import vatious modules
# 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, and similar issues
# associated with depreciation of equivalence of 1x1 matrix and scalar
%reset -f
import os
from datetime import date, datetime
from math import exp, pi, sin, cos, sqrt, floor, ceil, atan2, log
import numpy as np
import scipy.sparse.linalg as las
from scipy import sparse
import scipy.linalg as la
import scipy.signal as sg
import matplotlib
from matplotlib import pyplot as plt
from matplotlib import patches as pch
from matplotlib.colors import ListedColormap
from mpl_toolkits.mplot3d import Axes3D
import re
# to undo error demonstrtion in edapy_13_01
npsave = np;
pisave = pi;
# 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 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 [5]:
# edapy_13_01: exemplary errors related to redefining variables
pisave=pi;
npsave=np;
print('the variable pi: %.4f (before)' % pi)
pi=2;
print('the variable pi: %.4f (after)' % pi)
np = 3;
x = np.zeros((3,1));
the variable pi: 3.1416 (before) the variable pi: 2.0000 (after)
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[5], line 10 7 print('the variable pi: %.4f (after)' % pi) 9 np = 3; ---> 10 x = np.zeros((3,1)) AttributeError: 'int' object has no attribute 'zeros'
In [8]:
# edapy_13_02: time arithmetic
# restore numpy, which was destroyed in previous cell
pi=pisave;
np=npsave;
print('seconds DT between 2008/2/11 4:4:1 and 2008/2/11 4:4:0');
DELTA = datetime(2008,2,11,4,4,1)-datetime(2008,2,11,4,4,0);
print('DT: %.2f\n' % DELTA.total_seconds() );
DELTA = datetime(1975,1,1,0,0,1)-datetime(1974,12,31,23,59,59);
print('leap second checker');
print('seconds DT between 1975/1/1 0:0:1 and 1974/12/31 23:59:59');
print('right answer is 3, calculation gives DT: %.2f\n' % DELTA.total_seconds() );
seconds DT between 2008/2/11 4:4:1 and 2008/2/11 4:4:0 DT: 1.00 leap second checker seconds DT between 1975/1/1 0:0:1 and 1974/12/31 23:59:59 right answer is 3, calculation gives DT: 2.00
In [13]:
# edapy_13_03: reading idiosynchratic text file
# since size of file unknaon, allocate very big arrays
NBIG=1000000;
t = np.zeros((NBIG,1));
d = np.zeros((NBIG,1));
# names of months
moname = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'];
# open file
filename='../Data/brf_raw.txt';
fid = open(filename,"r");
N=0; # line number counter
while 1: # loop over lines in file
myline = fid.readline(); # read line
if not myline:
fid.close();
break;
# now process the line
myline = myline.replace("\t"," "); # tab to space
Ns = len(myline);
while 1: # replace multiple spaces with single space
myline = myline.replace(" "," ");
Nsnew = len(myline);
if( Nsnew==Ns ):
break;
Ns = Nsnew;
s = myline.split(); # split into words, result is list
s2 = s[0].split("-"); # spit first word into two
hrmn1 = s2[0]; # first hour-minute string
hrmn2 = s2[1]; # second hour-minute string
hr1 = int(hrmn1[0:2]); # first hour
hr2 = int(hrmn2[0:2]); # secong hour
mn1 = int(hrmn1[2:4]); # first minute
mn2 = int(hrmn2[2:4]); # second minute
sc1 = 0.0;
sc2 = 0.0;
# some end times listed as 24:00:00, change to 23:59:59
if( (hr2 == 24) and (mn2==0) ):
hr2 = 23;
mn2 = 59;
sc2 = 59;
da = int(s[1]);
yr = int(s[3]);
va = float(s[4])
# look-up month
mo=(-1);
for imo in range(12):
if( s[2]==moname[imo]):
mo=imo+1;
break;
# start time
t1=datetime(yr,mo,da,hr1,mn1,0);
# end time
t2=datetime(yr,mo,da,hr2,mn2,0);
if( N==0 ): # reference time is year of first sample
tref=datetime(yr,1,1,0,0,0);
DT1 = (t1-tref).total_seconds();
DT2 = (t2-tref).total_seconds();
t[N,0] = 0.5*(DT1+DT2)/86400; # time since reference time
d[N,0] = va; # data value
N=N+1; # increment line counter
# truncate to actual length
t = np.copy( t[0:N,0:1] );
d = np.copy( d[0:N,0:1] );
fig1 = plt.figure(1,figsize=(10,6));
plt.axis([0, np.max(t), np.min(d), np.max(d)]);
plt.plot( t, d, 'k-');
plt.xlabel('time (days');
plt.ylabel('temperature (degC)');
plt.title('Black Rock Forest');
In [14]:
# edapy_13_04, example of use of eda_draw() function
# define an examplary matrix, G
N=40;
v = eda_cvec(np.linspace(1,N,N));
G=np.matmul( v, v.T );
# define a exemplacy vector, m
m = eda_cvec(1.0*np.linspace(1,N,N));
# multiply
d = np.matmul( G, m);
# draw
eda_draw( 'title d', d, '=', 'title G', G, 'title m', m);
plt.show();
# add caption
print('Fig. The equation Gm=d, depicted using eda_draw()');
Fig. The equation Gm=d, depicted using eda_draw()
In [15]:
# edapy_13_05, simple examples of functions
# computes area, a, of circle of radius, r.
def areaofcircle(r):
a = pi*(r**2);
return a;
# computes circumference, c, and area, a, of
# a rectangle of length, l, and width, w
def CandAofRectangle(a,b):
circ = 2*(a+b);
area = a*b;
return circ, area;
# calculate area of a circle
radius=2;
area = areaofcircle(radius);
print('circle:');
print('radius: %.6f area: %.6f' % (radius,area) );
print(' ');
# calculate circumference and area of a rectangle
a=2;
b=4;
circ, area = CandAofRectangle(a,b);
print('rectangle;');
print('length: %.6f width: %.6f circumference: %.6f area: %.6f' % (a,b,circ,area) );
print(' ');
circle: radius: 2.000000 area: 12.566371 rectangle; length: 2.000000 width: 4.000000 circumference: 12.000000 area: 8.000000
In [16]:
# edapy_13_06, examples of reorganizing a matrix
# define a matrix
A = np.array( [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]] );
N,M = np.shape(A);
print('test matrix A');
print(A);
print(' ');
# create a column-vector, b, containing the elements of A,
# organized column-wise.
c = np.reshape(A,(N*M,1),order='F');
print('transpose of a column vector c formed from A');
print(c.T);
print(' ');
# create a row-vector, c, containing the elements of A, organized
# column-wise.
r = np.reshape(A,(1,N*M),order='F');
print('row vector formed from A');
print(r);
print(' ');
# the matrix reformed from the vector
B = np.reshape(c,(N,M),order='F');
print('the matrix reformed from the vector c');
print(B);
print(' ');
# when you 'unravel' an array with elements A[i,j] into
# column vector with elements b[k,0], you sometimes need
# to determine the correspondence beteen (i,j) and k.
# and vice versa:
i=1;
j=0;
print('for a (%d,%d) matrix M[i,j] unraveled into a vector c[k,0]' % (N,M) );
k = np.ravel_multi_index( (i,j), (N,M), order='F') ;
print('(i,j) of (%d,%d) corresponds to k of %d' % (i,j,k) );
print('A(i,j)=%.0f and c[k]=%.0f' % (A[i,j], c[k,0]) );
ii, jj = np.unravel_index( k, (N,M), order='F' ) ;
print('k of %d corresponds to (i,j) of (%d,%d)' % (k,i,j) );
test matrix A [[ 1 2 3 4] [ 5 6 7 8] [ 9 10 11 12]] transpose of a column vector c formed from A [[ 1 5 9 2 6 10 3 7 11 4 8 12]] row vector formed from A [[ 1 5 9 2 6 10 3 7 11 4 8 12]] the matrix reformed from the vector c [[ 1 2 3 4] [ 5 6 7 8] [ 9 10 11 12]] for a (3,4) matrix M[i,j] unraveled into a vector c[k,0] (i,j) of (1,0) corresponds to k of 1 A(i,j)=5 and c[k]=5 k of 1 corresponds to (i,j) of (1,0)
In [17]:
# edapy_13_07, function Phi used in Lagrange multiplier figure
# make an (x,y) grid
N=51;
M=51;
Dx=0.01;
Dy=0.01;
x = eda_cvec(Dx*np.linspace(0,N-1,N));
y = eda_cvec(Dx*np.linspace(0,M-1,M));
# a long-wavelength, 2D sinusoid
Phi= np.multiply( np.sin(x), (np.sin(y)).T );
eda_draw(Phi);
In [18]:
# edapy_13_08, supports Figure 13.3
# covariance corresponding to 2nd Derivative smoothness
# constants
s = 0.1; # scale factor
g2 = 1.0; # amplitude
# stamdard setup of x-axis
N2 = 100;
N = 2*N2+1;
Dx = 1.0;
x = eda_cvec(Dx*np.linspace(-N2,N2,N));
xmin = np.min(x);
xmax = np.max(x);
# autocovariance, with lag |x| = |xi - xj|
c = g2*np.multiply((1.0+s*abs(x)),np.exp(-s*np.abs(x)));
# Gussian autocovariance for comparison
cg = g2*np.exp(-0.5*s*s*np.power(x,2));
# plot c(x) vs x
fig1 = plt.figure(1,figsize=(8,4));
ax1 = plt.subplot(1,1,1);
plt.axis([xmin, xmax, -0.1, 1.1])
plt.plot(x,c,'k-',lw=3);
plt.plot(x,cg,'k:',lw=2);
plt.plot(x,0.0*c,'k-',lw=1);
plt.xlabel('x');
plt.ylabel('c(x)');
plt.show();
In [19]:
# edapy_13_09, computing ln det Cm
# model parameter axis, t, with sampling Dt
M=64; # number of model parameters
Dt=1.0; # time sampling
tmax=Dt*(M-1); # maximum time
t=eda_cvec(np.linspace(0,tmax,M)); # time
print('M: %d' % (M) );
# covariance matrix, of cos(w0t) exp(-st) form
n=6; # number of oscillations on (0,tmax) interal
w0 = n*2*pi/tmax; # angular frequency of oscilation
tmid=tmax/2.0;
s = Dt/60.0; # decay rate
gamma = 1.0;
gamma2 = gamma**2; # variance
# build covariances via outer product
Tm = np.abs(np.matmul(t,np.ones((1,M)))-np.matmul(np.ones((M,1)),(t.T)));
Cm = gamma2*np.multiply(np.cos(w0*Tm),np.exp(-s*Tm));
eda_draw('title Cm',Cm);
# Brute force method
D = la.det(Cm)
logdetCm_A = log(D);
print('Brute force method, log det Cm = %e' % (logdetCm_A) );
# Eigenvalue method
LAMBDA, V = la.eig(Cm);
LAMBDA = np.real(LAMBDA);
logdetCm_B = np.sum(np.log(LAMBDA));
print('Eigenvalue method, log det Cm = %e' % (logdetCm_B) );
# Cholesky method
R = la.cholesky(Cm);
logdetCm_C = 2.0*np.sum(np.log(np.diag(R)));
print('Cholesky method, log det Cm = %e' % (logdetCm_C) );
print(' ');
print("brute force works fine for M=%d, but fails for M=230 and above" %(M));
M: 64
Brute force method, log det Cm = -1.860150e+02 Eigenvalue method, log det Cm = -1.860150e+02 Cholesky method, log det Cm = -1.860150e+02 brute force works fine for M=64, but fails for M=230 and above
In [ ]: