edapy_10

In [2]:

    # edapy_10_00 clear all variables and import vatious modules
    %reset -f
    import os
    from datetime import date
    from math import exp, pi, sin, cos, tan, sqrt, floor, ceil, log
    import numpy as np
    import scipy.sparse.linalg as las
    import scipy.interpolate as ip
    import scipy.spatial as sp
    from scipy import sparse
    import scipy.linalg as la
    import scipy.signal as sg
    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;

    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;

    def GLSFilterMul(v):
        # this function is used by the bicongugate gradient solver to solve the
        # geneneralized least squares problem Fm=f with F=[G;H] and f=[d,h],
        # where G is a topplitz matrix built from the filter g
        # Note that "H" must be sparese and called "edaHsparse" and that
        # the filter "g" bust be called edafilterg and must be a column vector
        global edaHsparse;
        global edafilterg;
        # rearrange as column-vector
        s = np.shape(v);
        M = s[0];
        if(len(s)==1):
            vv = np.zeros((s[0],1));
            vv[:,0] = v;
        else:
            vv=v;
        N, i = np.shape(edafilterg);
        # implements (GT G + HT H) v
        temp1 = np.zeros((N+M-1,1));
        temp1[:,0] = np.convolve(edafilterg.ravel(),vv.ravel()); # G v is of length N
        a = np.zeros((N,1));
        a[:,0] = temp1[0:N,0];
        b = edaHsparse * v; # H v is of length K
        temp2 = np.zeros((N+N-1,1));
        temp2[:,0] = np.convolve((np.flipud(edafilterg)).ravel(),a.ravel()); # GT (G v) is of length M
        a2 = temp2[N-1:N+M-1,0];
        b2 = edaHsparse.transpose()*b; # HT (H v) is of length M
        # FT F v = GT G v + HT H v
        return (a2+b2);

    def GLSautoMul(v):
        # this function is used by the bicongugate gradient solver to solve the
        # geneneralized least squares problem Fm=f with F=[G;H] and f=[d,h],
        # where GTG is built from the autocorrelation of the filter g
        global edaHsparse;
        global edaautoc;
        # rearrange as column-vector
        s = np.shape(v);
        M = s[0];
        if(len(s)==1):
            vv = np.zeros((s[0],1));
            vv[:,0] = v;
        else:
            vv=v;
        # FT F v = GT G v + HT H v
        GTGv=np.zeros((M,1));
        for i in range(M):
            u = np.concatenate( (np.flipud(edaautoc[0:i+1,0]), edaautoc[1:M-i,0]), axis=0 );
            GTGv[i,0] =  np.matmul( u.T, vv );
        Hv = edaHsparse*vv;
        HTHv = edaHsparse.transpose()*Hv;
        return (GTGv + HTHv);

    def chebyshevfilt(d, Dt, flow, fhigh):
        #chebyshevfilt
        # chebyshev IIR bandpass filter
        # d - input array of data
        # Dt - sampling interval
        # flow - low pass frequency, Hz
        # fhigh - high pass frequency, Hz
        # dout - output array of data
        # u - the numerator filter
        # v - the denominator filter
        # these filters can be used again using dout=filter(u,v,din);

        # make sure input timeseries is a column vector
        s = np.shape(d);
        N = s[0];
        if(N==1):
            dd = np.zeros((N,1));
            dd[:,0] = d;
        else:
            dd=d;

        # sampling rate
        rate=1/Dt;

        # ripple parameter, set to ten percent
        ripple=0.1;

        # normalise frequency
        fl=2.0*flow/rate;
        fh=2.0*fhigh/rate;

        # center frequency
        cf = 4 * tan( (fl*pi/2) ) * tan( (fh*pi/2) );

        # bandwidth
        bw = 2 * ( tan( (fh*pi/2) ) - tan( (fl*pi/2) ) );

        # ripple parameter factor
        rpf = sqrt((sqrt((1.0+1.0/(ripple*ripple))) + 1.0/ripple));
        a = 0.5*(rpf-1.0/rpf);
        b = 0.5*(rpf+1.0/rpf);

        u=np.zeros((5,1));
        v=np.zeros((5,1));
        theta = 3*pi/4;
        sr = a * cos(theta);
        si = b * sin(theta);
        es = sqrt(sr*sr+si*si);
        tmp= 16.0 - 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw - 4.0*bw*cf*sr + cf*cf;
        v[0,0] = 1.0;
        v[1,0] = 4.0*(-16.0 + 8.0*bw*sr - 2.0*bw*cf*sr + cf*cf)/tmp;
        v[2,0] = (96.0 - 16.0*cf - 8.0*es*es*bw*bw + 6.0*cf*cf)/tmp;
        v[3,0] = (-64.0 - 32.0*bw*sr + 8.0*bw*cf*sr + 4.0*cf*cf)/tmp;
        v[4,0] = (16.0 + 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw + 4.0*bw*cf*sr + cf*cf)/tmp;
        tmp = 4.0*es*es*bw*bw/tmp;
        u[0,0] = tmp;
        u[1,0] = 0.0;
        u[2,0] = -2.0*tmp;
        u[3,0] = 0.0;
        u[4,0] = tmp;

        dout = sg.lfilter(u.ravel(),v.ravel(),dd.ravel());
        return (dout,u,v);

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

    def eda_sign(x):
    # sign function; apparently Python doesn't have one
        if x>0.0:
            s=1.0;
        elif x<0.0:
            s=(-1.0);
        else:
            s=0.0;
        return(s);

    # function to perform the multiplication y = (Ccc + sigmad2 I) v
    def eda_CpImul(v):
        global eda_N;
        global eda_td;
        global eda_sigmad2;
        global eda_gamma2;
        global eda_w0est;
        global eda_sest;

        # rearrange as column-vector
        s = np.shape(v);
        M = s[0];
        if(len(s)==1):
            vv = np.zeros((s[0],1));
            vv[:,0] = v;
        else:
            vv=v;
    # Note (C + sigmad2 I) v implies yi = SUMj Cij vj + sigmad2 vi
        y = np.zeros((eda_N,1));
        for i in range(eda_N):
            y[i,0]=eda_sigmad2*vv[i,0];
            for j in range(eda_N):
                Tabs = abs(eda_td[i,0]-eda_td[j,0]);
                Cij = eda_gamma2*cos(eda_w0est*Tabs)*exp(-eda_sest*Tabs);
                y[i,0]=y[i,0]+Cij*vv[j,0];
        return y;

In [2]:

    # edapy10_01, fitting a high-order polynomial
    # illustrate problem associated with interpolation by
    # fitting a high-order polynomial

    # times
    N=11;
    t = eda_cvec( [0.01, 0.11, 0.24, 0.38, 0.42, 0.57, 0.62, 0.77, 0.81, 0.91, 0.99] );

    # simple function
    d =  t/2.0 + eda_cvec([0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1]);

    # plot simple function
    fig1 = plt.figure(1,figsize=(10,6));
    plt.axis([0, 1, -5, 5]);
    plt.plot( t, d, 'ko');
    plt.xlabel('time t (s)');
    plt.ylabel('d(t)');

    # linear representation of polynomial relationship
    G=np.zeros((N,N));
    for i in range(N):
        G[:,i]=np.power(t,i).ravel();

    # solve for coefficients
    m = la.solve(G,d);

    # evaluate at equally-spaced values of t
    Np=101;
    tp = eda_cvec(np.linspace(0.0,1.0,Np));
    Gp = np.zeros((Np,N));
    for i in range(N):
        Gp[:,i] = np.power( tp, i ).ravel();
    dp = np.matmul(Gp,m);

    # plot simple function
    plt.plot(tp,dp,'k-');
    plt.show();

