edapy_05

In [1]:

    # edapy_05_00 clear all variables and import vatious modules
    # this cell must be exectuted first, before any of the others

    %reset -f
    import os
    from datetime import date
    from math import exp, pi, sin, sqrt, floor, ceil
    import numpy as np
    import scipy.linalg as la
    import scipy.sparse.linalg as las
    from scipy import sparse
    from matplotlib import pyplot as plt
    import matplotlib.cm as cm
    from matplotlib.colors import ListedColormap

    # 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 perform the multiplication FT (F v)
    def FTFmul(v):
        # this function is used by the bicongugate gradient solver to solve the geneneralized least
        # squares problem Fm=f.  Note that "F" must be called "edaFsparse".
        global edaFsparse;
        s = np.shape(v);
        if(len(s)==1):
            vv = np.zeros((s[0],1));
            vv[:,0] = v;
        else:
            vv=v;
        temp = edaFsparse*vv;
        return edaFsparse.transpose()*temp;

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

    # edapy_05_01: application of Bayes theorem to least squares

    # set up vectors m1 and m2
    L = 40;
    Dm = 1.0;
    m1 = eda_cvec( np.linspace(0,Dm*(L-1),L) );
    m2 = eda_cvec( np.linspace(0,Dm*(L-1),L) );

    # make normal distribution pp(m)
    m1bar=10;
    m2bar=20;
    Cpm=[[6,0],[0,10]];
    CpmI=la.inv(Cpm);
    norm=1/(2*pi*sqrt(la.det(Cpm)));
    Pp=np.zeros((L,L));
    dm = np.zeros((2,1));
    for i in range(L):
        for j in range(L):
            dm = eda_cvec( [m1[i,0]-m1bar, m2[j,0]-m2bar] );
            Pp[i,j]=norm*np.exp( -0.5 * np.matmul( np.matmul( dm.T, CpmI ), dm ));

    # search grid for maximum
    k = np.argmax( Pp, axis=0 );
    Pptemp = np.max( Pp, axis=0 );
    j = np.argmax(Pptemp);
    i=k[j];
    m1bar=m1[i];
    m2bar=m2[j];

    # make normal distribution pdgm(m)
    # with a one-row data kernel that the
    # difference of two model parameters is d1

    G = np.zeros((1,2));
    G[0,:] = [1, -1];
    d1 = 0;
    Cd = 3;
    CdI = 1/Cd;
    norm = 1/sqrt(2*pi*Cd);
    Pdgm = np.zeros((L,L));
    m = np.zeros((2,1));
    for i in range(L):
        for j in range(L):
            m[:,0] = [m1[i,0], m2[j,0]];
            dd = d1 - np.matmul( G, m );
            Pdgm[i,j]=norm*np.exp( -0.5 * CdI * np.matmul( dd.T,dd ) );

    # multiply the two to get Pmgd
    Pmgd = np.multiply( Pp , Pdgm);

    # search grid for maximum
    k = np.argmax( Pmgd, axis=0 );     # row indices of Pmgd with the largest value in each col
    Pmgdtemp = np.max( Pmgd, axis=0 ); # maximums in each column
    j = np.argmax(Pmgdtemp);           # indice to the column with the overall maximum
    i=k[j];                            # row of the overall maximum
    m1est=m1[i];                       # m1 position on the grid of the maximum
    m2est=m2[j];                       # m2 position on the grid of the maximum

    print('m1: prior', m1bar, 'estimate', m1est);
    print('m2: prior', m2bar, 'estimate', m2est);

    # plot result
    eda_draw('title pp(d)',Pp,'title p(d|m)',Pdgm,'title p(m|d)',Pmgd);

    m1: prior [10.] estimate [13.]
    m2: prior [20.] estimate [15.]

In [29]:

    # edapy_05_02: illustrate the product of twos Normal p.d.f.s

    # set up vectors m1 and m2
    L = 40;
    Dd = 1.0;
    d1 = eda_cvec( np.linspace(0,Dd*(L-1),L) );
    d2 = eda_cvec( np.linspace(0,Dd*(L-1),L) );

    # first Normal p.d.f.
    d1bar=15; # d1 mean
    d2bar=15; # d2 mean
    C=[[10,8],[8,10]]; # covariance
    CI=la.inv(C);      # its inverse
    norm=1/(2*pi*sqrt(la.det(C)));  # normalization factor for 2D Normal pdf
    pa=np.zeros((L,L)); # matrix to hold p.d.f.
    dd = np.zeros((2,1)); # vector to hold difference between value and mean
    for i in range(L):
        for j in range(L):
            dd[:,0] = [d1[i,0]-d1bar, d2[j,0]-d2bar];
            pa[i,j]=norm*np.exp( -0.5 * np.matmul( np.matmul( dd.T, CI ), dd ));

    # second Normal p.d.f.
    d1bar=15;
    d2bar=25;
    C=[[10,-8],[-8,10]];
    CI=la.inv(C);
    norm=1/(2*pi*sqrt(la.det(C)));
    pb=np.zeros((L,L));
    dd = np.zeros((2,1));
    for i in range(L):
        for j in range(L):
            dd[:,0] = [d1[i,0]-d1bar, d2[j,0]-d2bar];
            pb[i,j]=norm*np.exp( -0.5 * np.matmul( np.matmul( dd.T, CI ), dd ));

    # product pc = pa * pb
    pc=np.zeros((L,L));
    pc = np.multiply(pa, pb);

    # normalize to unit area
    norm = (Dd**2)*np.sum(pc);  # sum is approximation for integral
    pc = pc/norm;

    # plot results
    eda_draw('title pA',pa,'title pB',pb,'title pC', pc);

