{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eda07_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, sqrt, floor, ceil\n",
    "import numpy as np\n",
    "import scipy.sparse.linalg as las\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",
    "# 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",
    "\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",
    "# 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",
    "\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 19688 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24620 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda07_01\n",
    "# example of a filter problem\n",
    "# using impulse repsonse of a heat-producing layer\n",
    "# with solution by damped least squares\n",
    "# Note: symbol q for theta, temperture\n",
    "\n",
    "# generic time setup, except make two time arrays.\n",
    "# one in seconds for calulations, one in days for plotting\n",
    "M=1024;    # length of time series\n",
    "N=M;       # number of model parameters equal number of data\n",
    "Dtdays = 0.5; # sampling interval in days\n",
    "tdays = eda_cvec( Dtdays*np.linspace(0,Dtdays*(N-1),N) ); # time vector (days)\n",
    "Dtseconds = Dtdays*3600*24; # samping interval in seconds\n",
    "tseconds = Dtseconds*tdays;  # time vector (seconds)\n",
    "\n",
    "# materal constants\n",
    "rho = 1500; # density, kg/m3\n",
    "cp = 800; # heat capacity, J / kg-K\n",
    "kappa = 1.25e-6; # thermal diffusivity, m2/s\n",
    "z = 1.0; # depth of observation, m\n",
    "\n",
    "# impulse response, whcih is a analysic formula that's\n",
    "# fairly complicated, so break it into three pieces,\n",
    "# P1, P2 and P3\n",
    "# impulse respone is g, but time=0 will be handled separately\n",
    "tsecondswo0 = eda_cvec( tseconds[1:M,0] ); # time vector without zero\n",
    "P1 = (1.0/(rho*cp));\n",
    "P2 = (Dtseconds/sqrt(2*pi));\n",
    "P3 = np.reciprocal( np.sqrt(2*kappa*tsecondswo0));\n",
    "P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tsecondswo0) );\n",
    "g = eda_cvec( 0.0, P1 * P2 * np.multiply(P3 ,P4) );\n",
    "\n",
    "#true model parameters, heat production\n",
    "t1 = Dtdays*N/4;     # time of heat pulse 2\n",
    "t2 = Dtdays*7*N/16;  # time of heat pulse 2\n",
    "sd = Dtdays*N/256;   # variance\n",
    "Ah=10;               # amplitude\n",
    "P5 = np.exp(-np.power(tdays-t1,2) / (2*sd**2) ); # gaussian pulse 1\n",
    "P6 = np.exp(-np.power(tdays-t2,2) / (2*sd**2) ); # gaussian pulse 2\n",
    "htrue = Ah*P5 + 0.5*Ah*P6;  # sum pulses\n",
    "\n",
    "tmp = np.convolve(g.ravel(), htrue.ravel()); # convolve heat production with impulse response\n",
    "qtrue = eda_cvec(tmp[0:N]);  # keep only first N values\n",
    "\n",
    "# plot impulse response\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([tdays[0,0], tdays[M-1,0], 0, 1.1*np.max(g)] );\n",
    "plt.plot(tdays[0:M],g,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('g(t)');\n",
    "\n",
    "# plot true heat production\n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:M],htrue,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "\n",
    "# plot true temperature\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qtrue)] );\n",
    "plt.plot(tdays[0:N],qtrue,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('qtrue(t)');\n",
    "plt.show();\n",
    "\n",
    "# Toeplitz matrix version of impulse response\n",
    "Gcol = np.concatenate((g,np.zeros((N-M,1))),axis=0); # 1st col of Toeplitx matric\n",
    "v = np.zeros((1,1));\n",
    "v[0,0]=g[0,0];\n",
    "Grow = np.concatenate((v,np.zeros((1,M-1))),axis=1); # 1st row of Toeplitz matrix\n",
    "G = la.toeplitz( Gcol, Grow ); # impulse response as a Toeplitz matrix\n",
    "\n",
    "# temperature = G times heat_production\n",
    "qtrue2 = np.matmul(G,htrue);\n",
    "\n",
    "# create simulated observations, qobs = qtrue+noise\n",
    "sigmad=0.001*np.amax(qtrue);  # noise is 0.1% amplitude of signal\n",
    "qobs=qtrue+np.random.normal(0,sigmad,(N,1)); # synthetic observations\n",
    "\n",
    "# plot simulated observations\n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );\n",
    "plt.plot(tdays[0:N],qobs,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('qobs(t)');\n",
    "\n",
    "# solve for h given qobs and G using least squares\n",
    "GTG = np.matmul(G.T,G);  # G-transpose time G\n",
    "GTGmax = np.amax(GTG);   # maximum value\n",
    "\n",
    "# solve with damped least squares\n",
    "e2=1e-4*GTGmax;          # damping value as fraction of GTG\n",
    "hest1 = la.solve(GTG+e2*np.identity(M), np.matmul(G.T,qobs) ); \n",
    "\n",
    "# plot true heat production\n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],htrue,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "\n",
    "# plot estimated heat production\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],hest1,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('hest1(t)');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 40336 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 28916 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda07_02\n",
    "# example of a filter problem\n",
    "# using impulse repsonse of a heat-producing layer\n",
    "# with solution via generalized least squares\n",
    "# with the prior information of smoothness\n",
    "# Note: symbol q for theta, temperture.\n",
    "# This is the non-sparse implementation.\n",
    "\n",
    "# generic time setup, except make two time arrays.\n",
    "# one in seconds for calulations, one in days for plotting\n",
    "M=1024;    # length of time series\n",
    "N=M;       # number of model parameters equal number of data\n",
    "Dtdays = 0.5; # sampling interval in days\n",
    "tdays = eda_cvec( Dtdays*np.linspace(0,Dtdays*(N-1),N) ); # time vector (days)\n",
    "Dtseconds = Dtdays*3600*24; # samping interval in seconds\n",
    "tseconds = Dtseconds*tdays;  # time vector (seconds)\n",
    "\n",
    "# materal constants\n",
    "rho = 1500; # density, kg/m3\n",
    "cp = 800; # heat capacity, J / kg-K\n",
    "kappa = 1.25e-6; # thermal diffusivity, m2/s\n",
    "z = 1.0; # depth of observation, m\n",
    "\n",
    "# impulse response, whcih is a analysic formula that's\n",
    "# fairly complicated, so break it into three pieces,\n",
    "# P1, P2 and P3\n",
    "# impulse respone is g, but time=0 will be handled separately\n",
    "tsecondswo0 = eda_cvec( tseconds[1:M,0] ); # time vector without zero\n",
    "P1 = (1.0/(rho*cp));\n",
    "P2 = (Dtseconds/sqrt(2*pi));\n",
    "P3 = np.reciprocal( np.sqrt(2*kappa*tsecondswo0));\n",
    "P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tsecondswo0) );\n",
    "g = eda_cvec( 0.0, P1 * P2 * np.multiply(P3 ,P4) );\n",
    "\n",
    "#true model parameters, heat production\n",
    "t1 = Dtdays*N/4;     # time of heat pulse 2\n",
    "t2 = Dtdays*7*N/16;  # time of heat pulse 2\n",
    "sd = Dtdays*N/256;   # variance\n",
    "Ah=10;               # amplitude\n",
    "P5 = np.exp(-np.power(tdays-t1,2) / (2*sd**2) ); # gaussian pulse 1\n",
    "P6 = np.exp(-np.power(tdays-t2,2) / (2*sd**2) ); # gaussian pulse 2\n",
    "htrue = Ah*P5 + 0.5*Ah*P6;  # sum pulses\n",
    "\n",
    "# Toeplitz matrix version of impulse response\n",
    "Gcol = np.concatenate((g,np.zeros((N-M,1))),axis=0); # 1st row of Toeplitz matrix\n",
    "v = np.zeros((1,1));\n",
    "v[0,0]=g[0,0];\n",
    "Grow = np.concatenate((v,np.zeros((1,M-1))),axis=1); # 1st col of Toeplitx matrix\n",
    "G = la.toeplitz( Gcol, Grow );  # Toeplitz matrix version of impulse response\n",
    "\n",
    "# temperature = G times heat_production\n",
    "qtrue2 = np.matmul(G,htrue);\n",
    "\n",
    "# create simulated observations, qobs = qtrue+noise\n",
    "sigmad=0.01*np.amax(qtrue);  # noise leve is 1% of maximum temperature\n",
    "qobs=qtrue+np.random.normal(0,sigmad,(N,1)); # temperature = signal + noise\n",
    "\n",
    "# for comparison purposes, solve for h given qobs and g using least squares\n",
    "GTG = np.matmul(G.T,G);  # G_transpose time G\n",
    "GTGmax = np.amax(GTG);   # maximum value of GTG\n",
    "\n",
    "# solved by damped least squares\n",
    "e2=1e-4*GTGmax;  # damping is 0.01% maximum value of GTG\n",
    "# damped least squares solution\n",
    "hest2 = la.solve(GTG+e2*np.identity(M), np.matmul(G.T,qobs) );\n",
    "\n",
    "# generalized least squares with F = [G; H]  and  f = [qobs, 0]\n",
    "# second derivative smoothing in interior of interval; nothing on ends\n",
    "# H is a Toeplitz matrix with diaginals e, -2e and e, where e is a constant\n",
    "# that is proportional to 1/sqrt(variance_of_h).  It is being tuned\n",
    "# here by trail and error to achieve the desired degree of smoothness\n",
    "e=10*np.amax(G);# e is 10 times the maximum value of G\n",
    "# setup for a Toeplitx matrix\n",
    "Hcol = np.zeros((M-2,1));  # 1st col of H\n",
    "Hcol[0,0] = e;\n",
    "Hrow = np.zeros((1,M)); # first row of H\n",
    "Hrow[0,0] = e;\n",
    "Hrow[0,1] = -2*e;\n",
    "Hrow[0,2] = e;\n",
    "H = la.toeplitz( Hcol, Hrow ); # construct H\n",
    "h = np.zeros((M-2,1)); # since prior info is that 2nd derivative is zero, h=0\n",
    "F = np.concatenate( (G,H), axis=0 );  # combine G and H into F\n",
    "f = np.concatenate( (qobs,h), axis=0 ); # combine qobs and h into f\n",
    "FTF = np.matmul(F.T,F); # F_transpose times F\n",
    "FTf = np.matmul(F.T,f); # F_transpose times f\n",
    "hest3 = la.solve(FTF, FTf);  # generalized least squares solution\n",
    "\n",
    "# figure 1\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "# plot observations\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );\n",
    "plt.plot(tdays[0:N],qobs,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "plt.title('Damped Least Squares');\n",
    "\n",
    "# plot true heat production\n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],htrue,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "\n",
    "# plot estimated heat production\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],hest2,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('hest(t)');\n",
    "plt.show();\n",
    "\n",
    "# figure 2\n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "\n",
    "# plot observations\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );\n",
    "plt.plot(tdays[0:N],qobs,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "plt.title('Generalized Least Squares with 2nd Derivative Smoothing');\n",
    "\n",
    "# plot true heat production\n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],htrue,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "\n",
    "# plot estimated heat production\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],hest3,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('hest1(t)');\n",
    "plt.show();\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 39536 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 29324 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda07_03\n",
    "# example of a filter problem\n",
    "# using impulse repsonse of a heat-producing layer\n",
    "# with solution via generalized least squares\n",
    "# with the prior information of smoothness\n",
    "# Note: symbol q for theta, temperture\n",
    "# This is the sparse implementation.\n",
    "\n",
    "# generic time setup, except make two time arrays.\n",
    "# one in seconds for calulations, one in days for plotting\n",
    "M=1024;    # length of time series\n",
    "N=M;       # number of model parameters equal number of data\n",
    "Dtdays = 0.5; # sampling interval in days\n",
    "tdays = eda_cvec( Dtdays*np.linspace(0,Dtdays*(N-1),N) ); # time vector (days)\n",
    "Dtseconds = Dtdays*3600*24; # samping interval in seconds\n",
    "tseconds = Dtseconds*tdays;  # time vector (seconds)\n",
    "\n",
    "# materal constants\n",
    "rho = 1500; # density, kg/m3\n",
    "cp = 800; # heat capacity, J / kg-K\n",
    "kappa = 1.25e-6; # thermal diffusivity, m2/s\n",
    "z = 1.0; # depth of observation, m\n",
    "\n",
    "# impulse response, whcih is a analysic formula that's\n",
    "# fairly complicated, so break it into three pieces,\n",
    "# P1, P2 and P3\n",
    "# impulse respone is g, but time=0 will be handled separately\n",
    "tsecondswo0 = eda_cvec( tseconds[1:M,0] ); # time vector without zero\n",
    "P1 = (1.0/(rho*cp));\n",
    "P2 = (Dtseconds/sqrt(2*pi));\n",
    "P3 = np.reciprocal( np.sqrt(2*kappa*tsecondswo0));\n",
    "P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tsecondswo0) );\n",
    "g = eda_cvec( 0.0, P1 * P2 * np.multiply(P3 ,P4) );\n",
    "edafilterg = g;  # edafilterg is a global variable defined in GLSFilterMul()\n",
    "\n",
    "# true heat production, sum of two pulses\n",
    "t1 = Dtdays*N/4;      # time of first pulse\n",
    "t2 = Dtdays*7*N/16;   # time of second pulse\n",
    "sd = Dtdays*N/256;    # variance of pulse\n",
    "Ah=10;                # amplitude of pulse\n",
    "P5 = np.exp(-np.power(tdays-t1,2) / (2*sd**2) );\n",
    "P6 = np.exp(-np.power(tdays-t2,2) / (2*sd**2) );\n",
    "htrue = Ah*P5 + 0.5*Ah*P6;\n",
    "\n",
    "# predict true temperature by colvolving htrue with impulse response\n",
    "tmp = np.convolve(g.ravel(), htrue.ravel()); # colvolve\n",
    "qtrue = eda_cvec(tmp[0:N]); # keep only first N values\n",
    "\n",
    "# create simulated observations, qobs = qtrue+noise\n",
    "sigmad=0.01*np.amax(qtrue);  # stddev of noise is 1% of signal\n",
    "qobs=qtrue+np.random.normal(0,sigmad,(N,1));\n",
    "\n",
    "# for comparison purposes, solve for h given qobs and G using damped least squares\n",
    "# build Toeplitz matrix G equivalent to g\n",
    "Gcol = np.concatenate((g,np.zeros((N-M,1))),axis=0); # 1st col of Toeplitz matrix\n",
    "v = np.zeros((1,1));\n",
    "v[0,0]=g[0,0];\n",
    "Grow = np.concatenate((v,np.zeros((1,M-1))),axis=1); # 1st row of Toeplitz matrix\n",
    "G = la.toeplitz( Gcol, Grow );\n",
    "qtrue2 = np.matmul(G,htrue); # predict true temperature from true heat prodcution\n",
    "# solve with damped least squares\n",
    "GTG = np.matmul(G.T,G); # G_transpose times G\n",
    "GTGmax = np.amax(GTG);  # maximum value of GTG\n",
    "e2=1e-4*GTGmax; # damping is 0.01% of maximum value of GTG\n",
    "# damped least squares osolution\n",
    "hest2 = la.solve(GTG+e2*np.identity(M), np.matmul(G.T,qobs) );\n",
    "\n",
    "# generalized leasr squares with F = [G; H]  and  f = [qobs, 0]\n",
    "# second derivative smoothing in interior of interval; nothing on ends\n",
    "e=10.0*np.amax(G);\n",
    "edaHsparse = sparse.diags([e, -2*e, e], [0, 1, 2], shape=(M-2, M));\n",
    "h = np.zeros((M-2,1));\n",
    "\n",
    "# set up right hand side, F'f = GT qobs + HT h\n",
    "temp = eda_cvec( np.convolve(np.flipud(edafilterg).ravel(),qobs.ravel()) );\n",
    "FTfa = eda_cvec( temp[N-1:N+M-1,0] );\n",
    "FTfb = edaHsparse.transpose()*h;\n",
    "FTf=FTfa+FTfb;\n",
    "\n",
    "# define linear operator needed for conjugate gradienet solver\n",
    "LO=las.LinearOperator(shape=(M,M),matvec=GLSFilterMul,rmatvec=GLSFilterMul);\n",
    "\n",
    "# solve for estimated model parameters using congugate gradient algrorithm\n",
    "hest3 = np.zeros((M,1));\n",
    "q=las.cg(LO,FTf,tol=1e-6, maxiter=(3*(N+M)+100));\n",
    "# note that q is a tupple, the first element of which is the solution\n",
    "hest3 = eda_cvec(q[0]);\n",
    "\n",
    "# figure 1\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "# plot observations\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );\n",
    "plt.plot(tdays[0:N],qobs,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "plt.title('Damped Least Squares');\n",
    "\n",
    "# plot true heat production\n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],htrue,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "\n",
    "# plot estimated heat production\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],hest2,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('hest(t)');\n",
    "plt.show();\n",
    "\n",
    "# figure 2\n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "\n",
    "# plot observations\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );\n",
    "plt.plot(tdays[0:N],qobs,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "plt.title('Generalized Least Squares with 2nd Derivative Smoothing');\n",
    "\n",
    "# plot true heat production\n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],htrue,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "\n",
    "# plot estimated heat production\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],hest3,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('hest1(t)');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 31396 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24964 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda07_04\n",
    "# example of estimating a filter r that predicts river\n",
    "# discharge d from precipitation data g; that is, g*r=d\n",
    "\n",
    "# in this example, the filter length is chosen to be\n",
    "# M<<N, to embody the idea that the response of the river\n",
    "# (meaning the length M of the filter r) is short\n",
    "\n",
    "M=30; # length of filter r\n",
    "\n",
    "# load merged discharge and precipitation data\n",
    "# column 1: time in inter count of days starting 1/1/2002\n",
    "# column 2: Hudson River discharge in cfs\n",
    "# column 3: Albany NY precipitation in inches\n",
    "D = np.genfromtxt('precip_discharge_data.txt', delimiter='\\t')\n",
    "[Nlong, K]=D.shape;\n",
    "tlong = eda_cvec(D[:,0]); # time in days\n",
    "dlong = eda_cvec(D[:,1]/35.3146); # discharge in cubic meters per second\n",
    "glong = eda_cvec(D[:,2]*25.4) # precipitation, im millimeters pper day;\n",
    "Dt = tlong[1,0]-tlong[0,0]; # sampling interval\n",
    "\n",
    "# in this example, we use only a 101 day portion of the\n",
    "# dataset, mainly so that when we plot it, storm events\n",
    "# are clearly visible\n",
    "\n",
    "# select small segment of time, discharge and precip data\n",
    "N=101;\n",
    "istart = 900-1;\n",
    "t = eda_cvec( tlong[istart:istart+N,0] );\n",
    "d = eda_cvec( dlong[istart:istart+N,0] );\n",
    "g = eda_cvec( glong[istart:istart+N,0] );\n",
    "tmin = t[0,0];    # start time\n",
    "tmax = t[N-1,0];  # end time\n",
    "edafilterg = g; # edafilterg is a global variable defined in GLSFilterMul()\n",
    "\n",
    "# figure 1: discharge and precipitation data\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "# plot observations\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([tmin, tmax, -20, 60 ] );\n",
    "plt.plot(t,g,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('precipitation (mm/day)');\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([tmin, tmax, -1000, 2000 ] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('discharge (m3/s)');\n",
    "plt.show();\n",
    "\n",
    "# Given the filter equation g*r=d, solve for r given d and g\n",
    "# using genralized least squares\n",
    "\n",
    "# Note that the convolution relationship g*r=d is equivalent\n",
    "# to the matrix equation G r = d where G is the (g*) operator,\n",
    "# Since we are using filterfun(), we do not need to construct\n",
    "# G.  We do, however, need to implement the prior information\n",
    "# equation Hr=h.\n",
    "\n",
    "# This code can toggle between the prior information of r\n",
    "# being small, its first derivative being small and its second\n",
    "# derivative being small, depending on whether DERIVATIVE is 0,\n",
    "# 1 or 2. (However, the results turn out not to be especially\n",
    "# sensitive to the choice).\n",
    "DERIVATIVE=1;\n",
    "if( DERIVATIVE==2 ):\n",
    "    # The prior information equation H*r = h is implemented\n",
    "    # for the smoothness condition of the second derivative\n",
    "    # being close to zero\n",
    "    e=2*np.amax(g); # the coefficient was manually tuned\n",
    "                    # to produce a reasonably smooth filter\n",
    "    K=M;\n",
    "    L=N+K;\n",
    "    # second derivative operator; no end conditions\n",
    "    edaHsparse = sparse.diags([e, -2*e, e], [0, 1, 2], shape=(M-2, M));\n",
    "    h = np.zeros((M-2,1));\n",
    "elif (DERIVATIVE==1):\n",
    "    # The prior information equation H*r = h is implemented\n",
    "    # for the smoothness condition of the first derivative\n",
    "    # being close to zero\n",
    "    e=4.0*np.amax(g); # the coefficient was manually tuned\n",
    "                      # % to produce a reasonably smooth filter\n",
    "    edaHsparse = sparse.diags([-e, e], [0, 1], shape=(M-1, M));\n",
    "    h = np.zeros((M-1,1));\n",
    "else: # DERRIVATIVE == 0\n",
    "    # The prior information equation H*r = h is implemented\n",
    "    # for the smallness condition of the filter is close to zero\n",
    "    e=0.1*max(abs(g)); # the coefficient was manually tuned\n",
    "                       # to produce a reasonably small filter\n",
    "    edaHsparse = sparse.diags([e], [0], shape=(M, M));\n",
    "    h = np.zeros((M,1));\n",
    "    \n",
    "# set up right hand side, F'f = GT qobs + HT h\n",
    "temp = eda_cvec( np.convolve(np.flipud(edafilterg).ravel(),d.ravel()) );\n",
    "FTfa = eda_cvec( temp[N-1:N+M-1,0] );\n",
    "FTfb = edaHsparse.transpose()*h;\n",
    "FTf=FTfa+FTfb;\n",
    "\n",
    "# define linear operator needed for conjugate gradienet solver\n",
    "LO=las.LinearOperator(shape=(M,M),matvec=GLSFilterMul,rmatvec=GLSFilterMul);\n",
    "\n",
    "# solve for estimated model parameters using congugate gradient algrorithm\n",
    "# note las.cg() returns a tupple, the first element of which is the solution\n",
    "q=las.cg(LO,FTf,tol=1e-6, maxiter=(3*(N+M)+100));\n",
    "r = eda_cvec( q[0] );\n",
    "\n",
    "# filter time axis\n",
    "tr = eda_cvec( Dt*np.linspace(0,M-1,M) );\n",
    "\n",
    "# predicted discharge, convolve and truncate to length N\n",
    "temp = eda_cvec( np.convolve(g.ravel(),r.ravel()) );\n",
    "dpre = eda_cvec( temp[0:N,0] );\n",
    "\n",
    "# figure 2: filter and predicted discharge \n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "\n",
    "# plot observations\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([0, M, -5, 10 ] );\n",
    "plt.plot(tr,r,'k-');\n",
    "plt.plot([tr[0,0], tr[M-1,0]],[0,0],'k:');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('filter r(t)');\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([tmin, tmax, -1000, 2000 ] );\n",
    "plt.plot(t,d,'k:');\n",
    "plt.plot(t,dpre,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('discharge (m3/s)');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 20964 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 29260 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda07_05\n",
    "# prediction error filter (pef) for Neuse Rievr Hydrograph\n",
    "# this code uses generalized least squares to implement\n",
    "# the prior information that pef(1)=(-1).  Two solutions\n",
    "# are calculated, one via standard matrix algebra,\n",
    "# (FTF)\\(FTf) and the other via the biconjugate gradient method\n",
    "\n",
    "# read in data\n",
    "D = np.genfromtxt('neuse.txt', delimiter='\\t')\n",
    "[N, K]=D.shape;\n",
    "t = eda_cvec( D[:,0] ); # time in days\n",
    "d = eda_cvec( D[:,1] ); # discharge in cubic meters per second\n",
    "Dt = t[1,0]-t[0,0];     # sampling interval\n",
    "tmin = t[0,0];          # start time\n",
    "tmax = t[N-1,0];        # end time\n",
    "\n",
    "# sizes of various arrays\n",
    "M=100; # length of pef is user tunable\n",
    "K=1;   # one piece of prior information\n",
    "L=N+K; # total number of equations\n",
    "\n",
    "# for a prediction error filter, pef, the data is the filter, g\n",
    "g = d;\n",
    "edafilterg = d; # edafilterg is a global variable defined in GLSFilterMul()\n",
    "\n",
    "# prior information that first element of the pef is (-1)\n",
    "e=(10^1)*np.amax(g); # large weight\n",
    "H=np.zeros((K,M));   # the pror information equation, Hm=h, such that m(1)=(-1);\n",
    "H[0,0]=e;            \n",
    "h=eda_cvec( (-e) );\n",
    "\n",
    "# Data Kernel (needed only for standard matrix algebra calculation\n",
    "# data kernel is Toeoplitz, so use standard setup\n",
    "v = np.zeros((1,1));\n",
    "v[0,0]=g[0,0];\n",
    "Grow = np.concatenate((v,np.zeros((1,M-1))),axis=1);\n",
    "G = la.toeplitz( g, Grow );\n",
    "# for pef, Gm=rhs=0\n",
    "\n",
    "# F needed only for standard matrix algebra calculation\n",
    "F = np.concatenate( (G, H), axis=0 );\n",
    "rhs = np.zeros((N,1));\n",
    "f = np.concatenate( (rhs, h), axis=0 );\n",
    "\n",
    "# standard GLS solution\n",
    "FTF = np.matmul( F.T, F );  # F_transpose times F\n",
    "FTf = np.matmul( F.T, f );  # F_transpose times f\n",
    "pef1 = la.solve(FTF, FTf ); # least squares solution\n",
    "\n",
    "# now set up for biconjugate gradient method\n",
    "\n",
    "# Hm=h\n",
    "e=(10^1)*np.amax(g); # large weight\n",
    "# edaHsparse is global variable defined in GLSFilterMul()\n",
    "edaHsparse = sparse.diags([e], [0], shape=(1, M)); \n",
    "h=np.zeros((1,1));\n",
    "h[0,0]=(-e);\n",
    "\n",
    "# FT f can be done analytucally:\n",
    "# FTf = GT 0 + HT h = 0 + [e 0 0 ... 0]T [-e] = [-e^2 0 0 ... 0]T\n",
    "FTf = np.zeros((M,1));\n",
    "FTf[0,0] = -(e**2);\n",
    "           \n",
    "# define linear operator needed for conjugate gradienet solver\n",
    "LO=las.LinearOperator(shape=(M,M),matvec=GLSFilterMul,rmatvec=GLSFilterMul);\n",
    "\n",
    "# solve for estimated model parameters using congugate gradient algrorithm\n",
    "r = np.zeros((M,1));\n",
    "# conjugate gradient solver returns soln vector as first element of a tupple\n",
    "q=las.cg(LO,FTf,tol=1e-6, maxiter=(3*(N+M)+100));\n",
    "pef2 = eda_cvec(q[0]);\n",
    "\n",
    "# compute prediction error by colvolving g with pef, then truncate to length N\n",
    "temp = np.convolve(pef2.ravel(),g.ravel());\n",
    "perror = eda_cvec( temp[0:N] );\n",
    "           \n",
    "# short time axis\n",
    "tpef = eda_cvec( t[0:M,0] );\n",
    "           \n",
    "# plot pef\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "# plot prediction error filters\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([ 0, Dt*(M-1), np.min(pef1), np.max(pef1)] );\n",
    "plt.plot(tpef,pef1,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('pef(t) (std)');\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([ 0, Dt*(M-1), np.min(pef1), np.max(pef1)] );\n",
    "plt.plot(tpef,pef2,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('pef(t) (std)');\n",
    "plt.show();\n",
    "\n",
    "# plot hydrograph data and prediction error for first year\n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "\n",
    "Np = 365;\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([tmin, tmin+Dt*Np, -5000, 15000]);\n",
    "plt.plot(t[0:Np,0],d[0:Np,0],'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('d(t)');\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([tmin, tmin+Dt*Np, -5000, 15000]);\n",
    "plt.plot(t[0:Np,0],perror[0:Np,0],'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('e(t)');\n",
    "plt.show();        \n",
    "           "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30340 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 14240 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda07_06 Script supporting CribSheet 7.2\n",
    "# Least squares filter estimation of filter m\n",
    "# obeying g*m=d\n",
    "\n",
    "# definition of time\n",
    "N=101;  # number of data\n",
    "Dt = 1; # sampling interval\n",
    "t = eda_cvec( Dt*np.linspace(0,N-1,N) );\n",
    "tmin = t[0,0];   # start time\n",
    "tmax = t[N-1,0]; # end time\n",
    "\n",
    "# make some synthetic timeseries g\n",
    "# draw from an exponential distribution so its all positive\n",
    "g = np.random.exponential(1,(N,1));\n",
    "edafilterg = g;  # edafilterg is a global variable defined in GLSFilterMul()\n",
    "\n",
    "# exemplary true filter is decaying exponential\n",
    "M=10;    # length of filter\n",
    "c=5*Dt;  # decay constant\n",
    "mtrue = eda_cvec( 0.3*np.exp(-t[0:M,0]/c) );\n",
    "\n",
    "# the true data is g covolved with m; then truncate to length N\n",
    "temp = np.convolve(g.ravel(),mtrue.ravel());\n",
    "dtrue = eda_cvec( temp[0:N] );\n",
    "\n",
    "# the observed data is the true data plus noise\n",
    "# the noise is Normally distributed with zero mean and variance sigmad^2\n",
    "sigmad = 0.1;  # standard deviation of noise\n",
    "noise = np.random.normal(0,sigmad,(N,1)); # random noise\n",
    "d = dtrue + noise; # data = signal + noise\n",
    "\n",
    "# now invert for the filter using the observed data and the time series g\n",
    "\n",
    "# prior infoirmation of smoothness, with no end conditions\n",
    "epsilon = 0.001;       # 1 / (prior standard deviation) \n",
    "e = epsilon/(Dt**2);   # include factor of (Dt**(-2)) from 2nd derivative\n",
    "# make Toeplitz matrix of second derivatives\n",
    "edaHsparse = sparse.diags([e, -2*e, e], [0, 1, 2], shape=(M-2, M));\n",
    "h = np.zeros((M-2,1)); # prior info that second derivative is zero,\n",
    "\n",
    "# set up right hand side, F'f = GT qobs + HT h\n",
    "temp = np.convolve(np.flipud(edafilterg).ravel(),d.ravel());\n",
    "FTfa = eda_cvec(temp[N-1:N+M-1]);\n",
    "FTfb = edaHsparse.transpose()*h;\n",
    "FTf=FTfa+FTfb;\n",
    "\n",
    "# define linear operator needed for conjugate gradienet solver\n",
    "LO=las.LinearOperator(shape=(M,M),matvec=GLSFilterMul,rmatvec=GLSFilterMul);\n",
    "\n",
    "# solve for estimated model parameters using congugate gradient algrorithm\n",
    "# returned value 1 is a tuple whose first element is the solution vector\n",
    "q=las.cg(LO,FTf,tol=1e-6, maxiter=(3*(N+M)+100));\n",
    "mest = eda_cvec(q[0]);\n",
    "\n",
    "# predicted data by colvolving mest and g, then truncat to length N\n",
    "temp = np.convolve(g.ravel(),mest.ravel());\n",
    "dpre = eda_cvec( temp[0:N] );\n",
    "\n",
    "# compute the error\n",
    "e = d-dpre;             # individual errors\n",
    "E = np.matmul(e.T, e);  # total error\n",
    "\n",
    "# Figure 1\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([tmin, tmax, -5, 5] );\n",
    "plt.plot(t,g,'k-');\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('g(t)');\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([tmin, tmax, -5, 5] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.plot(t,dpre,'k:');\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('d(t)');\n",
    "plt.show();\n",
    "\n",
    "# Figure 2\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([t[0,0], t[M-1,0], 0, 0.5] );\n",
    "plt.plot( t[0:M], mtrue.ravel(),'k-');\n",
    "plt.plot( t[0:M], mest.ravel(), 'k-', color = (0.8,0.8,0.8));\n",
    "plt.xlabel('t');\n",
    "plt.ylabel('m(t)');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 roots not outside the unit circle\n",
      "error of filter rebuilt from its roots: 0.000000\n",
      "error of inverse filter: 0.000000\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 15764 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda07_07\n",
    "# use z-transforms to:\n",
    "#    check if a filter is minumum phase\n",
    "#    represent it as a cascade of length-2 filters\n",
    "#    compute the inverse filter\n",
    "\n",
    "# define a filter\n",
    "g = eda_cvec( [6, 5.5, 5.3, 5, 4, 3, 2] );\n",
    "N,K = np.shape(g);\n",
    "gN = g[N-1,0];\n",
    "\n",
    "# find roots of g\n",
    "pol = (np.flipud(g)).ravel(); # pol is a list of coefficients\n",
    "r = np.roots(pol);  # r is a list of roots\n",
    "\n",
    "# check size of roots\n",
    "badroots=0;\n",
    "for j in range(N-1):\n",
    "    if( abs(r[j]) <= 1 ):\n",
    "        badroots=badroots+1;\n",
    "\n",
    "printstr = \"%d roots not outside the unit circle\" % (badroots);\n",
    "print(printstr);\n",
    "\n",
    "# rebuild filter from its roots\n",
    "gp = np.zeros((1,1), dtype=np.complex);\n",
    "gp[0,0]=gN;\n",
    "for j in range(N-1):\n",
    "    temp = np.convolve(gp.ravel(),[-r[j],1]);\n",
    "    K = len(temp);\n",
    "    gp = np.zeros((K,1), dtype=np.complex);\n",
    "    gp[:,0] = temp;\n",
    "    \n",
    "# there should be no imiginary part\n",
    "gp = np.real(gp);\n",
    "\n",
    "# check that it worked by computing error\n",
    "E1 = np.matmul( (g-gp).T , (g-gp) );\n",
    "printstr = \"error of filter rebuilt from its roots: %f\" % (E1);\n",
    "print(printstr);\n",
    "\n",
    "# construct inverse filter, gi, of length Ni\n",
    "Ni = 50;\n",
    "gi = np.zeros((Ni,1),dtype=np.complex);\n",
    "gi[0,0]=1/gN;\n",
    "# filter cascade, one filter per loop\n",
    "for j in range(N-1):\n",
    "    # construct inverse filter of a length-2 filter\n",
    "    temp = np.zeros((Ni,1),dtype=np.complex);\n",
    "    temp[0,0] = 1;\n",
    "    for k  in range(1,Ni):\n",
    "        temp[k]=temp[k-1]/r[j];\n",
    "    temp = -temp/r[j];\n",
    "    temp2 = np.convolve(gi.ravel(),temp.ravel());\n",
    "    gi[0:Ni,0]=temp2[0:Ni];\n",
    "\n",
    "# delete imaginary part (which should be zero)\n",
    "gi = np.real(gi);\n",
    "\n",
    "# test quality of inverse filter; truncate to length Ni\n",
    "temp=np.convolve(g.ravel(),gi.ravel());\n",
    "gig = eda_cvec( temp[0:Ni] );\n",
    "\n",
    "# impulse function\n",
    "d = np.zeros((Ni,1));\n",
    "d[0]=1;\n",
    "E2 = np.matmul( (gig-d).T, (gig-d) );  # error\n",
    "printstr = \"error of inverse filter: %f\" % (E2);\n",
    "print(printstr);\n",
    "\n",
    "# plot results\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "# plot filter\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([0, Ni, np.min(g), np.max(g)] );\n",
    "gx = np.concatenate( (g, np.zeros((Ni-N,1))), axis=0);\n",
    "plt.plot( np.linspace(0,Ni-1,Ni), gx.ravel(),'k-');\n",
    "plt.xlabel('element j');\n",
    "plt.ylabel('g(j)');\n",
    "\n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis([0, Ni, np.min(gi), np.max(gi)] );\n",
    "plt.plot( np.linspace(0,Ni-1,Ni), gi.ravel(),'k-');\n",
    "plt.xlabel('element j');\n",
    "plt.ylabel('ginv(j)');\n",
    "\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis([0, Ni, -1.2, 1.2] );\n",
    "plt.plot( np.linspace(0,Ni-1,Ni), gig.ravel(),'k-');\n",
    "plt.xlabel('element j');\n",
    "plt.ylabel('[g*ginv](j)');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 24624 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sum of elements of the unit spike: 1.000000\n",
      "sum of elements of the filtered spike: 1.000000\n"
     ]
    }
   ],
   "source": [
    "# eda07_08\n",
    "# example of recursive smoothing filter\n",
    "# this version used for-loops, for demonstration purposes\n",
    "# the method in eda07_09 is prefered\n",
    "\n",
    "# create two synthetic datasets\n",
    "# h1, a unit spike\n",
    "# h2, random noise with zero mean and unit variance\n",
    "N=101;   # number of data\n",
    "Dt=1.0;  # sampling interval\n",
    "t = eda_cvec( Dt*np.linspace(0,N-1,N) );  # time vector\n",
    "h1=np.zeros((N,1)); # h1 is a time series with a spike in its middle\n",
    "No2 = floor(N/2);   # middle at N/2\n",
    "h1[No2,0]=1.0;      # spike has unit amplitude\n",
    "h2=np.random.normal(0.0,1.0,(N,1));  # h2 is a time series of random noise\n",
    "\n",
    "# apply recursive filter\n",
    "q1=np.zeros((N,1));\n",
    "q1[0,0]=0.5*h1[0,0];\n",
    "q2=np.zeros((N,1));\n",
    "q2[0,0]=0.5*h2[0,0];\n",
    "for j in range(2,N):\n",
    "    q1[j,0]=0.5*(h1[j,0]+q1[j-1,0]);\n",
    "    q2[j,0]=0.5*(h2[j,0]+q2[j-1,0]);\n",
    "    \n",
    "# plot results\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [t[0,0], t[N-1,0], -1.1, 1.1] );\n",
    "plt.plot( t, h1, 'k.');\n",
    "plt.plot( t, q1, 'k-', color=(0.7,0.7,0.7));\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('h1(t) and q1(t)');\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([t[0,0], t[N-1,0], -5, 5]);\n",
    "plt.plot( t, h2, 'k.');\n",
    "plt.plot( t, q2, 'k-', color=(0.7,0.7,0.7));\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('h2(t) and q2(t)');\n",
    "plt.show();\n",
    "\n",
    "displaystr = \"sum of elements of the unit spike: %f\" % (np.sum(h1));\n",
    "print(displaystr);\n",
    "displaystr = \"sum of elements of the filtered spike: %f\" % (np.sum(q1));\n",
    "print(displaystr);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 24152 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sum of elements of the unit spike: 1.000000\n",
      "sum of elements of the filtered spike: 1.000000\n"
     ]
    }
   ],
   "source": [
    "# eda07_09\n",
    "# example of recursive smoothing filter\n",
    "# this version used scipy.signal.lfilter()\n",
    "\n",
    "# create two synthetic datasets\n",
    "# h1, a unit spike\n",
    "# h2, random noise with zero mean and unit variance\n",
    "N=101;   # number of data\n",
    "Dt=1.0;  # sampling interval\n",
    "t = eda_cvec( Dt*np.linspace(0,N-1,N) ); # time vector\n",
    "h1=np.zeros((N,1)); # h1 is a time series with a spike in its middle\n",
    "No2 = floor(N/2);   # middle position\n",
    "h1[No2,0]=1.0;      # unit spike\n",
    "h2=np.random.normal(0.0,1.0,(N,1)); # h2 is a random time series\n",
    "\n",
    "Nu = 2;                      # length of filter u\n",
    "u = eda_cvec( [0.5, 0.0] );  # filter u\n",
    "\n",
    "Nv = 2;                      # length of filter v\n",
    "v = eda_cvec([1.0, -0.5] );  # filter v\n",
    "\n",
    "temp1 = sg.lfilter(u.ravel(),v.ravel(),h1.ravel()); # apply filter to h1\n",
    "q1 = eda_cvec( temp1[0:N] ); # truncate to length N\n",
    "\n",
    "temp2 = sg.lfilter(u.ravel(),v.ravel(),h2.ravel()); # apply filter to h1\n",
    "q2 = eda_cvec( temp2[0:N] ); # truncate to length N\n",
    "\n",
    "# plot results\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [t[0,0], t[N-1,0], -1.1, 1.1] );\n",
    "plt.plot( t, h1, 'k.');\n",
    "plt.plot( t, q1, 'k-', color=(0.7,0.7,0.7));\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('h1(t) and q1(t)');\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([t[0,0], t[N-1,0], -5, 5]);\n",
    "plt.plot( t, h2, 'k.');\n",
    "plt.plot( t, q2, 'k-', color=(0.7,0.7,0.7));\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('h2(t) and q2(t)');\n",
    "plt.show();\n",
    "\n",
    "displaystr = \"sum of elements of the unit spike: %f\" % (np.sum(h1));\n",
    "print(displaystr);\n",
    "displaystr = \"sum of elements of the filtered spike: %f\" % (np.sum(q1));\n",
    "print(displaystr);"
   ]
  },
  {
   "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
}