In [3]:

    # eda10_02, linear interpolation

    # times
    N=11;
    t = eda_cvec( [0.01, 0.11, 0.24, 0.38, 0.42, 0.57, 0.62, 0.77, 0.81, 0.91, 0.99] );

    # simple function
    d = t/2 + eda_cvec( [0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1] );

    # plot simple function
    fig1 = plt.figure(1,figsize=(10,6));
    plt.axis([0, 1, -5, 5]);
    plt.plot( t, d, 'ko');
    plt.xlabel('time t (s)');
    plt.ylabel('d(t)');

    # qually-spaced values of t
    Np=101;
    tp = eda_cvec( np.linspace(0.0,1.0,Np) );

    # linear interpolation
    f = ip.interp1d(t.ravel(),d.ravel(),axis=0,
           kind='linear',fill_value='extrapolate');
    dp = eda_cvec(f(tp));

    # plot simple function
    plt.plot(tp,dp,'k-');
    plt.show();

In [5]:

    # edapy_10_03, cubic interpolation

    # times
    N=11;
    t = eda_cvec( [0.01, 0.11, 0.24, 0.38, 0.42, 0.57, 0.62, 0.77, 0.81, 0.91, 0.99] );

    # simple function
    d = t/2.0 + eda_cvec( [0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1, -0.1, 0.1] );

    # plot simple function
    fig1 = plt.figure(1,figsize=(10,6));
    plt.axis([0, 1, -5, 5]);
    plt.plot( t, d, 'ko');
    plt.xlabel('time t (s)');
    plt.ylabel('d(t)');

    # qually-spaced values of t
    Np=101;
    tp = eda_cvec( np.linspace(0,1,Np) );

    # cubic interpolation
    f = ip.interp1d(t.ravel(),d.ravel(),axis=0,
            kind='cubic',fill_value='extrapolate');
    dp=eda_cvec(f(tp));

    # plot simple function
    plt.plot(tp,dp,'k-');
    plt.show();

In [12]:

    # edapy_10_04: make correlated random time series

    # tunable parameters
    GAUSSIAN = 0; # Gaussian on 1, Oscillating exponential on 0
    MPLOT = 5;    # scale of model parameter plot
    CPLOT = 1;    # scale of covariance plot

    # standard time axis
    N=256;  # must be even, best if a power of 2
    Dt=1.0;
    tmax=Dt*(N-1);
    t=eda_cvec(np.linspace(0, tmax, N));
    tmid=tmax/2.0;

    # exact gamma*cos(w0T)exp(-sT) covariance function
    n=6;  # number of oscillations on (0,tmax) interal
    w0 = n*2*pi/tmax;  # angular frequency of oscilation
    s = 20.0/tmax;     # decay rate
    gamma = 1.0;       # sqrt variance
    gamma2 = gamma**2; #variance
    if( GAUSSIAN ):
        c = gamma2*np.exp(-0.5*s*np.power((t-tmid),2));
    else:
        c = gamma2*np.multiply( np.cos(w0*(t-tmid)), np.exp(-s*abs(t-tmid)) );

    print("gamma: %.4f s: %.4f" % (gamma, s) );

    # realize a time series with a specified covariance
    ctilde = np.fft.fft(c,axis=0);  # Fourier transform of c
    casd = np.sqrt(np.abs(ctilde)); # amp. spec. density of time series
    # make random numbers
    d = np.random.normal(loc=0.0,scale=1.0,size=(N,1));
    d = d - np.mean(d);            # remove mean
    dtilde = np.fft.fft(d,axis=0); # Fourier transform of d
    dasd = np.abs(dtilde);         # a.s.d. of d
    # adjust so has a.s.d. of c
    dtilde = np.multiply(np.divide(casd,dasd),dtilde);
    d = np.real(np.fft.ifft(dtilde,axis=0)); # inverse Fourier transform
    vard = np.var(d);     # variance of d
    d = (gamma/sqrt(vard))*d; # normalize to variance gamma^2

    # estimated covariance
    cest = eda_cvec(np.correlate(d.ravel(),d.ravel(),mode='full') )/N;

    # plot this covariance
    fig1 = plt.figure(1,figsize=(10,6));
    plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);
    plt.xlabel('time lag t (s)');
    plt.ylabel('c(t)');
    plt.plot( t-tmid,c,'k-',lw=3);
    tc = Dt*eda_cvec(np.linspace(-N+1,N-1,2*N-1));
    plt.plot( tc,cest,'k:','LineWidth',3);

    # plot time series
    fig2 = plt.figure(2,figsize=(10,6));
    plt.axis([ 0, tmax, -MPLOT, MPLOT]);
    plt.xlabel('time t (s)');
    plt.ylabel('d(t)');
    plt.plot( t,d,'k-',lw=3);
    plt.show();

    gamma: 1.0000 s: 0.0784