In [11]:

    # edapy_05_03: straight-line fit with unequal variance
    # illustrate least squares fitting of a line
    # when the data have unequal variance
    # this version uses G and d

    # set of vector of x's
    N = 51;
    N2 = floor((N-1)/2);
    Dx = 2;
    xmin = 0.0;
    xmax = xmin+Dx*(N-1);
    x = eda_cvec( np.linspace(xmin, xmax, N) );

    # prepare some test data hat consists of
    # two sections, that differ slightly in
    # interrcept and slope
    M=2;
    a = 1.0;
    b = 2.0;
    dobs = a + b*x;
    Da = 1;
    Db = 0.3;
    dobs[N2:N,0:1] = dobs[N2:N,0:1] + Da + Db*x[N2:N,0:1];

    # Part 1: variance, high for thr first section low for the second section
    sd = np.zeros((N,1));
    sd[0:N2,0:1] = 20.0;
    sd[N2:N,0:1] = 2.0;
    sd2i = np.power(sd,-2);
    Cdi = np.diag(sd2i.ravel());

    # set up the data kernel
    G=np.zeros((N,2));
    G[:,0]=np.ones((N,1)).ravel();
    G[:,1]=x.ravel();

    # solve for estimated model parameters
    FTF = np.matmul( G.T , np.matmul( Cdi, G ) );
    FTf = np.matmul( G.T , np.matmul( Cdi, dobs ) );
    mest = la.solve(FTF,FTf);
    # predict data
    dpre = mest[0] + mest[1] * x;

    #  plot data
    fig1 = plt.figure(1,figsize=(8,6));
    ax1 = plt.subplot(2,1,1);
    plt.axis([x[0,0], x[N-1,0], 0, 300])
    plt.plot(x,dobs,'ko');
    # 95% error bars
    for p in range(N):
        plt.plot( [x[p], x[p]], [dobs[p]-2*sd[p], dobs[p]+2*sd[p]], 'k-', lw=2 );
    plt.plot(x,dpre,'k-');
    plt.xlabel('x');
    plt.ylabel('d');
    titlestr = "intercept %.2f slope %.2f" % (mest[0,0], mest[1,0]);
    plt.title(titlestr);

    # Part 2: variance, high for thr first section low for the second section
    sd = np.zeros((N,1));
    sd[0:N2,0:1] = 2.0;
    sd[N2:N,0:1] = 20.0;
    sd2i = np.power(sd,-2);
    Cdi = np.diag(sd2i.ravel());

    # set up the data kernel
    G[0:N,0:1]=np.ones((N,1));
    G[0:N,1:2]=x;

    # solve for estimated model parameters
    FTF = np.matmul( G.T , np.matmul( Cdi, G ) );
    FTf = np.matmul( G.T , np.matmul( Cdi, dobs ) );
    mest = la.solve(FTF,FTf);
    # predict data
    dpre = mest[0,0] + mest[1,0] * x;

    #  plot data
    ax1 = plt.subplot(2,1,2);
    plt.axis([x[0,0], x[N-1,0], 0, 300])
    plt.plot(x,dobs,'ko');
    for p in range(N):
        plt.plot( [x[p], x[p]], [dobs[p]-2*sd[p], dobs[p]+2*sd[p]], 'k-', lw=2 );
    plt.plot(x,dpre,'k-');
    plt.xlabel('x');
    plt.ylabel('d');
    titlestr = "intercept %.2f slope %.2f" % (mest[0,0], mest[1,0]);
    plt.title(titlestr);
    plt.show();

