{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eda10_00 clear all variables and import vatious modules\n",
    "%reset -f\n",
    "import os\n",
    "from datetime import date\n",
    "from math import exp, pi, sin, cos, tan, sqrt, floor, ceil, log\n",
    "import numpy as np\n",
    "import scipy.sparse.linalg as las\n",
    "import scipy.interpolate as ip\n",
    "import scipy.spatial as sp\n",
    "from scipy import sparse\n",
    "import scipy.linalg as la\n",
    "import scipy.signal as sg\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "\n",
    "# eda_draw function makes a \"pictorial matrix equation\"\n",
    "# arguments are vectors, matrices and strings\n",
    "# which are plotted in the order that the appear\n",
    "# except that strings starting with 'title ' are plotted\n",
    "# under the subseqeunt matrix or vector\n",
    "# always returns a status of 1\n",
    "def eda_draw(*argv):\n",
    "    bw = np.zeros((256,4));\n",
    "    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "    bw[:,0] = v;\n",
    "    bw[:,1] = v;\n",
    "    bw[:,2] = v;\n",
    "    bw[:,3] = np.ones(256);\n",
    "    bwcmap = ListedColormap(bw);\n",
    "    # size of plot\n",
    "    W = 16;\n",
    "    H = 4;\n",
    "    fig1 = plt.figure(1);\n",
    "    # figsize width and height in inches\n",
    "    fig1.set_size_inches(W,H);\n",
    "    ax1 = plt.subplot(1,1,1);\n",
    "    plt.axis([0, W, -H/2, H/2]);\n",
    "    plt.axis('off');\n",
    "    LM = W/6;    # matrix width and heoght\n",
    "    LV = W/40;   # vector width\n",
    "    FS = 0.12;    # character width\n",
    "    TO = 0.4;    # title vertical offset\n",
    "    SP = 0.2;    # space between objects\n",
    "    LS = 0.2;    # leading space\n",
    "    p = LS; # starting x-position\n",
    "    istitle=0; # flags presence of a title\n",
    "    for a in argv:\n",
    "        if isinstance(a,np.ndarray):\n",
    "            sh = np.shape(a);\n",
    "            if len(sh) == 1:  # conversion to nx1 array\n",
    "                n = sh[0];\n",
    "                m = 1;\n",
    "                ap = a;\n",
    "                a = np.zeros((n,1));\n",
    "                a[:,0] = ap;\n",
    "            else:\n",
    "                n = sh[0];\n",
    "                m = sh[1];\n",
    "            if m==1:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LV;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LV;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "            else:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LM;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LM;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "        elif isinstance(a,str):\n",
    "            ns = len(a);\n",
    "            istitle=0;\n",
    "            if( ns>=6 ):\n",
    "                if 'title ' in a[0:6]:\n",
    "                    istitle=1;\n",
    "                    titlestr=a[6:];\n",
    "            if( istitle != 1):\n",
    "                plt.text(p,0,a);\n",
    "                p = p + ns*FS + SP;\n",
    "    plt.show();\n",
    "    return 1;\n",
    "\n",
    "\n",
    "def FTFmul(v):\n",
    "    # this function is used by the bicongugate gradient solver to solve the geneneralized least\n",
    "    # squares problem Fm=f.  Note that \"F\" must be called \"edaFsparse\".\n",
    "    global edaFsparse;\n",
    "    s = np.shape(v);\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    temp = edaFsparse*vv;\n",
    "    return edaFsparse.transpose()*temp;\n",
    "\n",
    "\n",
    "def GLSFilterMul(v):\n",
    "    # this function is used by the bicongugate gradient solver to solve the\n",
    "    # geneneralized least squares problem Fm=f with F=[G;H] and f=[d,h],\n",
    "    # where G is a topplitz matrix built from the filter g\n",
    "    # Note that \"H\" must be sparese and called \"edaHsparse\" and that\n",
    "    # the filter \"g\" bust be called edafilterg and must be a column vector\n",
    "    global edaHsparse;\n",
    "    global edafilterg;\n",
    "    # rearrange as column-vector\n",
    "    s = np.shape(v);\n",
    "    M = s[0];\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    N, i = np.shape(edafilterg);\n",
    "    # implements (GT G + HT H) v\n",
    "    temp1 = np.zeros((N+M-1,1));\n",
    "    temp1[:,0] = np.convolve(edafilterg.ravel(),vv.ravel()); # G v is of length N\n",
    "    a = np.zeros((N,1));\n",
    "    a[:,0] = temp1[0:N,0];\n",
    "    b = edaHsparse * v; # H v is of length K\n",
    "    temp2 = np.zeros((N+N-1,1));\n",
    "    temp2[:,0] = np.convolve((np.flipud(edafilterg)).ravel(),a.ravel()); # GT (G v) is of length M\n",
    "    a2 = temp2[N-1:N+M-1,0];\n",
    "    b2 = edaHsparse.transpose()*b; # HT (H v) is of length M\n",
    "    # FT F v = GT G v + HT H v\n",
    "    return (a2+b2);\n",
    "\n",
    "\n",
    "def GLSautoMul(v):\n",
    "    # this function is used by the bicongugate gradient solver to solve the\n",
    "    # geneneralized least squares problem Fm=f with F=[G;H] and f=[d,h],\n",
    "    # where GTG is built from the autocorrelation of the filter g\n",
    "    global edaHsparse;\n",
    "    global edaautoc;\n",
    "    # rearrange as column-vector\n",
    "    s = np.shape(v);\n",
    "    M = s[0];\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    # FT F v = GT G v + HT H v\n",
    "    GTGv=np.zeros((M,1));\n",
    "    for i in range(M):\n",
    "        u = np.concatenate( (np.flipud(edaautoc[0:i+1,0]), edaautoc[1:M-i,0]), axis=0 );\n",
    "        GTGv[i,0] =  np.matmul( u.T, vv );\n",
    "    Hv = edaHsparse*vv;\n",
    "    HTHv = edaHsparse.transpose()*Hv;\n",
    "    return (GTGv + HTHv);\n",
    "\n",
    "\n",
    "def chebyshevfilt(d, Dt, flow, fhigh):\n",
    "    #chebyshevfilt\n",
    "    # chebyshev IIR bandpass filter\n",
    "    # d - input array of data\n",
    "    # Dt - sampling interval\n",
    "    # flow - low pass frequency, Hz\n",
    "    # fhigh - high pass frequency, Hz\n",
    "    # dout - output array of data\n",
    "    # u - the numerator filter\n",
    "    # v - the denominator filter\n",
    "    # these filters can be used again using dout=filter(u,v,din);\n",
    "\n",
    "    # make sure input timeseries is a column vector\n",
    "    s = np.shape(d);\n",
    "    N = s[0];\n",
    "    if(N==1):\n",
    "        dd = np.zeros((N,1));\n",
    "        dd[:,0] = d;\n",
    "    else:\n",
    "        dd=d;\n",
    "        \n",
    "    # sampling rate\n",
    "    rate=1/Dt;\n",
    "\n",
    "    # ripple parameter, set to ten percent\n",
    "    ripple=0.1;  \n",
    "\n",
    "    # normalise frequency\n",
    "    fl=2.0*flow/rate;\n",
    "    fh=2.0*fhigh/rate;\n",
    "\n",
    "    # center frequency \n",
    "    cf = 4 * tan( (fl*pi/2) ) * tan( (fh*pi/2) );\n",
    "\n",
    "    # bandwidth\n",
    "    bw = 2 * ( tan( (fh*pi/2) ) - tan( (fl*pi/2) ) );\n",
    "\n",
    "    # ripple parameter factor\n",
    "    rpf = sqrt((sqrt((1.0+1.0/(ripple*ripple))) + 1.0/ripple));\n",
    "    a = 0.5*(rpf-1.0/rpf);\n",
    "    b = 0.5*(rpf+1.0/rpf);\n",
    "\n",
    "    u=np.zeros((5,1));\n",
    "    v=np.zeros((5,1));\n",
    "    theta = 3*pi/4;\n",
    "    sr = a * cos(theta);\n",
    "    si = b * sin(theta);\n",
    "    es = sqrt(sr*sr+si*si);\n",
    "    tmp= 16.0 - 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw - 4.0*bw*cf*sr + cf*cf;\n",
    "    v[0,0] = 1.0;\n",
    "    v[1,0] = 4.0*(-16.0 + 8.0*bw*sr - 2.0*bw*cf*sr + cf*cf)/tmp;\n",
    "    v[2,0] = (96.0 - 16.0*cf - 8.0*es*es*bw*bw + 6.0*cf*cf)/tmp;\n",
    "    v[3,0] = (-64.0 - 32.0*bw*sr + 8.0*bw*cf*sr + 4.0*cf*cf)/tmp;\n",
    "    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;\n",
    "    tmp = 4.0*es*es*bw*bw/tmp;\n",
    "    u[0,0] = tmp;\n",
    "    u[1,0] = 0.0;\n",
    "    u[2,0] = -2.0*tmp;\n",
    "    u[3,0] = 0.0;\n",
    "    u[4,0] = tmp;\n",
    "\n",
    "    dout = sg.lfilter(u.ravel(),v.ravel(),dd.ravel());\n",
    "    return (dout,u,v);\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec1(v) from entity v that is array-like,\n",
    "# including a number, a list of numbers, a tuple\n",
    "# of numbers, an N-by-0 np array and a N-by1 np array.\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it.\n",
    "# This version takes just one input argment.\n",
    "# whereas c=eda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def eda_cvec1(v):\n",
    "    if isinstance(v, int):\n",
    "        w = np.zeros((1,1),dtype=int);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, float):\n",
    "        w = np.zeros((1,1),dtype=float);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, complex):\n",
    "        w = np.zeros((1,1),dtype=complex);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, np.ndarray):\n",
    "        s = np.shape(v);\n",
    "        if len(s) == 1:\n",
    "            return np.copy(np.reshape(v,(s[0],1)));\n",
    "        else:\n",
    "            [r,c]=s;\n",
    "            if( c==1 ):\n",
    "                return(np.copy(v));\n",
    "            elif(r==1):\n",
    "                return(np.copy(v.T));\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: %d by %d ndarray not allowed\" % (r, c));\n",
    "    elif isinstance(v, list):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, tuple):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"eda_cvec: %s not supported\" % type(v));\n",
    "        \n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec(v1, v2, ...) from a list of several\n",
    "# array-like entities v1, v2, including a number\n",
    "# a list of numbers, a tuple of numbers, an N-by-0 np array\n",
    "# and a N-by-1 np array. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context.\n",
    "# This version concatenates many argments,\n",
    "# whereas c=eda_cvec1(v1) takes just one argiment.\n",
    "# I recommend always using eda_cvec(v1, v2, ...)\n",
    "def eda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        Nt = Nt + N;\n",
    "        if (t==int) and isinstance(v[0,0],float):\n",
    "            t=float;\n",
    "        elif isinstance(v[0,0],complex):\n",
    "            t=complex;\n",
    "    w = np.zeros((Nt,1),dtype=t);\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 22120 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda10_01\n",
    "# illustrate problem associated with interpolation by\n",
    "# fitting a high-order polynomial\n",
    "\n",
    "# times\n",
    "N=11;\n",
    "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] );\n",
    "\n",
    "# simple function\n",
    "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]);\n",
    "\n",
    "# plot simple function\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis([0, 1, -5, 5]);\n",
    "plt.plot( t, d, 'ko');\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t)');\n",
    "\n",
    "# linear representation of polynomial relationship\n",
    "G=np.zeros((N,N));\n",
    "for i in range(N):\n",
    "    G[:,i]=np.power(t,i).ravel();\n",
    "\n",
    "# solve for coefficients\n",
    "m = la.solve(G,d);\n",
    "\n",
    "# evaluate at equally-spaced values of t\n",
    "Np=101;\n",
    "tp = eda_cvec(np.linspace(0.0,1.0,Np));\n",
    "Gp = np.zeros((Np,N));\n",
    "for i in range(N):\n",
    "    Gp[:,i] = np.power( tp, i ).ravel();\n",
    "dp = np.matmul(Gp,m);\n",
    "\n",
    "# plot simple function\n",
    "plt.plot(tp,dp,'k-');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 10272 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda10_02\n",
    "# example of linear interpolation\n",
    "\n",
    "# times\n",
    "N=11;\n",
    "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] );\n",
    "\n",
    "# simple function\n",
    "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] );\n",
    "\n",
    "# plot simple function\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis([0, 1, -5, 5]);\n",
    "plt.plot( t, d, 'ko');\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t)');\n",
    "\n",
    "# qually-spaced values of t\n",
    "Np=101;\n",
    "tp = eda_cvec( np.linspace(0.0,1.0,Np) );\n",
    "\n",
    "# linear interpolation\n",
    "f = ip.interp1d(t.ravel(),d.ravel(),axis=0,kind='linear',fill_value='extrapolate');\n",
    "dp = eda_cvec(f(tp));\n",
    "\n",
    "# plot simple function\n",
    "plt.plot(tp,dp,'k-');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12548 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda10_03\n",
    "# example of  ubic spline interpolation\n",
    "\n",
    "# times\n",
    "N=11;\n",
    "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] );\n",
    "\n",
    "# simple function\n",
    "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] );\n",
    "\n",
    "# plot simple function\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis([0, 1, -5, 5]);\n",
    "plt.plot( t, d, 'ko');\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t)');\n",
    "\n",
    "# qually-spaced values of t\n",
    "Np=101;\n",
    "tp = eda_cvec( np.linspace(0,1,Np) );\n",
    "\n",
    "# linear interpolation\n",
    "f = ip.interp1d(t.ravel(),d.ravel(),axis=0,kind='cubic',fill_value='extrapolate');\n",
    "dp=eda_cvec(f(tp));\n",
    "\n",
    "# plot simple function\n",
    "plt.plot(tp,dp,'k-');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 33980 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda10_04\n",
    "# simple example of kriging\n",
    "\n",
    "# the model parameters are the values of a sinusoidal\n",
    "# function evenly spaced along the t axis\n",
    "M=101;   # length of time series\n",
    "Dt=1.0;  # sampling interval\n",
    "tmin=0.0;  # start time\n",
    "t = eda_cvec( tmin+Dt*np.linspace(0,M-1,M) ); # time vector\n",
    "tmax=t[M-1,0]; # end time\n",
    "\n",
    "# the true model parameters are sinusoidally varying\n",
    "freq=0.02;\n",
    "mtrue = np.sin(2*pi*freq*t);\n",
    "\n",
    "# the data are the values the function measured on only a\n",
    "# subset of randomply chosen points along the x axis\n",
    "N=40;\n",
    "rowindex = np.sort( np.random.randint(M,size=N) );\n",
    "tobs = eda_cvec( t[rowindex,:] );\n",
    "dobs = eda_cvec( mtrue[rowindex,:] );\n",
    "\n",
    "# estimate the model parameters at these points\n",
    "test = t;\n",
    "\n",
    "# plot the observed data with circles\n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "plt.axis([tmin, tmax, -2, 2]);\n",
    "plt.plot( t, mtrue, 'k-');\n",
    "plt.plot( tobs, dobs, 'ko');\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t) and m(t)');\n",
    "\n",
    "# assume the autocorrelation function of the data\n",
    "# is a a Normal function with variance L^2.  The normalization\n",
    "# constant cancels, so just set to unity\n",
    "\n",
    "# set the constant L arbitrarily\n",
    "L2 = 20;\n",
    "\n",
    "# basics equations are Aw=v(i), dest(i)=dobs'*w(i);\n",
    "\n",
    "# set up the kriging matrix, A\n",
    "wN = np.ones((N,1));\n",
    "wM = np.ones((M,1));\n",
    "A = np.exp( -np.power(np.abs(np.matmul(tobs,wN.T)-np.matmul(wN,tobs.T)), 2) / (2*L2) );\n",
    "\n",
    "# set up kriging vector, v, but with vector for each interplation\n",
    "# point a column of the matrix, V\n",
    "V = np.exp( -np.power(np.abs(np.matmul(tobs,wM.T)-np.matmul(wN,test.T)), 2) / (2*L2) );\n",
    "\n",
    "# solve, but add some damping, just in case A is near singular.  This selects the smallest w\n",
    "# that solves the equation\n",
    "e = 1.0e-6; # damping constant, to insure aganst sigular matrix\n",
    "mest = eda_cvec( np.matmul( dobs.T, la.solve(A+e*np.identity(N),V) ) );\n",
    "    \n",
    "# plot estimated data with line\n",
    "plt.plot(test,mest,'k-',color=(0.7,0.7,0.7));\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 93516 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 12960 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 16400 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda10_05\n",
    "# 2D interpolation of simulated pressure data\n",
    "# this script\n",
    "#             creates simulated data (Figue 1a)\n",
    "#             creates Delaunay triangle mesh (Figure 1b)\n",
    "#             finds triangle enclosing an arbitrary point (Figure 2)\n",
    "#             linear interpolation (Figure 3)\n",
    "\n",
    "# rows of grid are x\n",
    "I = 41;      # number of x's in grid\n",
    "Dx = 1.0;    # sampling interval\n",
    "xmin = 0.0;  # start position\n",
    "xmax = 40.0; # end position\n",
    "# vector of x's\n",
    "x = xmin + (xmax-xmin)*eda_cvec(np.linspace(0,I-1,I)/(I-1));\n",
    "\n",
    "# columns of grid are y\n",
    "J = 41;      # number of y's in grid\n",
    "Dy = 1.0;    # sampling interval\n",
    "ymin = 0.0;  # start position \n",
    "ymax = 40.0; # end position\n",
    "# vector of y's\n",
    "y = ymin + (ymax-ymin)*eda_cvec(np.linspace(0,J-1,J)/(J-1));\n",
    "\n",
    "# grid 10% populated with data\n",
    "N0 = floor(I*J/10);  # number of data\n",
    "rowindex0=np.random.randint(I,size=(N0,1)); # row index of data\n",
    "colindex0=np.random.randint(J,size=(N0,1)); # col index of data\n",
    "\n",
    "# remove duplicate points, because the interpolation\n",
    "# function doesn't seem to like them, by checking\n",
    "# eack point against a matrix F of ones and zeros\n",
    "F = np.zeros((I,J),dtype=np.int);\n",
    "rowindex1=np.zeros((N0,1),dtype=np.int);\n",
    "colindex1=np.zeros((N0,1),dtype=np.int);\n",
    "N=0;\n",
    "for i in range(N0):\n",
    "    irow = rowindex0[i,0];\n",
    "    icol = colindex0[i,0];\n",
    "    if ( F[irow,icol] == 0 ):\n",
    "        rowindex1[N,0]=irow;\n",
    "        colindex1[N,0]=icol;\n",
    "        N=N+1;\n",
    "        F[irow,icol] = 1;               \n",
    "rowindex=eda_cvec(rowindex1[0:N,0]);\n",
    "colindex=eda_cvec(colindex1[0:N,0]);\n",
    "\n",
    "# corresponding x and y values of data\n",
    "xobs = x[rowindex,0];\n",
    "yobs = y[colindex,0];\n",
    "\n",
    "# create simulated pressure data (solution to Laplace's equation)\n",
    "kappa = 0.1;\n",
    "dtrue = 10.0*np.multiply( np.sin(kappa*xobs), np.exp(-kappa*yobs) );\n",
    "sigmad = 0.1;\n",
    "dobs = dtrue + np.random.normal(0.0,sigmad,(N,1) );\n",
    "\n",
    "# populate matrix A with observations\n",
    "A = np.zeros((I,J));\n",
    "for i in range(N):\n",
    "    A[rowindex[i,0],colindex[i,0]] = dobs[i,0];\n",
    "\n",
    "# B&W colormap\n",
    "bw = np.zeros((256,4));\n",
    "v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "bw[:,0] = v;\n",
    "bw[:,1] = v;\n",
    "bw[:,2] = v;\n",
    "bw[:,3] = np.ones(256);\n",
    "bwcmap = ListedColormap(bw);\n",
    "    \n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis([ymin, ymax, xmin, xmax]);\n",
    "left=ymin;\n",
    "right=ymax;\n",
    "bottom=xmin;\n",
    "top=xmax;\n",
    "plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );\n",
    "plt.xlabel('y');\n",
    "plt.ylabel('x');\n",
    "plt.title('data');\n",
    "\n",
    "# triangulation\n",
    "XY = np.concatenate( (xobs,yobs), axis=1 );\n",
    "mytri = sp.Delaunay( XY  );\n",
    "\n",
    "# plot triangles\n",
    "TRI = mytri.simplices;\n",
    "[NTRI, i] = np.shape(TRI);\n",
    "\n",
    "# plot\n",
    "ax1 = plt.subplot(1,2,2);\n",
    "plt.axis([ymin, ymax, xmin, xmax]);\n",
    "left=ymin;\n",
    "right=ymax;\n",
    "bottom=xmin;\n",
    "top=xmax;\n",
    "plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );\n",
    "for i in range(NTRI):\n",
    "    v1=TRI[i,0];\n",
    "    v2=TRI[i,1];\n",
    "    v3=TRI[i,2];\n",
    "    plt.plot( [XY[v1,1], XY[v2,1]], [XY[v1,0], XY[v2,0]], 'k-');\n",
    "    plt.plot( [XY[v2,1], XY[v3,1]], [XY[v2,0], XY[v3,0]], 'k-');\n",
    "    plt.plot( [XY[v3,1], XY[v1,1]], [XY[v3,0], XY[v1,0]], 'k-');\n",
    "plt.xlabel('y');\n",
    "plt.ylabel('x');\n",
    "plt.title('triangles');\n",
    "plt.show();\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(16,6));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([ymin, ymax, xmin, xmax]);\n",
    "left=ymin;\n",
    "right=ymax;\n",
    "bottom=xmin;\n",
    "top=xmax;\n",
    "plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );\n",
    "plt.xlabel('y');\n",
    "plt.ylabel('x');\n",
    "plt.title('one triangle enclosing a specified point');\n",
    "\n",
    "# find triangle enclosing exemplary point and plot them\n",
    "x0=[20.5,30.5];\n",
    "k = mytri.find_simplex( x0 );\n",
    "if( k != -1 ):\n",
    "    plt.plot( x0[1], x0[0], 'ko');\n",
    "    v1=TRI[k,0];\n",
    "    v2=TRI[k,1];\n",
    "    v3=TRI[k,2];\n",
    "    plt.plot( [XY[v1,1], XY[v2,1]], [XY[v1,0], XY[v2,0]], 'k-');\n",
    "    plt.plot( [XY[v2,1], XY[v3,1]], [XY[v2,0], XY[v3,0]], 'k-');\n",
    "    plt.plot( [XY[v3,1], XY[v1,1]], [XY[v3,0], XY[v1,0]], 'k-'); \n",
    "plt.show();\n",
    "\n",
    "# 2D interpolation\n",
    "XGRID = np.matmul( x, np.ones((1,J)) );\n",
    "YGRID = np.matmul( np.ones((I,1)), y.T );\n",
    "# note: could also use method='cubic' for smoother interpolation\n",
    "BB=ip.griddata( np.concatenate((xobs,yobs),axis=1),dobs,(XGRID,YGRID), method='linear', fill_value=0.0);\n",
    "B = BB[:,:,0];\n",
    "\n",
    "# plot\n",
    "fig3 = plt.figure(3,figsize=(16,6));\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis([ymin, ymax, xmin, xmax]);\n",
    "left=ymin;\n",
    "right=ymax;\n",
    "bottom=xmin;\n",
    "top=xmax;\n",
    "plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );\n",
    "plt.xlabel('y');\n",
    "plt.ylabel('x');\n",
    "plt.title('data');\n",
    "ax1 = plt.subplot(1,2,2);\n",
    "plt.axis([ymin, ymax, xmin, xmax]);\n",
    "left=ymin;\n",
    "right=ymax;\n",
    "bottom=xmin;\n",
    "top=xmax;\n",
    "plt.imshow( np.flipud(B), cmap=bwcmap, vmin=np.min(B), vmax=np.max(B), extent=(left,right,bottom,top) );\n",
    "plt.xlabel('y');\n",
    "plt.ylabel('x');\n",
    "plt.title('interpolated data');\n",
    "plt.show();\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 98808 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 13252 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 16080 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda10_06\n",
    "# 2D interpolation of simulated pressure data\n",
    "# this script\n",
    "#             creates simulated data (Figue 1a)\n",
    "#             creates Delaunay triangle mesh (Figure 1b)\n",
    "#             finds triangle enclosing an arbitrary point (Figure 2)\n",
    "#             cubic spline interpolation (Figure 3)\n",
    "\n",
    "# rows of grid are x\n",
    "I = 41;       # number of x's\n",
    "Dx = 1.0;     # sampling intervak\n",
    "xmin = 0.0;   # starting x\n",
    "xmax = 40.0;  # ending x\n",
    "# vector of x's\n",
    "x = xmin + (xmax-xmin)*eda_cvec( np.linspace(0,I-1,I)/(I-1) );\n",
    "\n",
    "# columns of grid are y\n",
    "J = 41;       # number of y's\n",
    "Dy = 1.0;     # sampling interval\n",
    "ymin = 0.0    # starting y\n",
    "ymax = 40.0;  # ending y\n",
    "# vector of y's\n",
    "y = ymin + (ymax-ymin)*eda_cvec(np.linspace(0,J-1,J)/(J-1));\n",
    "\n",
    "# grid 10% populated with data\n",
    "N0 = floor(I*J/10);\n",
    "\n",
    "# randomly chosen points on grid\n",
    "rowindex0=np.random.randint(I,size=(N0,1));\n",
    "colindex0=np.random.randint(J,size=(N0,1));\n",
    "\n",
    "# remove duplicate points, because the interpolation\n",
    "# function doesn't seem to like them, by checking\n",
    "# eack point against a matrix F of ones and zeros\n",
    "F = np.zeros((I,J),dtype=np.int);\n",
    "rowindex1=np.zeros((N0,1),dtype=np.int);\n",
    "colindex1=np.zeros((N0,1),dtype=np.int);\n",
    "N=0;\n",
    "for i in range(N0):\n",
    "    irow = rowindex0[i,0];\n",
    "    icol = colindex0[i,0];\n",
    "    if ( F[irow,icol] == 0 ):\n",
    "        rowindex1[N,0]=irow;\n",
    "        colindex1[N,0]=icol;\n",
    "        N=N+1;\n",
    "        F[irow,icol] = 1;              \n",
    "rowindex=eda_cvec(rowindex1[0:N,0]);\n",
    "colindex=eda_cvec(colindex1[0:N,0]);\n",
    "xobs = x[rowindex,0];\n",
    "yobs = y[colindex,0];\n",
    "\n",
    "# create simulated pressure data (solution to Laplace's equation)\n",
    "kappa = 0.1;\n",
    "dtrue = 10*np.multiply( np.sin(kappa*xobs), np.exp(-kappa*yobs) );\n",
    "sigmad = 0.1;\n",
    "dobs = dtrue + np.random.normal(0.0,sigmad,(N,1) );\n",
    "\n",
    "# populate matrix A with observations\n",
    "A = np.zeros((I,J));\n",
    "for i in range(N):\n",
    "    A[rowindex[i,0],colindex[i,0]] = dobs[i,0];\n",
    "\n",
    "# B&W colormap\n",
    "bw = np.zeros((256,4));\n",
    "v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "bw[:,0] = v;\n",
    "bw[:,1] = v;\n",
    "bw[:,2] = v;\n",
    "bw[:,3] = np.ones(256);\n",
    "bwcmap = ListedColormap(bw);\n",
    "    \n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis([ymin, ymax, xmin, xmax]);\n",
    "left=ymin;\n",
    "right=ymax;\n",
    "bottom=xmin;\n",
    "top=xmax;\n",
    "plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );\n",
    "plt.xlabel('y');\n",
    "plt.ylabel('x');\n",
    "plt.title('data');\n",
    "\n",
    "# triangulation\n",
    "XY = np.concatenate( (xobs,yobs), axis=1 );\n",
    "mytri = sp.Delaunay( XY  );\n",
    "\n",
    "# plot triangles\n",
    "TRI = mytri.simplices;\n",
    "[NTRI, i] = np.shape(TRI);\n",
    "\n",
    "# plot\n",
    "ax1 = plt.subplot(1,2,2);\n",
    "plt.axis([ymin, ymax, xmin, xmax]);\n",
    "left=ymin;\n",
    "right=ymax;\n",
    "bottom=xmin;\n",
    "top=xmax;\n",
    "plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );\n",
    "for i in range(NTRI):\n",
    "    v1=TRI[i,0];\n",
    "    v2=TRI[i,1];\n",
    "    v3=TRI[i,2];\n",
    "    plt.plot( [XY[v1,1], XY[v2,1]], [XY[v1,0], XY[v2,0]], 'k-');\n",
    "    plt.plot( [XY[v2,1], XY[v3,1]], [XY[v2,0], XY[v3,0]], 'k-');\n",
    "    plt.plot( [XY[v3,1], XY[v1,1]], [XY[v3,0], XY[v1,0]], 'k-');\n",
    "plt.xlabel('y');\n",
    "plt.ylabel('x');\n",
    "plt.title('triangles');\n",
    "plt.show();\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(16,6));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([ymin, ymax, xmin, xmax]);\n",
    "left=ymin;\n",
    "right=ymax;\n",
    "bottom=xmin;\n",
    "top=xmax;\n",
    "plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );\n",
    "plt.xlabel('y');\n",
    "plt.ylabel('x');\n",
    "plt.title('one triangle enclosing a specified point');\n",
    "\n",
    "# find triangle enclosing exemplary point and plot them\n",
    "x0=[20.5,30.5];\n",
    "k = mytri.find_simplex( x0 );\n",
    "if( k != -1 ):\n",
    "    plt.plot( x0[1], x0[0], 'ko');\n",
    "    v1=TRI[k,0];\n",
    "    v2=TRI[k,1];\n",
    "    v3=TRI[k,2];\n",
    "    plt.plot( [XY[v1,1], XY[v2,1]], [XY[v1,0], XY[v2,0]], 'k-');\n",
    "    plt.plot( [XY[v2,1], XY[v3,1]], [XY[v2,0], XY[v3,0]], 'k-');\n",
    "    plt.plot( [XY[v3,1], XY[v1,1]], [XY[v3,0], XY[v1,0]], 'k-'); \n",
    "plt.show();\n",
    "\n",
    "# 2D interpolation\n",
    "XGRID = np.matmul( x, np.ones((1,J)) );\n",
    "YGRID = np.matmul( np.ones((I,1)), y.T );\n",
    "# note: could also use method='linear' for less-smooth interpolation\n",
    "BB=ip.griddata( np.concatenate((xobs,yobs),axis=1),dobs,(XGRID,YGRID), method='cubic', fill_value=0.0);\n",
    "B = BB[:,:,0];\n",
    "\n",
    "# plot\n",
    "fig3 = plt.figure(3,figsize=(16,6));\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis([ymin, ymax, xmin, xmax]);\n",
    "left=ymin;\n",
    "right=ymax;\n",
    "bottom=xmin;\n",
    "top=xmax;\n",
    "plt.imshow( np.flipud(A), cmap=bwcmap, vmin=np.min(A), vmax=np.max(A), extent=(left,right,bottom,top) );\n",
    "plt.xlabel('y');\n",
    "plt.ylabel('x');\n",
    "plt.title('data');\n",
    "ax1 = plt.subplot(1,2,2);\n",
    "plt.axis([ymin, ymax, xmin, xmax]);\n",
    "left=ymin;\n",
    "right=ymax;\n",
    "bottom=xmin;\n",
    "top=xmax;\n",
    "plt.imshow( np.flipud(B), cmap=bwcmap, vmin=np.min(B), vmax=np.max(B), extent=(left,right,bottom,top) );\n",
    "plt.xlabel('y');\n",
    "plt.ylabel('x');\n",
    "plt.title('interpolated data');\n",
    "plt.show();\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13312 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Symmetry test:\n",
      "    error in transform 0.0000\n",
      "    error in function  0.0000\n"
     ]
    }
   ],
   "source": [
    "# eda10_07\n",
    "# illustrate two-dimensional fourier transform\n",
    "# Part 1, FFT of cosine wave, reordered so zero us in center\n",
    "# Part 2, Illustration of use of symmetry\n",
    "\n",
    "# standard distance/wavenumber set up\n",
    "Nx=32;               # numer of x's, must be even\n",
    "Nxo2 = floor(Nx/2);  # integer version of Nx/2\n",
    "Nkx=Nxo2+1;          # number of non-negative x-wavenumbers kx\n",
    "Dx=1.0;              # sampling interval in x\n",
    "x = eda_cvec(Dx*np.linspace(0,Nx-1,Nx)); # vector of x's\n",
    "xmax = Dx*(Nx-1);    # ending x\n",
    "kxmax=2*pi/(2*Dx);   # nyquist kx\n",
    "Dkx=kxmax/(Nx/2);    # sampling in kx\n",
    "# full (pos and neg) x-wavenumber vector\n",
    "kx=eda_cvec( Dkx * np.concatenate( ( np.linspace(0,Nxo2,Nkx), np.linspace(-Nxo2+1,-1,Nxo2-1) ), axis=0 ) );\n",
    "\n",
    "Ny=32;               # numer of y's, must be even\n",
    "Nyo2 = floor(Ny/2);  # integer version of Ny/2            \n",
    "Nky=Nyo2+1;          # number of non-negative y-wavenumbers ky\n",
    "Dy=1.0;              # sampling interval in y\n",
    "y = eda_cvec(Dy*np.linspace(0,Ny-1,Ny)); # vector of y's\n",
    "ymax = Dy*(Ny-1);    # ending y\n",
    "kymax=2*pi/(2*Dy);   # nyquist ky\n",
    "Dky=kymax/(Ny/2);    # sampling in ky\n",
    "# full (pos and neg) y-wavenumber vector\n",
    "ky=eda_cvec(Dky * np.concatenate( ( np.linspace(0,Nyo2,Nky), np.linspace(-Nyo2+1,-1,Nyo2-1) ), axis=0 ));\n",
    "\n",
    "# Part 1: cosine wave with radial wavenumber kr1 in the theta diection\n",
    "theta = (pi/180.0)*30.0;\n",
    "t1 = eda_cvec( [cos(theta), sin(theta)] );\n",
    "kr1 = kxmax/4.0;\n",
    "F = np.zeros((Nx,Ny));\n",
    "for n in range(Nx):\n",
    "    for m in range(Ny):\n",
    "        F[n,m] = cos( kr1*t1[0,0]*x[n,0] + kr1*t1[1,0]*y[m,0] );\n",
    "\n",
    "# amplitude spectral density\n",
    "Ftt = np.fft.fft2(F);\n",
    "norm = ((Dx**2)/xmax) * ((Dy**2)/ymax);\n",
    "Fasd = norm * np.sqrt( np.power( np.abs(Ftt), 2));\n",
    "\n",
    "# put into more reasonable order, with zero in the center\n",
    "Fasd2=np.zeros((Nx,Ny));\n",
    "Fasd2[Nx-Nkx:Nx,Ny-Nky:Ny] = Fasd[0:Nkx,0:Nky];\n",
    "Fasd2[0:Nkx-2,Ny-Nky:Ny] = Fasd[Nkx:Nx,0:Nky];\n",
    "Fasd2[0:Nkx-2,0:Nky-2] = Fasd[Nkx:Nx,Nky:Ny];\n",
    "Fasd2[Nkx-2:Nx,0:Nky-2] = Fasd[0:Nkx,Nky:Ny];\n",
    "\n",
    "# Part 2: example of the use of symmetry\n",
    "\n",
    "# new function F is noise, which has a complicated Fourier Transform\n",
    "F=np.random.normal(0.0,1.0,(Nx,Ny) );\n",
    "Ftt = np.fft.fft2(F);\n",
    "\n",
    "# for real data, use symmetries to build right hand columns of\n",
    "# 2D Fourier Transform from left hand columns\n",
    "# First copy left Ny/2+1 columns; they are treated as known\n",
    "Ftt2=np.zeros((Nx,Ny),dtype=np.complex);\n",
    "Ftt2[:,0:Nky] = Ftt[:,0:Nky];\n",
    "# then rebuild right from left columns\n",
    "for m in range(1,Nyo2):\n",
    "    mp=Ny-m;\n",
    "    Ftt2[0,mp] = np.conj(Ftt2[0,m]);\n",
    "    Ftt2[1:Nx,mp] = np.flipud(np.conj(Ftt2[1:Nx,m]));\n",
    "\n",
    "# invert back to F\n",
    "F2 = np.real(np.fft.ifft2(Ftt2));\n",
    "\n",
    "eda_draw('Part 1', 'title ASD', Fasd, 'title ASD rearranged', Fasd2, 'Part 2', 'title F2', F, 'title reconstructed F2', F2);\n",
    "\n",
    "# compute error, as a check\n",
    "Ett=np.sum(np.abs(Ftt2-Ftt));\n",
    "E = np.sum(np.abs(F2-F));\n",
    "print('Symmetry test:');\n",
    "print(\"    error in transform %.4f\" % (Ett) );\n",
    "print(\"    error in function  %.4f\" % (E)    );\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 49208 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda10_04_GPR\n",
    "# simple example of Gaussian Process Regression\n",
    "\n",
    "# define the time t axis\n",
    "M=101;      # number of pints\n",
    "Dt=1.0;     # sampling interval\n",
    "tmin=0.0;   # start time\n",
    "t = eda_cvec( tmin+Dt*np.linspace(0,M-1,M) ); # time vector\n",
    "tmax=t[M-1,0]; # ending time\n",
    "\n",
    "# the true model parameters are sinusoidally varying with frequency f0\n",
    "f0=0.02;\n",
    "mtrue = np.sin(2*pi*f0*t);\n",
    "\n",
    "# the data are the values the function measured on only a\n",
    "# subset of randomly chosen points along the time axis\n",
    "N=40;  # number of randomly chosen points\n",
    "rowindex = np.sort( np.random.randint(M,size=(N)));\n",
    "tobs = eda_cvec(t[rowindex,:]);\n",
    "dobs = eda_cvec(mtrue[rowindex,:]);\n",
    "\n",
    "# add noise to data\n",
    "sigmad = 0.2; # standard deviation of noise\n",
    "sigma2d = sigmad**2; # variance of noise\n",
    "nse = np.random.normal(0.0,sigmad,(N,1)); # noise vector\n",
    "dobs=dobs+nse; # observed = signal + noise\n",
    "\n",
    "# estimate the model parameters at these times\n",
    "test = t;\n",
    "\n",
    "# plot the observed data with circles\n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "plt.axis([tmin, tmax, -2, 2]);\n",
    "plt.plot( t, mtrue, 'k-');\n",
    "plt.plot( tobs, dobs, 'ko');\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t) and m(t)');\n",
    "\n",
    "# The auto-covariance function of the data is assumed to be\n",
    "# a Gaussian function with ampltude gamma2 and variance L^2.\n",
    "\n",
    "# set the constants\n",
    "gamma2 = 1.0;\n",
    "L = 5;\n",
    "L2 = L**2;\n",
    "\n",
    "# GPR equation is m = CpMN inv( sigma2d * eye(N,N) + CpNN ) dobs;\n",
    "# but this is implemented in serval steps\n",
    "# Q = sigma2d * eye(N,N) + CpNN\n",
    "# Q v = dobs\n",
    "# m = CpMN v\n",
    "\n",
    "# vectors of ones used to make auto-covariance matrices\n",
    "wN = np.ones((N,1));\n",
    "wM = np.ones((M,1));\n",
    "\n",
    "# set up the \"right\" autocovariance matrix, CpNN\n",
    "CpNN = np.exp( -np.power(np.abs(np.matmul(tobs,wN.T)-np.matmul(wN,tobs.T)), 2) / (2*L2) );\n",
    "\n",
    "# set up the \"left\" autocovariance matrix, CpMN\n",
    "CpMN = np.exp( -np.power(np.abs(np.matmul(wM,tobs.T)-np.matmul(test,wN.T)), 2) / (2*L2) );\n",
    "\n",
    "# GPR formula\n",
    "Q = sigma2d * np.identity(N) + CpNN;\n",
    "mest = np.matmul( CpMN, la.solve(Q,dobs) );\n",
    "    \n",
    "# plot estimated data with line\n",
    "plt.plot(test,mest,'k-',color=(0.7,0.7,0.7));\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