In [11]:

    # edapy_10_05, Gaussian Process Regression example
    # this example uses synthetic data
    # and a gamma2*cos(w0T)*exp(-sT) covariance function,
    # with T = |t(i)-t(j)|

    # tunable parameters
    MPLOT = 5;    # scale of model parameter plot
    CPLOT = 1;    # scale of covariance plot
    gamma = 1;    # sqrt(variance) of covariance
    sigmad = 0.2; # noise level of observed data
    wbias=1.0;    # compute GPR with wrong w0
    sbias=1.0;    # compute GPR with wrong s
    PLOTCONF=1;   # plot confidence intervals on 1
    USEFORRLOOP=0; # for loops on 1, outer products on 0

    # axes definitons

    # model parameter axis, tm, with sampling Dt
    M=256;  # must be even, best if a power of 2
    Dtm=1.0;
    tmax=Dtm*(M-1);
    tm=Dtm*eda_cvec(np.linspace(0,tmax,M));

    # finely sampled model parameter vector, with sampling
    # Dt/K, where K is a smallish power of 2, is used in
    # the initialization process, as described below
    K=8; # must be even, best if a power of 2
    Mf=M*K;
    Dtf=Dtm/K;
    tf=Dtf*eda_cvec(np.linspace(0,Mf-1,Mf));

    # Part 1: realization of a finely sampled time
    # % series, mf, with the desired covariance

    # PART 1: create true "fine" covariance function

    n=6;  # number of oscillations on (0,tmax) interal
    w0 = n*2*pi/tmax;  # angular frequency of oscilation
    tmid=tmax/2.0;
    s = 20.0/tmax;     # decay rate
    gamma2 = gamma**2; # variance
    cf = gamma2*np.multiply(np.cos(w0*(tf-tmid)),np.exp(-s*np.abs(tf-tmid)));

    print("gamma: %.4f w0: %.4f s: %.4f" % (gamma, w0, s) );

    # PART 2: create "fine" time series with thus
    # covariance using Fourier techniques

    # amplitude spectral density of true covariance
    cftilde = np.fft.fft(cf,axis=0);
    cfasd = np.sqrt(np.abs(cftilde));

    # To create a time series with a specified covariannce,
    # start with a times series of uniformly-distrbuted
    # random noise, which has a uniform a.s.d., Fourier
    # transorm it, multiply teh transform by the a.s.d
    # of the specified covariance, trasnform back, and adjust
    # the amplitude to the desired variance.
    # random numbers
    mf = np.random.normal(loc=0.0,scale=1.0,size=(Mf,1));
    mf = mf - np.mean(mf); # remove mean
    mftilde = np.fft.fft(mf,axis=0); # Fourier transform of mf
    mfasd = np.abs(mftilde);   # a.s.d. of mf
    # adjust so has a.s.d. of c
    mftilde = np.multiply(np.divide(cfasd,mfasd),mftilde);
    # inverse Fourier transform
    mf = np.real(np.fft.ifft(mftilde,axis=0));
    varmf = np.var(mf);      # variance of mf
    # normalize mf to variance gamma^2
    mf = (gamma/sqrt(varmf))*mf;

    # PART 3: Regularly sub-sample mf to get mtrue, irregularly
    # sub-sample mf to get true data dtrue, and add uncorrelated
    # noise to dtrue to get observed data dobs

    mtrue = np.copy(mf[0:Mf:K]);

    Ntarget = 50; # fewer observed data than true data
    # create list, j, of N unique, randomly chosen indices
    # that are not in the sequence 1, K+1, 2K+1, ...
    j=np.zeros((Ntarget,1),dtype=int);
    Nj=0;
    while (Nj<Ntarget):
        i=np.random.randint(low=0,high=Mf-1,size=(1));
        if( (((i-1)%K) != 0) and (not np.isin(i,j)) ):
            j[Nj,0]=i;
            Nj=Nj+1;
    # sort for plotting purposes
    j = eda_cvec(np.sort(j.ravel()));
    j = eda_cvec(np.unique(j.ravel()));
    N,i = np.shape(j);
    print('target N: %d, actual %d' % (Ntarget,N) );
    # times of observations
    td = np.copy(tf[j.ravel()]);
    # irregularly sub-sampled data
    dtrue = np.copy(mf[j.ravel()]);
    # observed data is true data + uncorrelated random noise
    sigmad2 = sigmad**2;
    dobs = dtrue+np.random.normal(loc=0,scale=sigmad,size=(N,1));

    # the data kernel, G, for the equation d=Gm
    G = np.concatenate( (np.zeros((N,M)), np.zeros((N,N))), axis=1 );

    # Since the true (w0, s) are in general unknown
    # this code allows them to be assigned differently
    # from the true values
    w0est = wbias*w0;  # w0est = 1.0*w0 is the true value
    sest = sbias*s;

    if( USEFORRLOOP ):
        # build covariances via for loop
        Ccc = np.zeros((N,N)); # covariance matrix, C(control,control)
        for i in range(N):
            for j in range(N):
                T = abs(td(i)-td(j));
                Ccc[i,j] = gamma2*cos(w0est*T)*exp(-sest*T);

        Ctc = np.zeros((M,N)); #  covariance matrix, C(target,control)
        for i in rane(M):
            for j in range(N):
                T = abs(tm(i)-td(j));
                Ctc[i,j] = gamma2*cos(w0est*T)*exp(-sest*T);
    else: # use outer product
        # build covariances via outer product
        Tdd = np.abs(np.matmul(td,np.ones((1,N)))-np.matmul(np.ones((N,1)),(td.T)));
        Ccc = gamma2*np.multiply(np.cos(w0est*Tdd),np.exp(-sest*Tdd));
        Tmd = np.abs(np.matmul(tm,np.ones((1,N)))-np.matmul(np.ones((M,1)),(td.T)));
        Ctc = gamma2*np.multiply(np.cos(w0est*Tmd),np.exp(-sest*Tmd));

    # solve GPE problem
    A = (Ccc + sigmad2*np.identity(N) );
    mest = np.matmul(Ctc,la.solve(A,dobs));

    # estimated covariance
    cest = eda_cvec(
        np.correlate(mest.ravel(),mest.ravel(),mode='full'))/M;

    # plot this covariance
    fig1 = plt.figure(1,figsize=(10,6));
    plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);
    plt.xlabel('time lag t (s)');
    plt.ylabel('c(t)');
    plt.plot( tf-tmid,cf,'k-',lw=3);
    tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));
    plt.plot(tc,cest,'k:',lw=3);
    plt.show();

    # time series figure
    fig2 = plt.figure(2,figsize=(10,6));
    plt.axis( [ 0, tmax, -MPLOT, MPLOT] );
    plt.xlabel('time t (s)');
    plt.ylabel('d(t) and  m(t)');
    plt.plot(tm,mtrue,'k-',lw=3);
    plt.plot(td,dtrue,'ko',lw=2);
    plt.plot(tm,mest,'-',lw=2,color=[0.8,0.8,0.8]);
    if( PLOTCONF ): # plot 95% confidence intervals
        # MATLAB: B = Ctc/A solves BA = Ctc
        # take transpose, so AT BT = CctT, then solve for BT
        # BT = solve(AT,CtcT) = solve(A,CtcT), since A=AT
        # MATLAB: B BT becomes BTT * BT
        BT = la.solve(A,Ctc.T)
        Cttest = sigmad2*np.matmul(BT.T,BT);
        sigmamest = np.sqrt(eda_cvec(np.diag(Cttest)));
        plt.plot(tm,mest-2*sigmamest,'k-',lw=1);
        plt.plot(tm,mest+2*sigmamest,'k-',lw=1);
    plt.show();

    gamma: 1.0000 w0: 0.1478 s: 0.0784
    target N: 50, actual 50