In [12]:

    # edapy_05_04: Straight-line fit with diagonal variance
    # illustrate least squares fitting of a line
    # when the data have unequal variance
    # this version uses F and f

    # set of vector of x's
    N = 51;
    N2 = floor((N-1)/2);
    Dx = 2;
    xmin = 0.0;
    xmax = xmin+Dx*(N-1);
    x = eda_cvec( np.linspace(xmin, xmax, N) );

    M=2; # number of model parametres

    # prepare some test data hat consists of
    # two sections, that differ slightly in
    # interrcept and slope
    a = 1.0;
    b = 2.0;
    dobs = a + b*x;
    Da = 1;
    Db = 0.3;
    dobs[N2:N,0] = np.copy( dobs[N2:N,0] + Da + Db*x[N2:N,0] );

    # Part 1: variance, high for thr first section low for the second section
    sd = np.zeros((N,1));
    sd[0:N2,0] = 20;
    sd[N2:N,0] = 2;

    # set up the data kernel
    G=np.zeros((N,2));
    G[0:N,0:1]=np.ones((N,1));
    G[0:N,1:2]=x;

    # set up F and f
    F = np.zeros((N,M));
    for i in range(N):
        F[i:i+1,0:M] = G[i:i+1,0:M] / sd[i,0];
    f = np.divide( dobs, sd );

    # solve for estimated model parameters
    FTF = np.matmul( F.T, F );
    FTf = np.matmul( F.T, f );
    mest = la.solve(FTF,FTf);

    # predict data
    dpre = mest[0] + mest[1] * x;

    #  plot data
    fig1 = plt.figure(1,figsize=(8,6));
    ax1 = plt.subplot(2,1,1);
    plt.axis([x[0,0], x[N-1,0], 0, 300])
    plt.plot(x,dobs,'ko');
    # 95% error barx
    for p in range(N):
        plt.plot( [x[p], x[p]], [dobs[p]-2*sd[p], dobs[p]+2*sd[p]], 'k-', 'LineWidth', 2 );
    plt.plot(x,dpre,'k-');
    plt.xlabel('x');
    plt.ylabel('d');
    titlestr = "intercept %.2f slope %.2f" % (mest[0,0], mest[1,0]);
    plt.title(titlestr);
    plt.show();

In [2]:

    # edapy_05_05: using smoothness to fill in data gaps
    # simple example of smoothness applied
    # to the problem of filling in data gaps
    # this version uses non-sparse arrays

    # the model parameters are the values of a sinusoidal
    # function evenly spaced along the x axis
    M=101;
    Dx=1.0;
    Dx2 = Dx**2;
    xmin=0.0;
    xmax = xmin+Dx*(M-1);
    x = eda_cvec( np.linspace(xmin, xmax, M) );

    # the data are the values the function measured on
    # only a subset of points along the x axis
    N=40;

    # randomly choose N of the M points to have data
    rowindex = np.sort( np.random.randint(M,size=N) );  # N random integers between 1 and M
    xd = eda_cvec( x[rowindex,:] );  # xd are the x's of the N data

    # data is sinusoidal
    freq=0.02;
    d = np.sin(2*pi*freq*xd);

    # prior variances
    # 1% error of signal level of 1 i 0.1
    sigmadi =  1.0/0.01;
    # maximum slope of about (1-0)/10 = 0.1
    sigmaD2i = 1.0/0.1;
    # maximum 2nd der about (0-2*1-0)/(20*20) = 0.005
    sigmaD1i = 1.0/0.005;

    # build F and f
    K=M;
    F=np.zeros((N+K,M));
    fobs=np.zeros((N+K,1));
    nextrow = 0;

    # the model parameter equals the observed data
    for p in range(N):
        F[nextrow,rowindex[p]] = sigmadi;
        fobs[nextrow,0]=sigmadi*d[p,0];
        nextrow=nextrow+1;

    # smoothness (small second derivative) in the interior
    for p in range(M-2):
        F[nextrow,p] = sigmaD2i/Dx2;
        F[nextrow,p+1] = -2*sigmaD2i/Dx2;
        F[nextrow,p+2] = sigmaD2i/Dx2;
        fobs[nextrow,0]=0.0;
        nextrow=nextrow+1;

    # flatness (small first derivative) at the left end
    F[nextrow,0]=-sigmaD1i*Dx;
    F[nextrow,1]=sigmaD1i*Dx;
    fobs[nextrow,0]=0;
    nextrow=nextrow+1;

    # flatness (small first derivative) at the right end
    F[nextrow,M-2]=-sigmaD1i*Dx;
    F[nextrow,M-1]=sigmaD1i*Dx;
    fobs[nextrow,0]=0;
    nextrow=nextrow+1;

    print("built %d rows out of %d" % (nextrow,N+K));

    # solve generalized least squares equation
    FTF = np.matmul( F.T, F );
    FTf = np.matmul( F.T, fobs );
    mest = np.zeros((M,1));
    mest[0:M,0:1] = la.solve(FTF,FTf);

    # error
    e = d-mest[rowindex];      # individual error
    E = np.matmul( e.T, e );   # overall error

    # plot data
    fig1 = plt.figure(1,figsize=(8,6));
    ax1 = plt.subplot(2,1,1);
    plt.axis([x[0,0], x[M-1,0], -1.5, 1.5])
    plt.plot(xd,d,'ko');
    plt.plot(x,mest,'k-');
    plt.xlabel('x');
    plt.ylabel('d');
    plt.show();

    built 141 rows out of 141

In [21]:

    # eda05_06: sparse matrix from list of non-zero elements

    # build lists
    L=7;
    M=6;
    Frow = eda_cvec( 0, 2, 3, 5);
    Fcol = eda_cvec( 0, 2, 0, 4);
    Fval = eda_cvec( 1.0, 2.0, 3.0, 4.0);
    KK, tmp = np.shape(Fval);

    # build sparse matrix
    edaFsparse = sparse.coo_matrix((Fval[0:KK,0],
                    (Frow[0:KK,0],Fcol[0:KK,0])), shape=(L,M) );

    M=sparse.csr_matrix.todense(edaFsparse);
    print(M);

    [[1. 0. 0. 0. 0. 0.]
     [0. 0. 0. 0. 0. 0.]
     [0. 0. 2. 0. 0. 0.]
     [3. 0. 0. 0. 0. 0.]
     [0. 0. 0. 0. 0. 0.]
     [0. 0. 0. 0. 4. 0.]
     [0. 0. 0. 0. 0. 0.]]

