In [2]:
# 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
# 2021/12/30 -  review, fixed some issues including np.where()
# 2022/06/26 -  checked with most recent Python & etc. software
# 2023/07/13 -  fixed some issues, incl colormap, loglof plt, eda_cvec()
# 2024/05/23 -  Patches to accomodate new verions of numpy, scipy, matplotlib
#               patched dot product to return scalar value
 
# A reset deletes any previously-created variables from memory
%reset -f

# import various packages

import numpy as np                    # matrices & etc
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
import matplotlib
import matplotlib.cm as cm
from matplotlib import pyplot as plt  # general plotting 
from matplotlib import patches        # plotting shapes

# 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_01_01 print the date
thedate = date.today();
print(thedate);
2024-05-25
In [5]:
# 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");
c:\bill\EDAWM-3rdEd\PythonPatch20240524\edapy\Notebooks
c:\bill\EDAWM-3rdEd\PythonPatch20240524\edapy\TestFolder
c:\bill\EDAWM-3rdEd\PythonPatch20240524\edapy\TestFolder
['TextDocument1.txt', 'TextDocument2.txt']
In [6]:
# edapy_01_03: simple algebra
a = 3.5;
b = 4.1;
c = a+b;
print(c);
7.6
In [7]:
# edapy_01_04: algebra with square root and power
a = 3.0;
b = 4.0;
c = sqrt(a**2 + b**2);
print(c);
5.0
In [8]:
# edapy_01_05: algebra with sine
n = 2; x = 3; x0 = 1; L = 5;
c = sin(n*pi*(x-x0)/L);
print(c);
0.5877852522924732
In [9]:
# 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);
List
mylist [1, 2, 4, 8]
mylist[0] 1
len(mylist) 4
Tuple
mytuple (1, 2, 4, 8)
mytuple[0] 1
len(mytuple) 4
4x0 np.array
[0. 0. 0. 0.]
1x4 np.array
[[0. 0. 0. 0.]]
4x1 np.array
[[0.]
 [0.]
 [0.]
 [0.]]
In [10]:
# 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);
c0
 [[0.]
 [0.]
 [0.]
 [0.]]
c1
 [[1.]
 [1.]
 [1.]
 [1.]]
c
 [[1.]
 [2.]
 [3.]
 [4.]]
r
 [[1. 2. 3. 4.]]
A0
 [[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
A1
 [[1. 1. 1.]
 [1. 1. 1.]
 [1. 1. 1.]]
A
 [[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]
shape of c:
 4 1
shape of M
 3 3
c
 [[1.]
 [2.]
 [3.]
 [4.]]
c
 [[1.]
 [2.]
 [3.]
 [4.]]
c
 [[1.]
 [2.]
 [3.]
 [4.]]
c
 [[1.]
 [2.]
 [3.]
 [4.]]
c
 [[1.]
 [2.]
 [3.]
 [4.]]
In [11]:
# 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)
c
 [[1.]
 [3.]
 [5.]]
r
 [[2. 4. 6.]]
A
 [[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]
In [12]:
# 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);
a
 [[1.]
 [3.]
 [5.]]
b
 [[2.]
 [4.]
 [6.]]
c
 [[ 3.]
 [ 7.]
 [11.]]
A
 [[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]
B
 [[1. 1. 1.]
 [1. 1. 1.]
 [1. 1. 1.]]
C
 [[ 2.  3.  4.]
 [ 5.  6.  7.]
 [ 8.  9. 10.]]
In [13]:
# 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
# 2024/05/2024 the new verion of numpy depreciated using 1x1 arrays
# as scalars (which worked fine in the past), so I have added s=s[0,0]
# which converts s to a scalar.
s = np.matmul(a.T,b); s=s[0,0];
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);
a
 [[1.]
 [3.]
 [5.]]
b
 [[2.]
 [4.]
 [6.]]
A
 [[1. 0. 2.]
 [0. 1. 0.]
 [2. 0. 1.]]
B
 [[ 1.  0. -1.]
 [ 0.  2.  0.]
 [-1.  0.  3.]]
s
 44.0
T
 [[ 2.  4.  6.]
 [ 6. 12. 18.]
 [10. 20. 30.]]
c
 [[11.]
 [ 3.]
 [ 7.]]
P
 [[-1.  0.  5.]
 [ 0.  2.  0.]
 [ 1.  0.  1.]]
d
 [[ 2.]
 [12.]
 [30.]]
In [15]:
# edapy_01_11: accessing elements of a matrix

# column vector a
a = eda_cvec([1.,2.,3.]);
print("a\n",a);

# square matrix M
B = np.array( [[1., 2., 3.],
               [4., 5., 6.],
               [7., 8., 9.]] );
print("B\n",B);

# one element of a column vector
s = np.copy(a[1,0]);
print("s\n",s);

# one element of a matrix
t = np.copy(B[1,2]);
print("t\n",t);

# column vector = one column of a matrix
b = np.copy( B[0:3,1:2] );
print("b\n",b);

# row vector = one row of a matrix
r = np.copy( B[1:2,0:3] );
print("r\n",r);

# column vector = one row of a matrix
c = np.copy( B[1:2,0:3].T );
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( B[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 =  B[1:2+1,1:2+1];
print("X\n",X);
B[1,1]=7;
print("X\n",X);
print("T\n",T);

# set first column of matrix B to a column vector b
b = eda_cvec([10.,11.,12.]);
B= np.array( [[1., 2., 3.],
               [4., 5., 6.],
               [7., 8., 9.]] );
print("B\n",B);
Bp = np.copy(B);
# 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
Bp[0:3,0:1] = b;
print("Bp\n",Bp);

# ravel a N by M matrix A into a NM by 0 numpy array
Blist = B.ravel();
print("B.ravel()",Blist);
a
 [[1.]
 [2.]
 [3.]]
B
 [[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]
s
 2.0
t
 6.0
b
 [[2.]
 [5.]
 [8.]]
r
 [[4. 5. 6.]]
c
 [[4.]
 [5.]
 [6.]]
T
 [[5. 6.]
 [8. 9.]]
X
 [[5. 6.]
 [8. 9.]]
X
 [[7. 6.]
 [8. 9.]]
T
 [[5. 6.]
 [8. 9.]]
B
 [[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]
Bp
 [[10.  2.  3.]
 [11.  5.  6.]
 [12.  8.  9.]]
B.ravel() [1. 2. 3. 4. 5. 6. 7. 8. 9.]
In [16]:
# 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,3));
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-', lw=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, '-', 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
s
 [[ 3.1286893 ]
 [ 3.05165059]
 [ 2.89947011]
 [ 2.67589505]
 [ 2.38643058]
 [ 2.03820426]
 [ 1.6397906 ]
 [ 1.20099984]
 [ 0.73263649]
 [ 0.24623319]
 [-0.24623319]
 [-0.73263649]
 [-1.20099984]
 [-1.6397906 ]
 [-2.03820426]
 [-2.38643058]
 [-2.67589505]
 [-2.89947011]
 [-3.05165059]
 [-3.1286893 ]]
length
 20
a
[[0.        ]
 [0.00782172]
 [0.02327257]
 [0.0459721 ]
 [0.07536136]
 [0.1107167 ]
 [0.15116755]
 [0.19571788]
 [0.2432707 ]
 [0.29265512]
 [0.34265512]
 [0.39203954]
 [0.43959236]
 [0.48414269]
 [0.52459354]
 [0.55994888]
 [0.58933814]
 [0.61203766]
 [0.62748851]
 [0.63531024]
 [0.63531024]]
atotal
 [0.63531024]
atotal
 0.6353102368087353
No description has been provided for this image
No description has been provided for this image
In [17]:
# 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);
A
 [[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]
a
 [[1.]
 [5.]
 [9.]]
B
 [[3. 2. 1.]
 [6. 5. 4.]
 [9. 8. 7.]]
c
 [[1. 1.]
 [2. 2.]
 [1. 1.]
 [4. 4.]
 [3. 3.]
 [2. 2.]
 [6. 6.]
 [4. 4.]
 [9. 6.]
 [2. 2.]
 [1. 1.]
 [4. 4.]]
In [18]:
# 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
r = np.where(a>=6.); # r is a data structure descibed below
k = r[0];            # k is a Nx0 list of indices
d[k,0:1]=6.;         # reset the rows of d that violate the condition

# r is a tupple
print("type r\n",type(r));
print("shape r\n",np.shape(r));
# k=r[0] gives the rows satisfying the a>=6 condition. It is a 3-by-0 np.array
print("type k=r[0]\n",type(k));
print("shape k\n",np.shape(k));
print("k\n",k);
print(" ");
print("d\n",d);
A
 [[1. 2. 3.]
 [4. 5. 6.]
 [7. 8. 9.]]
a
 [[1.]
 [5.]
 [9.]]
N
 [[3. 2. 1.]
 [6. 5. 4.]
 [9. 8. 7.]]
a
 [[ 1.]
 [ 2.]
 [10.]
 [ 4.]
 [ 3.]
 [ 2.]
 [ 6.]
 [ 4.]
 [ 9.]
 [ 2.]
 [ 1.]
 [ 4.]]
b
 [[1.]
 [2.]
 [6.]
 [4.]
 [3.]
 [2.]
 [6.]
 [4.]
 [6.]
 [2.]
 [1.]
 [4.]]
c
 [[1.]
 [2.]
 [6.]
 [4.]
 [3.]
 [2.]
 [6.]
 [4.]
 [6.]
 [2.]
 [1.]
 [4.]]
type r
 <class 'tuple'>
shape r
 (2, 3)
type k=r[0]
 <class 'numpy.ndarray'>
shape k
 (3,)
k
 [2 6 8]
 
d
 [[1.]
 [2.]
 [6.]
 [4.]
 [3.]
 [2.]
 [6.]
 [4.]
 [6.]
 [2.]
 [1.]
 [4.]]
In [19]:
# 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);
A:
 [[1.  0.5 0. ]
 [0.5 5.  0. ]
 [0.  0.  9. ]]
B:
 [[ 1.05263158 -0.10526316  0.        ]
 [-0.10526316  0.21052632  0.        ]
 [ 0.          0.          0.11111111]]
AB:
 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
BA:
 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
In [20]:
# 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);
b
 [[2.]
 [4.]
 [6.]]
c
 [[ 2.00000000e+00]
 [-1.24900090e-16]
 [ 4.16333634e-17]]
error
 [[-4.4408921e-16]
 [-8.8817842e-16]
 [-8.8817842e-16]]
D
 [[ 4.00000000e+00  3.50000000e+00  3.00000000e+00]
 [-1.66666667e-01 -8.33333333e-02 -1.87350135e-16]
 [-1.66666667e-01 -8.33333333e-02  6.24500451e-17]]
error
 [[-6.66133815e-16  0.00000000e+00  0.00000000e+00]
 [-1.77635684e-15 -8.88178420e-16  0.00000000e+00]
 [-2.66453526e-15  0.00000000e+00  0.00000000e+00]]
In [21]:
# 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 (days)");           # label time axis
plt.ylabel("discharge (cfs)");       # label data axis
plt.show();                           # display plot
No description has been provided for this image
In [22]:
# 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();
No description has been provided for this image
In [23]:
# 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-"); # plot data as lines
plt.loglog(E, r, "ko"); # 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();
No description has been provided for this image
In [24]:
# 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");
No description has been provided for this image
In [26]:
# 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)); E=E[0,0];
        # 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
mycmap=matplotlib.colormaps['jet']; # colormap
plt.imshow(Pp,cmap=mycmap);         # plot Pp
plt.show();                         # show image
No description has been provided for this image
In [28]:
# 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)); E=E[0,0];
        # 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
mycmap=matplotlib.colormaps['jet']; # colormap
plt.imshow(Pp,cmap=mycmap);         # plot Pp
plt.show();                         # show image
No description has been provided for this image
In [ ]: