{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numpy  - vectors and matrices\n",
    "\n",
    "Let's start with numpy and working with vectors and matrices, at first we \n",
    "have a look on numpy arrays and try the matrix-vector-multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.23, 2.23, 3.23, 4.23] <class 'list'>\n",
      "[[1.23 2.23 3.23 4.23]] <class 'numpy.ndarray'>\n",
      "2 (1, 4)\n",
      "[[1. 2. 3. 4.]\n",
      " [5. 6. 7. 8.]] <class 'numpy.ndarray'>\n",
      "2 (2, 4)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# have a list of values\n",
    "vl=[1.23,2.23,3.23,4.23]\n",
    "\n",
    "# and make it to a numpy array\n",
    "v=np.array([vl])\n",
    "\n",
    "print(vl,type(vl))\n",
    "print(v,type(v))\n",
    "print(v.ndim,v.shape)\n",
    "\n",
    "# create a matrix\n",
    "# see the value 1.0 instead of 1\n",
    "# what is the difference ? \n",
    "A=np.array( [[1.0,2,3,4],[5,6,7,8]] )\n",
    "print(A,type(A))\n",
    "print(A.ndim,A.shape)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we look at the matrix - vector multiplication, there is a pitfall for numpy-beginners, if you have a matrix A and a vector v and write ```b=A*v``` you get **not** what you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.23  4.46  9.69 16.92]\n",
      " [ 6.15 13.38 22.61 33.84]] <class 'numpy.ndarray'> 2 (2, 4)\n"
     ]
    }
   ],
   "source": [
    "# matrix-vector multiplication -- wrong way \n",
    "# the * is a simple array multiplication and NOT the matrix-multiplication\n",
    "b=A*v\n",
    "print(b,type(b),b.ndim,b.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next try, instead of ```b=A*v``` we do \n",
    "  * make v from a row vector to the correct column vector -- try out what happens if you skip this\n",
    "  * use ```b=A.dot(u)``` or ```b=A@u``` for the matrix vector product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.23]\n",
      " [2.23]\n",
      " [3.23]\n",
      " [4.23]] <class 'numpy.ndarray'> 2 (4, 1)\n",
      "[[32.3 ]\n",
      " [75.98]] <class 'numpy.ndarray'> 2 (2, 1)\n",
      "[[32.3 ]\n",
      " [75.98]] <class 'numpy.ndarray'> 2 (2, 1)\n"
     ]
    }
   ],
   "source": [
    "# matrix-vector multiplication again\n",
    "\n",
    "# make v from a row vector to a column vector\n",
    "# x.T gives in numpy the transposed of x\n",
    "u=v.T\n",
    "print(u,type(u),u.ndim,u.shape)\n",
    "\n",
    "b=A.dot(u)\n",
    "print(b,type(b),b.ndim,b.shape)\n",
    "\n",
    "#you can also use the @ operator for matrix multiplication\n",
    "c=A@u\n",
    "print(c,type(c),c.ndim,c.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task\n",
    "\n",
    "Play around with the next code blocks to understand the element access and data types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 2. 3. 4.]\n",
      " [5. 6. 7. 8.]]\n",
      "1.0\n",
      "7.0\n",
      "[[ 1.    2.    3.    4.  ]\n",
      " [ 5.    6.    7.   -3.14]]\n"
     ]
    }
   ],
   "source": [
    "# access to elements \n",
    "print(A)\n",
    "print(A[0,0])\n",
    "print(A[1,2])\n",
    "\n",
    "# try the next with 1.0 and 1 in the definition of A above \n",
    "# what happens ? \n",
    "A[1,3]=-3.14\n",
    "print(A)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 2 3 4]\n",
      "[1. 2. 3. 4.]\n",
      "[1 2 3 4]\n",
      "[1.5 2.  3.  4. ]\n"
     ]
    }
   ],
   "source": [
    "x=np.array([1,2,3,4])\n",
    "print(x)\n",
    "y=np.array([1.0,2.0,3.0,4.0])\n",
    "print(y)\n",
    "\n",
    "x[0]=1.5\n",
    "y[0]=1.5\n",
    "\n",
    "print(x)\n",
    "print(y)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task\n",
    "\n",
    "Write a function mymat(n) which takes a size argument n>0 and returns a regular nxn Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mymat(n):\n",
    "    A=np.zeros((n,n))\n",
    "    # TODO\n",
    "    return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0.]\n",
      " [0. 0. 0.]\n",
      " [0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "print (mymat(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task\n",
    "\n",
    "Now use your function to create for n=3,5,10 and 15 matrices A and right hand sides $$b = A \\cdot [1,1, ... 1]^T$$ and find out \n",
    "how you can use numpy.linalg.solve for solving linear systems Ax=b. Check if you get the correct solution $$[1,1,...1]^T$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "n=3\n",
    "A=mymat(n)\n",
    "# TODO\n",
    "x = 0 # replace with solution of Ax=b\n",
    "print(x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Toeplitz and circulant matrices \n",
    "\n",
    "We want to create a Toeplitz matrix and also a circulant matrix, both incredibly important in signal processing. \n",
    "\n",
    "Toeplitz matrices are constant along the diagonals from left to rigtht, the structure is\n",
    "$$T=\\left[ \\begin{matrix} a & b & c & d & e \\\\ f & a &  b &  c & d  \\\\ g & f & a &  b &  c \\\\ h& g & f & a &  b   \n",
    " \\\\ i & h & g & f & a  \\end{matrix} \\right] $$\n",
    "\n",
    "As you can see - a Toeplitz matrix is fully defined by setting the first row and the first column.\n",
    "\n",
    "# Task\n",
    "\n",
    "Write a function which takes 2 vectors **c** ( the first column )  and **r** (the first row) \n",
    "as numpy arrays and returns the Toeplitz matrix, defined by this vectors. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def toeplitz(c,r):\n",
    "    if(r[0] != c[0] ):\n",
    "        raise ValueError('Error c[0] must be equal r[0] but  %f %f ' % (c[0],r[0]))\n",
    "    \n",
    "    # TODO\n",
    "    T=0 # replace this\n",
    "    return T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    }
   ],
   "source": [
    "\n",
    "r=np.array([1,2,3,4])\n",
    "c=np.array([1,5,6,7])\n",
    "T=toeplitz(c,r)\n",
    "print(T)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Circulant matrices \n",
    "\n",
    "Circulant matrices are square matrices and fully defined by its first column **c** and have the structure\n",
    "\n",
    "$$C=\\left[ \\begin{matrix} \n",
    "a & e & d & c & b \\\\ \n",
    "b & a & e & d & c  \\\\ \n",
    "c & b & a & e & d \\\\ \n",
    "d & c & b & a & e    \\\\ \n",
    "e & d & c & b & a  \\end{matrix} \\right] $$\n",
    "\n",
    "You can see that Circulant matrices are a special kind of Toeplitz matrices.\n",
    "\n",
    "## Task\n",
    "\n",
    "Write a function which takes a vector **c** as 1 numpy array \n",
    "and returns the Circulant matrix, defined by this vector. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def circulant(c):\n",
    "    # TODO\n",
    "    C=np.zeros((4,4)) # remove this \n",
    "    return C\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "C=circulant(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0. 0.]\n",
      " [0. 0. 0. 0.]\n",
      " [0. 0. 0. 0.]\n",
      " [0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "print(C)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourier matrix\n",
    "\n",
    "The Fourier matrix $F$ of size $n \\times n$ is defined as\n",
    "$$\n",
    "F=\\frac{1}{\\sqrt{n}}\n",
    "\\left[\n",
    "\\begin{array}\n",
    "\\cdot & \\cdot & \\cdot \\\\\n",
    "\\cdot & \\omega^{jk} & \\cdot \\\\\n",
    "\\cdot & \\cdot & \\cdot \\\\\n",
    "\\end{array}\n",
    "\\right]_{j,k=0}^{n-1}\n",
    "$$\n",
    "where $\\omega$ is the n-th primitive root of unity, i.e. $$\\omega=e^{-2\\Pi i / n}$$\n",
    "\n",
    "## Task\n",
    "\n",
    "Write a function which computes for given size **n** the Fourier matrix $F$.\n",
    "\n",
    "(Later we will use scipy - there is the Fouriermatrix given as build in, but now we do it by hand)\n",
    "  * Hint 1 : The complex number **i** you can get in numpy as **1.j**\n",
    "  * Hint 2 : $\\Pi$ you get as np.pi\n",
    "  * Hint 3 : in complex computations initialize your data types correct, f.e. \n",
    "    ``` F=np.zeros((n,n),dtype=complex) ``` . Try out what happens if you do it as \n",
    "    ``` F=np.zeros((n,n)) ```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fourier(n):\n",
    "    # create a complex zero matrix\n",
    "    F=np.zeros((n,n),dtype=complex)\n",
    "    # TODO\n",
    "\n",
    "    return F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "F=fourier(4)\n",
    "print(F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diagonalization of a circulant matrix\n",
    "\n",
    "With the help of our Fouriermatrix $F$  we can diagonlize the circulant matrix $C$ with\n",
    "$$ D = \\bar{F}C F $$\n",
    "where $\\bar{F}$ is the complex conjugate matrix of $F$, let's try this out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
      " [0.+0.j 0.+0.j 0.+0.j 0.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "# matrix multiplication with the @ operator\n",
    "D=np.conjugate(F)@C@F\n",
    "print(D)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result $D$ is not very good readable because of the zeros written as small numbers in range $10^{-16}$. So we do some matrix cosmetic and cut off small numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0. 0.]\n",
      " [0. 0. 0. 0.]\n",
      " [0. 0. 0. 0.]\n",
      " [0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "tol=1e-12\n",
    "D.real[abs(D.real) < tol] = 0.0\n",
    "D.imag[abs(D.imag) < tol] = 0.0\n",
    "\n",
    "print(D.real)\n",
    "#print(D.imag)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now $D$ looks like a diagonal matrix :-) \n",
    "\n",
    "The construction of $F$ and $C$ can be avoided if the **Fast Fourier Transfor (FFT)** is used and the complexity reduces to something that is $\\mathcal{O}(n\\mathrm{log}(n))$ rather than $\\mathcal{O}(n^2)$\n",
    "\n",
    "[here](https://math.mit.edu/icg/resources/teaching/18.085-spring2015/toeplitz.pdf) you can find the information of how to embedd a Toeplitz matrix into a circulant matrix. Implement this function and verify with a DFT of the double size that this is diagonalized using $F$.\n",
    "\n",
    "\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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}