In [13]:

    # edapy_10_06, GPE using bicg()
    # this example uses synthetic data
    # and a gamma2*cos(w0T)*exp(-sT) covariance function,
    # with T = |t(i)-t(j)|

    # tunable parameters
    MPLOT = 5;    # scale of model parameter plot
    CPLOT = 1;    # scale of covariance plot
    gamma = 1;    # sqrt(variance) of covariance
    sigmad = 0.2; # noise level of observed data
    wbias=1.0;    # compute GPR with wrong w0
    sbias=1.0;    # compute GPR with wrong s
    PLOTCONF=1;   # plot confidence intervals on 1
    USEFORRLOOP=0; # for loops on 1, outer products on 0

    # axes definitons

    # model parameter axis, tm, with sampling Dt
    M=256;  # must be even, best if a power of 2
    Dtm=1.0;
    tmax=Dtm*(M-1);
    tm=Dtm*eda_cvec(np.linspace(0,tmax,M));

    # finely sampled model parameter vector, with sampling
    # Dt/K, where K is a smallish power of 2, is used in
    # the initialization process, as described below
    K=8; # must be even, best if a power of 2
    Mf=M*K;
    Dtf=Dtm/K;
    tf=Dtf*eda_cvec(np.linspace(0,Mf-1,Mf));

    # Part 1: realization of a finely sampled time
    # % series, mf, with the desired covariance

    # PART 1: create true "fine" covariance function

    n=6;  # number of oscillations on (0,tmax) interal
    w0 = n*2*pi/tmax;  # angular frequency of oscilation
    tmid=tmax/2.0;
    s = 20.0/tmax;     # decay rate
    gamma2 = gamma**2; # variance
    cf = gamma2*np.multiply(np.cos(w0*(tf-tmid)),np.exp(-s*np.abs(tf-tmid)));

    print("gamma: %.4f w0: %.4f s: %.4f" % (gamma, w0, s) );

    # PART 2: create "fine" time series with thus
    # covariance using Fourier techniques

    # amplitude spectral density of true time series
    cftilde = np.fft.fft(cf,axis=0);
    cfasd = np.sqrt(np.abs(cftilde));

    # To create a time series with a specified covariannce,
    # start with a times series of uniformly-distrbuted
    # random noise, which has a uniform a.s.d., Fourier
    # transorm it, multiply teh transform by the a.s.d
    # of the specified covariance, trasnform back, and adjust
    # the amplitude to the desired variance.
    # random numbers
    mf = np.random.normal(loc=0.0,scale=1.0,size=(Mf,1));
    mf = mf - np.mean(mf); # remove mean
    mftilde = np.fft.fft(mf,axis=0); # Fourier transform of mf
    mfasd = np.abs(mftilde);   # a.s.d. of mf
    # adjust so has a.s.d. of c
    mftilde = np.multiply(np.divide(cfasd,mfasd),mftilde);
    # inverse Fourier transform
    mf = np.real(np.fft.ifft(mftilde,axis=0));
    varmf = np.var(mf);      # variance of mf
    # normalize mf to variance gamma^2
    mf = (gamma/sqrt(varmf))*mf;

    # PART 3: Regularly sub-sample mf to get mtrue, irregularly
    # sub-sample mf to get true data dtrue, and add uncorrelated
    # noise to dtrue to get observed data dobs

    mtrue = np.copy(mf[0:Mf:K]);

    Ntarget = 50; # fewer observed data than true data
    # create list, j, of N unique, randomly chosen indices
    # that are not in the sequence 1, K+1, 2K+1, ...
    j=np.zeros((Ntarget,1),dtype=int);
    Nj=0;
    while (Nj<Ntarget):
        i=np.random.randint(low=0,high=Mf-1,size=(1));
        if( (((i-1)%K) != 0) and (not np.isin(i,j)) ):
            j[Nj,0]=i;
            Nj=Nj+1;
    # sort for plotting purposes
    j = eda_cvec(np.sort(j.ravel()));
    j = eda_cvec(np.unique(j.ravel()));
    N, i =np.shape(j);
    print('target N: %d, actual %d' % (Ntarget, N) );
    # times of observations
    td = np.copy(tf[j.ravel()]);
    # irregularly sub-sampled data
    dtrue = np.copy(mf[j.ravel()]);
    # observed data is true data + uncorrelated random noise
    sigmad2 = sigmad**2;
    dobs = dtrue+np.random.normal(loc=0,scale=sigmad,size=(N,1));

    # Since the true (w0, s) are in general unknown
    # this code allows them to be assigned differently
    # from the true values
    w0est = wbias*w0;  # w0est = 1.0*w0 is the true value
    sest = sbias*s;

    # set globals
    eda_N=N;
    eda_td=np.copy(td);
    eda_sigmad2=sigmad2;
    eda_gamma2=gamma2;
    eda_w0est=w0est;
    eda_sest=sest;

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

    # solve using congugate gradient algorithm
    u = np.zeros((N,1));
    q=las.bicg(LO,dobs,tol=1e-12,maxiter=(3*(N+M)+100));
    u[:,0] = q[0];
    mest = np.zeros((M,1));
    for i in range(M):
        mest[i,0]=0.0;
        for j in range(N):
            Tabs = abs(tm[i,0]-td[j,0]);
            Cij = gamma2*cos(w0est*Tabs)*exp(-sest*Tabs);
            mest[i,0]=mest[i,0]+Cij*u[j,0];

    dpre = np.zeros((N,1));
    for i in range(N):
        dpre[i,0]=0.0;
        for j in range(N):
            Tabs = abs(td[i,0]-td[j,0]);
            Cij = gamma2*cos(w0est*Tabs)*exp(-sest*Tabs);
            dpre[i,0]=dpre[i,0]+Cij*u[j,0];

    # estimated covariance
    cest = eda_cvec(
        np.correlate(mest.ravel(),mest.ravel(),mode='full'))/M;

    # plot this covariance
    fig1 = plt.figure(1,figsize=(10,6));
    plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);
    plt.xlabel('time lag t (s)');
    plt.ylabel('c(t)');
    plt.plot( tf-tmid,cf,'k-',lw=3);
    tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));
    plt.plot(tc,cest,'k:',lw=3);
    plt.show();

    # time series figure
    fig2 = plt.figure(2,figsize=(10,6));
    plt.axis( [ 0, tmax, -MPLOT, MPLOT] );
    plt.xlabel('time t (s)');
    plt.ylabel('d(t) and  m(t)');
    plt.plot(tm,mtrue,'k-',lw=3);
    plt.plot(td,dtrue,'ko',lw=2);
    plt.plot(tm,mest,'-',lw=2,color=[0.8,0.8,0.8]);
    plt.show();

    gamma: 1.0000 w0: 0.1478 s: 0.0784
    target N: 50, actual 50