In [19]:

    # edapy_05_07: sparse matrix implementation of Fm=f
    # simple example of smoothness applied
    # to the problem of filling in data gaps
    # this version uses sparse arrays

    # the model parameters are the values of a sinusoidal
    # function evenly spaced along the x axis
    M=101;
    Dx=1.0;
    Dx2 = Dx**2;
    xmin=0.0;
    xmax = xmin+Dx*(M-1);
    x = eda_cvec( np.linspace(xmin, xmax, M) );

    # the data are the values the function measured on
    # only a subset of points along the x axis
    N=40;

    # randomly choose N of the M points to have data
    rowindex = np.sort( np.random.randint(M,size=N) ); # N random integers between 1 and M
    xd = eda_cvec( x[rowindex,:] ); # x's of the data

    # data is sinusoidal
    freq=0.02;
    dobs = np.sin(2*pi*freq*xd);

    # prior variances
    # 1% error of signal level of 1 i 0.1
    sigmadi =  1.0/0.01;
    # maximum slope of about (1-0)/10 = 0.1
    sigmaD2i = 1.0/0.1;
    # maximum 2nd der about (0-2*1-0)/(20*20) = 0.005
    sigmaD1i = 1.0/0.005;

    # set up f
    L=N+M;  # rows in F, N data, and N priors
    fobs=np.zeros((L,1));

    # set up F.  In this impletation, the non-zero
    # elements of F are put first in a vector Fval, with
    # two auxilary arrays Frow and Fcol
    K = int;
    K = N + 3*(M-2) + 2*2 + 100;  # needs to be at least as large as the number of non-zero elments in F
    Frow = np.zeros((K,1),dtype=np.int32);  # the row index vector is integer
    Fcol = np.zeros((K,1),dtype=np.int32);  # the col index vector is integer
    Fval = np.zeros((K,1),dtype=np.double); # the value vector is double precision

    # a trick for keeping tracks of rows is to use a row counter
    # an increment it after each row is added.  The counters
    # are defined as integers.
    nextrow = int;  # next row in F; note is an integer
    nextrow = 0;    # starting value of zero
    KK = int;  # next element in Fval; note is an integer
    KK = 0;    # starting value of zero

    # the first N rows of F are that the model parameter
    # equals the observed data
    for p in range(N):
        # F[nextrow,rowindex[p]] = 1;
        Frow[KK,0]=nextrow;
        Fcol[KK,0]=rowindex[p];
        Fval[KK,0]=sigmadi;
        KK=KK+1;
        # f[nextrow]=d[p];
        fobs[nextrow,0]=sigmadi*dobs[p,0];
        nextrow=nextrow+1;

    # implement the smoothness constraints in the interior
    # of the x interval
    for p in range(M-2):
        # F[nextrow,p] = shi/Dx2;
        Frow[KK,0]=nextrow;
        Fcol[KK,0]=p;
        Fval[KK,0]=sigmaD2i/Dx2;
        KK=KK+1;
        # F[nextrow,p+1] = -2*shi/Dx2;
        Frow[KK,0]=nextrow;
        Fcol[KK,0]=p+1;
        Fval[KK,0]=-2*sigmaD2i/Dx2;
        KK=KK+1;
        # F[nextrow,p+2] = shi/Dx2;
        Frow[KK,0]=nextrow;
        Fcol[KK,0]=p+2;
        Fval[KK,0]=sigmaD2i/Dx2;
        KK=KK+1;
        # f[nextrow,0]=0.0;
        fobs[nextrow,0]=0.0;
        nextrow=nextrow+1;

    # implement flatness constraint at left end of the x interval
    # F[nextrow,0]=-shi*Dx;
    Frow[KK,0]=nextrow;
    Fcol[KK,0]=0;
    Fval[KK,0]=-sigmaD1i*Dx;
    KK=KK+1;
    # F[nextrow,1]=shi*Dx;
    Frow[KK,0]=nextrow;
    Fcol[KK,0]=1;
    Fval[KK,0]=sigmaD1i*Dx;
    KK=KK+1;
    # f[nextrow,0]=0;
    fobs[nextrow,0]=0;
    nextrow=nextrow+1;

    # implement flatness constraint at right end of the x interval
    # F[nextrow,M-2]=-shi*Dx;
    Frow[KK,0]=nextrow;
    Fcol[KK,0]=M-2;
    Fval[KK,0]=-sigmaD1i*Dx;
    KK=KK+1;
    # F[nextrow,M-1]=shi*Dx;
    Frow[KK,0]=nextrow;
    Fcol[KK,0]=M-1;
    Fval[KK,0]=sigmaD1i*Dx;
    KK=KK+1;
    # f[nextrow]=0;
    fobs[nextrow]=0;
    nextrow=nextrow+1;

    print("%d non-zero-elements, expected %d" % (KK, N+3*(M-2)+2*2) );
    print("%d rows, expected %d" % (nextrow, N+M) );
    # actual rows in F

    # the "F" matrix must be called "edaFsparse" in order to correctly interface
    # witk the FYFmul() function called by the biconjugate gradient solver
    # Note the reduced length of the arrays Fval, Frow, Fcol
    edaFsparse = sparse.coo_matrix((Fval[0:KK,0], (Frow[0:KK,0],
                                Fcol[0:KK,0])), shape=((N+M),M) );
    del Fval;  # delete no loger needed arrays
    del Frow;
    del Fcol;

    # define linear operator needed for conjugate gradienet solver
    LO=las.LinearOperator(shape=(M,M),
                          matvec=FTFmul,rmatvec=FTFmul);

    # solve for eusing congugate gradient algorithm
    mest = np.zeros((M,1));
    q=las.bicg(LO,edaFsparse.T*fobs,tol=1e-12,
               maxiter=(3*(N+M)+100));
    mest[:,0] = q[0];

    # error
    e = dobs-mest[rowindex];  # individual errors
    E = np.matmul( e.T, e );    # overall error
    print('E:',E);

    # plot data
    fig1 = plt.figure(1,figsize=(8,6));
    ax1 = plt.subplot(2,1,1);
    plt.axis([x[0,0], x[M-1,0], -1.5, 1.5])
    plt.plot(xd,dobs,'ko');
    plt.plot(x,mest,'k-');
    plt.xlabel('x');
    plt.ylabel('d');
    plt.show();

    341 non-zero-elements, expected 341
    141 rows, expected 141
    E: [[3.1975708e-07]]

