{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_10_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( N==0 ):\n",
    "            continue; # skip vector of zero length\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",
    "\n",
    "def eda_sign(x):\n",
    "# sign function; apparently Python doesn't have one\n",
    "    if x>0.0:\n",
    "        s=1.0;\n",
    "    elif x<0.0:\n",
    "        s=(-1.0);\n",
    "    else:\n",
    "        s=0.0;\n",
    "    return(s);\n",
    "\n",
    "# function to perform the multiplication y = (Ccc + sigmad2 I) v\n",
    "def eda_CpImul(v):\n",
    "    global eda_N;\n",
    "    global eda_td;\n",
    "    global eda_sigmad2;\n",
    "    global eda_gamma2;\n",
    "    global eda_w0est;\n",
    "    global eda_sest;\n",
    "\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",
    "# Note (C + sigmad2 I) v implies yi = SUMj Cij vj + sigmad2 vi\n",
    "    y = np.zeros((eda_N,1));\n",
    "    for i in range(eda_N):\n",
    "        y[i,0]=eda_sigmad2*vv[i,0];\n",
    "        for j in range(eda_N):\n",
    "            Tabs = abs(eda_td[i,0]-eda_td[j,0]);\n",
    "            Cij = eda_gamma2*cos(eda_w0est*Tabs)*exp(-eda_sest*Tabs);\n",
    "            y[i,0]=y[i,0]+Cij*vv[j,0];\n",
    "    return y;"
   ]
  },
  {
   "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": [
    "# edapy10_01, fitting a high-order polynomial\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": 3,
   "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, 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,\n",
    "       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": 5,
   "metadata": {
    "scrolled": true
   },
   "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": [
    "# edapy_10_03, cubic 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",
    "# cubic interpolation\n",
    "f = ip.interp1d(t.ravel(),d.ravel(),axis=0,\n",
    "        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": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gamma: 1.0000 s: 0.0784\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 30856 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 27512 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_10_04: make correlated random time series\n",
    "\n",
    "# tunable parameters\n",
    "MPLOT = 5;   # scale of model parameter plot\n",
    "CPLOT = 1;   # scale of covariance plot\n",
    "\n",
    "# standard time axis\n",
    "N=256;  # must be even, best if a power of 2\n",
    "Dt=1.0;\n",
    "tmax=Dt*(N-1); \n",
    "t=eda_cvec(np.linspace(0, tmax, N)); \n",
    "tmid=tmax/2.0;\n",
    "\n",
    "# exact gamma*cos(w0T)exp(-sT) covariance function\n",
    "n=6;  # number of oscillations on (0,tmax) interal\n",
    "w0 = n*2*pi/tmax;  # angular frequency of oscilation\n",
    "s = 20.0/tmax;     # decay rate\n",
    "gamma = 1.0;       # sqrt variance\n",
    "gamma2 = gamma**2; #variance\n",
    "c = gamma2*np.multiply( np.cos(w0*(t-tmid)), np.exp(-s*abs(t-tmid)) );\n",
    "\n",
    "print(\"gamma: %.4f s: %.4f\" % (gamma, s) );\n",
    "\n",
    "# realize a time series with a specified covariance\n",
    "ctilde = np.fft.fft(c,axis=0);  # Fourier transform of c\n",
    "casd = np.abs(ctilde);  # amplitude spectral density\n",
    "# make random numbers\n",
    "d = np.random.normal(loc=0.0,scale=1.0,size=(N,1)); \n",
    "d = d - np.mean(d);            # remove mean\n",
    "dtilde = np.fft.fft(d,axis=0); # Fourier transform of d\n",
    "dasd = np.abs(dtilde);         # a.s.d. of d\n",
    "# adjust so has a.s.d. of c\n",
    "dtilde = np.multiply(np.divide(casd,dasd),dtilde);\n",
    "d = np.real(np.fft.ifft(dtilde,axis=0)); # inverse Fourier transform\n",
    "vard = np.var(d);     # variance of d\n",
    "d = (gamma/sqrt(vard))*d; # normalize to variance gamma^2\n",
    "\n",
    "# estimated covariance\n",
    "cest = eda_cvec(np.correlate(d.ravel(),d.ravel(),mode='full') )/N; \n",
    "\n",
    "# plot this covariance\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);\n",
    "plt.xlabel('time lag t (s)');\n",
    "plt.ylabel('c(t)');\n",
    "plt.plot( t-tmid,c,'k-',lw=3);\n",
    "tc = Dt*eda_cvec(np.linspace(-N+1,N-1,2*N-1));\n",
    "plt.plot( tc,cest,'k:','LineWidth',3);\n",
    "\n",
    "# plot time series\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "plt.axis([ 0, tmax, -MPLOT, MPLOT]);\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t)');\n",
    "plt.plot( t,d,'k-',lw=3);\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gamma: 1.0000 w0: 0.1478 s: 0.0784\n",
      "target N: 50, actual 50\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 32952 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 59280 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_10_05, Gaussian Process Regression example\n",
    "# this example uses synthetic data\n",
    "# and a gamma2*cos(w0T)*exp(-sT) covariance function,\n",
    "# with T = |t(i)-t(j)|\n",
    "\n",
    "# tunable parameters\n",
    "MPLOT = 5;    # scale of model parameter plot\n",
    "CPLOT = 1;    # scale of covariance plot\n",
    "gamma = 1;    # sqrt(variance) of covariance\n",
    "sigmad = 0.2; # noise level of observed data\n",
    "wbias=1.0;    # compute GPR with wrong w0\n",
    "sbias=1.0;    # compute GPR with wrong s\n",
    "PLOTCONF=1;   # plot confidence intervals on 1\n",
    "USEFORRLOOP=0; # for loops on 1, outer products on 0\n",
    "\n",
    "# axes definitons\n",
    "\n",
    "# model parameter axis, tm, with sampling Dt\n",
    "M=256;  # must be even, best if a power of 2\n",
    "Dtm=1.0;\n",
    "tmax=Dtm*(M-1); \n",
    "tm=Dtm*eda_cvec(np.linspace(0,tmax,M)); \n",
    "\n",
    "# finely sampled model parameter vector, with sampling\n",
    "# Dt/K, where K is a smallish power of 2, is used in\n",
    "# the initialization process, as described below\n",
    "K=8; # must be even, best if a power of 2\n",
    "Mf=M*K;  \n",
    "Dtf=Dtm/K;\n",
    "tf=Dtf*eda_cvec(np.linspace(0,Mf-1,Mf));\n",
    "\n",
    "# Part 1: realization of a finely sampled time\n",
    "# % series, mf, with the desired covariance\n",
    "\n",
    "# PART 1: create true \"fine\" covariance function\n",
    "\n",
    "n=6;  # number of oscillations on (0,tmax) interal\n",
    "w0 = n*2*pi/tmax;  # angular frequency of oscilation\n",
    "tmid=tmax/2.0;\n",
    "s = 20.0/tmax;     # decay rate\n",
    "gamma2 = gamma**2; # variance\n",
    "cf = gamma2*np.multiply(np.cos(w0*(tf-tmid)),np.exp(-s*np.abs(tf-tmid)));\n",
    "\n",
    "print(\"gamma: %.4f w0: %.4f s: %.4f\" % (gamma, w0, s) );\n",
    "\n",
    "# PART 2: create \"fine\" time series with thus\n",
    "# covariance using Fourier techniques\n",
    "\n",
    "# amplitude spectral density of true covariance\n",
    "cftilde = np.fft.fft(cf,axis=0);\n",
    "cfasd = np.abs(cftilde);\n",
    "\n",
    "# To create a time series with a specified covariannce,\n",
    "# start with a times series of uniformly-distrbuted\n",
    "# random noise, which has a uniform a.s.d., Fourier\n",
    "# transorm it, multiply teh transform by the a.s.d\n",
    "# of the specified covariance, trasnform back, and adjust\n",
    "# the amplitude to the desired variance.\n",
    "# random numbers\n",
    "mf = np.random.normal(loc=0.0,scale=1.0,size=(Mf,1));\n",
    "mf = mf - np.mean(mf); # remove mean\n",
    "mftilde = np.fft.fft(mf,axis=0); # Fourier transform of mf\n",
    "mfasd = np.abs(mftilde);   # a.s.d. of mf\n",
    "# adjust so has a.s.d. of c\n",
    "mftilde = np.multiply(np.divide(cfasd,mfasd),mftilde);\n",
    "# inverse Fourier transform\n",
    "mf = np.real(np.fft.ifft(mftilde,axis=0));\n",
    "varmf = np.var(mf);      # variance of mf\n",
    "# normalize mf to variance gamma^2\n",
    "mf = (gamma/sqrt(varmf))*mf;\n",
    "\n",
    "# PART 3: Regularly sub-sample mf to get mtrue, irregularly\n",
    "# sub-sample mf to get true data dtrue, and add uncorrelated\n",
    "# noise to dtrue to get observed data dobs\n",
    "\n",
    "mtrue = np.copy(mf[0:Mf:K]);\n",
    "\n",
    "Ntarget = 50; # fewer observed data than true data\n",
    "# create list, j, of N unique, randomly chosen indices\n",
    "# that are not in the sequence 1, K+1, 2K+1, ...\n",
    "j=np.zeros((Ntarget,1),dtype=int);\n",
    "Nj=0;\n",
    "while (Nj<Ntarget):\n",
    "    i=np.random.randint(low=0,high=Mf-1,size=(1));\n",
    "    if( (((i-1)%K) != 0) and (not np.isin(i,j)) ):\n",
    "        j[Nj,0]=i;\n",
    "        Nj=Nj+1;\n",
    "# sort for plotting purposes\n",
    "j = eda_cvec(np.sort(j.ravel())); \n",
    "j = eda_cvec(np.unique(j.ravel())); \n",
    "N,i = np.shape(j);\n",
    "print('target N: %d, actual %d' % (Ntarget,N) );\n",
    "# times of observations\n",
    "td = np.copy(tf[j.ravel()]); \n",
    "# irregularly sub-sampled data\n",
    "dtrue = np.copy(mf[j.ravel()]);  \n",
    "# observed data is true data + uncorrelated random noise\n",
    "sigmad2 = sigmad**2;\n",
    "dobs = dtrue+np.random.normal(loc=0,scale=sigmad,size=(N,1));\n",
    "\n",
    "# the data kernel, G, for the equation d=Gm\n",
    "G = np.concatenate( (np.zeros((N,M)), np.zeros((N,N))), axis=1 );\n",
    "\n",
    "# Since the true (w0, s) are in general unknown\n",
    "# this code allows them to be assigned differently\n",
    "# from the true values\n",
    "w0est = wbias*w0;  # w0est = 1.0*w0 is the true value\n",
    "sest = sbias*s;\n",
    "\n",
    "if( USEFORRLOOP ):\n",
    "    # build covariances via for loop\n",
    "    Ccc = np.zeros((N,N)); # covariance matrix, C(control,control)\n",
    "    for i in range(N):\n",
    "        for j in range(N):\n",
    "            T = abs(td(i)-td(j));\n",
    "            Ccc[i,j] = gamma2*cos(w0est*T)*exp(-sest*T);\n",
    "    \n",
    "    Ctc = np.zeros((M,N)); #  covariance matrix, C(target,control)\n",
    "    for i in rane(M):\n",
    "        for j in range(N):\n",
    "            T = abs(tm(i)-td(j));\n",
    "            Ctc[i,j] = gamma2*cos(w0est*T)*exp(-sest*T);   \n",
    "else: # use outer product\n",
    "    # build covariances via outer product\n",
    "    Tdd = np.abs(np.matmul(td,np.ones((1,N)))-np.matmul(np.ones((N,1)),(td.T)));\n",
    "    Ccc = gamma2*np.multiply(np.cos(w0est*Tdd),np.exp(-sest*Tdd));\n",
    "    Tmd = np.abs(np.matmul(tm,np.ones((1,N)))-np.matmul(np.ones((M,1)),(td.T)));\n",
    "    Ctc = gamma2*np.multiply(np.cos(w0est*Tmd),np.exp(-sest*Tmd));\n",
    "\n",
    "# solve GPE problem\n",
    "A = (Ccc + sigmad2*np.identity(N) );\n",
    "mest = np.matmul(Ctc,la.solve(A,dobs));\n",
    "\n",
    "# estimated covariance\n",
    "cest = eda_cvec(\n",
    "    np.correlate(mest.ravel(),mest.ravel(),mode='full'))/M;\n",
    "\n",
    "# plot this covariance\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);\n",
    "plt.xlabel('time lag t (s)');\n",
    "plt.ylabel('c(t)');\n",
    "plt.plot( tf-tmid,cf,'k-',lw=3);\n",
    "tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));\n",
    "plt.plot(tc,cest,'k:',lw=3);\n",
    "plt.show();\n",
    "\n",
    "# time series figure\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "plt.axis( [ 0, tmax, -MPLOT, MPLOT] );\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t) and  m(t)');\n",
    "plt.plot(tm,mtrue,'k-',lw=3);\n",
    "plt.plot(td,dtrue,'ko',lw=2);\n",
    "plt.plot(tm,mest,'-',lw=2,color=[0.8,0.8,0.8]);\n",
    "if( PLOTCONF ): # plot 95% confidence intervals\n",
    "    # MATLAB: B = Ctc/A solves BA = Ctc\n",
    "    # take transpose, so AT BT = CctT, then solve for BT\n",
    "    # BT = solve(AT,CtcT) = solve(A,CtcT), since A=AT\n",
    "    # MATLAB: B BT becomes BTT * BT\n",
    "    BT = la.solve(A,Ctc.T)\n",
    "    Cttest = sigmad2*np.matmul(BT.T,BT);\n",
    "    sigmamest = np.sqrt(eda_cvec(np.diag(Cttest)));\n",
    "    plt.plot(tm,mest-2*sigmamest,'k-',lw=1);\n",
    "    plt.plot(tm,mest+2*sigmamest,'k-',lw=1);\n",
    "plt.show();\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gamma: 1.0000 w0: 0.1478 s: 0.0784\n",
      "target N: 50, actual 50\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 29492 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 39940 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_10_06, GPE using bicg()\n",
    "# this example uses synthetic data\n",
    "# and a gamma2*cos(w0T)*exp(-sT) covariance function,\n",
    "# with T = |t(i)-t(j)|\n",
    "\n",
    "# tunable parameters\n",
    "MPLOT = 5;    # scale of model parameter plot\n",
    "CPLOT = 1;    # scale of covariance plot\n",
    "gamma = 1;    # sqrt(variance) of covariance\n",
    "sigmad = 0.2; # noise level of observed data\n",
    "wbias=1.0;    # compute GPR with wrong w0\n",
    "sbias=1.0;    # compute GPR with wrong s\n",
    "PLOTCONF=1;   # plot confidence intervals on 1\n",
    "USEFORRLOOP=0; # for loops on 1, outer products on 0\n",
    "\n",
    "# axes definitons\n",
    "\n",
    "# model parameter axis, tm, with sampling Dt\n",
    "M=256;  # must be even, best if a power of 2\n",
    "Dtm=1.0;\n",
    "tmax=Dtm*(M-1); \n",
    "tm=Dtm*eda_cvec(np.linspace(0,tmax,M)); \n",
    "\n",
    "# finely sampled model parameter vector, with sampling\n",
    "# Dt/K, where K is a smallish power of 2, is used in\n",
    "# the initialization process, as described below\n",
    "K=8; # must be even, best if a power of 2\n",
    "Mf=M*K;  \n",
    "Dtf=Dtm/K;\n",
    "tf=Dtf*eda_cvec(np.linspace(0,Mf-1,Mf));\n",
    "\n",
    "# Part 1: realization of a finely sampled time\n",
    "# % series, mf, with the desired covariance\n",
    "\n",
    "# PART 1: create true \"fine\" covariance function\n",
    "\n",
    "n=6;  # number of oscillations on (0,tmax) interal\n",
    "w0 = n*2*pi/tmax;  # angular frequency of oscilation\n",
    "tmid=tmax/2.0;\n",
    "s = 20.0/tmax;     # decay rate\n",
    "gamma2 = gamma**2; # variance\n",
    "cf = gamma2*np.multiply(np.cos(w0*(tf-tmid)),np.exp(-s*np.abs(tf-tmid)));\n",
    "\n",
    "print(\"gamma: %.4f w0: %.4f s: %.4f\" % (gamma, w0, s) );\n",
    "\n",
    "# PART 2: create \"fine\" time series with thus\n",
    "# covariance using Fourier techniques\n",
    "\n",
    "# amplitude spectral density of true covariance\n",
    "cftilde = np.fft.fft(cf,axis=0);\n",
    "cfasd = np.abs(cftilde);\n",
    "\n",
    "# To create a time series with a specified covariannce,\n",
    "# start with a times series of uniformly-distrbuted\n",
    "# random noise, which has a uniform a.s.d., Fourier\n",
    "# transorm it, multiply teh transform by the a.s.d\n",
    "# of the specified covariance, trasnform back, and adjust\n",
    "# the amplitude to the desired variance.\n",
    "# random numbers\n",
    "mf = np.random.normal(loc=0.0,scale=1.0,size=(Mf,1));\n",
    "mf = mf - np.mean(mf); # remove mean\n",
    "mftilde = np.fft.fft(mf,axis=0); # Fourier transform of mf\n",
    "mfasd = np.abs(mftilde);   # a.s.d. of mf\n",
    "# adjust so has a.s.d. of c\n",
    "mftilde = np.multiply(np.divide(cfasd,mfasd),mftilde);\n",
    "# inverse Fourier transform\n",
    "mf = np.real(np.fft.ifft(mftilde,axis=0));\n",
    "varmf = np.var(mf);      # variance of mf\n",
    "# normalize mf to variance gamma^2\n",
    "mf = (gamma/sqrt(varmf))*mf;\n",
    "\n",
    "# PART 3: Regularly sub-sample mf to get mtrue, irregularly\n",
    "# sub-sample mf to get true data dtrue, and add uncorrelated\n",
    "# noise to dtrue to get observed data dobs\n",
    "\n",
    "mtrue = np.copy(mf[0:Mf:K]);\n",
    "\n",
    "Ntarget = 50; # fewer observed data than true data\n",
    "# create list, j, of N unique, randomly chosen indices\n",
    "# that are not in the sequence 1, K+1, 2K+1, ...\n",
    "j=np.zeros((Ntarget,1),dtype=int);\n",
    "Nj=0;\n",
    "while (Nj<Ntarget):\n",
    "    i=np.random.randint(low=0,high=Mf-1,size=(1));\n",
    "    if( (((i-1)%K) != 0) and (not np.isin(i,j)) ):\n",
    "        j[Nj,0]=i;\n",
    "        Nj=Nj+1;\n",
    "# sort for plotting purposes\n",
    "j = eda_cvec(np.sort(j.ravel())); \n",
    "j = eda_cvec(np.unique(j.ravel()));\n",
    "N, i =np.shape(j);\n",
    "print('target N: %d, actual %d' % (Ntarget, N) );\n",
    "# times of observations\n",
    "td = np.copy(tf[j.ravel()]); \n",
    "# irregularly sub-sampled data\n",
    "dtrue = np.copy(mf[j.ravel()]);  \n",
    "# observed data is true data + uncorrelated random noise\n",
    "sigmad2 = sigmad**2;\n",
    "dobs = dtrue+np.random.normal(loc=0,scale=sigmad,size=(N,1));\n",
    "\n",
    "# Since the true (w0, s) are in general unknown\n",
    "# this code allows them to be assigned differently\n",
    "# from the true values\n",
    "w0est = wbias*w0;  # w0est = 1.0*w0 is the true value\n",
    "sest = sbias*s;\n",
    "\n",
    "# set globals\n",
    "eda_N=N;\n",
    "eda_td=np.copy(td);\n",
    "eda_sigmad2=sigmad2;\n",
    "eda_gamma2=gamma2;\n",
    "eda_w0est=w0est;\n",
    "eda_sest=sest;\n",
    "\n",
    "# define linear operator needed for conjugate gradienet solver\n",
    "LO=las.LinearOperator(shape=(N,N),\n",
    "                      matvec=eda_CpImul,rmatvec=eda_CpImul);\n",
    "\n",
    "# solve using congugate gradient algorithm\n",
    "u = np.zeros((N,1));\n",
    "q=las.bicg(LO,dobs,tol=1e-12,maxiter=(3*(N+M)+100));\n",
    "u[:,0] = q[0];\n",
    "mest = np.zeros((M,1));\n",
    "for i in range(M):\n",
    "    mest[i,0]=0.0;\n",
    "    for j in range(N):\n",
    "        Tabs = abs(tm[i,0]-td[j,0]);\n",
    "        Cij = gamma2*cos(w0est*Tabs)*exp(-sest*Tabs);\n",
    "        mest[i,0]=mest[i,0]+Cij*u[j,0];\n",
    "\n",
    "dpre = np.zeros((N,1));\n",
    "for i in range(N):\n",
    "    dpre[i,0]=0.0;\n",
    "    for j in range(N):\n",
    "        Tabs = abs(td[i,0]-td[j,0]);\n",
    "        Cij = gamma2*cos(w0est*Tabs)*exp(-sest*Tabs);\n",
    "        dpre[i,0]=dpre[i,0]+Cij*u[j,0];\n",
    "\n",
    "# estimated covariance\n",
    "cest = eda_cvec(\n",
    "    np.correlate(mest.ravel(),mest.ravel(),mode='full'))/M;\n",
    "\n",
    "# plot this covariance\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);\n",
    "plt.xlabel('time lag t (s)');\n",
    "plt.ylabel('c(t)');\n",
    "plt.plot( tf-tmid,cf,'k-',lw=3);\n",
    "tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));\n",
    "plt.plot(tc,cest,'k:',lw=3);\n",
    "plt.show();\n",
    "\n",
    "# time series figure\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "plt.axis( [ 0, tmax, -MPLOT, MPLOT] );\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t) and  m(t)');\n",
    "plt.plot(tm,mtrue,'k-',lw=3);\n",
    "plt.plot(td,dtrue,'ko',lw=2);\n",
    "plt.plot(tm,mest,'-',lw=2,color=[0.8,0.8,0.8]);\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gamma: 1.0000 w0: 0.0738 s: 0.0167\n",
      "target N: 50, actual 50\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 35108 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 42704 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1: qg 0.2213 nu -1.00 alpha 1.000000e+00 PSI -1.811824e+03 PSIlast -1.811824e+03\n",
      "2: qg 0.1963 nu -1.00 alpha 1.000000e+00 PSI -1.825454e+03 PSilast -1.825454e+03\n",
      "2: qg 0.1713 nu -1.00 alpha 1.000000e+00 PSI -1.845824e+03 PSilast -1.845824e+03\n",
      "2: qg 0.1463 nu -1.00 alpha 1.000000e+00 PSI -1.872873e+03 PSilast -1.872873e+03\n",
      "2: qg 0.1213 nu -1.00 alpha 1.000000e+00 PSI -1.906843e+03 PSilast -1.906843e+03\n",
      "2: qg 0.0963 nu -1.00 alpha 1.000000e+00 PSI -1.943066e+03 PSilast -1.943066e+03\n",
      "2: qg 0.0713 nu -1.00 alpha 1.000000e+00 PSI -1.963933e+03 PSilast -1.963933e+03\n",
      "2: qg 0.0713 nu 1.00 alpha 5.000000e-01 PSI -1.963933e+03 PSIlast -1.963933e+03\n",
      "2: qg 0.0713 nu -1.00 alpha 2.500000e-01 PSI -1.963933e+03 PSIlast -1.963933e+03\n",
      "2: qg 0.0651 nu -1.00 alpha 2.500000e-01 PSI -1.965123e+03 PSilast -1.965123e+03\n",
      "2: qg 0.0651 nu 1.00 alpha 1.250000e-01 PSI -1.965123e+03 PSIlast -1.965123e+03\n",
      "2: qg 0.0651 nu -1.00 alpha 6.250000e-02 PSI -1.965123e+03 PSIlast -1.965123e+03\n",
      "2: qg 0.0635 nu -1.00 alpha 6.250000e-02 PSI -1.965189e+03 PSilast -1.965189e+03\n",
      "2: qg 0.0635 nu 1.00 alpha 3.125000e-02 PSI -1.965189e+03 PSIlast -1.965189e+03\n",
      "2: qg 0.0635 nu -1.00 alpha 1.562500e-02 PSI -1.965189e+03 PSIlast -1.965189e+03\n",
      "2: qg 0.0631 nu -1.00 alpha 1.562500e-02 PSI -1.965192e+03 PSilast -1.965192e+03\n",
      "2: qg 0.0631 nu 1.00 alpha 7.812500e-03 PSI -1.965192e+03 PSIlast -1.965192e+03\n",
      "2: qg 0.0631 nu -1.00 alpha 3.906250e-03 PSI -1.965192e+03 PSIlast -1.965192e+03\n",
      "2: qg 0.0631 nu 1.00 alpha 1.953125e-03 PSI -1.965192e+03 PSIlast -1.965192e+03\n",
      "2: qg 0.0631 nu -1.00 alpha 9.765625e-04 PSI -1.965192e+03 PSIlast -1.965192e+03\n",
      "2: qg 0.0631 nu -1.00 alpha 9.765625e-04 PSI -1.965192e+03 PSilast -1.965192e+03\n",
      "2: qg 0.0631 nu 1.00 alpha 4.882812e-04 PSI -1.965192e+03 PSIlast -1.965192e+03\n",
      "2: qg 0.0631 nu -1.00 alpha 2.441406e-04 PSI -1.965192e+03 PSIlast -1.965192e+03\n",
      "2: qg 0.0631 nu -1.00 alpha 2.441406e-04 PSI -1.965192e+03 PSilast -1.965192e+03\n",
      "2: qg 0.0631 nu 1.00 alpha 1.220703e-04 PSI -1.965192e+03 PSIlast -1.965192e+03\n",
      "2: qg 0.0631 nu -1.00 alpha 6.103516e-05 PSI -1.965192e+03 PSIlast -1.965192e+03\n",
      "2: qg 0.0631 nu -1.00 alpha 6.103516e-05 PSI -1.965192e+03 PSilast -1.965192e+03\n",
      "2: qg 0.0631 nu 1.00 alpha 3.051758e-05 PSI -1.965192e+03 PSIlast -1.965192e+03\n",
      "2: qg 0.0631 nu -1.00 alpha 1.525879e-05 PSI -1.965192e+03 PSIlast -1.965192e+03\n",
      "2: qg 0.0631 nu 1.00 alpha 1.525879e-05 PSI -1.965192e+03 PSilast -1.965192e+03\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 45208 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 38980 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "qtrue 7.377517e-02 qest 6.308996e-02 relating error 1.448348e-01\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 14292 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edama_10_07, GPR Tuning Example\n",
    "\n",
    "# this example uses synthetic data\n",
    "# and a gamma2*cos(w0T)*exp(-sT) covariance function,\n",
    "# with T = |t(i)-t(j)|\n",
    "\n",
    "# tunable parameters\n",
    "MPLOT = 5;    # scale of model parameter plot\n",
    "CPLOT = 1;    # scale of covariance plot\n",
    "gamma = 1;    # sqrt(variance) of covariance\n",
    "sigmad = 0.2; # noise level of observed data\n",
    "wbias=3.0;    # compute GPR with wrong w0\n",
    "sbias=1.0;    # compute GPR with wrong s\n",
    "\n",
    "# axes definitons\n",
    "\n",
    "# model parameter axis, tm, with sampling Dt\n",
    "M=512;  # must be even, best if a power of 2\n",
    "Dtm=1.0;\n",
    "tmax=Dtm*(M-1); \n",
    "tm=eda_cvec(np.linspace(0,tmax,M)); \n",
    "\n",
    "# finely sampled model parameter vector, with sampling\n",
    "# Dt/K, where K is a smallish power of 2, is used in\n",
    "# the initialization process, as described below\n",
    "K=8; # must be even, best if a power of 2\n",
    "Mf=M*K;  \n",
    "Dtf=Dtm/K;\n",
    "tf=eda_cvec(np.linspace(0,tmax,Mf));\n",
    "\n",
    "# Part 1: realization of a finely sampled time\n",
    "# % series, mf, with the desired covariance\n",
    "\n",
    "# PART 1: create true \"fine\" covariance function\n",
    "\n",
    "n=6;  # number of oscillations on (0,tmax) interal\n",
    "w0 = n*2*pi/tmax;  # angular frequency of oscilation\n",
    "tmid=tmax/2.0;     # midpoint\n",
    "s = Dtm/60.0;      # decay rate\n",
    "gamma2 = gamma**2; # variance\n",
    "cf = gamma2*np.multiply(np.cos(w0*(tf-tmid)),np.exp(-s*np.abs(tf-tmid)));\n",
    "\n",
    "print(\"gamma: %.4f w0: %.4f s: %.4f\" % (gamma, w0, s) );\n",
    "\n",
    "# PART 2: create \"fine\" time series with thus\n",
    "# covariance using Fourier techniques\n",
    "\n",
    "# amplitude spectral density of true covariance\n",
    "cftilde = np.fft.fft(cf,axis=0);\n",
    "cfasd = np.abs(cftilde);\n",
    "\n",
    "# To create a time series with a specified covariannce,\n",
    "# start with a times series of uniformly-distrbuted\n",
    "# random noise, which has a uniform a.s.d., Fourier\n",
    "# transorm it, multiply teh transform by the a.s.d\n",
    "# of the specified covariance, trasnform back, and adjust\n",
    "# the amplitude to the desired variance.\n",
    "# random numbers\n",
    "mf = np.random.normal(loc=0.0,scale=1.0,size=(Mf,1));\n",
    "mf = mf - np.mean(mf); # remove mean\n",
    "mftilde = np.fft.fft(mf,axis=0); # Fourier transform of mf\n",
    "mfasd = np.abs(mftilde);   # a.s.d. of mf\n",
    "# adjust so has a.s.d. of c\n",
    "mftilde = np.multiply(np.divide(cfasd,mfasd),mftilde);\n",
    "# inverse Fourier transform\n",
    "mf = np.real(np.fft.ifft(mftilde,axis=0));\n",
    "varmf = np.var(mf);      # variance of mf\n",
    "# normalize mf to variance gamma^2\n",
    "mf = (gamma/sqrt(varmf))*mf;\n",
    "\n",
    "# PART 3: Regularly sub-sample mf to get mtrue, irregularly\n",
    "# sub-sample mf to get true data dtrue, and add uncorrelated\n",
    "# noise to dtrue to get observed data dobs\n",
    "\n",
    "mtrue = np.copy(mf[0:Mf:K]);\n",
    "\n",
    "Ntarget = 50; # fewer observed data than true data\n",
    "# create list, j, of N unique, randomly chosen indices\n",
    "# that are not in the sequence 1, K+1, 2K+1, ...\n",
    "j=np.zeros((Ntarget,1),dtype=int);\n",
    "Nj=0;\n",
    "while (Nj<Ntarget):\n",
    "    i=np.random.randint(low=0,high=Mf-1,size=(1));\n",
    "    if( (((i-1)%K) != 0) and (not np.isin(i,j)) ):\n",
    "        j[Nj,0]=i;\n",
    "        Nj=Nj+1;\n",
    "# sort for plotting purposes\n",
    "j = eda_cvec(np.sort(j.ravel())); \n",
    "j = eda_cvec(np.unique(j));\n",
    "N, i =np.shape(j);\n",
    "print('target N: %d, actual %d' % (Ntarget, N) );\n",
    "# times of observations\n",
    "td = np.copy(tf[j.ravel()]); \n",
    "# irregularly sub-sampled data\n",
    "dtrue = np.copy(mf[j.ravel()]);  \n",
    "# observed data is true data + uncorrelated random noise\n",
    "sigmad2 = sigmad**2;\n",
    "sigmadr = 1.0/sigmad;\n",
    "dobs = dtrue+np.random.normal(loc=0,scale=sigmad,size=(N,1));\n",
    "\n",
    "# the data kernel, G, for the equation d=Gm\n",
    "G = np.concatenate( (np.zeros((N,M)), np.zeros((N,N))), axis=1 );\n",
    "\n",
    "# PART 4: Solve GPR problem for initial guess of parameter\n",
    "w0est = wbias*w0;  # w0est = 1.0*w0 is the true value\n",
    "sest = sbias*s;\n",
    " \n",
    "# build covariances via outer product\n",
    "Tdd = np.abs(np.matmul(td,np.ones((1,N)))-np.matmul(np.ones((N,1)),(td.T)));\n",
    "Ccc = gamma2*np.multiply(np.cos(w0est*Tdd),np.exp(-sest*Tdd));\n",
    "Tmd = np.abs(np.matmul(tm,np.ones((1,N)))-np.matmul(np.ones((M,1)),(td.T)));\n",
    "Ctc = gamma2*np.multiply(np.cos(w0est*Tmd),np.exp(-sest*Tmd));\n",
    "\n",
    "# solve GPE problem\n",
    "A = (Ccc + sigmad2*np.identity(N) );\n",
    "mest = np.matmul(Ctc,la.solve(A,dobs));\n",
    "\n",
    "# estimated covariance\n",
    "cest = eda_cvec(np.correlate(mest.ravel(),mest.ravel(),mode='full'))/M;\n",
    "\n",
    "# plot startng covariance\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);\n",
    "plt.xlabel('time lag t (s)');\n",
    "plt.ylabel('c(t)');\n",
    "plt.plot( tf-tmid,cf,'k-',lw=3);\n",
    "tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));\n",
    "plt.plot(tc,cest,'k:',lw=3);\n",
    "plt.show();\n",
    "\n",
    "# time series figure\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "plt.axis( [ 0, tmax, -MPLOT, MPLOT] );\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t) and  m(t)');\n",
    "plt.plot(tm,mtrue,'k-',lw=3);\n",
    "plt.plot(td,dtrue,'ko',lw=2);\n",
    "plt.plot(tm,mest,'-',lw=2,color=[0.8,0.8,0.8]);\n",
    "plt.show();\n",
    "\n",
    "# Part 5: Set-up for Tuning\n",
    "\n",
    "# Gradient Descent Method adjustment of w0, now called q\n",
    "qtrue=w0;     # true parameter, q\n",
    "qg = w0est;   # inital guess of parameter, q\n",
    "Dq = 0.025;   # initial step size for gradiet descent\n",
    "NITER = 30;   # maximum number of iterations\n",
    "qg = qg-Dq;   # so first iteration starts qg\n",
    "alpha = 1.0;  # step size reducer\n",
    "nu = 1.0;     # sign of gradient\n",
    "mtpri = np.zeros((M,1)); # prior model parameters are zero\n",
    "mcpri = np.zeros((N,1)); # prior model parameters are zero\n",
    "# concatenate into overall prior vector mpri and times tall\n",
    "mpri = np.concatenate((mtpri,mcpri),axis=0); # prior model pars\n",
    "tall = np.concatenate((tm,td),axis=0); # corresponding times\n",
    "qsave = np.zeros((NITER,1)); # save estimates, for plot\n",
    "\n",
    "# Part 6: Tuning\n",
    "\n",
    "for iter in range(NITER): # loop over descent increments\n",
    "    \n",
    "    qg = qg+nu*alpha*Dq; # new guess of parameter, q\n",
    "\n",
    "    # full covariance matrix and its derivative\n",
    "    Tall = np.matmul(tall,np.ones((1,M+N))) - np.matmul(np.ones((M+N,1)),tall.T);\n",
    "    Call = gamma2*np.multiply(np.cos(qg*Tall),np.exp(-sest*np.abs(Tall)));\n",
    "    dCalldq = -gamma2*np.multiply(\n",
    "              np.multiply(Tall,np.sin(qg*Tall)),np.exp(-sest*abs(Tall)));\n",
    "    \n",
    "    # break out four quandrants; note these are views, not arrays,\n",
    "    # so we need to be careful not to overwrite Ctt\n",
    "    Ctt = Call[0:M,0:M];\n",
    "    dCttdq = dCalldq[0:M,0:M];\n",
    "    Ctc = Call[0:M,M:M+N];\n",
    "    dCtcdq = dCalldq[0:M,M:M+N];\n",
    "    Cct = Ctc.T;\n",
    "    Ccc = Call[M:M+N,M:M+N];\n",
    "    dCccdq = dCalldq[M:M+N,M:M+N];\n",
    "    \n",
    "    # A = Ccc + sigmad^2 I\n",
    "    A = Ccc + sigmad2*np.identity(N);\n",
    "    dAdq = dCccdq;\n",
    "    \n",
    "    # various functions of the covariance matrix\n",
    "    CallI = la.inv(Call);\n",
    "    CallIsr = la.sqrtm(CallI);\n",
    "    dCallIdq = -np.matmul(np.matmul(CallI,dCalldq), CallI);\n",
    "    dCallIsrdq = la.solve_sylvester( CallIsr, CallIsr, dCallIdq );\n",
    "    R = la.cholesky(Call);\n",
    "    logdetCall = 2.0*np.sum(np.log(np.diag(R)));\n",
    "    dlogdetCalldq = np.trace(np.matmul(CallI,dCalldq));\n",
    "\n",
    "    # GPR solution\n",
    "    u=la.solve(A,dobs);\n",
    "    mtg = np.matmul(Ctc,u);\n",
    "    mcg = np.matmul(Ccc,u);\n",
    "    mg = np.concatenate( (mtg,mcg), axis=0 );\n",
    "    dg = eda_cvec(mg[M:M+N]);\n",
    "\n",
    "    # derivative of predicted data\n",
    "    v = la.solve( A, np.matmul(dAdq,u) );\n",
    "    dmestdq = np.concatenate((\n",
    "        np.matmul(dCtcdq,u) - np.matmul(Ctc,v),\n",
    "        np.matmul(dCccdq,u) - np.matmul(Ccc,v)),\n",
    "        axis=0 );\n",
    "    ddpredq = np.copy(dmestdq[M:N+M]);\n",
    "    detildedq = -sigmadr*ddpredq;\n",
    "    e = dobs-dg;\n",
    "    etilde = sigmadr*e;\n",
    "    \n",
    "    # derivative of total data error\n",
    "    E = np.matmul(etilde.T,etilde);\n",
    "    dEdp = 2.0*np.matmul(etilde.T,detildedq);\n",
    "    \n",
    "    # prior error and its derivative\n",
    "    el=(mpri-mg);\n",
    "    eltilde = np.matmul(CallIsr,el);\n",
    "    deltildedq = np.matmul(dCallIsrdq,el)-np.matmul(CallIsr,dmestdq);\n",
    "    \n",
    "    # total prior error and its derivative\n",
    "    L = np.matmul(eltilde.T,eltilde);\n",
    "    dLdp = 2.0*np.matmul(eltilde.T,deltildedq);\n",
    "    \n",
    "    # function, PSI being minimized\n",
    "    PSI = logdetCall + E + L;\n",
    "    # and its deriative\n",
    "    dPSIdp = dlogdetCalldq + dEdp + dLdp;\n",
    "    nu = -eda_sign(dPSIdp);  # downhill vector\n",
    "   \n",
    "    # gradient method method\n",
    "    if( iter==0 ):  # on the first iteration do nothing but save\n",
    "        qglast = qg;\n",
    "        PSIlast = PSI;\n",
    "        print(\"1: qg %.4f nu %.2f alpha %e PSI %e PSIlast %e\" % (qg,nu,alpha,PSI,PSIlast));\n",
    "    elif (PSI < PSIlast ):\n",
    "        qglast = qg;\n",
    "        PSIlast = PSI;\n",
    "        print(\"2: qg %.4f nu %.2f alpha %e PSI %e PSilast %e\" % (qg,nu,alpha,PSI,PSIlast));\n",
    "    else: # psi increased, so decrease alpha\n",
    "        qg = qglast;\n",
    "        PSI = PSIlast;\n",
    "        alpha = alpha/2.0;\n",
    "        print(\"2: qg %.4f nu %.2f alpha %e PSI %e PSIlast %e\" % (qg,nu,alpha,PSI,PSIlast));\n",
    "    qsave[iter,0] = qg;\n",
    "\n",
    "pest = qg;\n",
    "mest = mg;\n",
    "mtest = np.copy(mest[0:M,0]);\n",
    "mcest = np.copy(mest[M:M+N,0]);\n",
    "dpre = mcest;\n",
    "\n",
    "# estimated covariance\n",
    "cest = eda_cvec(\n",
    "    np.correlate(mtest.ravel(),mtest.ravel(),mode='full'))/M;\n",
    "\n",
    "# plot final estimated covariance\n",
    "fig3 = plt.figure(3,figsize=(10,6));\n",
    "plt.axis([ -tmid, tmid, -CPLOT, CPLOT]);\n",
    "plt.xlabel('time lag t (s)');\n",
    "plt.ylabel('c(t)');\n",
    "plt.plot( tf-tmid,cf,'k-',lw=3);\n",
    "tc = eda_cvec(Dtm*np.linspace(-M+1,M-1,2*M-1));\n",
    "plt.plot(tc,cest,'k:',lw=3);\n",
    "plt.show();\n",
    "\n",
    "# plot final time series \n",
    "fig4 = plt.figure(4,figsize=(10,6));\n",
    "plt.axis( [ 0, tmax, -MPLOT, MPLOT] );\n",
    "plt.xlabel('time t (s)');\n",
    "plt.ylabel('d(t) and  m(t)');\n",
    "plt.plot(tm,mtrue,'k-',lw=3);\n",
    "plt.plot(td,dtrue,'ko',lw=2);\n",
    "plt.plot(tm,mtest,'-',lw=2,color=[0.8,0.8,0.8]);\n",
    "plt.show();\n",
    "\n",
    "print(\"qtrue %e qest %e relating error %e\" % (qtrue, pest, (qtrue-pest)/qtrue) );\n",
    "\n",
    "# plot p as a function of iteration number\n",
    "fig5 = plt.figure(5,figsize=(10,6));\n",
    "plt.axis( [ 0, NITER, 0, 0.3] );\n",
    "plt.xlabel('iteration');\n",
    "plt.ylabel('w0');\n",
    "itvec = eda_cvec(np.linspace(1,NITER,NITER));\n",
    "plt.plot(itvec,qsave,'k-',lw=2);\n",
    "plt.plot(itvec,qsave,'ko',lw=2);\n",
    "plt.plot(itvec,w0*np.ones((NITER,1)),'k:',lw=2);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 99484 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 13244 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 16372 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda10_08, 2D interpolation of simulated pressure data\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",
    "DOCUBIC = 0;  # linear interpolation on 0, cubic on 1\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",
    "# Delaunay triangulation\n",
    "XY = np.concatenate( (xobs,yobs), axis=1 );\n",
    "# XY is a N by 2 array of (x,y) positions of vertices\n",
    "mytri = sp.Delaunay( XY  );\n",
    "# TRI is a NTRI by 3 array of the indices of\n",
    "# the vertices associated with each triangle\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  point x0\n",
    "x0=[20.5,30.5];\n",
    "tri0 = mytri.find_simplex( x0 );\n",
    "\n",
    "k = tri0;\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",
    "# grid of (x,y) values\n",
    "XGRID = np.matmul( x, np.ones((1,J)) );\n",
    "YGRID = np.matmul( np.ones((I,1)), y.T );\n",
    "\n",
    "# interpolation\n",
    "if( DOCUBIC):\n",
    "    # linear interpolation\n",
    "    BB=ip.griddata( np.concatenate((xobs,yobs),axis=1),\n",
    "                   dobs,(XGRID,YGRID), method='cubic', fill_value=0.0);\n",
    "    B = BB[:,:,0];\n",
    "else:\n",
    "    # cubic intepolation\n",
    "    CC=ip.griddata( np.concatenate((xobs,yobs),axis=1),\n",
    "                   dobs,(XGRID,YGRID), method='linear', fill_value=0.0);\n",
    "    C = BB[:,:,0];\n",
    "    B = C;\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_10_09, is not provided, since the MATHLAB\n",
    "# version edama_10_09 addresses an issue that does\n",
    "# not arise in Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 7756 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 13720 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_10, two-dimensional Fourier transform\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",
    "eda_draw('title f(x,y)', F);\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": 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
}
