edapy_13

In [65]:

    # eda13_00 clear all variables and import vatious modules
    %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
    from matplotlib import pyplot as plt
    from matplotlib import patches as pch
    import matplotlib.cm as cm
    from matplotlib.colors import ListedColormap
    from mpl_toolkits.mplot3d import Axes3D
    import re

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

    # function to make a numpy N-by-1 column vector
    # c=eda_cvec(v1, v2, ...) from a list of several
    # array-like entities v1, v2, including a number
    # a list of numbers, a tuple of numbers, an N-by-0 np array
    # and a N-by-1 np array. The function
    # also insures that, if v is an np array, that
    # c is a copy, as contrasted to a view, of it
    # It promotes integers to floats, and integers
    # and floats to complex, by context.
    # This version concatenates many argments,
    # whereas c=eda_cvec1(v1) takes just one argiment.
    # I recommend always using eda_cvec(v1, v2, ...)
    def eda_cvec(*argv):
        t = int;
        Nt = 0;
        for a in argv:
            v = eda_cvec1(a);
            N,M = np.shape(v);
            Nt = Nt + N;
            if( N==0 ):
                continue; # skip vector of zero length
            if (t==int) and isinstance(v[0,0],float):
                t=float;
            elif isinstance(v[0,0],complex):
                t=complex;
        w = np.zeros((Nt,1),dtype=t);
        Nt = 0;
        for a in argv:
            v = eda_cvec1(a);
            N,M = np.shape(v);
            w[Nt:Nt+N,0] = v[0:N,0];
            Nt = Nt + N;
        return w;

In [91]:

    # 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)
    <ipython-input-91-e52a44371377> in <module>
          8
          9 np = 3;
    ---> 10 x = np.zeros((3,1));

    AttributeError: 'int' object has no attribute 'zeros'

In [92]:

    # 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 [93]:

    # 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 [94]:

    # 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 [3]:

    # 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 [22]:

    # 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 [84]:

    # 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 [66]:

    # edapy_12_08, 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 [ ]: