In [1]:
# edapy_01_00 clears all variables and import various modules
# you must run this cell first, before running any of others
# History:
# 2021/08/25 - Edited by W.Menke, synchronizing Python/MATLAB versions
# 2021/06/19 - Edited by W. Menke, while writing corresponding chapter
# 2020/07/21 - Created by W.Menke, by translation of Matlab version
# 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, sqrt # math functions
import scipy.linalg as la # linear algebra functions
import os # operating system
# 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):
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):
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):
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 [2]:
# edapy_01_01 print the date
thedate = date.today();
print(thedate);
In [4]:
# edapy_01_02: directory manipulation
# get current directory
edapy_dir=os.getcwd();
print(edapy_dir);
# get current folder, with is presumably edapy/Notebooks
edapy_dir=os.getcwd();
# change to the TestFolder, up-and-over from Notebooks
os.chdir("../TestFolder");
# print current working directory (folder)
thisdir=os.getcwd();
print(thisdir);
# change back to original folder
os.chdir(edapy_dir);
# change to the TestFolder, which is up-and-over
os.chdir("../TestFolder");
# print current working directory (folder)
thisdir=os.getcwd();
print(thisdir);
# change to folder one level up
os.chdir("..");
# change to the Notebooks folder
os.chdir("Notebooks");
# change to TestFolder (because it only contains two files)
os.chdir("../TestFolder");
mydirlist=os.listdir();
print(mydirlist);
os.chdir("../Notebooks");
In [9]:
# edapy_01_03: simple algebra
a = 3.5;
b = 4.1;
c = a+b;
print(c);
In [5]:
# edapy_01_04: algebra with square root and power
a = 3.0;
b = 4.0;
c = sqrt(a**2 + b**2);
print(c);
In [6]:
# edapy_01_05: algebra with sine
n = 2; x = 3; x0 = 1; L = 5;
c = sin(n*pi*(x-x0)/L);
print(c);
In [8]:
# edapy_01_06 list, tupples, and other list-like objects
print("List");
mylist = [1, 2, 4, 8];
N=len(mylist);
print("mylist", mylist);
print("mylist[0]", mylist[0]);
print("len(mylist)",N);
print("Tuple");
mytuple = (1, 2, 4, 8);
N=len(mytuple);
print("mytuple", mytuple);
print("mytuple[0]", mytuple[0]);
print("len(mytuple)",N);
v4x0 = np.zeros((4));
print("4x0 np.array");
print(v4x0);
v1x4 = np.zeros((1,4));
print("1x4 np.array");
print(v1x4);
v4x1 = np.zeros((4,1));
print("4x1 np.array");
print(v4x1);
In [13]:
# edapy_01_07, creating vectors and matrices
# a length-4 column-vector c of zeros
c0 = np.zeros((4,1));
print('c0\n',c0);
# a length-4 column-vector c of ones
c1 = np.ones((4,1));
print('c1\n',c1);
# column vector c
# note that decimal points are needed,
# else you will get an integer vector
# which is usually not what you want
c = eda_cvec( [1.,2.,3.,4.] );
print('c\n',c);
# row vector r
c = eda_cvec([1.,2.,3.,4.]);
r=c.T;
print('r\n',r);
# 3x3 square matrix A0 of zeros, and
# 3x3 square matrix A1 of ones
A0 = np.zeros((3,3));
print("A0\n",A0);
A1 = np.ones((3,3));
print("A1\n",A1);
# square matrix A
# note that decimal points are needed,
# else you will get an integer vector
# which is usually not what you want
A = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("A\n",A);
# find shape of c
Nc, Mc = np.shape(c);
print('shape of c:\n',Nc,Mc);
# find shape of c
NA, MA = np.shape(A);
print('shape of M\n',NA,MA);
# note that all these calls to eda_cvec
# achieve the same result
c = eda_cvec( [1.,2.,3.,4.] );
print('c\n',c);
c = eda_cvec( 1.,2.,3.,4. );
print('c\n',c);
c = eda_cvec( (1.,2.,3.,4.) );
print('c\n',c);
c = eda_cvec( 1.,[2.,3.,4.] );
print('c\n',c);
c = eda_cvec( [1.,2.],[3.,4.] );
print('c\n',c);
In [14]:
# edapy_01_08: defining vectors and matrices
# column vector c
c = eda_cvec( [1.,3.,5.] );
print("c\n",c);
# column vector r
r = eda_cvec( [2.,4.,6.] ).T;
print("r\n",r);
# square matrix A
A = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("A\n",A)
In [9]:
# edapy_01_09: vector and matrix addition
# column vector c
a = eda_cvec( [1.,3.,5.] );
print("a\n",a);
# column vector b
b = eda_cvec( [2.,4.,6.] );
print("b\n",b);
c = a+b;
print("c\n",c);
# square matrix A
A = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("A\n",A);
# square matrix B
B = np.ones((3,3));
print("B\n",B);
C = A+B;
print("C\n",C);
In [2]:
# edapy_01_10: vector and matrix multiplication
# column vector c
a = eda_cvec( [1.,3.,5.] );
print("a\n",a);
# column vector r
b = eda_cvec( [2.,4.,6.] );
print("b\n",b);
# square matrix A
A = np.array( [[1., 0., 2.],
[0., 1., 0.],
[2., 0., 1.]] );
print("A\n",A);
# square matrix B
B = np.array( [[1., 0., -1.],
[0., 2., 0.],
[-1., 0., 3.]] );
print("B\n",B);
# inner product: s = a-transpose time b
s = np.matmul(a.T,b);
print("s\n",s);
# outer product: T = a times b-transpose
T = np.matmul(a,b.T);
print("T\n",T);
# matrix times a vector: P = A times B
c = np.matmul(A,a);
print("c\n",c);
# matrix times a matrix: P = A times B
P = np.matmul(A,B);
print("P\n",P);
# d = element-wise multiply of a and b
d = np.multiply(a,b);
print("d\n",d);
In [17]:
# edapy_01_11: accessing elements of a matrix
# column vector a
a = eda_cvec([1.,2.,3.]);
print("a\n",a);
# square matrix M
A = np.zeros((3,3));
A = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("M\n",A);
# one element of a column vector
s = a[1,0];
print("s\n",s);
# one element of a matrix
t = A[1,2];
print("t\n",t);
# column vector = one column of a matrix
b = np.copy( A[0:3,1:2] );
print("b\n",b);
# row vector = one row of a matrix
r = np.copy( A[1:2,0:3] );
print("r\n",r);
# row vector = one row of a matrix
c = eda_cvec( np.copy( A[1:2,0:3] ) );
print("c\n",c);
# small matrix = rectangular block of larger matrixx
# note that end of range is one larger than you want
T = np.copy( A[1:3,1:3] );
print("T\n",T);
# the purpose for the copy is to convert a view of a matrix
# into a separate matrix. Else changeing the original matrix
# will change the
X = A[1:2+1,1:2+1];
print("X\n",X);
A[1,1]=7;
print("X\n",X);
print("T\n",T);
# set first column of matrix A to a column vector b
b = eda_cvec([10.,11.,12.]);
A = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("A\n",A);
B = np.copy(A);
# the shape of the view on the left-hand-side of the equal sign
# must exactly match the shape of the quantity on the right hand side
B[0:3,0:1] = b;
print("B\n",B);
# ravel a N by M matrix A into a NM by 0 numpy array
Alist = A.ravel();
print("Alist",Alist);
In [32]:
# edapy_01_12: derivatives and integrals
# column vector t of N evenly-spaced points beterrn tmin and tmax
N=21; # length of t
tmin = 0.0; # minimum time
tmax = 1.0; # maximum time
Dt = (tmax-tmin)/(N-1); # time increment
t = eda_cvec( np.linspace( tmin, tmax, N ) );
# expemplary function
x=np.sin(pi*t);
# approximate slope; note one point shorter then x
# sapprox = eda_cvec( (x[1:N+1]-x[0:N-2+1])/Dt );
s = (1/Dt) * ( x[1:N,0:1]-x[0:N-1,0:1] );
Ns, M = np.shape(s);
print('s\n',s);
print('length\n',Ns);
# approximate integral; same length as x
# integral is approximated as a cumulative sum
a = eda_cvec( Dt*np.cumsum(x) );
print('a');
print(a);
# integral over all t; integral approximated as a sum
atotal = Dt*np.sum(x);
print('atotal\n', a[N-1] );
print('atotal\n', atotal );
# the rest of this section makes plots used in the text,
# but since plotting has not yet been discussed, readers
# may prefer to skip it
# true derivative and integral
strue = pi*np.cos(pi*t); # derivative
atrue = 1/pi-np.cos(pi*t)/pi; # integral, with integration constant
# chosen so integral is zero at t=0
# time series plot
fig1 = plt.figure(1,figsize=(8,4));
ax1 = plt.subplot(1, 1, 1);
plt.axis([0.0, 1.0, 0, 1.1]);
plt.plot(t,x,"-",lw=4,color=[0.8,0.8,0.8]);
plt.plot(t,x,"ko",lw=2);
plt.xlabel("t");
plt.ylabel("x(t)");
for i in range(N-1):
plt.plot( eda_cvec(t[i,0],t[i]), eda_cvec(0.0,0.1), 'k-', 'LineWidth', 2 );
plt.show();
# calculus plot
fig2 = plt.figure(2,figsize=(8,6));
ax1 = plt.subplot(3, 1, 1);
plt.axis([0.0, 1.0, 0, 1.1]);
plt.xlabel("t");
plt.ylabel("x(t)");
# plot exemplary rectangles, shaded
K=5;
for k in range(K):
if(x[k,0]>=0.0):
r=patches.Rectangle( (t[k,0],0.0), Dt, x[k,0], ec=[0.0,0.0,0.0], fc=[0.9,0.9,0.9] );
ax1.add_patch(r);
# plot function as continuous curve
plt.plot( t, x, lw=4, color=[0.9,0.9,0.9] );
# plot function as circles
plt.plot( t, x, 'ko', lw=2 );
# embolden one slope
plt.plot( eda_cvec(t[K,0],t[K+1,0]), eda_cvec(x[K,0],x[K+1,0]), '-', lw=6, color=[0,0,0]);
# plot trapezoids
for i in range(N-1):
tb = eda_cvec( t[i,0],t[i,0],t[i+1,0],t[i+1,0] );
xb = eda_cvec(0.0,x[i,0],x[i+1,0],0.0 );
plt.plot( tb, xb, 'k-', lw=2 );
# plot derivative
ax2 = plt.subplot(3, 1, 2);
plt.axis( [0, 1, -1.1*pi, 1.1*pi ] );
plt.xlabel("t");
plt.ylabel("s=dx/dx");
plt.plot( t, strue, '-', lw=4, Color=[0.8,0.8,0.8] );
plt.plot( t[0:N-1,0:1], s, 'k-', lw=2 );
# plot integral
ax3 = plt.subplot(3, 1, 3);
plt.axis( [0.0, 1.0, 0.0, 1.1*2/pi ] );
plt.xlabel("t");
plt.ylabel("integral of x");
plt.plot( t, atrue, 'k-', lw=4, color=[0.8,0.8,0.8] );
plt.plot( t, a, 'k-', lw=2 );
plt.plot( t[N-1,0], atotal, 'ko', lw=3 ); # last value
In [10]:
# edapy_01_13: for-loops
# square matrix A
A = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("A\n",A);
# for loop that extracts the main diagonal of a matrix
a = np.zeros((3,1));
for i in range(3):
a[i,0]=A[i,i];
print("a\n",a)
# two nested for loops used to rearrange (flip left-right)
# the elements of a matrix
B = np.zeros((3,3));
for i in range(3): # loop over all rows
for j in range(3): # loop over all cols
B[i,2-j]=A[i,j]; # left-right-flip
print("B\n",B);
# column vector of length 12
a = eda_cvec( [1.0, 2.0, 1.0, 4.0, 3.0, 2.0, 6.0, 4.0, 9.0, 2.0, 1.0, 4.0] );
N, M = np.shape(a);
# for loop that "clips" a vector
b = np.zeros((N,1));
for i in range(N):
if a[i,0] >= 6.:
b[i,0]=6.;
else:
b[i,0]=a[i,0];
# side-by-side concatenation vectors a and b for aeshetic printing
c = np.concatenate((a,b),axis=1);
print("c\n",c);
In [20]:
# edapy_01_14: alternatives to for-loops
# square matric A
A = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("A\n",A);
# extract diagonal of matcix into a vector
a = eda_cvec( np.diagonal(A) );
print("a\n",a);
# B is A flipped left to right
B = np.fliplr(A);
print("N\n",B);
# column vector a of lenth 12
a = eda_cvec( [1.0, 2.0, 10.0, 4.0, 3.0, 2.0, 6.0, 4.0, 9.0, 2.0, 1.0, 4.0] );
print("a\n",a);
# alternative method #1 of cliping a vector
b = np.multiply(a, a<6.) + np.multiply(6., a>=6.);
print("b\n",b);
# alternative method #2 of cliping a vector
c=np.copy(a); # start with c that is a copy of a
c[a>=6.]=6.;
print("c\n",c);
# alternative method #3 f cliping a vector
# we use the np.where() command to finds the
# rows of a that satisfy the condition
d=np.copy(a); # start with a d that is a copy of a
k = np.where(a>=6.); # k is a data structure descibed below
d[k]=6.; # reset the rows of d that violate the condition
# k is a tupple of length 2
print("type k\n",type(k));
print("len k\n",len(k));
# k[0] gives the rows satisfying the a>=6 condition. It is a 3-by-0 np.array
print("type k[0]\n",type(k[0]));
print("shape k[0]\n",np.shape(k[0]));
print("k[0]\n",k[0]);
print("d\n",d);
In [13]:
# edapy_01_15: matrix inverse
# square matrix A
A = np.array( [[1.0, 0.5, 0.0],
[0.5, 5.0, 0.0],
[0.0, 0.0, 9.0]] );
print("A:\n",A);
# B = inverse of A
B = np.linalg.inv(A);
print("B:\n",B);
# A times B
I1 = np.matmul(A,B);
print("AB:\n",I1);
# B times A
I2 = np.matmul(B,A);
print("BA:\n",I2);
In [17]:
# edapy_01_16: solving linear system
# square matrix A
A = np.array([[1.0, 5.0, 13.0],
[2.0, 7.0, 17.0],
[3.0, 11.0, 19.0]] );
# column vercor b
b = eda_cvec( [2., 4., 6.] );
print("b\n",b);
# when Ac=b
# alternative to c = A-inverse times b
c = np.linalg.solve(A,b);
print("c\n",c);
# check result is correct
error1=b-np.matmul(A,c);
print("error\n",error1);
B = np.array( [[1.0, 2.0, 3.0],
[4.0, 5.0, 6.0],
[7.0, 8.0, 9.0]] );
# when AD=B
# alternative to D = A-inverse B
D = np.linalg.solve(A,B);
print("D\n",D);
# check result is correct
error2=B-np.matmul(A,D);
print("error\n",error2);
In [20]:
# edapy_01_17: load text file and make very-quick plot
# this example uses the neuse.txt stream flow dataset
# load rectangular, tab-separated table of numbers into a N by M matrix D
D = np.genfromtxt("../Data/neuse.txt", delimiter="\t");
N, M = np.shape(D);
t = np.copy( D[0:N,0:1] ); # time in column 1
d = np.copy( D[0:N,1:2] ); # discharge in column 2
# plot the data
fig1 = plt.figure(1,figsize=(8,4)); # smallish figure
ax1 = plt.subplot(1, 1, 1); # only one subplot
plt.plot(t,d,"k-"); # plot d(t) with black line
plt.xlabel("time in days"); # label time axis
plt.ylabel("discharge in cfs"); # label data axis
plt.show(); # display plot
In [21]:
# edapy_01_18: load text file and mmke nicer plot
# this example uses the neuse.txt stream flow dataset
# load data from text fole
D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
[N, M]=D.shape;
# time in column 1, discharge in colun 2
t = eda_cvec( D[:,0] );
d = eda_cvec( D[:,1] );
# nicer plot
fig1 = plt.figure(1,figsize=(8,4));
ax1 = plt.subplot(1, 1, 1);
plt.axis([0, 4500, 0, 1.1*np.max(d)])
plt.plot(t,d,"k-");
plt.xlabel("time (days)");
plt.ylabel("discharge (cfs)");
plt.title("Neuse River Hydrograph");
plt.show();
In [22]:
# edapy_01_19: log-log plot
# get 10-year list of earthquake magnitudes from file
# note: since the file contains only one column of data, the
# np.genfromtxt puts the data into an N by 0 array. So we
# call eda_cvec to put it into a standard columnn vector
d = eda_cvec( np.genfromtxt('../Data/earthquake_magnitudes.txt', delimiter='\t') );
time = 10.0; # 10 years of data
# count up earthquakes with magnitudes in set of magnitude bins
# 51 bins from magnitude 5 to magnitude 10
c, e = np.histogram(d,51,(5,10));
Nc = len(c); # counts of earthquakes in each bin
Ne = len(e); # magnitudes at the edges of the bins
count = eda_cvec(c);
edges = eda_cvec(e);
# convert count to rate
r = eda_cvec( count[0:Nc-1,0]/time );
# compute center point of magnitude bins
m = eda_cvec( (edges[0:Ne-2,0]+edges[1:Ne-1,0])/2 );
# convert magntude to energy using the Gutenberg-Richter magnitude-energy relation
log10Ee = 1.5*m + 11.8; # in ergs
log10Ej = log10Ee - 7; # in joules
E = 1e-7 * np.power(10,log10Ej);
# log-log plot
fig2 = plt.figure(2,figsize=(8,4));
ax1 = plt.subplot(1, 1, 1);
plt.loglog(E, r, "k-", basex=10.0, basey=10.0); # plot data as lines
plt.loglog(E, r, "ko", basex=10.0, basey=10.0); # plot data as circles
plt.axis([10**5, 10**12, 10**-1, 10**3])
plt.xlabel("earthquake energy (joules)");
plt.ylabel("earthquakes rate (per year)");
plt.title("Earthquake rate vs energy");
plt.show();
In [25]:
# edapy_01_20: write matrix to file
# load the neuse.txt data set
# convert discharge from cfs to m3/s
# and write to a file neuse_metric.txt
# get matrix from file
D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
N, M = np.shape(D);
# time in column 0
t = D[0:N,0:1];
# dischange in cfs in column 1
d = D[0:N,1:2];
# conversion from cfs to m3/s
f = 35.3146;
dm = d/f;
# make simple plot to veryify result "by eye"
# plot the data
fig1 = plt.figure(1,figsize=(8,4)); # smallish figure
ax1 = plt.subplot(1, 1, 1); # only one subplot
plt.plot(t,dm,"k-"); # plot d(t) with black line
plt.xlabel("time in days"); # label time axis
plt.ylabel("discharge in m3/s"); # label data axis
plt.show(); # display plot
# concatenate time, discharge into a matrix C
C = np.concatenate( (t,dm),axis=1);
# save matrix to file of tab-separated values
np.savetxt("../Data/neuse_metric.txt",C, delimiter="\t");
In [29]:
# edapy_01_21: big-picture commenting (recommended)
# This script evaluates the 2D Normal p.d.f. p(d1,d2) with
# mean (d1bar,d2bar) and covariance C on a L-by-L grid
# grid size
L=101;
# d1 axis
d1min = 0.0;
d1max = 10.0;
d1 = eda_cvec( np.linspace( d1min, d1max, L ) );
# d2 axis
d2min = 0.0;
d2max = 10.0;
d2 = eda_cvec( np.linspace( d2min, d2max, L ) );
# mean
d1bar = 4.0;
d2bar = 6.0;
# covariance matrix
C = np.array( [[1.0, 0.5],
[0.5, 2.0]] );
# Evaluate Normal p.d.f. Pp
# with mean dbar and covariance C
dd = np.zeros((2,1));
CI = np.linalg.inv(C);
detC = np.linalg.det(C);
norm = 1/(2*pi*sqrt(detC));
Pp = np.zeros((L,L));
for i in range(L): # loop over d1
for j in range(L): # loop over d2
dd[0,0] = d1[i,0]-d1bar;
dd[1,0] = d2[j,0]-d2bar;
E = np.matmul(dd.T,np.matmul(CI,dd));
# Normal p.d.f.
Pp[i,j] = norm*exp(-0.5* E);
fig1 = plt.figure(1,figsize=(8,4)); # figure 1
ax1 = plt.subplot(1, 1, 1); # single subplot
plt.title("normal p.d.f."); # title
plt.imshow(Pp); # plot Pp
plt.show(); # show image
In [30]:
# edapy_01_22 Comments that focus on minutia (not recommended)
# create and plot the array p(d1,d2)
# define L
L=101;
# create d1
d1min = 0.0;
d1max = 10.0;
d1 = eda_cvec( np.linspace( d1min, d1max,L ) );
# create d2
d2min = 0.0;
d2max = 10.0;
d2 = eda_cvec( np.linspace( d2min, d2max, L ) );
# define d1bar and d2bar
d1bar = 4.;
d2bar = 6.;
# define C
C = np.zeros((2,2));
C = np.array( [[1.0, 0.5],
[0.5, 2.0]] );
# define dd
dd = np.zeros((2,1));
# take inverse
CI = np.linalg.inv(C);
# take determinant
detC = np.linalg.det(C);
# define normalization
norm = 1/(2*pi*sqrt(detC));
Pp = np.zeros((L,L)); # initialize Pp
for i in range(L): # loop over i
for j in range(L): # loop over j
# det up dd
dd[:,0] = [d1[i,0]-d1bar, d2[j,0]-d2bar];
# compute E
E = np.matmul(dd.T,np.matmul(CI,dd));
# peform exponential
Pp[i,j] = norm*exp(-0.5 * E );
fig1 = plt.figure(1,figsize=(8,4)); # figure 1
ax1 = plt.subplot(1, 1, 1); # single subplot
plt.title("normal p.d.f."); # title
plt.imshow(Pp); # plot Pp
plt.show(); # show image