In [4]:

    # edama_10_07, GPR Tuning Example

    # this example uses synthetic data
    # and a gamma2*cos(w0T)*exp(-sT) covariance function,
    # with T = |t(i)-t(j)|

    # tunable parameters
    MPLOT = 5;    # scale of model parameter plot
    CPLOT = 1;    # scale of covariance plot
    gamma = 1;    # sqrt(variance) of covariance
    sigmad = 0.2; # noise level of observed data
    wbias=3.0;    # compute GPR with wrong w0
    sbias=1.0;    # compute GPR with wrong s

    # axes definitons

    # model parameter axis, tm, with sampling Dt
    M=512;  # must be even, best if a power of 2
    Dtm=1.0;
    tmax=Dtm*(M-1);
    tm=eda_cvec(np.linspace(0,tmax,M));

    # finely sampled model parameter vector, with sampling
    # Dt/K, where K is a smallish power of 2, is used in
    # the initialization process, as described below
    K=8; # must be even, best if a power of 2
    Mf=M*K;
    Dtf=Dtm/K;
    tf=eda_cvec(np.linspace(0,tmax,Mf));

    # Part 1: realization of a finely sampled time
    # % series, mf, with the desired covariance

    # PART 1: create true "fine" covariance function

    n=6;  # number of oscillations on (0,tmax) interal
    w0 = n*2*pi/tmax;  # angular frequency of oscilation
    tmid=tmax/2.0;     # midpoint
    s = Dtm/60.0;      # decay rate
    gamma2 = gamma**2; # variance
    cf = gamma2*np.multiply(np.cos(w0*(tf-tmid)),np.exp(-s*np.abs(tf-tmid)));

    print("gamma: %.4f w0: %.4f s: %.4f" % (gamma, w0, s) );

    # PART 2: create "fine" time series with thus
    # covariance using Fourier techniques

    # amplitude spectral density of true time series
    cftilde = np.fft.fft(cf,axis=0);
    cfasd = np.sqrt(np.abs(cftilde));

    # To create a time series with a specified covariannce,
    # start with a times series of uniformly-distrbuted
    # random noise, which has a uniform a.s.d., Fourier
    # transorm it, multiply teh transform by the a.s.d
    # of the specified covariance, trasnform back, and adjust
    # the amplitude to the desired variance.
    # random numbers
    mf = np.random.normal(loc=0.0,scale=1.0,size=(Mf,1));
    mf = mf - np.mean(mf); # remove mean
    mftilde = np.fft.fft(mf,axis=0); # Fourier transform of mf
    mfasd = np.abs(mftilde);   # a.s.d. of mf
    # adjust so has a.s.d. of c
    mftilde = np.multiply(np.divide(cfasd,mfasd),mftilde);
    # inverse Fourier transform
    mf = np.real(np.fft.ifft(mftilde,axis=0));
    varmf = np.var(mf);      # variance of mf
    # normalize mf to variance gamma^2
    mf = (gamma/sqrt(varmf))*mf;

    # PART 3: Regularly sub-sample mf to get mtrue, irregularly
    # sub-sample mf to get true data dtrue, and add uncorrelated
    # noise to dtrue to get observed data dobs

    mtrue = np.copy(mf[0:Mf:K]);

    Ntarget = 150; # fewer observed data than true data
    # create list, j, of N unique, randomly chosen indices
    # that are not in the sequence 1, K+1, 2K+1, ...
    j=np.zeros((Ntarget,1),dtype=int);
    Nj=0;
    while (Nj<Ntarget):
        i=np.random.randint(low=0,high=Mf-1,size=(1));
        if( ((i%K) != 0) and (not np.isin(i,j)) ):
            j[Nj,0]=i;
            Nj=Nj+1;
    # sort for plotting purposes
    j = eda_cvec(np.sort(j.ravel()));
    j = eda_cvec(np.unique(j));
    N, i =np.shape(j);
    print('target N: %d, actual %d' % (Ntarget, N) );
    # times of observations
    td = np.copy(tf[j.ravel()]);
    # irregularly sub-sampled data
    dtrue = np.copy(mf[j.ravel()]);
    # observed data is true data + uncorrelated random noise
    sigmad2 = sigmad**2;
    sigmadr = 1.0/sigmad;
    dobs = dtrue+np.random.normal(loc=0,scale=sigmad,size=(N,1));

    # the data kernel, G, for the equation d=Gm
    G = np.concatenate( (np.zeros((N,M)), np.zeros((N,N))), axis=1 );

    # PART 4: Solve GPR problem for initial guess of parameter
    w0est = wbias*w0;  # w0est = 1.0*w0 is the true value
    sest = sbias*s;

    # build covariances via outer product
    Tdd = np.abs(np.matmul(td,np.ones((1,N)))-np.matmul(np.ones((N,1)),(td.T)));
    Ccc = gamma2*np.multiply(np.cos(w0est*Tdd),np.exp(-sest*Tdd));
    Tmd = np.abs(np.matmul(tm,np.ones((1,N)))-np.matmul(np.ones((M,1)),(td.T)));
    Ctc = gamma2*np.multiply(np.cos(w0est*Tmd),np.exp(-sest*Tmd));

    # solve GPE problem
    A = (Ccc + sigmad2*np.identity(N) );
    mest = np.matmul(Ctc,la.solve(A,dobs));

    # estimated covariance
    cest = eda_cvec(np.correlate(mest.ravel(),mest.ravel(),mode='full'))/M;

    # plot startng covariance
    fig1 = plt.figure(1,figsize=(10,6));
    plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);
    plt.xlabel('time lag t (s)');
    plt.ylabel('c(t)');
    plt.plot( tf-tmid,cf,'k-',lw=3);
    tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));
    plt.plot(tc,cest,'k:',lw=3);
    plt.show();

    # time series figure
    fig2 = plt.figure(2,figsize=(10,6));
    plt.axis( [ 0, tmax, -MPLOT, MPLOT] );
    plt.xlabel('time t (s)');
    plt.ylabel('d(t) and  m(t)');
    plt.plot(tm,mtrue,'k-',lw=3);
    plt.plot(td,dtrue,'ko',lw=2);
    plt.plot(tm,mest,'-',lw=2,color=[0.8,0.8,0.8]);
    plt.show();

    # Part 5: Set-up for Tuning

    # Gradient Descent Method adjustment of w0, now called q
    qtrue=w0;     # true parameter, q
    qg = w0est;   # inital guess of parameter, q
    Dq = 0.025;   # initial step size for gradiet descent
    NITER = 30;   # maximum number of iterations
    qg = qg-Dq;   # so first iteration starts qg
    alpha = 1.0;  # step size reducer
    nu = 1.0;     # sign of gradient
    mtpri = np.zeros((M,1)); # prior model parameters are zero
    mcpri = np.zeros((N,1)); # prior model parameters are zero
    # concatenate into overall prior vector mpri and times tall
    mpri = np.concatenate((mtpri,mcpri),axis=0); # prior model pars
    tall = np.concatenate((tm,td),axis=0); # corresponding times
    qsave = np.zeros((NITER,1)); # save estimates, for plot

    # Part 6: Tuning

    for iter in range(NITER): # loop over descent increments

        qg = qg+nu*alpha*Dq; # new guess of parameter, q

        # full covariance matrix and its derivative
        Tall = np.matmul(tall,np.ones((1,M+N))) - np.matmul(np.ones((M+N,1)),tall.T);
        Call = gamma2*np.multiply(np.cos(qg*Tall),np.exp(-sest*np.abs(Tall)));
        dCalldq = -gamma2*np.multiply(
                  np.multiply(Tall,np.sin(qg*Tall)),np.exp(-sest*abs(Tall)));

        # break out four quandrants; note these are views, not arrays,
        # so we need to be careful not to overwrite Ctt
        Ctt = Call[0:M,0:M];
        dCttdq = dCalldq[0:M,0:M];
        Ctc = Call[0:M,M:M+N];
        dCtcdq = dCalldq[0:M,M:M+N];
        Cct = Ctc.T;
        Ccc = Call[M:M+N,M:M+N];
        dCccdq = dCalldq[M:M+N,M:M+N];

        # A = Ccc + sigmad^2 I
        A = Ccc + sigmad2*np.identity(N);
        dAdq = dCccdq;

        # various functions of the covariance matrix
        CallI = la.inv(Call);
        CallIsr = la.sqrtm(CallI);
        dCallIdq = -np.matmul(np.matmul(CallI,dCalldq), CallI);
        dCallIsrdq = la.solve_sylvester( CallIsr, CallIsr, dCallIdq );
        R = la.cholesky(Call);
        logdetCall = 2.0*np.sum(np.log(np.diag(R)));
        dlogdetCalldq = np.trace(np.matmul(CallI,dCalldq));

        # GPR solution
        u=la.solve(A,dobs);
        mtg = np.matmul(Ctc,u);
        mcg = np.matmul(Ccc,u);
        mg = np.concatenate( (mtg,mcg), axis=0 );
        dg = eda_cvec(mg[M:M+N]);

        # derivative of predicted data
        v = la.solve( A, np.matmul(dAdq,u) );
        dmestdq = np.concatenate((
            np.matmul(dCtcdq,u) - np.matmul(Ctc,v),
            np.matmul(dCccdq,u) - np.matmul(Ccc,v)),
            axis=0 );
        ddpredq = np.copy(dmestdq[M:N+M]);
        detildedq = -sigmadr*ddpredq;
        e = dobs-dg;
        etilde = sigmadr*e;

        # derivative of total data error
        E = np.matmul(etilde.T,etilde);
        dEdp = 2.0*np.matmul(etilde.T,detildedq);

        # prior error and its derivative
        el=(mpri-mg);
        eltilde = np.matmul(CallIsr,el);
        deltildedq = np.matmul(dCallIsrdq,el)-np.matmul(CallIsr,dmestdq);

        # total prior error and its derivative
        L = np.matmul(eltilde.T,eltilde);
        dLdp = 2.0*np.matmul(eltilde.T,deltildedq);

        # function, PSI being minimized
        PSI = logdetCall + E + L;
        # and its deriative
        dPSIdp = dlogdetCalldq + dEdp + dLdp;
        nu = -eda_sign(dPSIdp);  # downhill vector

        # gradient method method
        if( iter==0 ):  # on the first iteration do nothing but save
            qglast = qg;
            PSIlast = PSI;
            print("1: qg %.4f nu %.2f alpha %e PSI %e PSIlast %e" % (qg,nu,alpha,PSI,PSIlast));
        elif (PSI < PSIlast ):
            qglast = qg;
            PSIlast = PSI;
            print("2: qg %.4f nu %.2f alpha %e PSI %e PSilast %e" % (qg,nu,alpha,PSI,PSIlast));
        else: # psi increased, so decrease alpha
            qg = qglast;
            PSI = PSIlast;
            alpha = alpha/2.0;
            print("2: qg %.4f nu %.2f alpha %e PSI %e PSIlast %e" % (qg,nu,alpha,PSI,PSIlast));
        qsave[iter,0] = qg;

    pest = qg;
    mest = mg;
    mtest = np.copy(mest[0:M,0]);
    mcest = np.copy(mest[M:M+N,0]);
    dpre = mcest;

    # estimated covariance
    cest = eda_cvec(
        np.correlate(mtest.ravel(),mtest.ravel(),mode='full'))/M;

    # plot final estimated covariance
    fig3 = plt.figure(3,figsize=(10,6));
    plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);
    plt.xlabel('time lag t (s)');
    plt.ylabel('c(t)');
    plt.plot( tf-tmid,cf,'k-',lw=3);
    tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));
    plt.plot(tc,cest,'k:',lw=3);
    plt.show();

    # plot final time series
    fig4 = plt.figure(4,figsize=(10,6));
    plt.axis( [ 0, tmax, -MPLOT, MPLOT] );
    plt.xlabel('time t (s)');
    plt.ylabel('d(t) and  m(t)');
    plt.plot(tm,mtrue,'k-',lw=3);
    plt.plot(td,dtrue,'ko',lw=2);
    plt.plot(tm,mtest,'-',lw=2,color=[0.8,0.8,0.8]);
    plt.show();

    print("qtrue %e qest %e relating error %e" % (qtrue, pest, (qtrue-pest)/qtrue) );

    # plot p as a function of iteration number
    fig5 = plt.figure(5,figsize=(10,6));
    plt.axis( [ 0, NITER, 0, 0.3] );
    plt.xlabel('iteration');
    plt.ylabel('w0');
    itvec = eda_cvec(np.linspace(1,NITER,NITER));
    plt.plot(itvec,qsave,'k-',lw=2);
    plt.plot(itvec,qsave,'ko',lw=2);
    plt.plot(itvec,w0*np.ones((NITER,1)),'k:',lw=2);

    gamma: 1.0000 w0: 0.0738 s: 0.0167
    target N: 150, actual 150

    1: qg 0.2213 nu -1.00 alpha 1.000000e+00 PSI -2.027642e+03 PSIlast -2.027642e+03
    2: qg 0.1963 nu -1.00 alpha 1.000000e+00 PSI -2.114795e+03 PSilast -2.114795e+03
    2: qg 0.1713 nu -1.00 alpha 1.000000e+00 PSI -2.195783e+03 PSilast -2.195783e+03
    2: qg 0.1463 nu -1.00 alpha 1.000000e+00 PSI -2.266060e+03 PSilast -2.266060e+03
    2: qg 0.1213 nu -1.00 alpha 1.000000e+00 PSI -2.321390e+03 PSilast -2.321390e+03
    2: qg 0.0963 nu -1.00 alpha 1.000000e+00 PSI -2.357814e+03 PSilast -2.357814e+03
    2: qg 0.0713 nu -1.00 alpha 1.000000e+00 PSI -2.371202e+03 PSilast -2.371202e+03
    2: qg 0.0713 nu 1.00 alpha 5.000000e-01 PSI -2.371202e+03 PSIlast -2.371202e+03
    2: qg 0.0713 nu -1.00 alpha 2.500000e-01 PSI -2.371202e+03 PSIlast -2.371202e+03
    2: qg 0.0713 nu 1.00 alpha 1.250000e-01 PSI -2.371202e+03 PSIlast -2.371202e+03
    2: qg 0.0713 nu -1.00 alpha 6.250000e-02 PSI -2.371202e+03 PSIlast -2.371202e+03
    2: qg 0.0698 nu 1.00 alpha 6.250000e-02 PSI -2.371216e+03 PSilast -2.371216e+03
    2: qg 0.0698 nu -1.00 alpha 3.125000e-02 PSI -2.371216e+03 PSIlast -2.371216e+03
    2: qg 0.0698 nu 1.00 alpha 1.562500e-02 PSI -2.371216e+03 PSIlast -2.371216e+03
    2: qg 0.0702 nu 1.00 alpha 1.562500e-02 PSI -2.371222e+03 PSilast -2.371222e+03
    2: qg 0.0702 nu -1.00 alpha 7.812500e-03 PSI -2.371222e+03 PSIlast -2.371222e+03
    2: qg 0.0702 nu 1.00 alpha 3.906250e-03 PSI -2.371222e+03 PSIlast -2.371222e+03
    2: qg 0.0703 nu 1.00 alpha 3.906250e-03 PSI -2.371222e+03 PSilast -2.371222e+03
    2: qg 0.0703 nu -1.00 alpha 3.906250e-03 PSI -2.371222e+03 PSilast -2.371222e+03
    2: qg 0.0703 nu 1.00 alpha 1.953125e-03 PSI -2.371222e+03 PSIlast -2.371222e+03
    2: qg 0.0703 nu -1.00 alpha 9.765625e-04 PSI -2.371222e+03 PSIlast -2.371222e+03
    2: qg 0.0703 nu -1.00 alpha 9.765625e-04 PSI -2.371222e+03 PSilast -2.371222e+03
    2: qg 0.0703 nu 1.00 alpha 4.882812e-04 PSI -2.371222e+03 PSIlast -2.371222e+03
    2: qg 0.0703 nu -1.00 alpha 2.441406e-04 PSI -2.371222e+03 PSIlast -2.371222e+03
    2: qg 0.0703 nu -1.00 alpha 2.441406e-04 PSI -2.371222e+03 PSilast -2.371222e+03
    2: qg 0.0703 nu 1.00 alpha 2.441406e-04 PSI -2.371222e+03 PSilast -2.371222e+03
    2: qg 0.0703 nu -1.00 alpha 1.220703e-04 PSI -2.371222e+03 PSIlast -2.371222e+03
    2: qg 0.0703 nu 1.00 alpha 6.103516e-05 PSI -2.371222e+03 PSIlast -2.371222e+03
    2: qg 0.0703 nu 1.00 alpha 6.103516e-05 PSI -2.371222e+03 PSilast -2.371222e+03
    2: qg 0.0703 nu -1.00 alpha 6.103516e-05 PSI -2.371222e+03 PSilast -2.371222e+03

    qtrue 7.377517e-02 qest 7.031538e-02 relating error 4.689643e-02