In [5]:

    # edapy_05_08: folding out a matrix into a vector

    I = 4;   # number of rows of a matrix A
    J = 3;   # number of coluns of a matrix A
    M = I*J; # total number of elements

    # make a matrix A
    A=np.zeros((I,J));
    A[:,:] = [ [1, 5, 9], [2, 6, 10], [3, 7, 11], [4, 8, 12]];
    print('A\n', A);

    # make a Mx1 vector m
    m=np.zeros((M,1));

    # PART 1: using a formula
    print('\nPart 1: with a formula');

    # copy elements of A[i,j] into m[k]
    for i in range(I):
        for j in range(J):
            k = i*J+j; # row-wise unfold
            m[k,0] = A[i,j];
    print('m\n', m);

    # copy the elements of m[k] into a matrix B[i,j]
    B=np.zeros((I,J));
    for k in range(M):
        i = floor(k/J); # row-wise fold
        j = k-i*J;
        B[i,j]=m[k];
    print('B\n',B);

    # PART 2: with index arrays
    print('\nPart 2: with index arrays')

    # the idea is to define a set of "pointer arrays"
    # that translate between (i,j) and k.  The arrays
    # are built during one traverse of (i,j).  One
    # advanage is that no "formulas" are needed.
    i_of_k  = np.zeros((M,1), dtype=np.int32);  # i given k
    j_of_k =  np.zeros((M,1), dtype=np.int32);  # j given k
    k_of_ij = np.zeros((I,J), dtype=np.int32);  # k given (i,j)
    k=0; # initialize counter
    for i in range(I):
        for j in range(J):
            i_of_k[k,0]=i;
            j_of_k[k,0]=j;
            k_of_ij[i,j]=k;
            k=k+1;
    M=k;

    # copy A[i,j] to m[k]
    m=np.zeros((M,1));
    for i in range(I):
        for j in range(J):
            k = k_of_ij[i,j];
            m[k,0] = A[i,j];
    print('m\n', m);

    # copy m[k] to B[i,j]
    B=np.zeros((I,J));
    for k in range(M):
        i = i_of_k[k,0];
        j = j_of_k[k,0];
        B[i,j]=m[k,0];
    print('B\n',B);

    A
     [[ 1.  5.  9.]
     [ 2.  6. 10.]
     [ 3.  7. 11.]
     [ 4.  8. 12.]]

    Part 1: with a formula
    m
     [[ 1.]
     [ 5.]
     [ 9.]
     [ 2.]
     [ 6.]
     [10.]
     [ 3.]
     [ 7.]
     [11.]
     [ 4.]
     [ 8.]
     [12.]]
    B
     [[ 1.  5.  9.]
     [ 2.  6. 10.]
     [ 3.  7. 11.]
     [ 4.  8. 12.]]

    Part 2: with index arrays
    m
     [[ 1.]
     [ 5.]
     [ 9.]
     [ 2.]
     [ 6.]
     [10.]
     [ 3.]
     [ 7.]
     [11.]
     [ 4.]
     [ 8.]
     [12.]]
    B
     [[ 1.  5.  9.]
     [ 2.  6. 10.]
     [ 3.  7. 11.]
     [ 4.  8. 12.]]

