{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eda13_00 clear all variables and import vatious modules\n",
    "%reset -f\n",
    "import os\n",
    "from datetime import date, datetime\n",
    "from math import exp, pi, sin, cos, sqrt, floor, ceil, atan2, log\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",
    "from matplotlib import patches as pch\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import re\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the variable pi: 3.1416 (before)\n",
      "the variable pi: 2.0000 (after)\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'int' object has no attribute 'zeros'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-6-3fd9ab34cd5b>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      8\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      9\u001b[0m \u001b[0mnp\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m3\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 10\u001b[1;33m \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m;\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     11\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     12\u001b[0m \u001b[1;31m# at this point you must rerun the first cell,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAttributeError\u001b[0m: 'int' object has no attribute 'zeros'"
     ]
    }
   ],
   "source": [
    "# edapy_13_01: exemplary errors related to redefining variables\n",
    "pisave=pi;\n",
    "npsave=np;\n",
    "\n",
    "print('the variable pi: %.4f (before)' % pi)\n",
    "pi=2;\n",
    "print('the variable pi: %.4f (after)' % pi)\n",
    "\n",
    "np = 3;\n",
    "x = np.zeros((3,1));\n",
    "\n",
    "# at this point you must rerun the first cell,\n",
    "# to restore numpy\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seconds DT between 2008/2/11 4:4:1 and 2008/2/11 4:4:0\n",
      "DT: 1.00\n",
      "\n",
      "leap second checker\n",
      "seconds DT between 1975/1/1 0:0:1 and 1974/12/31 23:59:59\n",
      "right answer is 3, calculation gives DT: 2.00\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# edapy_13_02: time arithmetic\n",
    "\n",
    "# restore numpy, which was destroyed in previous cell\n",
    "pi=pisave;\n",
    "np=npsave;\n",
    "\n",
    "print('seconds DT between 2008/2/11 4:4:1 and 2008/2/11 4:4:0');\n",
    "DELTA = datetime(2008,2,11,4,4,1)-datetime(2008,2,11,4,4,0);\n",
    "print('DT: %.2f\\n' % DELTA.total_seconds() );\n",
    "\n",
    "DELTA = datetime(1975,1,1,0,0,1)-datetime(1974,12,31,23,59,59);\n",
    "print('leap second checker');\n",
    "print('seconds DT between 1975/1/1 0:0:1 and 1974/12/31 23:59:59');\n",
    "print('right answer is 3, calculation gives DT: %.2f\\n' % DELTA.total_seconds() );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35784 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_13_03: reading idiosynchratic text file\n",
    "\n",
    "# since size of file unknaon, allocate very big arrays\n",
    "NBIG=1000000;\n",
    "t = np.zeros((NBIG,1));\n",
    "d = np.zeros((NBIG,1));\n",
    "\n",
    "# names of months\n",
    "moname = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'];\n",
    "\n",
    "# open file\n",
    "filename='../Data/brf_raw.txt';\n",
    "fid = open(filename,\"r\");\n",
    "N=0; # line number counter\n",
    "while 1: # loop over lines in file\n",
    "    myline = fid.readline(); # read line\n",
    "    if not myline:\n",
    "        fid.close();\n",
    "        break;\n",
    "    # now process the line\n",
    "    myline = myline.replace(\"\\t\",\" \"); # tab to space\n",
    "    Ns = len(myline);\n",
    "    while 1: # replace multiple spaces with single space\n",
    "        myline = myline.replace(\"  \",\" \");\n",
    "        Nsnew = len(myline);\n",
    "        if( Nsnew==Ns ):\n",
    "            break;\n",
    "        Ns = Nsnew;\n",
    "    s = myline.split(); # split into words, result is list\n",
    "    s2 = s[0].split(\"-\"); # spit first word into two\n",
    "    hrmn1 = s2[0];  # first hour-minute string\n",
    "    hrmn2 = s2[1];  # second hour-minute string\n",
    "    hr1 = int(hrmn1[0:2]); # first hour\n",
    "    hr2 = int(hrmn2[0:2]); # secong hour\n",
    "    mn1 = int(hrmn1[2:4]); # first minute\n",
    "    mn2 = int(hrmn2[2:4]); # second minute\n",
    "    sc1 = 0.0;\n",
    "    sc2 = 0.0;\n",
    "    # some end times listed as 24:00:00, change to 23:59:59\n",
    "    if( (hr2 == 24) and (mn2==0) ):\n",
    "        hr2 = 23;\n",
    "        mn2 = 59;\n",
    "        sc2 = 59;\n",
    "    da = int(s[1]);\n",
    "    yr = int(s[3]);\n",
    "    va = float(s[4])\n",
    "    # look-up month\n",
    "    mo=(-1);\n",
    "    for imo in range(12):\n",
    "        if( s[2]==moname[imo]):\n",
    "            mo=imo+1;\n",
    "            break;\n",
    "    # start time\n",
    "    t1=datetime(yr,mo,da,hr1,mn1,0);\n",
    "    # end time\n",
    "    t2=datetime(yr,mo,da,hr2,mn2,0);\n",
    "    if( N==0 ): # reference time is year of first sample\n",
    "        tref=datetime(yr,1,1,0,0,0);\n",
    "    DT1 = (t1-tref).total_seconds();\n",
    "    DT2 = (t2-tref).total_seconds();\n",
    "    t[N,0] = 0.5*(DT1+DT2)/86400; # time since reference time\n",
    "    d[N,0] = va; # data value\n",
    "    N=N+1; # increment line counter\n",
    "\n",
    "# truncate to actual length\n",
    "t = np.copy( t[0:N,0:1] );\n",
    "d = np.copy( d[0:N,0:1] );\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "plt.axis([0, np.max(t), np.min(d), np.max(d)]);\n",
    "plt.plot( t, d, 'k-');\n",
    "plt.xlabel('time (days');\n",
    "plt.ylabel('temperature (degC)');\n",
    "plt.title('Black Rock Forest');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 5076 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fig. The equation Gm=d, depicted using eda_draw()\n"
     ]
    }
   ],
   "source": [
    "# edapy_13_04, example of use of eda_draw() function\n",
    "\n",
    "# define an examplary matrix, G\n",
    "N=40;\n",
    "v = eda_cvec(np.linspace(1,N,N));\n",
    "G=np.matmul( v, v.T );\n",
    "\n",
    "# define a exemplacy vector, m\n",
    "m = eda_cvec(1.0*np.linspace(1,N,N));\n",
    "\n",
    "# multiply\n",
    "d = np.matmul( G, m);\n",
    "\n",
    "# draw\n",
    "eda_draw( 'title d', d, '=', 'title G', G, 'title m', m);\n",
    "plt.show();\n",
    "\n",
    "# add caption\n",
    "print('Fig. The equation Gm=d, depicted using eda_draw()');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "circle:\n",
      "radius: 2.000000 area: 12.566371\n",
      " \n",
      "rectangle;\n",
      "length: 2.000000 width: 4.000000 circumference: 12.000000 area: 8.000000\n",
      " \n"
     ]
    }
   ],
   "source": [
    "# edapy_13_05, simple examples of functions\n",
    "\n",
    "# computes area, a, of circle of radius, r.\n",
    "def areaofcircle(r):\n",
    "    a = pi*(r**2);\n",
    "    return a;\n",
    "\n",
    "# computes circumference, c, and area, a, of\n",
    "# a rectangle of length, l, and width, w\n",
    "def CandAofRectangle(a,b):\n",
    "    circ = 2*(a+b);\n",
    "    area = a*b;\n",
    "    return circ, area;\n",
    "\n",
    "# calculate area of a circle \n",
    "radius=2;\n",
    "area = areaofcircle(radius);\n",
    "print('circle:');\n",
    "print('radius: %.6f area: %.6f' % (radius,area) );\n",
    "print(' ');\n",
    "\n",
    "# calculate circumference and area of a rectangle\n",
    "a=2;\n",
    "b=4;\n",
    "circ, area = CandAofRectangle(a,b);\n",
    "print('rectangle;');\n",
    "print('length: %.6f width: %.6f circumference: %.6f area: %.6f' % (a,b,circ,area) );\n",
    "print(' ');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test matrix A\n",
      "[[ 1  2  3  4]\n",
      " [ 5  6  7  8]\n",
      " [ 9 10 11 12]]\n",
      " \n",
      "transpose of a column vector c formed from A\n",
      "[[ 1  5  9  2  6 10  3  7 11  4  8 12]]\n",
      " \n",
      "row vector formed from A\n",
      "[[ 1  5  9  2  6 10  3  7 11  4  8 12]]\n",
      " \n",
      "the matrix reformed from the vector c\n",
      "[[ 1  2  3  4]\n",
      " [ 5  6  7  8]\n",
      " [ 9 10 11 12]]\n",
      " \n",
      "for a (3,4) matrix M[i,j] unraveled into a vector c[k,0]\n",
      "(i,j) of (1,0) corresponds to k of 1\n",
      "A(i,j)=5 and c[k]=5\n",
      "k of 1 corresponds to (i,j) of (1,0)\n"
     ]
    }
   ],
   "source": [
    "# edapy_13_06, examples of reorganizing a matrix\n",
    "\n",
    "# define a matrix\n",
    "A = np.array( [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]] );\n",
    "N,M = np.shape(A);\n",
    "print('test matrix A');\n",
    "print(A);\n",
    "print(' ');\n",
    "\n",
    "# create a column-vector, b, containing the elements of A,\n",
    "# organized column-wise.\n",
    "c = np.reshape(A,(N*M,1),order='F');\n",
    "print('transpose of a column vector c formed from A');\n",
    "print(c.T);\n",
    "print(' ');\n",
    "\n",
    "# create a row-vector, c, containing the elements of A, organized\n",
    "# column-wise.\n",
    "r = np.reshape(A,(1,N*M),order='F');\n",
    "print('row vector formed from A');\n",
    "print(r);\n",
    "print(' ');\n",
    "\n",
    "# the matrix reformed from the vector\n",
    "B = np.reshape(c,(N,M),order='F');\n",
    "print('the matrix reformed from the vector c');\n",
    "print(B);\n",
    "print(' ');\n",
    "\n",
    "# when you 'unravel' an array with elements A[i,j] into\n",
    "# column vector with elements b[k,0], you sometimes need\n",
    "# to determine the correspondence beteen (i,j) and k.\n",
    "# and vice versa:\n",
    "i=1;\n",
    "j=0;\n",
    "print('for a (%d,%d) matrix M[i,j] unraveled into a vector c[k,0]' % (N,M) );\n",
    "k = np.ravel_multi_index( (i,j), (N,M), order='F') ;\n",
    "print('(i,j) of (%d,%d) corresponds to k of %d' % (i,j,k) );\n",
    "print('A(i,j)=%.0f and c[k]=%.0f' % (A[i,j], c[k,0]) );\n",
    "ii, jj = np.unravel_index( k, (N,M), order='F' ) ;\n",
    "print('k of %d corresponds to (i,j) of (%d,%d)' % (k,i,j) );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 9056 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_13_07, function Phi used in Lagrange multiplier figure\n",
    "\n",
    "# make an (x,y) grid\n",
    "N=51;\n",
    "M=51;\n",
    "Dx=0.01;\n",
    "Dy=0.01;\n",
    "x = eda_cvec(Dx*np.linspace(0,N-1,N));\n",
    "y = eda_cvec(Dx*np.linspace(0,M-1,M));\n",
    "\n",
    "# a long-wavelength, 2D sinusoid\n",
    "Phi= np.multiply( np.sin(x), (np.sin(y)).T );\n",
    "\n",
    "eda_draw(Phi);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 21980 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_13_08, supports Figure 13.3\n",
    "# covariance corresponding to 2nd Derivative smoothness\n",
    "\n",
    "# constants\n",
    "s = 0.1;  # scale factor\n",
    "g2 = 1.0;  # amplitude\n",
    "\n",
    "# stamdard setup of x-axis\n",
    "N2 = 100;\n",
    "N = 2*N2+1;\n",
    "Dx = 1.0;\n",
    "x = eda_cvec(Dx*np.linspace(-N2,N2,N));\n",
    "xmin = np.min(x);\n",
    "xmax = np.max(x);\n",
    "\n",
    "# autocovariance, with lag |x| = |xi - xj|\n",
    "c = g2*np.multiply((1.0+s*abs(x)),np.exp(-s*np.abs(x)));\n",
    "\n",
    "# Gussian autocovariance for comparison\n",
    "cg = g2*np.exp(-0.5*s*s*np.power(x,2));\n",
    "\n",
    "# plot c(x) vs x\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([xmin, xmax, -0.1, 1.1])\n",
    "plt.plot(x,c,'k-',lw=3);\n",
    "plt.plot(x,cg,'k:',lw=2);\n",
    "plt.plot(x,0.0*c,'k-',lw=1);\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('c(x)');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M: 64\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 38396 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Brute force method, log det Cm = -1.860150e+02\n",
      "Eigenvalue  method, log det Cm = -1.860150e+02\n",
      "Cholesky    method, log det Cm = -1.860150e+02\n",
      " \n",
      "brute force works fine for M=64, but fails for M=230 and above\n"
     ]
    }
   ],
   "source": [
    "# edapy_13_09, computing ln det Cm\n",
    "\n",
    "# model parameter axis, t, with sampling Dt\n",
    "M=64;  # number of model parameters\n",
    "Dt=1.0;  # time sampling\n",
    "tmax=Dt*(M-1);  # maximum time\n",
    "t=eda_cvec(np.linspace(0,tmax,M)); # time\n",
    "\n",
    "print('M: %d' % (M) );\n",
    "\n",
    "# covariance matrix, of cos(w0t) exp(-st) form\n",
    "n=6;  # number of oscillations on (0,tmax) interal\n",
    "w0 = n*2*pi/tmax;  # angular frequency of oscilation\n",
    "tmid=tmax/2.0;\n",
    "s = Dt/60.0;      # decay rate\n",
    "gamma = 1.0;\n",
    "gamma2 = gamma**2; # variance\n",
    "\n",
    "# build covariances via outer product\n",
    "Tm = np.abs(np.matmul(t,np.ones((1,M)))-np.matmul(np.ones((M,1)),(t.T)));\n",
    "Cm = gamma2*np.multiply(np.cos(w0*Tm),np.exp(-s*Tm));\n",
    "\n",
    "eda_draw('title Cm',Cm);\n",
    "\n",
    "# Brute force method\n",
    "D = la.det(Cm)\n",
    "logdetCm_A = log(D);\n",
    "print('Brute force method, log det Cm = %e' % (logdetCm_A) );\n",
    "\n",
    "# Eigenvalue method\n",
    "LAMBDA, V = la.eig(Cm);\n",
    "LAMBDA = np.real(LAMBDA);\n",
    "logdetCm_B = np.sum(np.log(LAMBDA));\n",
    "print('Eigenvalue  method, log det Cm = %e' % (logdetCm_B) );\n",
    "\n",
    "# Cholesky method\n",
    "R = la.cholesky(Cm);\n",
    "logdetCm_C = 2.0*np.sum(np.log(np.diag(R)));\n",
    "print('Cholesky    method, log det Cm = %e' % (logdetCm_C) );\n",
    "print(' ');\n",
    "print(\"brute force works fine for M=%d, but fails for M=230 and above\" %(M));\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