In [27]:

    # eda10_08, 2D interpolation of simulated pressure data
    # 2D interpolation of simulated pressure data
    # this script
    #             creates simulated data (Figue 1a)
    #             creates Delaunay triangle mesh (Figure 1b)
    #             finds triangle enclosing an arbitrary point (Figure 2)
    #             linear interpolation (Figure 3)

    DOCUBIC = 0;  # linear interpolation on 0, cubic on 1

    # rows of grid are x
    I = 41;      # number of x's in grid
    Dx = 1.0;    # sampling interval
    xmin = 0.0;  # start position
    xmax = 40.0; # end position
    # vector of x's
    x = xmin + (xmax-xmin)*eda_cvec(np.linspace(0,I-1,I)/(I-1));

    # columns of grid are y
    J = 41;      # number of y's in grid
    Dy = 1.0;    # sampling interval
    ymin = 0.0;  # start position
    ymax = 40.0; # end position
    # vector of y's
    y = ymin + (ymax-ymin)*eda_cvec(np.linspace(0,J-1,J)/(J-1));

    # grid 10% populated with data
    N0 = floor(I*J/10);  # number of data
    rowindex0=np.random.randint(I,size=(N0,1)); # row index of data
    colindex0=np.random.randint(J,size=(N0,1)); # col index of data

    # remove duplicate points, because the interpolation
    # function doesn't seem to like them, by checking
    # eack point against a matrix F of ones and zeros
    F = np.zeros((I,J),dtype=np.int);
    rowindex1=np.zeros((N0,1),dtype=np.int);
    colindex1=np.zeros((N0,1),dtype=np.int);
    N=0;
    for i in range(N0):
        irow = rowindex0[i,0];
        icol = colindex0[i,0];
        if ( F[irow,icol] == 0 ):
            rowindex1[N,0]=irow;
            colindex1[N,0]=icol;
            N=N+1;
            F[irow,icol] = 1;
    rowindex=eda_cvec(rowindex1[0:N,0]);
    colindex=eda_cvec(colindex1[0:N,0]);

    # corresponding x and y values of data
    xobs = x[rowindex,0];
    yobs = y[colindex,0];

    # create simulated pressure data (solution to Laplace's equation)
    kappa = 0.1;
    dtrue = 10.0*np.multiply( np.sin(kappa*xobs), np.exp(-kappa*yobs) );
    sigmad = 0.1;
    dobs = dtrue + np.random.normal(0.0,sigmad,(N,1) );

    # populate matrix A with observations
    A = np.zeros((I,J));
    for i in range(N):
        A[rowindex[i,0],colindex[i,0]] = dobs[i,0];

    # B&W colormap
    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);

    fig1 = plt.figure(1,figsize=(16,6));
    ax1 = plt.subplot(1,2,1);
    plt.axis([ymin, ymax, xmin, xmax]);
    left=ymin;
    right=ymax;
    bottom=xmin;
    top=xmax;
    plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );
    plt.xlabel('y');
    plt.ylabel('x');
    plt.title('data');

    # Delaunay triangulation
    XY = np.concatenate( (xobs,yobs), axis=1 );
    # XY is a N by 2 array of (x,y) positions of vertices
    mytri = sp.Delaunay( XY  );
    # TRI is a NTRI by 3 array of the indices of
    # the vertices associated with each triangle
    TRI = mytri.simplices;
    [NTRI, i] = np.shape(TRI);

    # plot
    ax1 = plt.subplot(1,2,2);
    plt.axis([ymin, ymax, xmin, xmax]);
    left=ymin;
    right=ymax;
    bottom=xmin;
    top=xmax;
    plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );
    for i in range(NTRI):
        v1=TRI[i,0];
        v2=TRI[i,1];
        v3=TRI[i,2];
        plt.plot( [XY[v1,1], XY[v2,1]], [XY[v1,0], XY[v2,0]], 'k-');
        plt.plot( [XY[v2,1], XY[v3,1]], [XY[v2,0], XY[v3,0]], 'k-');
        plt.plot( [XY[v3,1], XY[v1,1]], [XY[v3,0], XY[v1,0]], 'k-');
    plt.xlabel('y');
    plt.ylabel('x');
    plt.title('triangles');
    plt.show();

    fig2 = plt.figure(2,figsize=(16,6));
    ax1 = plt.subplot(1,1,1);
    plt.axis([ymin, ymax, xmin, xmax]);
    left=ymin;
    right=ymax;
    bottom=xmin;
    top=xmax;
    plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );
    plt.xlabel('y');
    plt.ylabel('x');
    plt.title('one triangle enclosing a specified point');

    # find triangle enclosing  point x0
    x0=[20.5,30.5];
    tri0 = mytri.find_simplex( x0 );

    k = tri0;
    if( k != -1 ):
        plt.plot( x0[1], x0[0], 'ko');
        v1=TRI[k,0];
        v2=TRI[k,1];
        v3=TRI[k,2];
        plt.plot( [XY[v1,1], XY[v2,1]], [XY[v1,0], XY[v2,0]], 'k-');
        plt.plot( [XY[v2,1], XY[v3,1]], [XY[v2,0], XY[v3,0]], 'k-');
        plt.plot( [XY[v3,1], XY[v1,1]], [XY[v3,0], XY[v1,0]], 'k-');
    plt.show();

    # grid of (x,y) values
    XGRID = np.matmul( x, np.ones((1,J)) );
    YGRID = np.matmul( np.ones((I,1)), y.T );

    # interpolation
    if( DOCUBIC):
        # linear interpolation
        BB=ip.griddata( np.concatenate((xobs,yobs),axis=1),
                       dobs,(XGRID,YGRID), method='cubic', fill_value=0.0);
        B = BB[:,:,0];
    else:
        # cubic intepolation
        CC=ip.griddata( np.concatenate((xobs,yobs),axis=1),
                       dobs,(XGRID,YGRID), method='linear', fill_value=0.0);
        C = BB[:,:,0];
        B = C;

    # plot
    fig3 = plt.figure(3,figsize=(16,6));
    ax1 = plt.subplot(1,2,1);
    plt.axis([ymin, ymax, xmin, xmax]);
    left=ymin;
    right=ymax;
    bottom=xmin;
    top=xmax;
    plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );
    plt.xlabel('y');
    plt.ylabel('x');
    plt.title('data');
    ax1 = plt.subplot(1,2,2);
    plt.axis([ymin, ymax, xmin, xmax]);
    left=ymin;
    right=ymax;
    bottom=xmin;
    top=xmax;
    plt.imshow( np.flipud(B), cmap=bwcmap, vmin=np.min(B), vmax=np.max(B), extent=(left,right,bottom,top) );
    plt.xlabel('y');
    plt.ylabel('x');
    plt.title('interpolated data');
    plt.show();