In [25]:

    # edapy_05_09: fill in data gaps in 2D using diff eqn
    # use generealized least square to fill in
    # data gaps for data that satisfy Laplace's
    # equation in the in the interior of a grid,
    # and the normal derivative is zero on the
    # edges of the grid

    # load the data.  These data are only 'simulated'
    # (also called "synthetic") and were created by the
    # script eda05_09

    # load (x, y, pressure) data from file
    D = np.genfromtxt('../data/pressure.txt', delimiter='\t')
    [N, K]=D.shape;
    xobs = eda_cvec( D[:,0] );
    yobs = eda_cvec( D[:,1] );
    dobs = eda_cvec( D[:,2] );

    # assign variance of data
    sigmad = 0.1;
    rsigmad = 1/sigmad;
    rsigmad2 = rsigmad**2;

    # rows of grid are x
    I = 41;
    Dx = 1;
    xmin = 0;
    xmax = Dx*(I-1);
    x = eda_cvec( np.linspace(xmin, xmax, I) );
    rDx = 1/Dx;
    rDx2 = rDx**2;

    # columns of grid are y
    J = 41;
    Dy = 1;
    ymin = 0;
    ymax = Dy*(J-1);
    y = eda_cvec( np.linspace(ymin, ymax, J) );
    rDy = 1/Dy;
    rDy2 = rDy**2;

    # associate observations with points on the grid
    # by scaling the range (xmin, xmax) to (1, N)
    rowindex = np.zeros( (N,1), dtype=np.int32 );
    rowindex[:,0] = np.floor(I*(xobs-xmin)/(xmax-xmin)+0.5).ravel();
    colindex = np.zeros( (N,1), dtype=np.int32 );
    colindex[:,0] = np.floor(J*(yobs-ymin)/(ymax-ymin)+0.5).ravel();

    # don't let any data fall outside the grid.  In this
    # case, we just put it on the edge of the grid.  The
    # alternative is to throw it away.
    for i in range(N):
        if( rowindex[i,0]<0 ):
            rowindex[i,0]=1;
        elif( rowindex[i,0]>=I ):
            rowindex[i,0]=I-1;
        if( colindex[i,0]<0 ):
            colindex[i,0]=0;
        elif( colindex[i,0]>=J ):
            colindex[i,0]=J-1;

    # copy the data to the grid A for later use
    Aobs = np.zeros((I,J));
    for i in range(N):
        Aobs[rowindex[i,0],colindex[i,0]] = dobs[i,0];

    # one model parameter per grid point
    M = I*J;

    # Fm=f is LxM with N rows of data plus M rows of prior information
    L = N+M;

    # F will have no more than K elements
    K = N + 5*M + 100;

    # r.h.s. of Fm=F
    fobs = np.zeros((L,1));

    # reminders on how to move from grid to vector & vice-versa
    # m(k) = A(i,j);
    # k = i*J+j;
    # i = floor(k/J);
    # j = k-i*J;
    # B(i,j)=m(k);

    # in this implementation, I build vectors on zon-zero elements and
    # then change them into a sparse matrix at the end
    Frow = np.zeros((K,1),dtype=np.int32);
    Fcol = np.zeros((K,1),dtype=np.int32);
    Fval = np.zeros((K,1),dtype=np.double);

    # build F and f

    # nextrow is used to keep track of the current row of F and f
    nextrow = int;
    nextrow = 0;

    # KK is used to keep track of the current row of Fvec, Fcol, Frow
    KK = int;
    KK = 0;

    # k=dindex(i): data i corresponds to row k of m
    dindex = np.zeros((N,1),dtype=np.int32);

    # first part of Fm=f: the data
    for i in range(N):
        k = rowindex[i]*J+colindex[i];
        # F[q,k] = rsigmad;
        Frow[KK,0] = nextrow;
        Fcol[KK,0] = k;
        Fval[KK,0] = rsigmad;
        KK = KK+1;
        fobs[nextrow,0] = rsigmad*dobs[i,0];
        dindex[i:0] = k;
        nextrow = nextrow+1;

    # error of prior information
    sigmah = 1; # was 1000
    rsigmah = 1/sigmah;

    # second part, laplaces equation in interior
    # second derivative is a (1,-2,1)/Dx^2 in x plus a (1,-2,1)/Dy^2 in y
    for i in range(1,I-1): # Matlab [2:I-1]
        for j in range(1,J-1): # Matlab [2:J-1]
            # k = (i-1)*J+j; and F(q,k) = rsigmah*(-2*rDx2-2*rDy2);
            k = i*J+j;
            Frow[KK,0] = nextrow;
            Fcol[KK,0] = k;
            Fval[KK,0] = rsigmah*(-2*rDx2-2*rDy2);
            KK=KK+1;
            # k = ((i-1)-1)*J+j; and F(q,k) = rsigmah*rDx2;
            k = (i-1)*J+j;
            Frow[KK,0] = nextrow;
            Fcol[KK,0] = k;
            Fval[KK,0] = rsigmah*rDx2
            KK=KK+1;
            # k = ((i+1)-1)*J+j; and F(q,k) = rsigmah*rDx2;
            k = (i+1)*J+j;
            Frow[KK,0] = nextrow;
            Fcol[KK,0] = k;
            Fval[KK,0] = rsigmah*rDx2;
            KK=KK+1;
            #  k = (i-1)*J+(j-1); and F(q,k) = rsigmah*rDy2;
            k = i*J+(j-1);
            Frow[KK,0] = nextrow;
            Fcol[KK,0] = k;
            Fval[KK,0] = rsigmah*rDy2;
            KK=KK+1;
            # k = (i-1)*J+(j+1); and F(q,k) = rsigmah*rDy2;
            k = i*J+(j+1);
            Frow[KK,0] = nextrow;
            Fcol[KK,0] = k;
            Fval[KK,0] = rsigmah*rDy2;
            KK=KK+1;
            fobs[nextrow,0] = 0.0;
            nextrow = nextrow+1;

    # third part: top edge
    for j in range(1,J-1): # MATLAB [2:J-1]
        i = 0;
        # k = (i-1)*J+j; and F(q,k) = -rsigmah*rDx;
        k = i*J+j;
        Frow[KK,0] = nextrow;
        Fcol[KK,0] = k;
        Fval[KK,0] = -rsigmah*rDx;
        KK=KK+1;
        # k = ((i+1)-1)*J+j; and F(q,k) = rsigmah*rDx;
        k = (i+1)*J+j;
        Frow[KK,0] = nextrow;
        Fcol[KK,0] = k;
        Fval[KK,0] = rsigmah*rDx;
        KK=KK+1;
        fobs[nextrow,0] = 0.0;
        nextrow = nextrow+1;

    # fourth part: bottom edge
    for j in range(1,J-1): # MATLAB [2:J-1]
        i = I-1;
        # k = (i-1)*J+j; and F(q,k) = rsigmah*rDx;
        k = i*J+j;
        Frow[KK,0] = nextrow;
        Fcol[KK,0] = k;
        Fval[KK,0] = rsigmah*rDx;
        KK=KK+1;
        # k = ((i-1)-1)*J+j; and F(q,k) = -rsigmah*rDx;
        k = (i-1)*J+j;
        Frow[KK,0] = nextrow;
        Fcol[KK,0] = k;
        Fval[KK,0] = -rsigmah*rDx;
        KK=KK+1;
        fobs[nextrow,0] = 0.0;
        nextrow = nextrow+1;

    # fifth part: left edge
    for i in range(1,I-1): # MATLAB [2:I-1]:
        j = 0;
        # k = (i-1)*J+j; and F(q,k) = -rsigmah*rDy;
        k = i*J+j;
        Frow[KK,0] = nextrow;
        Fcol[KK,0] = k;
        Fval[KK,0] = -rsigmah*rDy;
        KK=KK+1;
        # k = (i-1)*J+(j+1); and F(q,k) = rsigmah*rDy;
        k = i*J+(j+1);
        Frow[KK,0] = nextrow;
        Fcol[KK,0] = k;
        Fval[KK,0] = rsigmah*rDy;
        KK=KK+1;
        fobs[nextrow,0] = 0.0;
        nextrow = nextrow+1;

    # sixth part: right edge
    for i in range(1,I-1): # MATLAB [2:I-1]
        j = J-1;
        # k = (i-1)*J+j; and F(q,k) = rsigmah*rDy;
        k = i*J+j;
        Frow[KK,0] = nextrow;
        Fcol[KK,0] = k;
        Fval[KK,0] = rsigmah*rDy;
        KK=KK+1;
        # k = (i-1)*J+(j-1); and F(q,k) = -rsigmah*rDy;
        k = i*J+(j-1);
        Frow[KK,0] = nextrow;
        Fcol[KK,0] = k;
        Fval[KK,0] = -rsigmah*rDy;
        KK=KK+1;
        fobs[nextrow,0] = 0.0;
        nextrow = nextrow+1;

    # seventh part: four corners
    # top left
    rDz = 1/sqrt((Dx^2)+(Dy^2));
    i=0;
    j=0;
    # k = (i-1)*J+j; and F(q,k) = -rsigmah*rDz;
    k = i*J+j;
    Frow[KK,0] = nextrow;
    Fcol[KK,0] = k;
    Fval[KK,0] = -rsigmah*rDz;
    KK=KK+1;
    # k = ((i+1)-1)*J+(j+1); F(q,k) = rsigmah*rDz;
    k=(i+1)*J+(j+1);
    Frow[KK,0] = nextrow;
    Fcol[KK,0] = k;
    Fval[KK,0] = rsigmah*rDz;
    KK=KK+1;
    fobs[nextrow,0]=0.0;
    nextrow=nextrow+1;
    # top right
    i=0;
    j=J-1;
    # k = (i-1)*J+j; and F(q,k) = rsigmah*rDz;
    k = i*J+j;
    Frow[KK,0] = nextrow;
    Fcol[KK,0] = k;
    Fval[KK,0] = rsigmah*rDz;
    KK=KK+1;
    # k = ((i+1)-1)*J+(j-1); F(q,k) = -rsigmah*rDz;
    k=(i+1)*J+(j-1);
    Frow[KK,0] = nextrow;
    Fcol[KK,0] = k;
    Fval[KK,0] = -rsigmah*rDz;
    KK=KK+1;
    fobs[nextrow,0]=0.0;
    nextrow=nextrow+1;
    # bottom left
    i=I-1;
    j=0;
    # k = (i-1)*J+j; and F(q,k) = -rsigmah*rDz;
    k = i*J+j;
    Frow[KK,0] = nextrow;
    Fcol[KK,0] = k;
    Fval[KK,0] = -rsigmah*rDz;
    KK=KK+1;
    # k = ((i-1)-1)*J+(j+1); and F(q,k) = rsigmah*rDz;
    k = (i-1)*J+(j+1);
    Frow[KK,0] = nextrow;
    Fcol[KK,0] = k;
    Fval[KK,0] = rsigmah*rDz;
    KK=KK+1;
    fobs[nextrow,0]=0.0;
    nextrow=nextrow+1;
    # bottom right
    i=I-1;
    j=0;
    # k = (i-1)*J+j; and F(q,k) = rsigmah*rDz;
    k = i*J+j;
    Frow[KK,0] = nextrow;
    Fcol[KK,0] = k;
    Fval[KK,0] = rsigmah*rDz;
    KK=KK+1;
    # k = ((i-1)-1)*J+(j-1); and F(q,k) = -rsigmah*rDz;
    k = (i-1)*J+(j-1);
    Frow[KK,0] = nextrow;
    Fcol[KK,0] = k;
    Fval[KK,0] = -rsigmah*rDz;
    KK=KK+1;
    fobs[nextrow,0]=0.0;
    nextrow=nextrow+1;

    # actual rows in F
    LL = int;
    LL = nextrow;

    print("%d non-zero elements" % (KK));
    print("%d rows, expected %d" % (LL, (N+M)));

    # the "F" matrix must be called "edaFsparse" in order to correctly interface
    # witk the FYFmul() function called by the biconjugate gradient solver
    # Note the reduced length of the arrays Fval, Frow, Fcol
    edaFsparse = sparse.coo_matrix((Fval[0:KK,0], (Frow[0:KK,0], Fcol[0:KK,0])), shape=(LL,M) );
    del Fval;  # delete the arrays used to build the sparse matrix
    del Frow;
    del Fcol;

    # define linear operator needed for conjugate gradient solver
    LO=las.LinearOperator(shape=(M,M),matvec=FTFmul,rmatvec=FTFmul);

    # solve for estimated model parameters using congugate gradient algrorithm
    mest = np.zeros((M,1));
    q=las.bicg(LO,edaFsparse.transpose()*fobs,tol=1e-6, maxiter=(3*(N+M)+100));
    mest[:,0] = q[0];

    # reorganize estimated model parameters back into grid
    Apre = np.zeros((I,J));
    for k in range(M):
        i = floor(k/J);
        j = k-i*J;
        Apre[i,j] = mest[k];

    # plot results
    eda_draw('title p-obs',Aobs,'title p-pre',Apre);

    # compute weighted prediction error Ed
    dpre=np.zeros((N,1));
    for i in range(N):
        k=dindex[i];
        dpre[i,0] = mest[k,0];
    e = dobs-dpre;
    Ed = np.matmul( e.T, e );
    sigmadpos = sqrt(Ed/N);
    print('posterior errror in data',sigmadpos);

    8093 non-zero elements
    1849 rows, expected 1849

    posterior errror in data 2.3415950502294525

