{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_04_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.linalg as la\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",
    "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( 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  0.]\n",
      " [ 1.  1.]\n",
      " [ 1.  2.]\n",
      " [ 1.  3.]\n",
      " [ 1.  4.]\n",
      " [ 1.  5.]\n",
      " [ 1.  6.]\n",
      " [ 1.  7.]\n",
      " [ 1.  8.]\n",
      " [ 1.  9.]\n",
      " [ 1. 10.]]\n"
     ]
    }
   ],
   "source": [
    "# eda04_01: data kernel for straight-line linear model\n",
    "\n",
    "# make the auxillary variable x\n",
    "N=11;\n",
    "x = eda_cvec( np.linspace(0,N-1,N) );\n",
    "\n",
    "# manualy set the columns of G\n",
    "G = np.zeros((N,2));\n",
    "# Note: In order to set a column of a matrix\n",
    "# to a Nx1 column vector, the view of the matrix\n",
    "# must have a shape of Nx1, e.g. G[0:N,0:1]\n",
    "G[0:N,0:1] = np.ones((N,1));\n",
    "G[0:N,1:2] = x;\n",
    "\n",
    "print(G);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.   0.   0.]\n",
      " [  1.   1.   1.]\n",
      " [  1.   2.   4.]\n",
      " [  1.   3.   9.]\n",
      " [  1.   4.  16.]\n",
      " [  1.   5.  25.]\n",
      " [  1.   6.  36.]\n",
      " [  1.   7.  49.]\n",
      " [  1.   8.  64.]\n",
      " [  1.   9.  81.]\n",
      " [  1.  10. 100.]]\n"
     ]
    }
   ],
   "source": [
    "# eda04_02: data kernel for quadratic-curve linear model\n",
    "\n",
    "# make the auxillary variable x\n",
    "N=11;\n",
    "x = eda_cvec( np.linspace(0,N-1,N) );\n",
    "\n",
    "# manualy set the columns of G\n",
    "G = np.zeros((N,3));\n",
    "# Note: In order to set a column of a matrix\n",
    "# to a Nx1 column vector, the view of the matrix\n",
    "# must have a shape of Nx1, e.g. G[0:N,0:1]\n",
    "G[0:N,0:1] = np.ones((N,1));\n",
    "G[0:N,1:2] = x;\n",
    "G[0:N,2:3] = np.power(x,2);\n",
    "\n",
    "print(G);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAADnCAYAAABWihTLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAJ90lEQVR4nO3dPY5cTRUG4Kr56ZmxnSJkicgpGUgmJPhWQMYWEHtgC0RIrACJjBUQkBAgkAhICZ2QeOzxTM9P9/0iJLpu2XNpdd/2236erEpn7pQd+VX5nKrDMBQAAADynBz6AAAAAGxHoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQSqADAAAIJdABAACEEugAAABCCXQAAAChBDoAAIBQAh0AAEAogQ4AACCUQAcAABBKoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQSqADAAAIdbarD9Vah3bv3bt3u/o834DXr1/XQ58BAACSuKEDAAAIJdABAACEEugAAABCCXQAAAChdjYUpef+/n5jfXFxsc9fBwAA8E1xQwcAABBKoAMAAAgl0AEAAISatYdusViMamr1ljQAAMA23NABAACEEugAAABCCXQAAAChBDoAAIBQex2KslwuN9avXr0a1RiKAgAAsB03dAAAAKEEOgAAgFACHQAAQKhZe+jW6/Wo5uREpgQAANiGNAUAABBKoAMAAAgl0AEAAIQS6AAAAELtdSjK/f39xnoYhlFNu+ehcQAAgGnc0AEAAIQS6AAAAEIJdAAAAKFm7aFbrVbjA5zt9QgAAABHyw0dAABAKIEOAAAglEAHAAAQSqADAAAINetQlKenp1HNYrHYWHtYHAAAYBo3dAAAAKEEOgAAgFACHQAAQKi99tAtl8uNde9h8fV6vbHu9dDpqwMAABhzQwcAABBKoAMAAAgl0AEAAIQS6AAAAELNOhSl97D4MAz7PAIAAMDRckMHAAAQSqADAAAIJdABAACE2msP3f39/ca610PXPize66nzsDgAAMCYGzoAAIBQAh0AAEAogQ4AACCUQAcAABBq1qEoDw8Po5rVarWxPj09HdW0Q1EMSQEAAHBDBwAAEEugAwAACCXQAQAAhNprD91yudxYPz4+jmraHrrew+Ltnh46AAAAN3QAAACxBDoAAIBQAh0AAEAogQ4AACDUrA+Lt+tSthuK0qsxKAUAAPjWuKEDAAAIJdABAACEEugAAABCzfqw+MPDw6imfWy87akrpZSTk83cqYcOAADADR0AAEAsgQ4AACCUQAcAABBKoAMAAAg168PiU4aiPD09jWrOzjaPOeXxcUNSAACAY+eGDgAAIJRABwAAEEqgAwAACDXrw+LtupRxX13vYfF2r9cf1+ura+mrAwAAjokbOgAAgFACHQAAQCiBDgAAIJRABwAAEGrWh8XbdSnjoSjtQ+OllHJ+fr6xPj09HdVMGYrSMiQFAABI5oYOAAAglEAHAAAQSqADAAAIJdABAACEmnUoyt3d3bM1vaEoq9Xqi+uek5NxVm2HoPQGqRiUAgAApHBDBwAAEEqgAwAACCXQAQAAhPrqeujah8Z7e1P646bofaftq9NTBwAAfK3c0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFB7HYrSDkHpDUVp95bL5ajm4uJiY312Nj52O+CkN8ykHXiyXq9HNVOGoBiUAgAAfA3c0AEAAIQS6AAAAEIJdAAAAKFmfVj89vZ2VNPu9XroLi8vN9anp6ejmt4j4a22967XCzflO1PoswMAAPbNDR0AAEAogQ4AACCUQAcAABBKoAMAAAh18KEonz59eramfVi8NxRlmyEkve9M+a7BKQAAwNfADR0AAEAogQ4AACCUQAcAABBq1h66tl+ulFJubm6erdlVD90wDF9c977d65drf25XPXU9+uwAAIDPcUMHAAAQSqADAAAIJdABAACEEugAAABC7XUoynK53Fi3A1B6ex8/fhzVnJ+fb6ynDApZr9fP7i0Wi1FNOxTl7Gz8V9QOQekNV2nP2Dvz1L1tGKYCAADHzw0dAABAKIEOAAAglEAHAAAQ6uAPi19fX2+sr66uRjVtH9uU/rDVajXaa3voejXt72r790oZ99n1Hjpv++x6j49PeZB8Sp9dr6bt69u2p04vHgAAfL3c0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFCzPizeezS8HYpyeXk5qmmHh/QeDX96etpYPzw8jGravYuLi1FN+/t7D4u3g1J6Ne2glN4AlCnDVHpDSabUtLZ9xHxXNQAAwO65oQMAAAgl0AEAAIQS6AAAAELN+rB4r4fu/fv3mwfq9KO1PXOPj4/P/q67u7tRzYsXLzbWvX69dq/XZ9f20PUeH2//HFP67Hp7u+qz6/XwTXmgfErNcz+zbQ0AAPBlbugAAABCCXQAAAChBDoAAIBQAh0AAECoWR8Wv7m5GdUsFouN9TAMo5r2QfDb29tRTTtw5eXLl6OadijK1dXVqGbKUJR2r/0zlLLd4JRe3baDU6Y8bD5lcMq+hqv0vgMAAPx//KsaAAAglEAHAAAQSqADAAAINevD4tfX16Oatmeu/ZlSxr132/bHtXvbPiy+TQ/dlJre3rZ9du3PbVvT9rrts1/vzZs3oz0AAODz3NABAACEEugAAABCCXQAAAChBDoAAIBQsw5FmVLz4cOHUc2UISTtXm+YyZTHvtvvTBlcsqvhJr263s9NqWkHlUz5zq4eMe/9rinDVd6+fTvaAwAAPs8NHQAAQCiBDgAAIJRABwAAEKq2D3tv/aFad/MhvlnDMNRDnwEAAJK4oQMAAAgl0AEAAIQS6AAAAEIJdAAAAKF2NhQFAACAebmhAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQSqADAAAIJdABAACEEugAAABCCXQAAAChBDoAAIBQAt2e1Vp/WGv9Q63137XWv9da/1pr/cWhzwUAAOQT6Pao1lpLKX8qpfxlGIY3wzD8tJTyy1LKjw57MgAA4BjUYRgOfYajVWv9rpTym2EYfn7oswAAAMfHDd1+/biU8o9DHwIAADhOAt2Maq2/q7X+s9b6t0OfBQAAyCfQ7de/Sik/+e9iGIZfl1K+K6X84GAnAgAAjoZAt19/LqVc1lp/9T97Lw51GAAA4LgYirJntdbXpZTfllJ+Vkr5TynlUynl98Mw/PGgBwMAAOIJdAAAAKH8l0sAAIBQAh0AAEAogQ4AACCUQAcAABBKoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQ6nvSIM38CIdJAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_03: data kernel for polynomial linear model\n",
    "\n",
    "# make the auxillary variable x on varange (0,1)\n",
    "N=32;\n",
    "x = eda_cvec( np.linspace(0,1,N) );\n",
    "Dx = x[1,0]-x[0,0];\n",
    "\n",
    "M=N;  # square G\n",
    "\n",
    "G = np.zeros((N,M)); # create a G containing zeroes\n",
    "# Note: In order to set a column of a matrix\n",
    "# to a Nx1 column vector, the view of the matrix\n",
    "# must have a shape of Nx1, e.g. G[0:N,0:1]\n",
    "G[0:N,0:1] = np.ones((N,1)); # handle 1st column individually\n",
    "for i in range(1,M): # loop over remaining columns\n",
    "    G[0:N,i:i+1] = np.power(x,i); # x raised to (i-1) power\n",
    "\n",
    "eda_draw('title G',G);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 8188 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_04: data kernel for Fourier sum linear model\n",
    "\n",
    "# make the auxillary variable x on varange (0,1)\n",
    "N=32;\n",
    "x = eda_cvec( np.linspace(0,1,N) );\n",
    "Dx = x[1,0]-x[0,0];\n",
    "xmax = x[N-1,0];\n",
    "\n",
    "M=N-2;  # G is not square, since lambda=infinity case is omitted\n",
    "\n",
    "# wavelengths lambda\n",
    "Mo2 = floor(M/2);\n",
    "lam= eda_cvec( 2*xmax * np.reciprocal( np.linspace(1,Mo2,Mo2) ) );\n",
    "\n",
    "G = np.zeros((N,M));   # create matrix of zeros\n",
    "for i in range(0,Mo2): # loop over (cos,sin) pairs\n",
    "    icos = 2*i;   # column of the cos in th i-th pair\n",
    "    isin = 2*i+1; # column of the sin in th i-th pair\n",
    "    G[0:N,icos:icos+1] = np.cos(2*pi*x/lam[i,0]);\n",
    "    G[0:N,isin:isin+1] = np.sin(2*pi*x/lam[i,0]);\n",
    "\n",
    "eda_draw('title G',G);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 9804 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_05: data kernel viewed as concatenation of columns\n",
    "\n",
    "# make the auxillary variable x on varange (0,1)\n",
    "N=32;\n",
    "x = eda_cvec( np.linspace(0,1,N) );\n",
    "Dx = x[1,0]-x[0,0];\n",
    "xmax = x[N-1,0];\n",
    "\n",
    "M=N-2;  # G no suqare, since lambda=infinity case is omitted\n",
    "Mo2 = floor(M/2);\n",
    "\n",
    "# wavelengths lambda\n",
    "lam = eda_cvec( 2*xmax * np.reciprocal( np.linspace(1,Mo2,Mo2) ) );\n",
    "\n",
    "G = np.zeros((N,M));\n",
    "# Note: In order to set a column of a matrix\n",
    "# to a Nx1 column vector, the view of the matrix\n",
    "# must have a shape of Nx1, e.g. G[0:N,0:1]\n",
    "for i in range(1,Mo2): # loop over (cos,sin) pairs\n",
    "    icos = 2*i;   # column of the cos in th i-th pair\n",
    "    isin = 2*i+1; # column of the sin in th i-th pair\n",
    "    G[0:N,icos:icos+1] = np.cos(2*pi*x/lam[i,0]);\n",
    "    G[0:N,isin:isin+1] = np.sin(2*pi*x/lam[i,0]);\n",
    "    \n",
    "g1 = eda_cvec( G[:,0] );\n",
    "g2 = eda_cvec( G[:,1] );\n",
    "g3 = eda_cvec( G[:,2] );\n",
    "g4 = eda_cvec( G[:,3] );\n",
    "gM = eda_cvec( G[:,M-1] );\n",
    "\n",
    "eda_draw('title G',G,'=','title g1',g1,'title g2',g2,'title g3',g3,'title g4',g4,'. . .','title gM',gM);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAADnCAYAAABWihTLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAKKUlEQVR4nO3dwatcdxUH8HPUgrYURFF5oqU7ESSUIgVRIgUxKMmi7ly5C2rBFvIXuBSSTTdKwYXgwoXiIhWLXVShQm0Vk1it4CYgdFFEixYXUvm5yFTa6TR9M+/emfed+Xzg8WZe7p13Fz8CX8475/QYowAAAMjzrl0/AAAAAJsR6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQSqADAAAIJdABAACEEugAAABCCXQAAAChBDoAAIBQAh0AAEAogQ4AACCUQAcAABBKoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0AAAAoQQ6AACAUAIdAABAKIEOAAAglEAHAAAQSqADAAAIJdABAACEEugAAABCvWeqD+rusfyzS5cu3fY9vNHR0VEf99pV5+3ChQtven/+/Pm33Ld8DYfrpOftzJkzt32/6TXsp3XOGwCsQ4UOAAAglEAHAAAQSqADAAAIJdABAACEmmwoyqqBJ1euXNnoPngnq4abXL16dZLPgWWrBpfcuHFja78LAODtqNABAACEEugAAABCCXQAAAChZu2hW6anjqmsWhq+bJOeuip9dbzVcfra5uqpO+7vBwAOkwodAABAKIEOAAAglEAHAAAQSqADAAAINdlQlFUMSmEumw4usXycTWw6lMTycQBgbip0AAAAoQQ6AACAUAIdAABAqFl76JbpqWNOm/S6WT7OcUzVs2b5OAAwNRU6AACAUAIdAABAKIEOAAAglEAHAAAQaqtDUVYxKIW5WD7OVOYcOGL5OABwEip0AAAAoQQ6AACAUAIdAABAqJ330C3TU8ecLB9nKnP1qFk+DgCsQ4UOAAAglEAHAAAQSqADAAAIJdABAACEOnVDUVYxKIW5WD7OVCwfBwB2QYUOAAAglEAHAAAQSqADAAAIFdFDt0xPHXOyfJypWD4OAMxNhQ4AACCUQAcAABBKoAMAAAgl0AEAAISKHIqyikEpzMXycaZi+TgAMDUVOgAAgFACHQAAQCiBDgAAINTe9NAt01PHnCwfZyqWjwMAJ6FCBwAAEEqgAwAACCXQAQAAhBLoAAAAQu3tUJRVDEphLpaPMxXLxwGAdajQAQAAhBLoAAAAQgl0AAAAoQ6qh26ZnjrmZPk4U7F8HAB4Oyp0AAAAoQQ6AACAUAIdAABAKIEOAAAg1EEPRVnFoBTmYvk4U0lcPl5lUAoAzEGFDgAAIJRABwAAEEqgAwAACKWH7h3oqWNOlo8zFcvHAeAwqdABAACEEugAAABCCXQAAAChBDoAAIBQhqJswKAU5rLN5eMn+X2cfpaPA8BhUKEDAAAIJdABAACEEugAAABC6aGbwFQ9dcf9LA6L5eNMxfJxANg/KnQAAAChBDoAAIBQAh0AAEAogQ4AACCUoSgzsXycuVg+zlQsHweAfCp0AAAAoQQ6AACAUAIdAABAKIEOAAAglKEoW3Lc4SYGpbCJbQ5KMSRlvyUOSjEkBYBDpkIHAAAQSqADAAAIJdABAACE0kO3Q5aPMxfLx5lKYk9dlb46AA6HCh0AAEAogQ4AACCUQAcAABBKoAMAAAhlKMopY1AKc7F8nKkkDkoxJAWAfaVCBwAAEEqgAwAACCXQAQAAhNJDd8rpqWNOm/S6WT7OssSeuip9dQDsBxU6AACAUAIdAABAKIEOAAAglEAHAAAQylCUQAalMBfLx5lK4qAUQ1IASKRCBwAAEEqgAwAACCXQAQAAhNJDtwf01DEny8eZylw9apaPA3DIVOgAAABCCXQAAAChBDoAAIBQAh0AAEAoQ1H2lEEpzMXycaZi+TgAnJwKHQAAQCiBDgAAIJRABwAAEEoP3YHQU8ecLB9nKpaPA8B6VOgAAABCCXQAAAChBDoAAIBQAh0AAEAoQ1EOmEEpzMXycaZi+TgA3J4KHQAAQCiBDgAAIJRABwAAEEoPHf+np445WT7OVBKXjx8dHc322QAcNhU6AACAUAIdAABAKIEOAAAglEAHAAAQylAUbsugFOZi+ThTSVg+fu7cuUk+BwCWqdABAACEEugAAABCCXQAAACh9NCxFj11zGmby8cvXry40X1kSFw+DgCbUKEDAAAIJdABAACEEugAAABCCXQAAAChDEXhxKYalHL58uUpHoc9MufycUNRDkvC8nEA2IQKHQAAQCiBDgAAIJRABwAAEEoPHZPbtKdODx3Hsc3l4+w3y8cB2AcqdAAAAKEEOgAAgFACHQAAQCiBDgAAIFSPMXb9DAAAAGxAhQ4AACCUQAcAABBKoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAIJRABwAAEEqgAwAACCXQAQAAhBLoAAAAQgl0E+ju93b3c919vbv/2N3fXvP+R7v7zrmej/3T3Te7+w/dfa27f7vmvc4ba+nuR7r7hcX/b4+uea/zBgAz6jHGrp8hXnd3Vd01xni1u++oqmeq6pExxrPHvP9mVX16jPG3GR+TPXKSM+O8sY7u/lRV/aiqHqiq/1TVk1X1jTHGX455/81y3gBgNip0Exi3vLp4e8fi601Jubvv7e4/d/cPuvtGd/+4u+/s7m9V1Uer6unufnrLj86ect6Y0Cer6tkxxr/HGK9V1a+q6qE3XuC8AcDuCHQT6e53d/e1qnq5qp4aY/xmxWWfqKrHxxhnquqfVfXNMcZjVfVSVT04xnhwe09MuFFVv+ju33X3xbe5xnljCi9U1dnu/uDiTye/XFUfX3Gd8wYAOyDQTWSM8d8xxn1V9bGqemDxZ0rL/jrG+PXi9Q+r6nNbe0D2zWfHGPdX1Zeq6uHuPrviGueNExtjvFhV36mqp+rWn1ter6rXVlzqvAHADgh0ExtjvFJVv6yqhxYDK65199df/+fly7f6cOyNMcZLi+8vV9VPq+rzzhtzGWN8f4xx/xjjbFX9vapecd4A4HQQ6CbQ3R/q7vcvXr+vqr5QVb8fY9y3+Pre4tJ7uvszi9dfrVvDU6qq/lVVd2/1oYnV3Xd1992vv66qL1bV884bc+nuDy++31NVX6mqx5w3ADgdBLppHNWtpv8bVfV83eqhe2LFdS9W1dcW132gqr67+PnjVfVzQwM4po9U1TPdfb2qnquqn40xnlxxnfPGVH7S3X+qqqtV9fAY4x8rrnHeAGAHrC3Yku6+t6qeGGOs6q2DSTlvbJPzBgC7o0IHAAAQSoUOAAAglAodAABAKIEOAAAglEAHAAAQSqADAAAIJdABAACEEugAAABCCXQAAAChBDoAAIBQAh0AAEAogQ4AACCUQAcAABBKoAMAAAgl0AEAAIQS6AAAAEIJdAAAAKEEOgAAgFACHQAAQCiBDgAAINT/AAourln0eqXJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_06: data kernel for the weighted average linear model\n",
    "# note: this implementation requires an odd number of weights\n",
    "# but does not require them to be symmetric\n",
    "\n",
    "N=32;\n",
    "M=N; # G is square\n",
    "\n",
    "# case 1: 3-point average\n",
    "w  = eda_cvec( 1, 2, 1 ); # unnormalized weights\n",
    "La, i = np.shape(w);      # length of average\n",
    "Lw = floor(La/2);         # position of central weight\n",
    "n = np.sum(w);            # sum of weights\n",
    "w = w/n;                  # normalized weights\n",
    "wf = np.flipud(w);        # weights in reversed order\n",
    "c = np.zeros((N,1));      # initialize left column of G\n",
    "r = np.zeros((M,1));      # initialize top row of G\n",
    "c[0:Lw+1,0:1]=wf[Lw:La,0:1]; # copy weights to c\n",
    "r[0:Lw+1,0:1]=w[Lw:La,0:1];  # copy weights to r\n",
    "G = la.toeplitz(c,r.T);      # data kernel G is toeplitz\n",
    "\n",
    "G3 = np.copy(G);  # copy for display purposes\n",
    "\n",
    "# case 2: 5-point average\n",
    "w  = eda_cvec( 1, 2, 3, 2, 1 );   # unnormalized weights\n",
    "La, i = np.shape(w);       # length of average\n",
    "Lw = floor(La/2);          # position of central weight\n",
    "n = np.sum(w);             # sum of weights\n",
    "w = w/n;                   # normalized weights\n",
    "wf = np.flipud(w);         # weights in reversed order\n",
    "c = np.zeros((N,1));       # initialize left column of G\n",
    "r = np.zeros((M,1));       # initialize top row of G\n",
    "c[0:Lw+1,0:1]=wf[Lw:La,0:1]; # copy weights to c\n",
    "r[0:Lw+1,0:1]=w[Lw:La,0:1];  # copy weights to r\n",
    "G = la.toeplitz(c,r.T);      # data kernel G is toeplitz\n",
    "\n",
    "G5 = np.copy(G);  # copy for display purposes\n",
    "\n",
    "# case 3: 9-point average\n",
    "w  = eda_cvec( 1, 2, 3, 4, 5, 4, 3, 2, 1 );   # unnormalized weights\n",
    "La, i = np.shape(w);        # length of average\n",
    "Lw = floor(La/2);           # length up to central weight\n",
    "n = np.sum(w);              # sum of weights\n",
    "w = w/n;                    # normalized weights\n",
    "wf = np.flipud(w);          # weights in reversed order\n",
    "c = np.zeros((N,1));        # initialize left column of G\n",
    "r = np.zeros((M,1));        # initialize left column of G\n",
    "c[0:Lw+1,0:1]=wf[Lw:La,0:1]; # only the first Lw element of the col are non-zero\n",
    "r[0:Lw+1,0:1]=w[Lw:La,0:1];  # only the first Lw element of the row are non-zero\n",
    "G = la.toeplitz(c,r.T);      # data kernel G is toeplitz\n",
    "\n",
    "G9 = np.copy(G);  # copy for display purposes\n",
    "\n",
    "eda_draw('title 3-pt',G3,'title 5-pt',G5,'title 9-pt',G9);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 24008 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_07, prediction error\n",
    "# reads in noisy, linear data from file\n",
    "# and computes and plots prediction error for a particular\n",
    "# choice of intercept and slope, m1 and m2\n",
    "\n",
    "# model paramerters are intercept and slope\n",
    "M=2;\n",
    "mest = eda_cvec( [1.1, 2.3] );  # model estimate\n",
    "\n",
    "# load data from file\n",
    "T = np.genfromtxt('../data/linedata01.txt', delimiter='\\t')\n",
    "[N, K]=T.shape;\n",
    "x = eda_cvec( T[:,0] );    # auxiliary varible\n",
    "dobs = eda_cvec( T[:,1] ); # observed data\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[0:N,0:1] = np.ones((N,1));\n",
    "G[0:N,1:2] = x;\n",
    "\n",
    "dpre = np.matmul(G,mest);  # predicted data\n",
    "e = dobs-dpre;             # individual errors\n",
    "E = np.matmul(e.T,e);      # total error\n",
    "\n",
    "#  plot data\n",
    "fig1 = plt.figure(1,figsize=(8,8));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.plot(x,dobs,'ko');\n",
    "plt.plot(x,dpre,'k-')\n",
    "for i in range(N):\n",
    "    plt.plot( [x[i,0], x[i,0]], [dobs[i,0], dpre[i,0]],'k:');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('d');\n",
    "titlestr = \"intercept %.2f slope %.2f\" % (mest[0,0], mest[1,0]);\n",
    "plt.title(titlestr);\n",
    "\n",
    "#  plot errors\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([-6.0, 6.0, -5, 5])\n",
    "plt.plot([-6.0,6.0],[0.0,0.0],'k-');\n",
    "plt.plot(x,e,'ko');\n",
    "for i in range(N):\n",
    "    plt.plot( [x[i,0], x[i,0]], [0.0, e[i,0]],'k:');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('e');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 10420 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Emin: 23.50\n",
      "m1true:  1.00   m2true:   2.00\n",
      "m1est:   0.96   m2est:    2.00\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 22524 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11\n"
     ]
    }
   ],
   "source": [
    "# eda04_08: grid search for intercept and slope\n",
    "# of trend line fitted to noisy data\n",
    "\n",
    "# model paramerters\n",
    "M=2;\n",
    "m = np.zeros((M,1));\n",
    "\n",
    "# load data from file\n",
    "T = np.genfromtxt('../data/linedata01.txt', delimiter='\\t')\n",
    "[N, K]=T.shape;\n",
    "x = eda_cvec( T[:,0] );    # auxially variable\n",
    "dobs = eda_cvec( T[:,1] ); # observed data\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[0:N,0:1] = np.ones((N,1));\n",
    "G[0:N,1:2] = x;\n",
    "\n",
    "# define the (m1, m2) grid\n",
    "L1 = 101;\n",
    "L2 = 101;\n",
    "m1min=0;\n",
    "m1max=4;\n",
    "m2min=0;\n",
    "m2max=4;\n",
    "m1 = eda_cvec( np.linspace(m1min,m1max,L1) );\n",
    "m2 = eda_cvec( np.linspace(m2min,m2max,L2) );\n",
    "\n",
    "# evaluate the error at each grid point\n",
    "E = np.zeros((L1,L2));\n",
    "for i in range(L1):\n",
    "    for j in range(L2):\n",
    "        m = eda_cvec([m1[i,0],m2[j,0]]); # put m's into vector form\n",
    "        dpre = np.matmul(G,m);     # predicted data\n",
    "        e = dobs - dpre;           # individual errors\n",
    "        E[i,j] = np.matmul(e.T,e); # total error\n",
    "\n",
    "# best estimate is the one with smallest error\n",
    "# find the (i,j) for which E(i,j) is minimum\n",
    "k = np.argmin( E, axis=0 );\n",
    "Etemp = np.min( E, axis=0 );\n",
    "j = np.argmin(Etemp);\n",
    "i=k[j];\n",
    "Emin = E[i,j];  # minimum error\n",
    "# estmated model parameters\n",
    "mbest = eda_cvec( [ m1[i,0],  m2[j,0] ] );\n",
    "dpre = np.matmul(G,mbest); # predicted data\n",
    "e = dobs-dpre; # individual errors\n",
    "\n",
    "eda_draw('title E(m1,m2)',E);\n",
    "\n",
    "print(\"Emin: %.2f\" % (Emin) );\n",
    "print(\"m1true:  %.2f   m2true:   %.2f\" % (1, 2) );\n",
    "print(\"m1est:   %.2f   m2est:    %.2f\" % (mbest[0,0], mbest[1,0]) );\n",
    "\n",
    "#  plot data\n",
    "fig2 = plt.figure(2,figsize=(8,8));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([-6.0, 6.0, -15.0, 15.0]);\n",
    "plt.plot(x,dobs,'ko');\n",
    "plt.plot(x,dpre,'k-');\n",
    "plt.xlabel('x)');\n",
    "plt.ylabel('d');\n",
    "titlestr = \"intercept %.2f slope %.2f\" % (mbest[0,0], mbest[1,0]);\n",
    "plt.title(titlestr);\n",
    "\n",
    "#  plot errors\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([-6.0, 6.0, -5, 5])\n",
    "plt.plot([-6.0,6.0],[0.0,0.0],'k-');\n",
    "plt.plot(x,e,'ko');\n",
    "for i in range(N):\n",
    "    plt.plot( [x[i,0], x[i,0]], [0.0, e[i,0]],'k:');\n",
    "plt.xlabel('x)');\n",
    "plt.ylabel('d');\n",
    "plt.show();\n",
    "\n",
    "print(N);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "left  shift: intercept 1.00  slope 2.00\n",
      "no    shift: intercept 1.00  slope 2.00\n",
      "right shift: intercept 1.00  slope 2.00\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 27844 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_09: grid search for intercept and slope\n",
    "# for problem of fittng straight line to noisy data\n",
    "\n",
    "# model paramerters\n",
    "M=2\n",
    "m = np.zeros((M,1));\n",
    "\n",
    "# load data from file\n",
    "T = np.genfromtxt('../data/linedata01.txt', delimiter='\\t')\n",
    "[N, K]=T.shape;\n",
    "x = eda_cvec( T[:,0] );    # auxiallary variable\n",
    "dobs = eda_cvec( T[:,1] ); # observed data\n",
    "\n",
    "# Part 1: shift left\n",
    "xshift=-2;\n",
    "xs = x+xshift;\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = xs.ravel();\n",
    "dpre = np.matmul(G,m);\n",
    "\n",
    "# define grid of intercepts and slopes\n",
    "L1 = 101;\n",
    "L2 = 101;\n",
    "m1min=-10;\n",
    "m1max=10;\n",
    "m2min=-10;\n",
    "m2max=10;\n",
    "m1 = eda_cvec( np.linspace(m1min,m1max,L1) );\n",
    "m2 = eda_cvec( np.linspace(m2min,m2max,L2) );\n",
    "e = np.zeros((N,1));\n",
    "E = np.zeros((L1,L2));\n",
    "\n",
    "for i in range(L1):\n",
    "    for j in range(L2):\n",
    "        m[:,0] = [m1[i], m2[j]];\n",
    "        dpre = np.matmul(G,m);\n",
    "        e = dobs - dpre;\n",
    "        E0 = np.matmul(e.T,e);\n",
    "        E[i,j] = E0;\n",
    "    \n",
    "# best estimate is the one with smallest error\n",
    "k = np.argmin( E, axis=0 );\n",
    "Etemp = np.min( E, axis=0 );\n",
    "j = np.argmin(Etemp);\n",
    "i=k[j];\n",
    "Ebest = E[i,j];\n",
    "mbest = eda_cvec( [ m1[i,0],  m2[j,0] ] );\n",
    "dpre = np.matmul(G,mbest);\n",
    "Esl = E;\n",
    "print(\"left  shift: intercept %.2f  slope %.2f\" % (mbest[0,0]-2*mbest[1,0], mbest[1,0]));\n",
    "\n",
    "\n",
    "# Part 2: no shift \n",
    "xshift=0;\n",
    "xs = x+xshift;\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = xs.ravel();\n",
    "dpre = np.matmul(G,m);\n",
    "\n",
    "# define grid of intercepts and slopes\n",
    "L1 = 101;\n",
    "L2 = 101;\n",
    "m1min=-10;\n",
    "m1max=10;\n",
    "m2min=-10;\n",
    "m2max=10;\n",
    "m1 = eda_cvec( np.linspace(m1min,m1max,L1) );\n",
    "m2 = eda_cvec( np.linspace(m2min,m2max,L2) );\n",
    "e = np.zeros((N,1));\n",
    "E = np.zeros((L1,L2));\n",
    "\n",
    "for i in range(L1):\n",
    "    for j in range(L2):\n",
    "        m = eda_cvec( [m1[i,0], m2[j,0]] );\n",
    "        dpre = np.matmul(G,m);\n",
    "        e = dobs-dpre;\n",
    "        E0 = np.matmul(e.T,e);\n",
    "        E[i,j] = E0;\n",
    "    \n",
    "# best estimate is the one with smallest error\n",
    "k = np.argmin( E, axis=0 );\n",
    "Etemp = np.min( E, axis=0 );\n",
    "j = np.argmin(Etemp);\n",
    "i=k[j];\n",
    "Ebest = E[i,j];\n",
    "mbest = eda_cvec( [ m1[i,0],  m2[j,0] ] );\n",
    "dpre = np.matmul(G,mbest);\n",
    "Ens = E;\n",
    "print(\"no    shift: intercept %.2f  slope %.2f\" % (mbest[0,0], mbest[1,0]));\n",
    "\n",
    "# Part 3: shift right\n",
    "xshift=2;\n",
    "xs = x+xshift;\n",
    "\n",
    "# linear data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "G[:,1] = xs.ravel();\n",
    "dpre = np.matmul(G,m);\n",
    "\n",
    "# define grid of intercepts and slopes\n",
    "L1 = 101;\n",
    "L2 = 101;\n",
    "m1min=-10;\n",
    "m1max=10;\n",
    "m2min=-10;\n",
    "m2max=10;\n",
    "m1 = eda_cvec( np.linspace(m1min,m1max,L1) );\n",
    "m2 = eda_cvec( np.linspace(m2min,m2max,L2) );\n",
    "e = np.zeros((N,1));\n",
    "E = np.zeros((L1,L2));\n",
    "\n",
    "for i in range(L1):\n",
    "    for j in range(L2):\n",
    "        m[:,0] = [m1[i], m2[j]];\n",
    "        dpre = np.matmul(G,m);\n",
    "        e = dobs-dpre;\n",
    "        E0 = np.matmul(e.T,e);\n",
    "        E[i,j] = E0;\n",
    "    \n",
    "# best estimate is the one with smallest error\n",
    "k = np.argmin( E, axis=0 );\n",
    "Etemp = np.min( E, axis=0 );\n",
    "j = np.argmin(Etemp);\n",
    "i=k[j];\n",
    "Ebest = E[i,j];\n",
    "mbest = eda_cvec( [ m1[i,0],  m2[j,0] ] );\n",
    "dpre = np.matmul(G,mbest);\n",
    "Esr = E;\n",
    "print(\"right shift: intercept %.2f  slope %.2f\" % (mbest[0,0]+2*mbest[1,0], mbest[1,0]));\n",
    "\n",
    "# lake log of error to reduce range of variation for plotting purposes\n",
    "eda_draw('title E shift left',np.log(Esl),'title E no shift',np.log(Ens),'title E shift right',np.log(Esr));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total error: 23.46\n",
      "intercept 0.94 +/- 0.97 (%95)\n",
      "slope     2.02 +/- 0.31 (%95)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 18684 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_10: least-squares solution\n",
    "# for straight line linear model\n",
    "\n",
    "# reads in noisy, linear data from file\n",
    "# and computes and plots predicted data\n",
    "# using least-squares estimate of intercept\n",
    "# and slope, m1 and m2\n",
    "\n",
    "# model paramerters\n",
    "M=2;\n",
    "mest = np.zeros((M,1));\n",
    "\n",
    "# load data from file\n",
    "T = np.genfromtxt('../data/linedata01.txt', delimiter='\\t')\n",
    "[N, K]=T.shape;\n",
    "x = eda_cvec( T[:,0] );\n",
    "dobs = eda_cvec( T[:,1] );\n",
    "\n",
    "# data kernel\n",
    "G = np.zeros((N,2));\n",
    "G[0:N,0:1] = np.ones((N,1));\n",
    "G[0:N,1:2] = x;\n",
    "\n",
    "# least squares solution\n",
    "GTG = np.matmul(G.T,G);\n",
    "mest = la.solve( GTG, np.matmul(G.T, dobs)  );\n",
    "\n",
    "# eda04_10: variance calculation\n",
    "dpre = np.matmul(G, mest); # predicted data\n",
    "e = dobs - dpre;           # error\n",
    "E = np.matmul(e.T,e);      # total error\n",
    "sigmad2 = E/(N-M);         # posterior variance of data\n",
    "Cm = sigmad2*la.inv(GTG);  # variance of model parameters\n",
    "sigmam1 = sqrt( Cm[0,0] ); # sqrt variance of model parameter 1\n",
    "sigmam2 = sqrt( Cm[1,1] ); # sqrt variance of model parameter 2\n",
    "\n",
    "# print model parameters and theoir 95% confidence intervals\n",
    "print(\"Total error: %.2f\" % E);\n",
    "str1 = \"intercept %.2f +/- %.2f (%%95)\" % (mest[0,0],2*sigmam1);\n",
    "str2 = \"slope     %.2f +/- %.2f (%%95)\" % (mest[1,0],2*sigmam2);\n",
    "print(str1);\n",
    "print(str2);\n",
    "\n",
    "#  plot data\n",
    "fig2 = plt.figure(2,figsize=(8,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([-6, 6, -15, 15])\n",
    "plt.plot(x,dobs,'ko');\n",
    "plt.plot(x,dpre,'k-');\n",
    "for i in range(N):\n",
    "    plt.plot([x[i,0],x[i,0]], [dobs[i,0],dpre[i,0]],'k-');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('d');\n",
    "titlestr = \"intercept %.2f slope %.2f\" % (mest[0,0], mest[1,0]);\n",
    "plt.title(titlestr);\n",
    "\n",
    "#  plot error, with different vertical scale\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis([-6, 6, -5, 5]);\n",
    "plt.plot(x,e,'ko');\n",
    "plt.plot([-6,6],[0,0],'k-');\n",
    "for i in range(N):\n",
    "    plt.plot([x[i,0],x[i,0]], [0,e[i,0]],'k:');\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('e');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "using prior     error: slope -3.137505e-02 +/- 1.626961e-05 deg C per year (95%)\n",
      "using posterior error: slope -3.137505e-02 +/- 9.144082e-03 deg C per year (95%)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 73304 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda04_11: trend of Black Rock Forest temperature data\n",
    "\n",
    "# load data from file\n",
    "D = np.genfromtxt('../data/brf_temp.txt', delimiter='\\t')\n",
    "[Nraw, K]=np.shape(D);\n",
    "traw = eda_cvec( D[:,0] );  # time (auxiliary variable)\n",
    "draw = eda_cvec( D[:,1] );  # temperature (data)\n",
    "\n",
    "# select only good data\n",
    "# r=where(); k=r[0]; computes an Nx0 array of indices that pass the test\n",
    "r = np.where( (draw[:,0]!=0)&(draw[:,0]>-40)&(draw[:,0]<38) );\n",
    "k = r[0];\n",
    "# select these indices from the raw data\n",
    "t = traw[k,0:1];\n",
    "dobs = draw[k,0:1];\n",
    "\n",
    "(N,temp)=np.shape(dobs); # length of good data\n",
    "\n",
    "# set up the data kernel\n",
    "Ty = 365.25;  # number of days in a year\n",
    "M=4;          # numnber of model parameters\n",
    "G = np.zeros((N,M)); # initialize data kernel to zero\n",
    "G[0:N,0:1] = np.ones((N,1)); # constant in time, m0\n",
    "G[0:N,1:2] = t;              # linear in time, m1\n",
    "# cosine in time, period of 1 year, m2\n",
    "G[0:N,2:3] = np.cos(2*pi*t/Ty);\n",
    "# sine in time, period of 1 year, m3\n",
    "G[0:N,3:4] = np.sin(2*pi*t/Ty);\n",
    "\n",
    "# simple least squares solution\n",
    "GTG = np.matmul(G.T,G);\n",
    "mest = la.solve( GTG, np.matmul(G.T, dobs)  );\n",
    "slope = Ty*mest[1,0]; # slope in deg C per year\n",
    "\n",
    "# variance calculation\n",
    "dpre = np.matmul(G, mest); # predicted data\n",
    "e = dobs - dpre;           # error\n",
    "E = np.matmul(e.T,e);      # total error\n",
    "sigmad2A = 0.01**2;        # prior error, 0.1 deg C\n",
    "sigmad2B = E/(N-M);        # posterior error , from fit\n",
    "GTGI = la.inv(GTG);        # inverse of GTG\n",
    "CmA = sigmad2A*GTGI;       # covariance of m using prior variance of data\n",
    "CmB = sigmad2B*GTGI;       # covariance of m using posterior variance of data\n",
    "errslopeA = Ty*sqrt(CmA[1,1]); # error in slope (prior case), with units of deg C per year\n",
    "errslopeB = Ty*sqrt(CmB[1,1]); # error in slope (prior case), with units of deg C per year\n",
    "\n",
    "# convert to per year\n",
    "slopestr = \"using prior     error: slope %e +/- %e deg C per year (95%%)\" % (slope, 2*errslopeA);\n",
    "print(slopestr);\n",
    "slopestr = \"using posterior error: slope %e +/- %e deg C per year (95%%)\" % (slope, 2*errslopeB);\n",
    "print(slopestr);\n",
    "\n",
    "#  plot data\n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([t[0,0], t[N-1,0], -40, 40])\n",
    "plt.plot(t,dobs,'k.');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('temperature (C)');\n",
    "ax2 = plt.subplot(3,1,2);\n",
    "plt.axis([t[0,0], t[N-1,0], -40, 40])\n",
    "plt.plot(t,dpre,'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('temperature (C)');\n",
    "ax3 = plt.subplot(3,1,3);\n",
    "plt.axis([t[0,0], t[N-1,0], -40, 40])\n",
    "plt.plot(t,e,'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('error (C)');\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
}
