edapy_01

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);

    2021-06-19

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");

    c:\bill\EDAWM-3rdEd\edapy\Notebooks
    c:\bill\EDAWM-3rdEd\edapy\TestFolder
    c:\bill\EDAWM-3rdEd\edapy\TestFolder
    ['TextDocument1.txt', 'TextDocument2.txt']

In [9]:

    # edapy_01_03: simple algebra
    a = 3.5;
    b = 4.1;
    c = a+b;
    print(c);

    7.6

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);

    5.0

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);

    0.5877852522924732

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);

    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 [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);

    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 [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)

    c
     [[1.]
     [3.]
     [5.]]
    r
     [[2. 4. 6.]]
    A
     [[1. 2. 3.]
     [4. 5. 6.]
     [7. 8. 9.]]

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);

    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 [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);

    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.]]
    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 [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);

    a
     [[1.]
     [2.]
     [3.]]
    M
     [[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.]]
    A
     [[1. 2. 3.]
     [4. 5. 6.]
     [7. 8. 9.]]
    B
     [[10.  2.  3.]
     [11.  5.  6.]
     [12.  8.  9.]]
    Alist [1. 2. 3. 4. 5. 6. 7. 8. 9.]

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

    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

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);

    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 [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);

    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 k
     <class 'tuple'>
    len k
     2
    type k[0]
     <class 'numpy.ndarray'>
    shape k[0]
     (3,)
    k[0]
     [2 6 8]
    d
     [[1.]
     [2.]
     [6.]
     [4.]
     [3.]
     [2.]
     [6.]
     [4.]
     [6.]
     [2.]
     [1.]
     [4.]]

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);

    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 [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);

    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 [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