In [28]:

    # edapy_05_10: make simulated (or "synthetic") pressure data

    # note: in order not to accidentally overwrite the data file
    # in the distribution, this code only writes the data to
    # mypressure.txt not pressure.txt; you must rename it to
    # have is used by eda05_08

    # rows of grid are x
    I = 41;
    Dx = 1;
    xmin = 0;
    xmax = Dx*(I-1);
    xobs = eda_cvec( np.linspace(xmin, xmax, I) );

    # columns of grid are y
    J = 41;
    Dy = 1;
    ymin = 0;
    ymax = Dy*(J-1);
    yobs = eda_cvec( np.linspace(ymin, ymax, J) );

    # randomly choose N of the I*J points to have data
    N = floor(I*J/10);
    rowindex = np.random.randint(I,size=N);
    xd = eda_cvec( xobs[rowindex,:] ); # x coords of data
    colindex = np.random.randint(J,size=N);
    yd = eda_cvec( yobs[colindex,:] ); # y coords of data

    # synthetoc data is chosen to be a solution to Laplace's equation
    kappa = 0.1;
    dtrue = eda_cvec(
        10*np.multiply(np.sin(kappa*xd),np.exp(-kappa*yd))  );

    # synthetic observed data are true data plus random noise
    sigmad = 0.1;
    noise = np.random.normal(loc=0.0, scale=sigmad, size=(N,1));
    dobs = dtrue + noise;

    # populate a matrix A with data and display it
    A = np.zeros((I,J));
    for i in range(N):
        A[rowindex[i],colindex[i]] = dobs[i];
    eda_draw('title pressure',A);

    # save matrix to file of tab-separated values
    # note that the filename has been changed to "pressure2.txt"
    # so that the original file "pressure.txr" is not overwritten
    # If you want to investigate the performance of eda05_08 with
    # different choices of data, change the filename in eda05_08 to
    # "pressure2.txt"
    DD = np.concatenate( (xd,yd,dobs),axis=1);  # concatenate (x,y,data) vectors into matrix
    np.savetxt("../data/mypressure.txt",DD, delimiter="\t"); # write the matrix to a file

In [ ]: