{
 "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\n",
    "from math import exp, pi, sin, cos, sqrt, floor, ceil, atan2\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",
    "\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",
    "    return 1;\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec1(v) from entity v that is array-like,\n",
    "# including a number, a list of numbers, a tuple\n",
    "# of numbers, an N-by-0 np array and a N-by1 np array.\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it.\n",
    "# This version takes just one input argment.\n",
    "# whereas c=eda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def eda_cvec1(v):\n",
    "    if isinstance(v, int):\n",
    "        w = np.zeros((1,1),dtype=int);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, float):\n",
    "        w = np.zeros((1,1),dtype=float);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, complex):\n",
    "        w = np.zeros((1,1),dtype=complex);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, np.ndarray):\n",
    "        s = np.shape(v);\n",
    "        if len(s) == 1:\n",
    "            return np.copy(np.reshape(v,(s[0],1)));\n",
    "        else:\n",
    "            [r,c]=s;\n",
    "            if( c==1 ):\n",
    "                return(np.copy(v));\n",
    "            elif(r==1):\n",
    "                return(np.copy(v.T));\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: %d by %d ndarray not allowed\" % (r, c));\n",
    "    elif isinstance(v, list):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, tuple):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"eda_cvec: %s not supported\" % type(v));\n",
    "        \n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec(v1, v2, ...) from a list of several\n",
    "# array-like entities v1, v2, including a number\n",
    "# a list of numbers, a tuple of numbers, an N-by-0 np array\n",
    "# and a N-by-1 np array. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context.\n",
    "# This version concatenates many argments,\n",
    "# whereas c=eda_cvec1(v1) takes just one argiment.\n",
    "# I recommend always using eda_cvec(v1, v2, ...)\n",
    "def eda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        Nt = Nt + N;\n",
    "        if (t==int) and isinstance(v[0,0],float):\n",
    "            t=float;\n",
    "        elif isinstance(v[0,0],complex):\n",
    "            t=complex;\n",
    "    w = np.zeros((Nt,1),dtype=t);\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": [
    "# eda13_01\n",
    "# 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": [
    "# eda13_02\n",
    "# simple examples of functions\n",
    "\n",
    "# definition of area of circle function\n",
    "def areaofcircle(radius):\n",
    "    area = pi*(radius**2);\n",
    "    return area;\n",
    "\n",
    "# definition of circumfernce and area of rectangle function\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 circumfernece 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": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test matrix A\n",
      "[[1 2 3]\n",
      " [4 5 6]\n",
      " [7 8 9]]\n",
      " \n",
      "column vector b formed from A\n",
      "[[1]\n",
      " [2]\n",
      " [3]\n",
      " [4]\n",
      " [5]\n",
      " [6]\n",
      " [7]\n",
      " [8]\n",
      " [9]]\n",
      " \n",
      "row vector formed from A\n",
      "[[1 2 3 4 5 6 7 8 9]]\n",
      " \n",
      "the matrix reformed from the vector b\n",
      "[[1 2 3]\n",
      " [4 5 6]\n",
      " [7 8 9]]\n",
      " \n",
      "for a (3,3) matrix M[i,j] unraveled into a vector b[k]\n",
      "(i,j) of (1,0) corresponds to k of 3\n",
      "k of 3 corresponds to (i,j) of (1,0)\n"
     ]
    }
   ],
   "source": [
    "# eda13_03\n",
    "# examples of reorganizing a matrix\n",
    "\n",
    "# define a matrix\n",
    "N=3;\n",
    "M=3;\n",
    "A = np.array( [[1, 2, 3],\n",
    "               [4, 5, 6],\n",
    "               [7, 8, 9]] );\n",
    "print('test matrix A');\n",
    "print(A);\n",
    "print(' ');\n",
    "\n",
    "# note that a matrix can be referred to using only a single\n",
    "# index, in which case, it acts like a vector constructed\n",
    "# by concatenating the columns of the matrix, one on top\n",
    "# of the other\n",
    "# disp('A(4) is A(1,2)');\n",
    "#(4)\n",
    "\n",
    "# create a column-vector, b, containing the elements of A, organized\n",
    "# column-wide.\n",
    "b = np.reshape(A,(N*M,1));\n",
    "print('column vector b formed from A');\n",
    "print(b);\n",
    "print(' ');\n",
    "\n",
    "# create a row-vector, c, containing the elements of A, organized\n",
    "# column-wise.\n",
    "c = np.reshape(A,(1,N*M));\n",
    "print('row vector formed from A');\n",
    "print(c);\n",
    "print(' ');\n",
    "\n",
    "# the matrix reformed from the vector\n",
    "B = np.reshape(b,(N,M));\n",
    "print('the matrix reformed from the vector b');\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 b[k]' % (N,M) );\n",
    "k = np.ravel_multi_index( (i,j), (N,M) ) ;\n",
    "print('(i,j) of (%d,%d) corresponds to k of %d' % (i,j,k) );\n",
    "ii, jj = np.unravel_index( k, (N,M) ) ;\n",
    "print('k of %d corresponds to (i,j) of (%d,%d)' % (k,i,j) );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "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": [
    "# eda12_04\n",
    "# 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": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "argument is a floating point numbers\n",
      "eda_cvec(1.0):\n",
      "[[1.]]\n",
      " \n",
      "argument are several floating point numbers\n",
      "eda_cvec(1.0,2.0,3.0):\n",
      "[[1.]\n",
      " [2.]\n",
      " [3.]]\n",
      " \n",
      "argument is a list of floating point numbers\n",
      "c=eda_cvec([1.0,2.0,3.0]):\n",
      "[[1.]\n",
      " [2.]\n",
      " [3.]]\n",
      " \n",
      "argument is a tupple of floating point numbers\n",
      "d=eda_cvec((4.0,5.0,6.0)):\n",
      "[[4.]\n",
      " [5.]\n",
      " [6.]]\n",
      " \n",
      "arguments are N-by-1 np arrays\n",
      "eda_cvec(c,d):\n",
      "[[1.]\n",
      " [2.]\n",
      " [3.]\n",
      " [4.]\n",
      " [5.]\n",
      " [6.]]\n",
      " \n",
      "argument is a list of integers\n",
      "eda_cvec([1,2,3]):\n",
      "[[1]\n",
      " [2]\n",
      " [3]]\n",
      " \n",
      "argument is a list containing a complex\n",
      "eda_cvec([complex(0,0,1.0),2,3]):\n",
      "[[0.+1.j]\n",
      " [2.+0.j]\n",
      " [3.+0.j]]\n",
      " \n"
     ]
    }
   ],
   "source": [
    "# eda12_05\n",
    "# example of use of eda_cvec() function\n",
    "\n",
    "# the rules:\n",
    "#   any number of arguments\n",
    "#   an argument can be:\n",
    "#       a number\n",
    "#       a list of numbers\n",
    "#       a tupple of numbers\n",
    "#       an N-by-0 np array\n",
    "#       an N-by-1 np array\n",
    "#       a 1-by-N np array\n",
    "#   all arguments are concatenated into a single N-by-1 np array\n",
    "#   if at least one argument is complex, everything becomes complex\n",
    "#   if there are no complex, but at least one argument is float, everything becomes float\n",
    "#   if all arguments are integer, its stays integer\n",
    "\n",
    "print(\"argument is a floating point numbers\");\n",
    "print(\"eda_cvec(1.0):\");\n",
    "a = eda_cvec(1.0);\n",
    "print(a);\n",
    "print(\" \");\n",
    "\n",
    "print(\"argument are several floating point numbers\");\n",
    "print(\"eda_cvec(1.0,2.0,3.0):\");\n",
    "b = eda_cvec(1.0,2.0,3.0);\n",
    "print(b);\n",
    "print(\" \");\n",
    "\n",
    "print(\"argument is a list of floating point numbers\");\n",
    "print(\"c=eda_cvec([1.0,2.0,3.0]):\");\n",
    "c = eda_cvec([1.0,2.0,3.0]);\n",
    "print(c);\n",
    "print(\" \");\n",
    "\n",
    "print(\"argument is a tupple of floating point numbers\");\n",
    "print(\"d=eda_cvec((4.0,5.0,6.0)):\");\n",
    "d = eda_cvec((4.0,5.0,6.0));\n",
    "print(d);\n",
    "print(\" \");\n",
    "\n",
    "print(\"arguments are N-by-1 np arrays\");\n",
    "print(\"eda_cvec(c,d):\");\n",
    "e = eda_cvec(c,d);\n",
    "print(e);\n",
    "print(\" \");\n",
    "\n",
    "print(\"argument is a list of integers\");\n",
    "print(\"eda_cvec([1,2,3]):\");\n",
    "e = eda_cvec([1,2,3]);\n",
    "print(e);\n",
    "print(\" \");\n",
    "\n",
    "print(\"argument is a list containing a complex\");\n",
    "print(\"eda_cvec([complex(0,0,1.0),2,3]):\");\n",
    "e = eda_cvec([complex(0.0,1.0),2.0,3]);\n",
    "print(e);\n",
    "print(\" \");"
   ]
  },
  {
   "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
}