In [ ]:

    # edapy_10_09, is not provided, since the MATHLAB
    # version edama_10_09 addresses an issue that does
    # not arise in Python

In [5]:

    # eda10_10, two-dimensional Fourier transform
    # illustrate two-dimensional fourier transform
    # Part 1, FFT of cosine wave, reordered so zero us in center
    # Part 2, Illustration of use of symmetry

    # standard distance/wavenumber set up
    Nx=32;               # number of x's, must be even
    Nxo2 = floor(Nx/2);  # integer version of Nx/2
    Nkx=Nxo2+1;          # number of non-negative x-wavenumbers kx
    Dx=1.0;              # sampling interval in x
    x = eda_cvec(Dx*np.linspace(0,Nx-1,Nx)); # vector of x's
    xmax = Dx*(Nx-1);    # ending x
    kxmax=2*pi/(2*Dx);   # nyquist kx
    Dkx=kxmax/(Nx/2);    # sampling in kx
    # full (pos and neg) x-wavenumber vector
    kx=eda_cvec( Dkx * np.concatenate( ( np.linspace(0,Nxo2,Nkx), np.linspace(-Nxo2+1,-1,Nxo2-1) ), axis=0 ) );

    Ny=32;               # number of y's, must be even
    Nyo2 = floor(Ny/2);  # integer version of Ny/2
    Nky=Nyo2+1;          # number of non-negative y-wavenumbers ky
    Dy=1.0;              # sampling interval in y
    y = eda_cvec(Dy*np.linspace(0,Ny-1,Ny)); # vector of y's
    ymax = Dy*(Ny-1);    # ending y
    kymax=2*pi/(2*Dy);   # nyquist ky
    Dky=kymax/(Ny/2);    # sampling in ky
    # full (pos and neg) y-wavenumber vector
    ky=eda_cvec(Dky * np.concatenate( ( np.linspace(0,Nyo2,Nky), np.linspace(-Nyo2+1,-1,Nyo2-1) ), axis=0 ));

    # Part 1: cosine wave with radial wavenumber kr1 in the theta diection
    theta = (pi/180.0)*30.0;
    t1 = eda_cvec( [cos(theta), sin(theta)] );
    kr1 = kxmax/4.0;
    F = np.zeros((Nx,Ny));
    for n in range(Nx):
        for m in range(Ny):
            F[n,m] = cos( kr1*t1[0,0]*x[n,0] + kr1*t1[1,0]*y[m,0] );

    eda_draw('title f(x,y)', F);

    # amplitude spectral density
    Ftt = np.fft.fft2(F);
    norm = ((Dx**2)/xmax) * ((Dy**2)/ymax);
    Fasd = norm * np.sqrt( np.power( np.abs(Ftt), 2));

    # put into more reasonable order, with zero in the center
    Fasd2=np.zeros((Nx,Ny));
    Fasd2[Nx-Nkx:Nx,Ny-Nky:Ny] = Fasd[0:Nkx,0:Nky];
    Fasd2[0:Nkx-2,Ny-Nky:Ny] = Fasd[Nkx:Nx,0:Nky];
    Fasd2[0:Nkx-2,0:Nky-2] = Fasd[Nkx:Nx,Nky:Ny];
    Fasd2[Nkx-2:Nx,0:Nky-2] = Fasd[0:Nkx,Nky:Ny];

    # Part 2: example of the use of symmetry

    # new function F is noise, which has a complicated Fourier Transform
    F=np.random.normal(0.0,1.0,(Nx,Ny) );
    Ftt = np.fft.fft2(F);

    # for real data, use symmetries to build right hand columns of
    # 2D Fourier Transform from left hand columns
    # First copy left Ny/2+1 columns; they are treated as known
    Ftt2=np.zeros((Nx,Ny),dtype=np.complex);
    Ftt2[:,0:Nky] = Ftt[:,0:Nky];
    # then rebuild right from left columns
    for m in range(1,Nyo2):
        mp=Ny-m;
        Ftt2[0,mp] = np.conj(Ftt2[0,m]);
        Ftt2[1:Nx,mp] = np.flipud(np.conj(Ftt2[1:Nx,m]));

    # invert back to F
    F2 = np.real(np.fft.ifft2(Ftt2));

    eda_draw('Part 1', 'title ASD', Fasd, 'title ASD rearranged', Fasd2, 'Part 2', 'title F2', F, 'title reconstructed F2', F2);

    # compute error, as a check
    Ett=np.sum(np.abs(Ftt2-Ftt));
    E = np.sum(np.abs(F2-F));
    print('Symmetry test:');
    print("    error in transform %.4f" % (Ett) );
    print("    error in function  %.4f" % (E)    );

    Symmetry test:
        error in transform 0.0000
        error in function  0.0000

