{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# edapy_01_00 clears all variables and import various modules\n",
    "# you must run this cell first, before running any of others\n",
    "\n",
    "# History:\n",
    "# 2021/08/25 -  Edited by W.Menke, synchronizing Python/MATLAB versions\n",
    "# 2021/06/19 -  Edited by W. Menke, while writing corresponding chapter\n",
    "# 2020/07/21 -  Created by W.Menke, by translation of Matlab version\n",
    "# 2021/12/30 -  review, fixed some issues including np.where()\n",
    " \n",
    "# A reset deletes any previously-created variables from memory\n",
    "%reset -f\n",
    "\n",
    "# import various packages\n",
    "\n",
    "import numpy as np                    # matrices & etc\n",
    "from matplotlib import pyplot as plt  # general plotting \n",
    "from matplotlib import patches        # plotting shapes\n",
    "from datetime import date             # date and time arithmetic\n",
    "from math import exp, pi, sin, sqrt   # math functions\n",
    "import scipy.linalg as la             # linear algebra functions\n",
    "import os                             # operating system\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",
    "\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2021-12-30\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_01 print the date\n",
    "thedate = date.today();\n",
    "print(thedate);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c:\\bill\\EDAWM-3rdEd\\edapy\\Notebooks\n",
      "c:\\bill\\EDAWM-3rdEd\\edapy\\TestFolder\n",
      "c:\\bill\\EDAWM-3rdEd\\edapy\\TestFolder\n",
      "['TextDocument1.txt', 'TextDocument2.txt']\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_02: directory manipulation\n",
    "\n",
    "# get current directory\n",
    "edapy_dir=os.getcwd();\n",
    "print(edapy_dir);\n",
    "\n",
    "# get current folder, with is presumably edapy/Notebooks\n",
    "edapy_dir=os.getcwd();\n",
    "\n",
    "# change to the TestFolder, up-and-over from Notebooks\n",
    "os.chdir(\"../TestFolder\");\n",
    "\n",
    "# print current working directory (folder)\n",
    "thisdir=os.getcwd();\n",
    "print(thisdir);\n",
    "\n",
    "# change back to original folder\n",
    "os.chdir(edapy_dir);\n",
    "\n",
    "# change to the TestFolder, which is up-and-over\n",
    "os.chdir(\"../TestFolder\");\n",
    "\n",
    "# print current working directory (folder)\n",
    "thisdir=os.getcwd();\n",
    "print(thisdir);\n",
    "\n",
    "# change to folder one level up\n",
    "os.chdir(\"..\");\n",
    "\n",
    "# change to the Notebooks folder\n",
    "os.chdir(\"Notebooks\");\n",
    "\n",
    "# change to TestFolder (because it only contains two files)\n",
    "os.chdir(\"../TestFolder\");\n",
    "\n",
    "mydirlist=os.listdir();\n",
    "print(mydirlist);\n",
    "\n",
    "os.chdir(\"../Notebooks\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.6\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_03: simple algebra\n",
    "a = 3.5;\n",
    "b = 4.1;\n",
    "c = a+b;\n",
    "print(c);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.0\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_04: algebra with square root and power\n",
    "a = 3.0;\n",
    "b = 4.0;\n",
    "c = sqrt(a**2 + b**2);\n",
    "print(c);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5877852522924732\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_05: algebra with sine\n",
    "n = 2; x = 3; x0 = 1; L = 5;\n",
    "c = sin(n*pi*(x-x0)/L);\n",
    "print(c);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "List\n",
      "mylist [1, 2, 4, 8]\n",
      "mylist[0] 1\n",
      "len(mylist) 4\n",
      "Tuple\n",
      "mytuple (1, 2, 4, 8)\n",
      "mytuple[0] 1\n",
      "len(mytuple) 4\n",
      "4x0 np.array\n",
      "[0. 0. 0. 0.]\n",
      "1x4 np.array\n",
      "[[0. 0. 0. 0.]]\n",
      "4x1 np.array\n",
      "[[0.]\n",
      " [0.]\n",
      " [0.]\n",
      " [0.]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_06 list, tupples, and other list-like objects\n",
    "\n",
    "print(\"List\");\n",
    "mylist = [1, 2, 4, 8];\n",
    "N=len(mylist);\n",
    "print(\"mylist\", mylist);\n",
    "print(\"mylist[0]\", mylist[0]);\n",
    "print(\"len(mylist)\",N);\n",
    "\n",
    "print(\"Tuple\");\n",
    "mytuple = (1, 2, 4, 8);\n",
    "N=len(mytuple);\n",
    "print(\"mytuple\", mytuple);\n",
    "print(\"mytuple[0]\", mytuple[0]);\n",
    "print(\"len(mytuple)\",N);\n",
    "\n",
    "v4x0 = np.zeros((4));\n",
    "print(\"4x0 np.array\");\n",
    "print(v4x0);\n",
    "v1x4 = np.zeros((1,4));\n",
    "print(\"1x4 np.array\");\n",
    "print(v1x4);\n",
    "v4x1 = np.zeros((4,1));\n",
    "print(\"4x1 np.array\");\n",
    "print(v4x1);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c0\n",
      " [[0.]\n",
      " [0.]\n",
      " [0.]\n",
      " [0.]]\n",
      "c1\n",
      " [[1.]\n",
      " [1.]\n",
      " [1.]\n",
      " [1.]]\n",
      "c\n",
      " [[1.]\n",
      " [2.]\n",
      " [3.]\n",
      " [4.]]\n",
      "r\n",
      " [[1. 2. 3. 4.]]\n",
      "A0\n",
      " [[0. 0. 0.]\n",
      " [0. 0. 0.]\n",
      " [0. 0. 0.]]\n",
      "A1\n",
      " [[1. 1. 1.]\n",
      " [1. 1. 1.]\n",
      " [1. 1. 1.]]\n",
      "A\n",
      " [[1. 2. 3.]\n",
      " [4. 5. 6.]\n",
      " [7. 8. 9.]]\n",
      "shape of c:\n",
      " 4 1\n",
      "shape of M\n",
      " 3 3\n",
      "c\n",
      " [[1.]\n",
      " [2.]\n",
      " [3.]\n",
      " [4.]]\n",
      "c\n",
      " [[1.]\n",
      " [2.]\n",
      " [3.]\n",
      " [4.]]\n",
      "c\n",
      " [[1.]\n",
      " [2.]\n",
      " [3.]\n",
      " [4.]]\n",
      "c\n",
      " [[1.]\n",
      " [2.]\n",
      " [3.]\n",
      " [4.]]\n",
      "c\n",
      " [[1.]\n",
      " [2.]\n",
      " [3.]\n",
      " [4.]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_07, creating vectors and matrices\n",
    "\n",
    "# a length-4 column-vector c of zeros\n",
    "c0 = np.zeros((4,1));\n",
    "print('c0\\n',c0);\n",
    "\n",
    "# a length-4 column-vector c of ones\n",
    "c1 = np.ones((4,1));\n",
    "print('c1\\n',c1);\n",
    "\n",
    "# column vector c\n",
    "# note that decimal points are needed,\n",
    "# else you will get an integer vector\n",
    "# which is usually not what you want\n",
    "c = eda_cvec( [1.,2.,3.,4.] );\n",
    "print('c\\n',c);\n",
    "\n",
    "# row vector r\n",
    "c = eda_cvec([1.,2.,3.,4.]);\n",
    "r=c.T;\n",
    "print('r\\n',r);\n",
    "\n",
    "# 3x3 square matrix A0 of zeros, and\n",
    "# 3x3 square matrix A1 of ones\n",
    "A0 = np.zeros((3,3));\n",
    "print(\"A0\\n\",A0);\n",
    "A1 = np.ones((3,3));\n",
    "print(\"A1\\n\",A1);\n",
    "\n",
    "# square matrix A\n",
    "# note that decimal points are needed,\n",
    "# else you will get an integer vector\n",
    "# which is usually not what you want\n",
    "A = np.array( [[1., 2., 3.],\n",
    "               [4., 5., 6.],\n",
    "               [7., 8., 9.]] );\n",
    "print(\"A\\n\",A);\n",
    "\n",
    "# find shape of c\n",
    "Nc, Mc = np.shape(c);\n",
    "print('shape of c:\\n',Nc,Mc);\n",
    "\n",
    "# find shape of c\n",
    "NA, MA = np.shape(A);\n",
    "print('shape of M\\n',NA,MA);\n",
    "\n",
    "# note that all these calls to eda_cvec\n",
    "# achieve the same result\n",
    "c = eda_cvec( [1.,2.,3.,4.] );\n",
    "print('c\\n',c);\n",
    "c = eda_cvec( 1.,2.,3.,4. );\n",
    "print('c\\n',c);\n",
    "c = eda_cvec( (1.,2.,3.,4.) );\n",
    "print('c\\n',c);\n",
    "c = eda_cvec( 1.,[2.,3.,4.] );\n",
    "print('c\\n',c);\n",
    "c = eda_cvec( [1.,2.],[3.,4.] );\n",
    "print('c\\n',c);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c\n",
      " [[1.]\n",
      " [3.]\n",
      " [5.]]\n",
      "r\n",
      " [[2. 4. 6.]]\n",
      "A\n",
      " [[1. 2. 3.]\n",
      " [4. 5. 6.]\n",
      " [7. 8. 9.]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_08: defining vectors and matrices\n",
    "\n",
    "# column vector c\n",
    "c = eda_cvec( [1.,3.,5.] );\n",
    "print(\"c\\n\",c);\n",
    "\n",
    "# column vector r\n",
    "r = eda_cvec( [2.,4.,6.] ).T;\n",
    "print(\"r\\n\",r);\n",
    "\n",
    "# square matrix A\n",
    "A = np.array( [[1., 2., 3.],\n",
    "               [4., 5., 6.],\n",
    "               [7., 8., 9.]] );\n",
    "print(\"A\\n\",A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a\n",
      " [[1.]\n",
      " [3.]\n",
      " [5.]]\n",
      "b\n",
      " [[2.]\n",
      " [4.]\n",
      " [6.]]\n",
      "c\n",
      " [[ 3.]\n",
      " [ 7.]\n",
      " [11.]]\n",
      "A\n",
      " [[1. 2. 3.]\n",
      " [4. 5. 6.]\n",
      " [7. 8. 9.]]\n",
      "B\n",
      " [[1. 1. 1.]\n",
      " [1. 1. 1.]\n",
      " [1. 1. 1.]]\n",
      "C\n",
      " [[ 2.  3.  4.]\n",
      " [ 5.  6.  7.]\n",
      " [ 8.  9. 10.]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_09: vector and matrix addition\n",
    "\n",
    "# column vector c\n",
    "a = eda_cvec( [1.,3.,5.] );\n",
    "print(\"a\\n\",a);\n",
    "\n",
    "# column vector b\n",
    "b = eda_cvec( [2.,4.,6.] );\n",
    "print(\"b\\n\",b);\n",
    "\n",
    "c = a+b;\n",
    "print(\"c\\n\",c);\n",
    "\n",
    "# square matrix A\n",
    "A = np.array( [[1., 2., 3.],\n",
    "               [4., 5., 6.],\n",
    "               [7., 8., 9.]] );\n",
    "print(\"A\\n\",A);\n",
    "\n",
    "# square matrix B\n",
    "B = np.ones((3,3));\n",
    "print(\"B\\n\",B);\n",
    "\n",
    "C = A+B;\n",
    "print(\"C\\n\",C);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a\n",
      " [[1.]\n",
      " [3.]\n",
      " [5.]]\n",
      "b\n",
      " [[2.]\n",
      " [4.]\n",
      " [6.]]\n",
      "A\n",
      " [[1. 0. 2.]\n",
      " [0. 1. 0.]\n",
      " [2. 0. 1.]]\n",
      "B\n",
      " [[ 1.  0. -1.]\n",
      " [ 0.  2.  0.]\n",
      " [-1.  0.  3.]]\n",
      "s\n",
      " [[44.]]\n",
      "T\n",
      " [[ 2.  4.  6.]\n",
      " [ 6. 12. 18.]\n",
      " [10. 20. 30.]]\n",
      "c\n",
      " [[11.]\n",
      " [ 3.]\n",
      " [ 7.]]\n",
      "P\n",
      " [[-1.  0.  5.]\n",
      " [ 0.  2.  0.]\n",
      " [ 1.  0.  1.]]\n",
      "d\n",
      " [[ 2.]\n",
      " [12.]\n",
      " [30.]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_10: vector and matrix multiplication\n",
    "\n",
    "# column vector c\n",
    "a = eda_cvec( [1.,3.,5.] );\n",
    "print(\"a\\n\",a);\n",
    "\n",
    "# column vector r\n",
    "b = eda_cvec( [2.,4.,6.] );\n",
    "print(\"b\\n\",b);\n",
    "\n",
    "# square matrix A\n",
    "A = np.array( [[1., 0., 2.],\n",
    "               [0., 1., 0.],\n",
    "               [2., 0., 1.]] );\n",
    "print(\"A\\n\",A);\n",
    "\n",
    "# square matrix B\n",
    "B = np.array( [[1., 0., -1.],\n",
    "               [0., 2.,  0.],\n",
    "               [-1., 0., 3.]] );\n",
    "print(\"B\\n\",B);\n",
    "\n",
    "# inner product: s = a-transpose time b\n",
    "s = np.matmul(a.T,b);\n",
    "print(\"s\\n\",s);\n",
    "\n",
    "# outer product: T = a times b-transpose\n",
    "T = np.matmul(a,b.T);\n",
    "print(\"T\\n\",T);\n",
    "\n",
    "# matrix times a vector: P = A times B\n",
    "c = np.matmul(A,a);\n",
    "print(\"c\\n\",c);\n",
    "\n",
    "# matrix times a matrix: P = A times B\n",
    "P = np.matmul(A,B);\n",
    "print(\"P\\n\",P);\n",
    "\n",
    "# d = element-wise multiply of a and b\n",
    "d = np.multiply(a,b);\n",
    "print(\"d\\n\",d);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a\n",
      " [[1.]\n",
      " [2.]\n",
      " [3.]]\n",
      "B\n",
      " [[1. 2. 3.]\n",
      " [4. 5. 6.]\n",
      " [7. 8. 9.]]\n",
      "s\n",
      " 2.0\n",
      "t\n",
      " 6.0\n",
      "b\n",
      " [[2.]\n",
      " [5.]\n",
      " [8.]]\n",
      "r\n",
      " [[4. 5. 6.]]\n",
      "c\n",
      " [[4.]\n",
      " [5.]\n",
      " [6.]]\n",
      "T\n",
      " [[5. 6.]\n",
      " [8. 9.]]\n",
      "X\n",
      " [[5. 6.]\n",
      " [8. 9.]]\n",
      "X\n",
      " [[7. 6.]\n",
      " [8. 9.]]\n",
      "T\n",
      " [[5. 6.]\n",
      " [8. 9.]]\n",
      "B\n",
      " [[1. 2. 3.]\n",
      " [4. 5. 6.]\n",
      " [7. 8. 9.]]\n",
      "Bp\n",
      " [[10.  2.  3.]\n",
      " [11.  5.  6.]\n",
      " [12.  8.  9.]]\n",
      "B.ravel() [1. 2. 3. 4. 5. 6. 7. 8. 9.]\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_11: accessing elements of a matrix\n",
    "\n",
    "# column vector a\n",
    "a = eda_cvec([1.,2.,3.]);\n",
    "print(\"a\\n\",a);\n",
    "\n",
    "# square matrix M\n",
    "B = np.array( [[1., 2., 3.],\n",
    "               [4., 5., 6.],\n",
    "               [7., 8., 9.]] );\n",
    "print(\"B\\n\",B);\n",
    "\n",
    "# one element of a column vector\n",
    "s = np.copy(a[1,0]);\n",
    "print(\"s\\n\",s);\n",
    "\n",
    "# one element of a matrix\n",
    "t = np.copy(B[1,2]);\n",
    "print(\"t\\n\",t);\n",
    "\n",
    "# column vector = one column of a matrix\n",
    "b = np.copy( B[0:3,1:2] );\n",
    "print(\"b\\n\",b);\n",
    "\n",
    "# row vector = one row of a matrix\n",
    "r = np.copy( B[1:2,0:3] );\n",
    "print(\"r\\n\",r);\n",
    "\n",
    "# column vector = one row of a matrix\n",
    "c = np.copy( B[1:2,0:3].T );\n",
    "print(\"c\\n\",c);\n",
    "\n",
    "# small matrix = rectangular block of larger matrixx\n",
    "# note that end of range is one larger than you want\n",
    "T = np.copy( B[1:3,1:3] );\n",
    "print(\"T\\n\",T);\n",
    "\n",
    "# the purpose for the copy is to convert a view of a matrix\n",
    "# into a separate matrix.  Else changeing the original matrix\n",
    "# will change the \n",
    "X =  B[1:2+1,1:2+1];\n",
    "print(\"X\\n\",X);\n",
    "B[1,1]=7;\n",
    "print(\"X\\n\",X);\n",
    "print(\"T\\n\",T);\n",
    "\n",
    "# set first column of matrix B to a column vector b\n",
    "b = eda_cvec([10.,11.,12.]);\n",
    "B= np.array( [[1., 2., 3.],\n",
    "               [4., 5., 6.],\n",
    "               [7., 8., 9.]] );\n",
    "print(\"B\\n\",B);\n",
    "Bp = np.copy(B);\n",
    "# the shape of the view on the left-hand-side of the equal sign\n",
    "# must exactly match the shape of the quantity on the right hand side\n",
    "Bp[0:3,0:1] = b;\n",
    "print(\"Bp\\n\",Bp);\n",
    "\n",
    "# ravel a N by M matrix A into a NM by 0 numpy array\n",
    "Blist = B.ravel();\n",
    "print(\"B.ravel()\",Blist);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "s\n",
      " [[ 3.1286893 ]\n",
      " [ 3.05165059]\n",
      " [ 2.89947011]\n",
      " [ 2.67589505]\n",
      " [ 2.38643058]\n",
      " [ 2.03820426]\n",
      " [ 1.6397906 ]\n",
      " [ 1.20099984]\n",
      " [ 0.73263649]\n",
      " [ 0.24623319]\n",
      " [-0.24623319]\n",
      " [-0.73263649]\n",
      " [-1.20099984]\n",
      " [-1.6397906 ]\n",
      " [-2.03820426]\n",
      " [-2.38643058]\n",
      " [-2.67589505]\n",
      " [-2.89947011]\n",
      " [-3.05165059]\n",
      " [-3.1286893 ]]\n",
      "length\n",
      " 20\n",
      "a\n",
      "[[0.        ]\n",
      " [0.00782172]\n",
      " [0.02327257]\n",
      " [0.0459721 ]\n",
      " [0.07536136]\n",
      " [0.1107167 ]\n",
      " [0.15116755]\n",
      " [0.19571788]\n",
      " [0.2432707 ]\n",
      " [0.29265512]\n",
      " [0.34265512]\n",
      " [0.39203954]\n",
      " [0.43959236]\n",
      " [0.48414269]\n",
      " [0.52459354]\n",
      " [0.55994888]\n",
      " [0.58933814]\n",
      " [0.61203766]\n",
      " [0.62748851]\n",
      " [0.63531024]\n",
      " [0.63531024]]\n",
      "atotal\n",
      " [0.63531024]\n",
      "atotal\n",
      " 0.6353102368087353\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 14452 characters]\n",
      "text/plain": [
       "<Figure size 576x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 35588 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_01_12: derivatives and integrals\n",
    "\n",
    "# column vector t of N evenly-spaced points beterrn tmin and tmax\n",
    "N=21;       # length of t\n",
    "tmin = 0.0;   # minimum time\n",
    "tmax = 1.0;   # maximum time\n",
    "Dt = (tmax-tmin)/(N-1);  # time increment\n",
    "t = eda_cvec( np.linspace( tmin, tmax, N ) );\n",
    "\n",
    "# expemplary function\n",
    "x=np.sin(pi*t);\n",
    "\n",
    "# approximate slope; note one point shorter then x\n",
    "# sapprox = eda_cvec( (x[1:N+1]-x[0:N-2+1])/Dt );\n",
    "s = (1/Dt) * ( x[1:N,0:1]-x[0:N-1,0:1] );\n",
    "Ns, M = np.shape(s);\n",
    "print('s\\n',s);\n",
    "print('length\\n',Ns);\n",
    "\n",
    "# approximate integral; same length as x\n",
    "# integral is approximated as a cumulative sum\n",
    "a = eda_cvec( Dt*np.cumsum(x) );\n",
    "print('a');\n",
    "print(a);\n",
    "\n",
    "# integral over all t; integral approximated as a sum\n",
    "atotal = Dt*np.sum(x);\n",
    "print('atotal\\n', a[N-1] );\n",
    "print('atotal\\n', atotal );\n",
    "\n",
    "# the rest of this section makes plots used in the text,\n",
    "# but since plotting has not yet been discussed, readers\n",
    "# may prefer to skip it\n",
    "\n",
    "# true derivative and integral\n",
    "strue = pi*np.cos(pi*t);      # derivative\n",
    "atrue = 1/pi-np.cos(pi*t)/pi; # integral, with integration constant\n",
    "                              # chosen so integral is zero at t=0\n",
    "    \n",
    "# time series plot\n",
    "fig1 = plt.figure(1,figsize=(8,3));\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis([0.0, 1.0, 0, 1.1]);\n",
    "plt.plot(t,x,\"-\",lw=4,color=[0.8,0.8,0.8]);\n",
    "plt.plot(t,x,\"ko\",lw=2);\n",
    "plt.xlabel(\"t\");\n",
    "plt.ylabel(\"x(t)\");\n",
    "for i in range(N-1):\n",
    "    plt.plot( eda_cvec(t[i,0],t[i]), eda_cvec(0.0,0.1), 'k-', lw=2 );\n",
    "plt.show();\n",
    "\n",
    "\n",
    "# calculus plot\n",
    "fig2 = plt.figure(2,figsize=(8,6));\n",
    "ax1 = plt.subplot(3, 1, 1);\n",
    "plt.axis([0.0, 1.0, 0, 1.1]);\n",
    "plt.xlabel(\"t\");\n",
    "plt.ylabel(\"x(t)\");\n",
    "# plot exemplary rectangles, shaded\n",
    "K=5;\n",
    "for k in range(K):\n",
    "    if(x[k,0]>=0.0):\n",
    "        r=patches.Rectangle( (t[k,0],0.0), Dt, x[k,0], ec=[0.0,0.0,0.0], fc=[0.9,0.9,0.9] );\n",
    "        ax1.add_patch(r);\n",
    "    # plot function as continuous curve\n",
    "plt.plot( t, x, lw=4, color=[0.9,0.9,0.9] );\n",
    "# plot function as circles\n",
    "plt.plot( t, x, 'ko', lw=2 );\n",
    "# embolden one slope\n",
    "plt.plot( eda_cvec(t[K,0],t[K+1,0]), eda_cvec(x[K,0],x[K+1,0]), '-', lw=6, color=[0,0,0]);\n",
    "# plot trapezoids\n",
    "for i in range(N-1):\n",
    "    tb = eda_cvec( t[i,0],t[i,0],t[i+1,0],t[i+1,0] );\n",
    "    xb = eda_cvec(0.0,x[i,0],x[i+1,0],0.0 );\n",
    "    plt.plot( tb, xb, 'k-', lw=2 );\n",
    "     \n",
    "# plot derivative\n",
    "ax2 = plt.subplot(3, 1, 2);\n",
    "plt.axis( [0, 1, -1.1*pi, 1.1*pi ] );\n",
    "plt.xlabel(\"t\");\n",
    "plt.ylabel(\"s=dx/dx\");\n",
    "plt.plot( t, strue, '-', lw=4, Color=[0.8,0.8,0.8] );\n",
    "plt.plot( t[0:N-1,0:1], s, 'k-', lw=2 );\n",
    "     \n",
    "# plot integral\n",
    "ax3 = plt.subplot(3, 1, 3);\n",
    "plt.axis( [0.0, 1.0, 0.0, 1.1*2/pi ] );\n",
    "plt.xlabel(\"t\");\n",
    "plt.ylabel(\"integral of x\");\n",
    "plt.plot( t, atrue, 'k-', lw=4, color=[0.8,0.8,0.8] );\n",
    "plt.plot( t, a, 'k-', lw=2 );\n",
    "plt.plot( t[N-1,0], atotal, 'ko', lw=3 ); # last value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A\n",
      " [[1. 2. 3.]\n",
      " [4. 5. 6.]\n",
      " [7. 8. 9.]]\n",
      "a\n",
      " [[1.]\n",
      " [5.]\n",
      " [9.]]\n",
      "B\n",
      " [[3. 2. 1.]\n",
      " [6. 5. 4.]\n",
      " [9. 8. 7.]]\n",
      "c\n",
      " [[1. 1.]\n",
      " [2. 2.]\n",
      " [1. 1.]\n",
      " [4. 4.]\n",
      " [3. 3.]\n",
      " [2. 2.]\n",
      " [6. 6.]\n",
      " [4. 4.]\n",
      " [9. 6.]\n",
      " [2. 2.]\n",
      " [1. 1.]\n",
      " [4. 4.]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_13: for-loops\n",
    "\n",
    "# square matrix A\n",
    "A = np.array( [[1., 2., 3.],\n",
    "               [4., 5., 6.],\n",
    "               [7., 8., 9.]] );\n",
    "print(\"A\\n\",A);\n",
    "\n",
    "# for loop that extracts the main diagonal of a matrix\n",
    "a = np.zeros((3,1));\n",
    "for i in range(3):\n",
    "    a[i,0]=A[i,i];\n",
    "print(\"a\\n\",a)\n",
    "\n",
    "# two nested for loops used to rearrange (flip left-right)\n",
    "# the elements of a matrix\n",
    "B = np.zeros((3,3));\n",
    "for i in range(3):     # loop over all rows\n",
    "    for j in range(3): # loop over all cols\n",
    "        B[i,2-j]=A[i,j]; # left-right-flip\n",
    "print(\"B\\n\",B);    \n",
    "\n",
    "# column vector of length 12\n",
    "\n",
    "a = eda_cvec( [1.0, 2.0, 1.0, 4.0, 3.0, 2.0, 6.0, 4.0, 9.0, 2.0, 1.0, 4.0] );\n",
    "N, M = np.shape(a);\n",
    "\n",
    "# for loop that \"clips\" a vector\n",
    "b = np.zeros((N,1));\n",
    "for i in range(N):\n",
    "    if a[i,0] >= 6.:\n",
    "        b[i,0]=6.;\n",
    "    else:\n",
    "        b[i,0]=a[i,0];\n",
    "\n",
    "# side-by-side concatenation vectors a and b for aeshetic printing\n",
    "c = np.concatenate((a,b),axis=1);\n",
    "print(\"c\\n\",c);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A\n",
      " [[1. 2. 3.]\n",
      " [4. 5. 6.]\n",
      " [7. 8. 9.]]\n",
      "a\n",
      " [[1.]\n",
      " [5.]\n",
      " [9.]]\n",
      "N\n",
      " [[3. 2. 1.]\n",
      " [6. 5. 4.]\n",
      " [9. 8. 7.]]\n",
      "a\n",
      " [[ 1.]\n",
      " [ 2.]\n",
      " [10.]\n",
      " [ 4.]\n",
      " [ 3.]\n",
      " [ 2.]\n",
      " [ 6.]\n",
      " [ 4.]\n",
      " [ 9.]\n",
      " [ 2.]\n",
      " [ 1.]\n",
      " [ 4.]]\n",
      "b\n",
      " [[1.]\n",
      " [2.]\n",
      " [6.]\n",
      " [4.]\n",
      " [3.]\n",
      " [2.]\n",
      " [6.]\n",
      " [4.]\n",
      " [6.]\n",
      " [2.]\n",
      " [1.]\n",
      " [4.]]\n",
      "c\n",
      " [[1.]\n",
      " [2.]\n",
      " [6.]\n",
      " [4.]\n",
      " [3.]\n",
      " [2.]\n",
      " [6.]\n",
      " [4.]\n",
      " [6.]\n",
      " [2.]\n",
      " [1.]\n",
      " [4.]]\n",
      "type r\n",
      " <class 'tuple'>\n",
      "shape r\n",
      " (2, 3)\n",
      "type k=r[0]\n",
      " <class 'numpy.ndarray'>\n",
      "shape k\n",
      " (3,)\n",
      "k\n",
      " [2 6 8]\n",
      " \n",
      "d\n",
      " [[1.]\n",
      " [2.]\n",
      " [6.]\n",
      " [4.]\n",
      " [3.]\n",
      " [2.]\n",
      " [6.]\n",
      " [4.]\n",
      " [6.]\n",
      " [2.]\n",
      " [1.]\n",
      " [4.]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_14: alternatives to for-loops\n",
    "\n",
    "# square matric A\n",
    "A = np.array( [[1., 2., 3.],\n",
    "               [4., 5., 6.],\n",
    "               [7., 8., 9.]] );\n",
    "print(\"A\\n\",A);\n",
    "\n",
    "# extract diagonal of matcix into a vector\n",
    "a = eda_cvec( np.diagonal(A) );\n",
    "print(\"a\\n\",a);\n",
    "\n",
    "# B is A flipped left to right\n",
    "B = np.fliplr(A);\n",
    "print(\"N\\n\",B);\n",
    "\n",
    "# column vector a of lenth 12\n",
    "a = eda_cvec( [1.0, 2.0, 10.0, 4.0, 3.0, 2.0, 6.0, 4.0, 9.0, 2.0, 1.0, 4.0] );\n",
    "print(\"a\\n\",a);\n",
    "\n",
    "# alternative method #1 of cliping a vector\n",
    "b = np.multiply(a, a<6.) + np.multiply(6., a>=6.);\n",
    "print(\"b\\n\",b);\n",
    "\n",
    "# alternative method #2 of cliping a vector\n",
    "c=np.copy(a); # start with c that is a copy of a\n",
    "c[a>=6.]=6.;\n",
    "print(\"c\\n\",c);\n",
    "\n",
    "# alternative method #3 f cliping a vector\n",
    "# we use the np.where() command to finds the\n",
    "# rows of a that satisfy the condition\n",
    "d=np.copy(a);        # start with a d that is a copy of a\n",
    "r = np.where(a>=6.); # r is a data structure descibed below\n",
    "k = r[0];            # k is a Nx0 list of indices\n",
    "d[k,0:1]=6.;         # reset the rows of d that violate the condition\n",
    "\n",
    "# r is a tupple\n",
    "print(\"type r\\n\",type(r));\n",
    "print(\"shape r\\n\",np.shape(r));\n",
    "# k=r[0] gives the rows satisfying the a>=6 condition. It is a 3-by-0 np.array\n",
    "print(\"type k=r[0]\\n\",type(k));\n",
    "print(\"shape k\\n\",np.shape(k));\n",
    "print(\"k\\n\",k);\n",
    "print(\" \");\n",
    "print(\"d\\n\",d);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A:\n",
      " [[1.  0.5 0. ]\n",
      " [0.5 5.  0. ]\n",
      " [0.  0.  9. ]]\n",
      "B:\n",
      " [[ 1.05263158 -0.10526316  0.        ]\n",
      " [-0.10526316  0.21052632  0.        ]\n",
      " [ 0.          0.          0.11111111]]\n",
      "AB:\n",
      " [[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n",
      "BA:\n",
      " [[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_15: matrix inverse\n",
    "\n",
    "# square matrix A\n",
    "A = np.array( [[1.0, 0.5, 0.0],\n",
    "               [0.5, 5.0, 0.0],\n",
    "               [0.0, 0.0, 9.0]] );\n",
    "print(\"A:\\n\",A);\n",
    "\n",
    "# B = inverse of A\n",
    "B = np.linalg.inv(A);\n",
    "print(\"B:\\n\",B);\n",
    "\n",
    "# A times B\n",
    "I1 = np.matmul(A,B);\n",
    "print(\"AB:\\n\",I1);\n",
    "\n",
    "# B times A\n",
    "I2 = np.matmul(B,A);\n",
    "print(\"BA:\\n\",I2);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b\n",
      " [[2.]\n",
      " [4.]\n",
      " [6.]]\n",
      "c\n",
      " [[ 2.00000000e+00]\n",
      " [-1.24900090e-16]\n",
      " [ 4.16333634e-17]]\n",
      "error\n",
      " [[-4.4408921e-16]\n",
      " [-8.8817842e-16]\n",
      " [-8.8817842e-16]]\n",
      "D\n",
      " [[ 4.00000000e+00  3.50000000e+00  3.00000000e+00]\n",
      " [-1.66666667e-01 -8.33333333e-02 -1.87350135e-16]\n",
      " [-1.66666667e-01 -8.33333333e-02  6.24500451e-17]]\n",
      "error\n",
      " [[-6.66133815e-16  0.00000000e+00  0.00000000e+00]\n",
      " [-1.77635684e-15 -8.88178420e-16  0.00000000e+00]\n",
      " [-2.66453526e-15  0.00000000e+00  0.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "# edapy_01_16: solving linear system\n",
    "\n",
    "# square matrix A\n",
    "A = np.array([[1.0, 5.0,  13.0],\n",
    "              [2.0, 7.0,  17.0],\n",
    "              [3.0, 11.0, 19.0]] );\n",
    "\n",
    "\n",
    "# column vercor b\n",
    "b = eda_cvec( [2., 4., 6.] );\n",
    "print(\"b\\n\",b);\n",
    "\n",
    "# when Ac=b\n",
    "# alternative to c = A-inverse times b\n",
    "c = np.linalg.solve(A,b);\n",
    "print(\"c\\n\",c);\n",
    "\n",
    "# check result is correct\n",
    "error1=b-np.matmul(A,c);\n",
    "print(\"error\\n\",error1);\n",
    "\n",
    "B = np.array( [[1.0, 2.0, 3.0],\n",
    "               [4.0, 5.0, 6.0],\n",
    "               [7.0, 8.0, 9.0]] );\n",
    "\n",
    "# when AD=B\n",
    "# alternative to D = A-inverse B\n",
    "D = np.linalg.solve(A,B);\n",
    "print(\"D\\n\",D);\n",
    "\n",
    "# check result is correct\n",
    "error2=B-np.matmul(A,D);\n",
    "print(\"error\\n\",error2);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 32616 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_01_17: load text file and make very-quick plot\n",
    "# this example uses the neuse.txt stream flow dataset\n",
    "\n",
    "# load rectangular, tab-separated table of numbers into a N by M matrix D\n",
    "D = np.genfromtxt(\"../Data/neuse.txt\", delimiter=\"\\t\");\n",
    "N, M = np.shape(D);\n",
    "\n",
    "t = np.copy( D[0:N,0:1] );  # time in column 1\n",
    "d = np.copy( D[0:N,1:2] );  # discharge in column 2\n",
    "\n",
    "# plot the data\n",
    "fig1 = plt.figure(1,figsize=(8,4));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);           # only one subplot\n",
    "plt.plot(t,d,\"k-\");                   # plot d(t) with black line\n",
    "plt.xlabel(\"time (days)\");           # label time axis\n",
    "plt.ylabel(\"discharge (cfs)\");       # label data axis\n",
    "plt.show();                           # display plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 36248 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_01_18: load text file and mmke nicer plot\n",
    "# this example uses the neuse.txt stream flow dataset\n",
    "\n",
    "# load data from text fole\n",
    "D = np.genfromtxt('../Data/neuse.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "\n",
    "# time in column 1, discharge in colun 2\n",
    "t = eda_cvec( D[:,0] );\n",
    "d = eda_cvec( D[:,1] );\n",
    "\n",
    "# nicer plot\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.axis([0, 4500, 0, 1.1*np.max(d)])\n",
    "plt.plot(t,d,\"k-\");\n",
    "plt.xlabel(\"time (days)\");\n",
    "plt.ylabel(\"discharge (cfs)\");\n",
    "plt.title(\"Neuse River Hydrograph\");\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 19728 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_01_19: log-log plot\n",
    "\n",
    "# get 10-year list of earthquake magnitudes from file\n",
    "# note: since the file contains only one column of data, the\n",
    "# np.genfromtxt puts the data into an N by 0 array.  So we\n",
    "# call eda_cvec to put it into a standard columnn vector\n",
    "d = eda_cvec( np.genfromtxt('../Data/earthquake_magnitudes.txt', delimiter='\\t') );\n",
    "\n",
    "time = 10.0; # 10 years of data\n",
    "\n",
    "# count up earthquakes with magnitudes in set of magnitude bins\n",
    "# 51 bins from magnitude 5 to magnitude 10\n",
    "c, e = np.histogram(d,51,(5,10));\n",
    "Nc = len(c);  # counts of earthquakes in each bin\n",
    "Ne = len(e);  # magnitudes at the edges of the bins\n",
    "count = eda_cvec(c);\n",
    "edges = eda_cvec(e);\n",
    "\n",
    "# convert count to rate\n",
    "r = eda_cvec( count[0:Nc-1,0]/time );\n",
    "\n",
    "# compute center point of magnitude bins\n",
    "m = eda_cvec( (edges[0:Ne-2,0]+edges[1:Ne-1,0])/2 );\n",
    "\n",
    "# convert magntude to energy using the Gutenberg-Richter magnitude-energy relation\n",
    "log10Ee = 1.5*m + 11.8; # in ergs\n",
    "log10Ej = log10Ee - 7;  # in joules\n",
    "E = 1e-7 * np.power(10,log10Ej);\n",
    "\n",
    "# log-log plot\n",
    "fig2 = plt.figure(2,figsize=(8,4));\n",
    "ax1 = plt.subplot(1, 1, 1);\n",
    "plt.loglog(E, r, \"k-\", basex=10.0, basey=10.0); # plot data as lines\n",
    "plt.loglog(E, r, \"ko\", basex=10.0, basey=10.0); # plot data as circles\n",
    "plt.axis([10**5, 10**12, 10**-1, 10**3])\n",
    "plt.xlabel(\"earthquake energy (joules)\");\n",
    "plt.ylabel(\"earthquakes rate (per year)\");\n",
    "plt.title(\"Earthquake rate vs energy\");\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 32716 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_01_20: write matrix to file\n",
    "\n",
    "# load the neuse.txt data set\n",
    "# convert discharge from cfs to m3/s\n",
    "# and write to a file neuse_metric.txt\n",
    "\n",
    "# get matrix from file\n",
    "D = np.genfromtxt('../Data/neuse.txt', delimiter='\\t')\n",
    "N, M = np.shape(D);\n",
    "\n",
    "# time in column 0\n",
    "t = D[0:N,0:1];\n",
    "\n",
    "# dischange in cfs in column 1\n",
    "d = D[0:N,1:2];\n",
    "\n",
    "# conversion from cfs to m3/s\n",
    "f = 35.3146;\n",
    "dm = d/f;\n",
    "\n",
    "# make simple plot to veryify result \"by eye\"\n",
    "# plot the data\n",
    "fig1 = plt.figure(1,figsize=(8,4));   # smallish figure\n",
    "ax1 = plt.subplot(1, 1, 1);           # only one subplot\n",
    "plt.plot(t,dm,\"k-\");                  # plot d(t) with black line\n",
    "plt.xlabel(\"time in days\");           # label time axis\n",
    "plt.ylabel(\"discharge in m3/s\");      # label data axis\n",
    "plt.show();                           # display plot\n",
    "\n",
    "# concatenate time, discharge into a matrix C\n",
    "C = np.concatenate( (t,dm),axis=1); \n",
    "# save matrix to file of tab-separated values\n",
    "np.savetxt(\"../Data/neuse_metric.txt\",C, delimiter=\"\\t\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 26352 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_01_21: big-picture commenting (recommended)\n",
    "\n",
    "# This script evaluates the 2D Normal p.d.f. p(d1,d2) with\n",
    "# mean (d1bar,d2bar) and covariance C on a L-by-L grid\n",
    "\n",
    "# grid size\n",
    "L=101;\n",
    "\n",
    "# d1 axis\n",
    "d1min = 0.0;\n",
    "d1max = 10.0;\n",
    "d1 = eda_cvec( np.linspace( d1min, d1max, L ) );\n",
    "\n",
    "# d2 axis\n",
    "d2min = 0.0;\n",
    "d2max = 10.0;\n",
    "d2 = eda_cvec( np.linspace( d2min, d2max, L ) );\n",
    "\n",
    "# mean\n",
    "d1bar = 4.0;\n",
    "d2bar = 6.0;\n",
    "\n",
    "# covariance matrix\n",
    "C = np.array( [[1.0, 0.5],\n",
    "               [0.5, 2.0]] );\n",
    "\n",
    "# Evaluate Normal p.d.f. Pp\n",
    "# with mean dbar and covariance C\n",
    "dd = np.zeros((2,1));     \n",
    "CI = np.linalg.inv(C);\n",
    "detC = np.linalg.det(C);\n",
    "norm = 1/(2*pi*sqrt(detC));\n",
    "Pp = np.zeros((L,L));\n",
    "for i in range(L):     # loop over d1\n",
    "    for j in range(L): # loop over d2\n",
    "        dd[0,0] = d1[i,0]-d1bar;\n",
    "        dd[1,0] = d2[j,0]-d2bar;\n",
    "        E = np.matmul(dd.T,np.matmul(CI,dd));\n",
    "        # Normal p.d.f.\n",
    "        Pp[i,j] = norm*exp(-0.5* E);\n",
    "        \n",
    "fig1 = plt.figure(1,figsize=(8,4)); # figure 1\n",
    "ax1 = plt.subplot(1, 1, 1);         # single subplot\n",
    "plt.title(\"normal p.d.f.\");         # title\n",
    "plt.imshow(Pp);                     # plot Pp\n",
    "plt.show();                         # show image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 26352 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# edapy_01_22 Comments that focus on minutia (not recommended)\n",
    "\n",
    "# create and plot the array p(d1,d2)\n",
    "\n",
    "# define L\n",
    "L=101;\n",
    "\n",
    "# create d1\n",
    "d1min = 0.0;\n",
    "d1max = 10.0;\n",
    "d1 = eda_cvec( np.linspace( d1min, d1max,L ) );\n",
    "\n",
    "# create d2\n",
    "d2min = 0.0;\n",
    "d2max = 10.0;\n",
    "d2 = eda_cvec( np.linspace( d2min, d2max, L ) );\n",
    "\n",
    "# define d1bar and d2bar\n",
    "d1bar = 4.;\n",
    "d2bar = 6.;\n",
    "\n",
    "# define C\n",
    "C = np.zeros((2,2));\n",
    "C = np.array( [[1.0, 0.5],\n",
    "               [0.5, 2.0]] );\n",
    "\n",
    "# define dd\n",
    "dd = np.zeros((2,1));  \n",
    "# take inverse\n",
    "CI = np.linalg.inv(C);  \n",
    "# take determinant\n",
    "detC = np.linalg.det(C);  \n",
    "# define normalization\n",
    "norm = 1/(2*pi*sqrt(detC));\n",
    "Pp = np.zeros((L,L));   # initialize Pp\n",
    "for i in range(L):      # loop over i\n",
    "    for j in range(L):  # loop over j\n",
    "        # det up dd\n",
    "        dd[:,0] = [d1[i,0]-d1bar, d2[j,0]-d2bar];\n",
    "        # compute E\n",
    "        E = np.matmul(dd.T,np.matmul(CI,dd));\n",
    "        # peform exponential\n",
    "        Pp[i,j] = norm*exp(-0.5 * E );\n",
    "        \n",
    "fig1 = plt.figure(1,figsize=(8,4)); # figure 1\n",
    "ax1 = plt.subplot(1, 1, 1);         # single subplot\n",
    "plt.title(\"normal p.d.f.\");         # title\n",
    "plt.imshow(Pp);                     # plot Pp\n",
    "plt.show();                         # show image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  }
 ],
 "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
}