In [34]:

    # edam_10_11: Filling in gaps with fft using POCS

    # Part 1, generate true synthetic data set
    # with Gaussian autocorrelation function
    # that simulates a crudely-layerd medium

    # standard distance/wavenumber set up
    Nx=128;              # number of x's, must be even
    Nxo2 = floor(Nx/2);  # integer version of Nx/2
    Nkx=Nxo2+1;          # number of non-negative x-wavenumbers kx
    Dx=1.0;              # sampling interval in x
    x = eda_cvec(Dx*np.linspace(0,Nx-1,Nx)); # vector of x's
    xmax = Dx*(Nx-1);    # ending x
    kxmax=2*pi/(2*Dx);   # nyquist kx
    Dkx=kxmax/(Nx/2);    # sampling in kx
    # full (pos and neg) x-wavenumber vector
    kx=eda_cvec( Dkx * np.concatenate( ( np.linspace(0,Nxo2,Nkx), np.linspace(-Nxo2+1,-1,Nxo2-1) ), axis=0 ) );

    Ny=128;               # numer of y's, must be even
    Nyo2 = floor(Ny/2);  # integer version of Ny/2
    Nky=Nyo2+1;          # number of non-negative y-wavenumbers ky
    Dy=1.0;              # sampling interval in y
    y = eda_cvec(Dy*np.linspace(0,Ny-1,Ny)); # vector of y's
    ymax = Dy*(Ny-1);    # ending y
    kymax=2*pi/(2*Dy);   # nyquist ky
    Dky=kymax/(Ny/2);    # sampling in ky
    # full (pos and neg) y-wavenumber vector
    ky=eda_cvec(Dky * np.concatenate( ( np.linspace(0,Nyo2,Nky), np.linspace(-Nyo2+1,-1,Nyo2-1) ), axis=0 ));

    # wavenumber variances
    # they scale a reciprocal of spatial variances
    # values are chosen to make quasi-layers
    sx = 5.0*Dkx;
    sy = 1.0*Dky;
    sx2 = sx**2;
    sy2 = sy**2;
    asd = np.zeros((Nx,Ny)); # target amplitude spectral density
    for ikx in range(Nx):
        for iky in range(Ny):
            psd=exp(-0.5*(kx[ikx,0]**2)/sx2-0.5*(ky[iky,0]**2)/sy2);
            asd[ikx,iky] = sqrt(psd);

    # start with random data and
    mtrue = np.random.normal(loc=0.0,scale=1.0,size=(Nx,Ny));
    Ftt = np.fft.fft2(mtrue);   # take its fft
    Ftt = np.multiply(asd,Ftt); # shape to target
    mtrue = np.real(np.fft.ifft2(Ftt)); # inverse fft
    sigmamtrue = sqrt(np.var(mtrue)); # true variance

    # normalization ratio, gives target asd the same
    # variance as the fft of the data
    R=sqrt(np.var(np.abs(Ftt))) / sqrt(np.var(np.abs(asd)));
    Rasd=R*asd;

    # sample mtrue along vertical profiles
    j = np.floor( eda_cvec( 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 )*Ny/8.0).astype(np.int);
    Nd,i = np.shape(j);
    d = np.zeros((Nx,Nd));
    for i in range(Nd):
        k = j[i,0];
        d[0:Nx,i:i+1]=mtrue[0:Nx,k:k+1];

    print('profiles through true model');
    eda_draw(d[0:Nx,0:1],d[0:Nx,1:2],d[0:Nx,2:3],d[0:Nx,3:4],d[0:Nx,4:5],d[0:Nx,5:6],d[0:Nx,6:7]);

    # Part 2, POCS method to fill in data gaps

    # make starting model that has the right a.s.d. and variance
    mest = np.random.normal(loc=0.0,scale=1.0,size=(Nx,Ny)); # random data
    Ftt = np.fft.fft2(mest); # take fft
    Ftt=np.multiply(asd,Ftt); # shape to target
    mest = np.real(np.fft.ifft2(Ftt)); # take inverse fft
    sigmamest = sqrt(np.var(mest)); # estimated variance
    mest = (sigmamtrue/sigmamest)*mest; # scale to right variance
    mstarting = np.copy(mest); # save starting model

    # POCS algorithm
    Niter=50;
    E=np.zeros((Niter,1));
    for iter in range(Niter): # loop over iterations
        # enforce data are satisfied
        for i in range(Nd):
            k = j[i,0];
            mest[0:Nx,k:k+1]=d[0:Nx,i:i+1];
        # enforce a.s.d. is not greater than target
        Ftt = np.fft.fft2(mest);
        for ikx in range(Nx): # loop over wavenumbers
            for iky in range(Ny):
                a=Rasd[ikx,iky];
                b=abs(Ftt[ikx,iky]);
                # if asd too big,
                # make it smaller, w/o altering phase
                if( b>a ):
                        Ftt[ikx,iky]=(a/b)*Ftt[ikx,iky];
        mest = np.real(np.fft.ifft2(Ftt));
        # model error
        e = eda_cvec( mtrue.ravel()-mest.ravel() );
        E[iter,0] = np.matmul(e.T,e);

    sigmamtrue=sqrt(np.var(mtrue));
    sigmamstarting=sqrt(np.var(mstarting));
    sigmamest=sqrt(np.var(mest));
    print('sqrt variances: starting %e, estimated %e, true %e' % (sigmamtrue, sigmamstarting, sigmamest));
    print('starting, estimated, true model');
    eda_draw('title starting',mstarting, 'title est', mest,'title true', mtrue);

    fig2 = plt.figure(2,figsize=(16,6));
    ax1 = plt.subplot(1,1,1);
    plt.axis([0, Niter, 0, E[0,0] ]);
    plt.plot( eda_cvec(np.linspace(1,Niter,Niter)), E, 'k-', lw=2);
    plt.xlabel('iteration, i');
    plt.ylabel('model error, E')
    plt.show();

    profiles through true model

    sqrt variances: starting 4.498101e-02, estimated 4.498101e-02, true 3.422609e-02
    starting, estimated, true model

In [ ]: