C1W2_3DIVIDE_N_CONQUER_ALGORITHMS_4THE_MASTER_METHOD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Divide and Conquer algorithms\n",
    "In this notebook we will code and analyse the following three divide and conquer algorithms.\n",
    "\n",
    "1. Counting the number of inversions\n",
    "2. Matrix multiplication using Strassen's algorithm\n",
    "3. Computing the closest points together.\n",
    "\n",
    "## Counting Inversions\n",
    "Counting the number of inversions is a measure of how similar (or not similar) are two arrays of numbers. This algorithm is used in recommender systems where the system recommends something (say a movie) to the consumer based on how similar their rating is to other people who rated similarly.\n",
    "\n",
    "Collaborative Filtering is a technique where we try to find similarities between likes of people. Suppose we want to find similarities of the movie rankings of two people, we will sort movies by ranking of one person's preference and then put ranks given by the second person for those movies.\n",
    "\n",
    "For example if your ranked 3 movies A, B and C as your favorite movies and your friend ranked their favorite movies as B A C then between you and the your friend the ranking list becomes 2 1 3 because A is your favorite movie but is second favorite for B. Similarly B is your second favorite movie but your friend's favorite. Both you and your friend's choice match for the third best movie.\n",
    "\n",
    "From the above example we see that if two people have similar tastes, both the arrays would be identical and there would be no inversions. Higher the number of inversions more are the differences in the preferences.\n",
    "\n",
    "We will use divide and conquer algorithm like Merge sort to find the number of inversions in an array. The number of inversions in an array A are the number of pairs of indices $(i, j)$ where $i < j$ and $A[i] < A[j]$\n",
    "\n",
    "An array that is sorted has no inversions, the converse is also true, that is, an array with no inversions is sorted and is not sorted if it has atleast one inversion.\n",
    "\n",
    "For example, consider the following array\n",
    "\n",
    "1 3 5 2 4 6\n",
    "\n",
    "The number of inversions in this array are the following pair of numbers (3, 2), (5, 2) and (5, 4)\n",
    "\n",
    "For an array of n numbers, we can have a maxumum of (n - 1) + (n - 2) + .. + 1 number of inversions = $n(n - 1) / 2$ number of inversions.\n",
    "\n",
    "Let us implement counting inversion in two different ways. First one is the Bruteforce approach and the second one is using divide and conquer approach that sorts the numbers using Merge Sort along with counting the number of inversions.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Brute force approach for Counting Inversions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "3\n",
      "15\n",
      "3\n",
      "28\n"
     ]
    }
   ],
   "source": [
    "def count_inversions_brute_force(arr):\n",
    "    count_inv = 0\n",
    "    for i in range(len(arr)-1):\n",
    "        for j in range(i+1, len(arr)):\n",
    "            if arr[i] > arr[j]:\n",
    "                count_inv +=1\n",
    "    print(count_inv)\n",
    "\n",
    "count_inversions_brute_force([1,2,3,4,5,6,7,8,9])\n",
    "count_inversions_brute_force([1,2,3,6,5,4,7,8,9])\n",
    "count_inversions_brute_force([6,5,4,3,2,1])\n",
    "count_inversions_brute_force([1, 3, 5, 2, 4, 6])\n",
    "count_inversions_brute_force([ 8, 7, 6, 5, 4, 3 ,2, 1 ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we will now load the two test case files provided <a href=http://algorithmsilluminated.org/ > here</a> and test our implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 315,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28\n"
     ]
    }
   ],
   "source": [
    "import urllib3\n",
    "http = urllib3.PoolManager()\n",
    "\n",
    "# Test case\n",
    "r1 = http.request('GET', \"http://algorithmsilluminated.org/datasets/problem3.5test.txt\")\n",
    "problem35teststr = r1.data.split('\\r\\n')\n",
    "del problem35teststr[-1]\n",
    "problem35test = [int(i) for i in problem35teststr]\n",
    "count_inversions_brute_force(problem35test)\n",
    "\n",
    "# Challenge data set\n",
    "###  Brute force approach not work  efficiently on larger problems for counting the inversions.\n",
    "#r2 = http.request('GET', \"http://algorithmsilluminated.org/datasets/problem3.5.txt\")\n",
    "#problem35str = r2.data.split('\\r\\n')\n",
    "#del problem35str[-1]\n",
    "#problem35 = [int(i) for i in problem35str]\n",
    "#count_inversions_brute_force(problem35)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Divide and Conquer approach for Counting Inversions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The brute force approach seems to be working fine for small inputs, it will however not work for efficiently on the larger 100000 numbers for counting the inversions. We will now implement the inversion counting piggy backed on merge sort."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 457,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n",
      "25\n",
      "0\n",
      "3\n",
      "15\n",
      "3\n",
      "28\n"
     ]
    }
   ],
   "source": [
    "def sort_and_count_inversions(arr, count_inversions):\n",
    "    n = len(arr)\n",
    "    if n == 1:\n",
    "        return arr, count_inversions\n",
    "\n",
    "    if n >= 1:\n",
    "        left_tree, count_inversions = sort_and_count_inversions(arr[ : n/2], count_inversions)   \n",
    "        right_tree, count_inversions = sort_and_count_inversions(arr[n/2 : ], count_inversions)\n",
    "        sorted_arr = [None] * n\n",
    "        i = 0\n",
    "        j = 0\n",
    "        for k in range(n):\n",
    "            if(i < len(left_tree) and j < len(right_tree)):\n",
    "                if left_tree[i] <= right_tree[j]:\n",
    "                    sorted_arr[k] = left_tree[i]\n",
    "                    i+=1 \n",
    "                elif left_tree[i] > right_tree[j]:\n",
    "                    sorted_arr[k] = right_tree[j]\n",
    "                    j+=1 \n",
    "                    count_inversions += len(left_tree) - i\n",
    "                \n",
    "            elif i >= len(left_tree) and j < len(right_tree):\n",
    "                sorted_arr[k] = right_tree[j]\n",
    "                j += 1\n",
    "            elif j >= len(right_tree) and i < len(left_tree):\n",
    "                sorted_arr[k] = left_tree[i]\n",
    "                i += 1\n",
    "    \n",
    "    return  sorted_arr, count_inversions\n",
    "    \n",
    "\n",
    "print(sort_and_count_inversions([4, 3, 2] , count_inversions = 0)[1])    \n",
    "print(sort_and_count_inversions([4, 3, 2, 10, 12, 1, 5, 6, 24, 33, 23,54, 12, 6 ] , count_inversions = 0)[1])    \n",
    "print(sort_and_count_inversions([1,2,3,4,5,6,7,8,9] , count_inversions = 0)[1])\n",
    "print(sort_and_count_inversions([1,2,3,6,5,4,7,8,9] , count_inversions = 0)[1])\n",
    "print(sort_and_count_inversions([6,5,4,3,2,1] , count_inversions = 0)[1])\n",
    "print(sort_and_count_inversions([1, 3, 5, 2, 4, 6] , count_inversions = 0)[1])\n",
    "print(sort_and_count_inversions([ 8, 7, 6, 5, 4, 3 ,2, 1 ] , count_inversions = 0)[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of splits in test array are  28\n",
      "Number of splits in the challenge data set are  2407905288\n"
     ]
    }
   ],
   "source": [
    "import urllib3\n",
    "http = urllib3.PoolManager()\n",
    "\n",
    "# Test case\n",
    "r1 = http.request('GET', \"http://algorithmsilluminated.org/datasets/problem3.5test.txt\")\n",
    "problem35teststr = r1.data.split('\\r\\n')\n",
    "del problem35teststr[-1]\n",
    "problem35test = [int(i) for i in problem35teststr]\n",
    "print(\"Number of splits in test array are  {}\".format(sort_and_count_inversions(problem35test , count_inversions = 0)[1]))\n",
    "\n",
    "# Challenge data set\n",
    "###**** Divide and Conquer approach works efficiently on larger problems for counting the inversions. ********************\n",
    "r2 = http.request('GET', \"http://algorithmsilluminated.org/datasets/problem3.5.txt\")\n",
    "problem35str = r2.data.split('\\r\\n')\n",
    "del problem35str[-1]\n",
    "problem35 = [int(i) for i in problem35str]\n",
    "print(\"Number of splits in the challenge data set are  {}\".format(sort_and_count_inversions(problem35 , count_inversions = 0)[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The divide and conquer method sort_and_count_inversions is trivial. All we do is to count the number of inversions on the left and right and get the corresponding halves sorted. We then find the split inversions and also merge the two sorted arrays. The total inversions are the number of inversions on the left plus the ones on the right plus the number of split inversions.\n",
    "\n",
    "The ingenuity lies in the count_inversions_and_sort function. This function piggy backs on the merge sort function and counts the number of inversions along with sorting the array. The function defined above is pretty straight forward and has comments inline giving explanation.\n",
    "\n",
    "If we have no inversions in a array A and we receive two halves of an array, then the elements in first half are strictly less than the elements in second half.\n",
    "\n",
    "The sort_and_count_inversions splits the input array in two and recursively sorts and counts inversions on the left and right half. The number of tasks at level n doubles than that of level n - 1 and input size given to each of task at level n is half of the input given to a unit on n - 1 level. This is similar to merge sort and given that the routine to sort and count inversions execute in linear time, the count sort_and_count_inversions also runs in $O(nlogn)$ time\n",
    "\n",
    "***\n",
    "\n",
    "#### Strassen's Matrix multiplication Algorithm\n",
    "Stratten's Matrix multiplication algorithm is a clever divide and conquer approach for multiplying two matrix. Before we see the details of matrix multiplication using Strassen's approach, let us implement matrix multiplication using simple, straight forward approach.\n",
    "\n",
    "For simplicity of implementing divide and conquer algorithms, we will assume that the matrix square and each side is a power of 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matmul_simple(A, B):\n",
    "    dim = len(A[0])\n",
    "    Z = [[0 for _ in range(dim)] for _ in range(dim)]\n",
    "    for i in range(dim):\n",
    "        for j in range(dim):\n",
    "            Z[i][j] = sum(A[i][k] * B[k][j] for k in range(dim))\n",
    "    return Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[90, 100, 110, 120],\n",
       " [202, 228, 254, 280],\n",
       " [314, 356, 398, 440],\n",
       " [426, 484, 542, 600]]"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "matmul_simple([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]], \n",
    "                  [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since there are three loops in the above multiplication (with variables i, j and k for each loop counter). The time complexity of a matrix multiplication is $O(n^3)$.\n",
    "\n",
    "Before we go ahead and implement a recursive divide and conquer algorithm, let us verify our result are correct by using a reference implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 316,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 90, 100, 110, 120],\n",
       "        [202, 228, 254, 280],\n",
       "        [314, 356, 398, 440],\n",
       "        [426, 484, 542, 600]])"
      ]
     },
     "execution_count": 316,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "A = np.matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]])\n",
    "A * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The divide and conquer approach will break split an $n \\times n$ matrix into 4 $\\frac{n}{2} \\times \\frac{n}{2}$ matrix and then multiple then recursively. Let us illustrate then as below\n",
    "\n",
    "Suppose the bigger matrix X and Y are split in the following two matrices. Then the resulting matrix multiplication is formed by multiplying and these 8 small matrices. This divide and conquer can continue until we get a matrix of size 1 where the recursive divide and conquer stops. This divide and conquer is summarized in the following picture where A, B, C, D and E, F, G, H are smaller matrices of the bigger matrices X and Y respectively.\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix} \n",
    "A & B \\\\\n",
    "C & D \n",
    "\\end{pmatrix}\n",
    "\\times\n",
    "\\begin{pmatrix} \n",
    "E & F \\\\\n",
    "G & H \n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix} \n",
    "AE + BG  & AF + BH \\\\\n",
    "CE + DG & CF + DH \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[176, 172, 168, 164, 160, 156, 152, 148], [148, 152, 156, 160, 164, 168, 172, 176], [176, 172, 168, 164, 160, 156, 152, 148], [148, 152, 156, 160, 164, 168, 172, 176], [176, 172, 168, 164, 160, 156, 152, 148], [148, 152, 156, 160, 164, 168, 172, 176], [176, 172, 168, 164, 160, 156, 152, 148], [148, 152, 156, 160, 164, 168, 172, 176]]\n"
     ]
    }
   ],
   "source": [
    "def slicemat(mat, r_from, r_to, c_from, c_to):\n",
    "    return [[mat[r][c] for c in range(c_from, c_to)] for r in range(r_from, r_to)]\n",
    "\n",
    "def hcat(U, V):\n",
    "    #print(\"U: {}\\n V: {}\".format(U,V))\n",
    "    return [U[0] + V[0]] if len(U) == 1 else [U[r] + V[r] for r in range(len(U))]\n",
    "\n",
    "def add_matrices(W, Z):\n",
    "    d = len(W)\n",
    "    #print(\"W: {}\\n Z: {}\".format(W,Z))\n",
    "    return [[W[i][j] + Z[i][j] for j in range(d)] for i in range(d)]\n",
    "\n",
    "    \n",
    "def matmul_recursively_8(X, Y):\n",
    "    #print(\"X: {}\\nY: {}\".format(X, Y))\n",
    "    dim = len(X[0])\n",
    "    if dim == 1:\n",
    "        return [[X[0][0] * Y[0][0]]]\n",
    "    else:\n",
    "        A = slicemat(X, 0, dim//2, 0, dim//2)\n",
    "        B = slicemat(X, 0, dim//2, dim//2, dim)\n",
    "        C = slicemat(X, dim//2, dim, 0, dim//2)\n",
    "        D = slicemat(X, dim//2, dim, dim//2, dim)\n",
    "        E = slicemat(Y, 0, dim//2, 0, dim//2)\n",
    "        F = slicemat(Y, 0, dim//2, dim//2, dim)\n",
    "        G = slicemat(Y, dim//2, dim, 0, dim//2)\n",
    "        H = slicemat(Y, dim//2, dim, dim//2, dim)\n",
    "        AE = matmul_recursively_8(A, E)\n",
    "        BG = matmul_recursively_8(B, G)\n",
    "        AF = matmul_recursively_8(A, F)\n",
    "        BH = matmul_recursively_8(B, H)\n",
    "        CE = matmul_recursively_8(C, E)\n",
    "        DG = matmul_recursively_8(D, G)\n",
    "        CF = matmul_recursively_8(C, F)\n",
    "        DH = matmul_recursively_8(D, H)\n",
    "        #print(\"AE: {}\\nBG: {}\".format(AE, BG))\n",
    "        return hcat(add_matrices(AE , BG) , add_matrices(AF , BH)) + hcat(add_matrices(CE , DG) , add_matrices(CF , DH))\n",
    "\n",
    "#X = [[1, 2], [8, 7]]\n",
    "#X = [[1, 2, 3, 4], [8, 7, 6, 4], [1, 2, 3, 4], [8, 7, 6, 4]]\n",
    "#X = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]\n",
    "X = [[1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1]]\n",
    "print(matmul_recursively_8(X, X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16], [1, 2, 3, 4]]\n",
      "4\n",
      "[]\n",
      "[1, 2]\n",
      "[[1, 2], [5, 6]]\n",
      "1\n",
      "2\n",
      "3\n"
     ]
    }
   ],
   "source": [
    "# Tidbits\n",
    "\n",
    "#List not slicable\n",
    "T = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]\n",
    "U = [[1, 2, 3, 4]]\n",
    "\n",
    "print(T+U)   #  +  --> cancels '][' and replaces it with ','\n",
    "print(len(T[0]))\n",
    "print(T[0:(d/2)][(d/2):])\n",
    "print(T[0:d/2][0][:d/2])\n",
    "print([T[0:d/2][0][:d/2]]+[T[0:d/2][1][:d/2]])\n",
    "print(len([[[5]]]))\n",
    "print(len([[3, 2], [4, 5]]))  # same dimension\n",
    "print(len([[3, 2], [4, 5], [5, 7]]))  # max dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[176, 172, 168, 164, 160, 156, 152, 148],\n",
       "        [148, 152, 156, 160, 164, 168, 172, 176],\n",
       "        [176, 172, 168, 164, 160, 156, 152, 148],\n",
       "        [148, 152, 156, 160, 164, 168, 172, 176],\n",
       "        [176, 172, 168, 164, 160, 156, 152, 148],\n",
       "        [148, 152, 156, 160, 164, 168, 172, 176],\n",
       "        [176, 172, 168, 164, 160, 156, 152, 148],\n",
       "        [148, 152, 156, 160, 164, 168, 172, 176]])"
      ]
     },
     "execution_count": 183,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.matrix([[1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1]])\n",
    "A * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scary looking function shown above is implementation of the same algorithm described earlier. We split the input array into 8 smaller arrays and perform these 8 matrix multiplication recursively before we combine the results into a bigger martix. Without looking at how we compute the time compexity of the algorithm it still is $O(n^3)$, no different than the simple approach tried earlier. The derivation of this time complexity will be done later when we study the Master method.\n",
    "\n",
    "***\n",
    "\n",
    "### Strassen's Subcubic Matrix Multiplication Algorithm\n",
    "\n",
    "The recursive approach above looks similar to the divide and conquer approach we saw for multiplication of two numbers. A clever trick used in Karatsuba multiplication reduced the number of recursive calls thus reducing the time complexity below $n^2$. Can we do something similar? The Strassen's approach reduces the recursive calls from 8 to 7 thus making the algorithm more efficient. How these 7 multiplications to be performed were discovered are not known and are as follows.\n",
    "\n",
    "$$ P_1 = A\\cdot (F - H)\\\\ P_2 = (A + B)\\cdot H\\\\ P_3 = (C + D)\\cdot E\\\\ P_4 = D\\cdot (G - E)\\\\ P_5 = (A + D)\\cdot (E + H)\\\\ P_6 = (B - D)\\cdot (G + H)\\\\ P_7 = (A - C)\\cdot (E + F)\\\\ $$$$ X\\cdot T = \\begin{pmatrix} A & B \\\\ C & D \\end{pmatrix} \\times \\begin{pmatrix} E & F \\\\ G & H \\end{pmatrix} = \\begin{pmatrix} AE + BG & AF + BH \\\\ CE + DG & CF + DH \\end{pmatrix} = \\begin{pmatrix} P_5 + P_4 - P_2 + P_6 & P_1 + P_2\\\\ P_3 + P_4 & P_1 + P_5 - P_3 - P_7 \\end{pmatrix} $$\n",
    "\n",
    "<br>\n",
    "Following code snippet is an implementation of strassen's matrix multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[176, 172, 168, 164, 160, 156, 152, 148], [148, 152, 156, 160, 164, 168, 172, 176], [176, 172, 168, 164, 160, 156, 152, 148], [148, 152, 156, 160, 164, 168, 172, 176], [176, 172, 168, 164, 160, 156, 152, 148], [148, 152, 156, 160, 164, 168, 172, 176], [176, 172, 168, 164, 160, 156, 152, 148], [148, 152, 156, 160, 164, 168, 172, 176]]\n"
     ]
    }
   ],
   "source": [
    "def sub_matrices(W, Z):\n",
    "    d = len(W)\n",
    "    #print(\"W: {}\\n Z: {}\".format(W,Z))\n",
    "    return [[W[i][j] - Z[i][j] for j in range(d)] for i in range(d)]\n",
    "\n",
    "\n",
    "def Strassen_Matrix_Multiplication(X, Y):\n",
    "    #print(\"X: {}\\nY: {}\".format(X, Y))\n",
    "    dim = len(X[0])\n",
    "    if dim == 1:\n",
    "        return [[X[0][0] * Y[0][0]]]\n",
    "    else:\n",
    "        A = slicemat(X, 0, dim//2, 0, dim//2)\n",
    "        B = slicemat(X, 0, dim//2, dim//2, dim)\n",
    "        C = slicemat(X, dim//2, dim, 0, dim//2)\n",
    "        D = slicemat(X, dim//2, dim, dim//2, dim)\n",
    "        E = slicemat(Y, 0, dim//2, 0, dim//2)\n",
    "        F = slicemat(Y, 0, dim//2, dim//2, dim)\n",
    "        G = slicemat(Y, dim//2, dim, 0, dim//2)\n",
    "        H = slicemat(Y, dim//2, dim, dim//2, dim)\n",
    "        P1 = Strassen_Matrix_Multiplication(A, sub_matrices(F, H))\n",
    "        P2 = Strassen_Matrix_Multiplication(add_matrices(A, B), H)\n",
    "        P3 = Strassen_Matrix_Multiplication(add_matrices(C, D), E)\n",
    "        P4 = Strassen_Matrix_Multiplication(D, sub_matrices(G, E))\n",
    "        P5 = Strassen_Matrix_Multiplication(add_matrices(A, D), add_matrices(E, H))\n",
    "        P6 = Strassen_Matrix_Multiplication(sub_matrices(B, D), add_matrices(G, H))\n",
    "        P7 = Strassen_Matrix_Multiplication(sub_matrices(A, C), add_matrices(E, F))\n",
    "        #print(\"AE: {}\\nBG: {}\".format(AE, BG))\n",
    "        return hcat(add_matrices(sub_matrices(add_matrices(P5 , P4), P2), P6) , add_matrices(P1 , P2)) + hcat(add_matrices(P3 , P4) , sub_matrices(sub_matrices(add_matrices(P1 , P5), P3), P7))\n",
    "\n",
    "#X = [[1, 2], [8, 7]]\n",
    "#X = [[1, 2, 3, 4], [8, 7, 6, 4], [1, 2, 3, 4], [8, 7, 6, 4]]\n",
    "#X = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]\n",
    "X = [[1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1], [1, 2, 3, 4, 5, 6, 7, 8], [8, 7, 6, 5, 4, 3, 2, 1]]\n",
    "print(matmul_recursively_8(X, X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time complexity analysis of Strassen's method will be done later when we explore the Master method later.\n",
    "\n",
    "***\n",
    "\n",
    "### Find the closest pair\n",
    "The final divide and conquer algorithm would be to find the closest pairs in $O(nlogn)$ time. For the problem we are given n points in a plane and the goal is to find two points with minimum euclidean distance between them.\n",
    "\n",
    "For two points $p_1 = (x_1, y_1)$ and $p_2 = (x_2, y_2)$ the euclidean distance $d(p_1, p_2)$ can be found using $\\sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}$\n",
    "\n",
    "Starting with the Brute force approach, we know that the Naive algorithm will compute the distance of all points with all other points and then finds two with minimum distance between them. This algorithm has a time complexity of $O(n^2)$. Quadratic time is certainly not the best algorithm in this case and we should be looking at something better.\n",
    "\n",
    "Let us start with points in 1 dimension. The naive approach still gives us the closest points in 1-D in $O(n^2)$ time. However, just by sorting these points in $O(nlogn)$ time using something like Merge Sort and then doing a linear scan in linear time to find the a pair of points closest to each other is easily achievable. Thus we have a better algorithm than in 1-D then quadratic time to find the closest pair.\n",
    "\n",
    "Let us now see how to find the closest points in 2-D. Let us consider the following points, (0, 0), (1, 1), (0.5, 4), (4, 0.5) and (-2, 2) and plot them to get a visual on them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAEY1JREFUeJzt3XGIZWd5x/HvbzdbdbA1YAYM2Z2dgqHUisb2EiL5J0QLqw0JbSNEtlZby0BRGqnQ1i5YDOwfUtDSpjRMGzG2F41E265WsREVW6jRu+kmNV0tS+luFgOZGk0M06asPv3j3rS7s7N7z525s3fvm+8HLuec9zxzznNI8tuTM+fum6pCktSWXbNuQJI0fYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUFXzOrEV111VS0vL8/q9JI0l44ePfqfVbU4rm5m4b68vMxgMJjV6SVpLiU52aXOxzKS1CDDXZIaZLhLUoMMd0lqkOEuSQ3qHO5Jdif55ySf3WTfi5Lcn+REkoeSLE+zSYl+H5aXYdeu4bLfn3VH0mVtkjv3O4HjF9j3TuB7VfVK4MPAB7fbmPR/+n1YWYGTJ6FquFxZMeCli+gU7kn2Ar8A/MUFSm4D7hutPwC8IUm2354EHDoE6+vnjq2vD8clbarrnfsfAb8D/OgC+68BHgeoqjPA08DLNxYlWUkySDJYW1vbQrt6QTp1arJxSePDPcktwJNVdfRiZZuMnTfzdlWtVlWvqnqLi2O/PSsNLS1NNi6p0537jcCtSf4D+ARwc5K/2lBzGtgHkOQK4GXAU1PsUy9khw/DwsK5YwsLw3FJmxob7lX1vqraW1XLwB3Al6rqVzaUHQHePlq/fVRz3p27tCUHD8LqKuzfD8lwubo6HJe0qS3/xWFJ7gIGVXUEuBf4yyQnGN6x3zGl/qShgwcNc2kCE4V7VX0F+Mpo/f1njf838JZpNiZJ2jq/oSpJDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDusyh+uIkX0/ySJLHknxgk5p3JFlLcmz0+Y2daVeS1EWXyTqeA26uqmeT7AH+Mcnnq+prG+rur6p3T79FSdKkxob7aC7UZ0ebe0Yf50eVpMtYp2fuSXYnOQY8CTxYVQ9tUvbLSR5N8kCSfVPtUpI0kU7hXlU/rKrrgL3A9UlevaHkM8ByVb0G+CJw32bHSbKSZJBksLa2tp2+JUkXMdHbMlX1fYYTZB/YMP7dqnputPnnwM9d4OdXq6pXVb3FxcUttCtJ6qLL2zKLSa4crb8EeCPwrQ01V5+1eStwfJpNSpIm0+VtmauB+5LsZviHwSer6rNJ7gIGVXUE+K0ktwJngKeAd+xUw5Kk8TJ8GebS6/V6NRgMZnJuSZpXSY5WVW9cnd9QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1qMs0ey9O8vUkjyR5LMkHNql5UZL7k5xI8lCS5Z1oVpLUTZc79+eAm6vqtcB1wIEkN2yoeSfwvap6JfBh4IPTbVOSNImx4V5Dz44294w+G+fmuw24b7T+APCGJJlal5KkiXR65p5kd5JjwJPAg1X10IaSa4DHAarqDPA08PJpNipJ6q5TuFfVD6vqOmAvcH2SV28o2ewu/byZt5OsJBkkGaytrU3erSSpk4nelqmq7wNfAQ5s2HUa2AeQ5ArgZcBTm/z8alX1qqq3uLi4pYYlSeN1eVtmMcmVo/WXAG8EvrWh7Ajw9tH67cCXquq8O3dJ0qVxRYeaq4H7kuxm+IfBJ6vqs0nuAgZVdQS4F/jLJCcY3rHfsWMdS5LGGhvuVfUo8LpNxt9/1vp/A2+ZbmuSpK3yG6qS1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAZ1mWZvX5IvJzme5LEkd25Sc1OSp5McG33ev9mxJEmXRpdp9s4A762qh5P8OHA0yYNV9a8b6v6hqm6ZfouSpEmNvXOvqieq6uHR+g+A48A1O92YJGnrJnrmnmSZ4XyqD22y+/VJHkny+SQ/c4GfX0kySDJYW1ubuFlJUjedwz3JS4FPAe+pqmc27H4Y2F9VrwX+BPibzY5RVatV1auq3uLi4lZ7liSN0Snck+xhGOz9qvr0xv1V9UxVPTta/xywJ8lVU+1UktRZl7dlAtwLHK+qD12g5hWjOpJcPzrud6fZqCSpuy5vy9wIvA34lyTHRmO/DywBVNU9wO3AbyY5A/wXcEdV1Q70K0nqYGy4V9U/AhlTczdw97SakiRtj99QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUFdZmLal+TLSY4neSzJnZvUJMkfJzmR5NEkP7sz7UqSuuhy534GeG9V/TRwA/CuJK/aUPMm4NrRZwX4s6l2+bx+H5aXYdeu4bLf35HTSNK8GxvuVfVEVT08Wv8BcBy4ZkPZbcDHauhrwJVJrp5qp/0+rKzAyZNQNVyurBjwkrSJiZ65J1kGXgc8tGHXNcDjZ22f5vw/ALbn0CFYXz93bH19OC5JOkfncE/yUuBTwHuq6pmNuzf5kfMmyE6ykmSQZLC2tjZZp6dOTTYuSS9gncI9yR6Gwd6vqk9vUnIa2HfW9l7gOxuLqmq1qnpV1VtcXJys06WlycYl6QWsy9syAe4FjlfVhy5QdgT41dFbMzcAT1fVE1PsEw4fhoWFc8cWFobjkqRzXNGh5kbgbcC/JDk2Gvt9YAmgqu4BPge8GTgBrAO/NvVODx4cLg8dGj6KWVoaBvvz45Kk/5Oq8x6NXxK9Xq8Gg8FMzi1J8yrJ0arqjavzG6qS1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAZ1mWbvI0meTPLNC+y/KcnTSY6NPu+ffpuSpEl0mWbvo8DdwMcuUvMPVXXLVDqSJG3b2Dv3qvoq8NQl6EWSNCXTeub++iSPJPl8kp+Z0jElSVvU5bHMOA8D+6vq2SRvBv4GuHazwiQrwArA0tLSFE4tSdrMtu/cq+qZqnp2tP45YE+Sqy5Qu1pVvarqLS4ubvfUkqQL2Ha4J3lFkozWrx8d87vbPa4kaevGPpZJ8nHgJuCqJKeBPwD2AFTVPcDtwG8mOQP8F3BHVdWOdSxJGmtsuFfVW8fsv5vhq5KSpMuE31CVpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDVobLgn+UiSJ5N88wL7k+SPk5xI8miSn51+m9Il0u/D8jLs2jVc9vuz7kjaki537h8FDlxk/5uAa0efFeDPtt+WNAP9PqyswMmTUDVcrqwY8JpLY8O9qr4KPHWRktuAj9XQ14Ark1w9rQalS+bQIVhfP3dsfX04Ls2ZaTxzvwZ4/Kzt06Ox8yRZSTJIMlhbW5vCqaUpOnVqsnHpMjaNcM8mY7VZYVWtVlWvqnqLi4tTOLU0RUtLk41Ll7FphPtpYN9Z23uB70zhuNKldfgwLCycO7awMByX5sw0wv0I8Kujt2ZuAJ6uqiemcFzp0jp4EFZXYf9+SIbL1dXhuDRnrhhXkOTjwE3AVUlOA38A7AGoqnuAzwFvBk4A68Cv7VSz0o47eNAwVxPGhntVvXXM/gLeNbWOJEnb5jdUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkN6hTuSQ4k+XaSE0l+b5P970iyluTY6PMb029VkuZcvw/Ly7Br13DZ7+/YqbpMs7cb+FPg5xlOhv2NJEeq6l83lN5fVe/egR4laf71+7CyAuvrw+2TJ4fbsCNTO3a5c78eOFFV/15V/wN8Arht6p1IUssOHfr/YH/e+vpwfAd0CfdrgMfP2j49Gtvol5M8muSBJPs2O1CSlSSDJIO1tbUttCtJc+rUqcnGt6lLuGeTsdqw/RlguapeA3wRuG+zA1XValX1qqq3uLg4WaeSNM+WliYb36Yu4X4aOPtOfC/wnbMLquq7VfXcaPPPgZ+bTnuS1IjDh2Fh4dyxhYXh+A7oEu7fAK5N8pNJfgy4AzhydkGSq8/avBU4Pr0WJakBBw/C6irs3w/JcLm6uiO/TIUOb8tU1Zkk7wa+AOwGPlJVjyW5CxhU1RHgt5LcCpwBngLesSPdStI8O3hwx8J8o1RtfHx+afR6vRoMBjM5tyTNqyRHq6o3rs5vqEpSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktSgTuGe5ECSbyc5keT3Ntn/oiT3j/Y/lGR52o1KkrobG+5JdgN/CrwJeBXw1iSv2lD2TuB7VfVK4MPAB6fdqC5j/T4sL8OuXcNlvz/rjqQXvC537tcDJ6rq36vqf4BPALdtqLkNuG+0/gDwhiSZXpu6bPX7sLICJ09C1XC5smLASzPWJdyvAR4/a/v0aGzTmqo6AzwNvHwaDeoyd+gQrK+fO7a+PhyXNDNdwn2zO/CNE692qSHJSpJBksHa2lqX/nS5O3VqsnFJl0SXcD8N7Dtrey/wnQvVJLkCeBnw1MYDVdVqVfWqqre4uLi1jnV5WVqabFzSJdEl3L8BXJvkJ5P8GHAHcGRDzRHg7aP124EvVdV5d+5q0OHDsLBw7tjCwnBc0syMDffRM/R3A18AjgOfrKrHktyV5NZR2b3Ay5OcAH4bOO91STXq4EFYXYX9+yEZLldXh+OSZiazusHu9Xo1GAxmcm5JmldJjlZVb1yd31CVpAYZ7pLUIMNdkhpkuEtSgwx3SWrQzN6WSbIGnNzij18F/OcU25klr+Xy1Mq1tHId4LU8b39Vjf0W6MzCfTuSDLq8CjQPvJbLUyvX0sp1gNcyKR/LSFKDDHdJatC8hvvqrBuYIq/l8tTKtbRyHeC1TGQun7lLki5uXu/cJUkXMbfhnuQPk3wryaNJ/jrJlbPuaauSvCXJY0l+lGTu3gYYN4H6PEnykSRPJvnmrHvZjiT7knw5yfHRv1t3zrqnrUry4iRfT/LI6Fo+MOuetiPJ7iT/nOSzO3meuQ134EHg1VX1GuDfgPfNuJ/t+CbwS8BXZ93IpDpOoD5PPgocmHUTU3AGeG9V/TRwA/CuOf7n8hxwc1W9FrgOOJDkhhn3tB13Mvzr03fU3IZ7Vf396O+aB/gawxmi5lJVHa+qb8+6jy3qMoH63Kiqr7LJLGLzpqqeqKqHR+s/YBgmG+c+ngs19Oxoc8/oM5e/LEyyF/gF4C92+lxzG+4b/Drw+Vk38QLVZQJ1zVCSZeB1wEOz7WTrRo8yjgFPAg9W1bxeyx8BvwP8aKdPdMVOn2A7knwReMUmuw5V1d+Oag4x/F/Q/qXsbVJdrmVOdZocXbOR5KXAp4D3VNUzs+5nq6rqh8B1o9+t/XWSV1fVXP1eJMktwJNVdTTJTTt9vss63KvqjRfbn+TtwC3AGy73OVvHXcsc6zKBumYgyR6Gwd6vqk/Pup9pqKrvJ/kKw9+LzFW4AzcCtyZ5M/Bi4CeS/FVV/cpOnGxuH8skOQD8LnBrVa3Pup8XsC4TqOsSSxKGcxsfr6oPzbqf7Uiy+PzbcEleArwR+NZsu5pcVb2vqvZW1TLD/06+tFPBDnMc7sDdwI8DDyY5luSeWTe0VUl+Mclp4PXA3yX5wqx76upCE6jPtqutS/Jx4J+An0pyOsk7Z93TFt0IvA24efTfx7HRHeM8uhr4cpJHGd5MPFhVO/oaYQv8hqokNWie79wlSRdguEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KD/BUrp6PuANc/TAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bdda750>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "points = [(0, 0), (1, 1), (0.5, 4), (4, 0.5), (-2, 2)] \n",
    "plt.plot([x[0] for x in points], [x[1] for x in points], 'ro')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visually the points (0, 0) and (1, 1) appear to be closest. Just like in case of 1-D let us sort these points by x axis and then by y axis. Thus we get a couple of copies of the original array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('px:', [(-2, 2), (0, 0), (0.5, 4), (1, 1), (4, 0.5)], 'py:', [(0, 0), (4, 0.5), (1, 1), (-2, 2), (0.5, 4)])\n"
     ]
    }
   ],
   "source": [
    "px = sorted(points, key = lambda x: x[0])\n",
    "py = sorted(points, key = lambda x: x[1])\n",
    "print('px:', px, 'py:', py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the above output, we see that the in px and py, the points sorted by x and y corrdinates, the closest pair of points (0, 0) and (1, 1) do not appear consecutively.\n",
    "\n",
    "Our divide and conquer algorithm in 2-D space essetntally does the following\n",
    "\n",
    "1. Define a function ClosestPair that accepts two inputs, $p_x$ and $p_y$, the points sorted by x and y coordinates respectively.\n",
    "2. Define the following four arrays\n",
    "- $l_x$: The Left half of $p_x$\n",
    "- $r_x$: The Right half of $p_x$\n",
    "- $l_y$: The Left half of $p_x$, sorted by their occurance in $p_y$\n",
    "- $r_y$: The Right half of $p_x$, sorted by their occurance in $p_y$\n",
    "3. Recursively call ClosestPair twice by passing $l_x, l_y$ and $r_x, r_y$ to find two pair of points $(l_1, l_2)$ and $(r_1,r_2)$, the closest points on the left and right half respectively.\n",
    "4. We then implement a fuction to find the closes points split in left and right. This is the most challenging bit which we will see later. This function called ClosestSplitPair returns the pair of points closest to each other but split in left and right halves. We call these points $(s_1, s_2)$\n",
    "5. Return the closest from these three pair of points $(l_1, l_2), (r_1, r_2), (s_1, s_2)$\n",
    "Now assuming that we implement the function CloseseSplitPair (point 4 above) in linear time, we can implement this entire algorithm in $O(nlogn)$ time as the analysis is very similar to the analysis of Merge sort where we split the input in two and make two recursive calls. Then perform linear work outside the recursive calls.\n",
    "\n",
    "Implementing ClosestSplitPair is not very straightforward and require us to make a subtle change. The relaxed criteria now is that this function can return garbage if the overall closest split pair is present entirely either on the left or right and needs to return the correct answer only if the overall closest pair is split. With this relaxed criteria we now have\n",
    "\n",
    "$$\\delta := min\\{d(l_1, d_2), d(r_1, r_2)\\}\\\\\n",
    "(s_1, s_2) = ClosestSplitPair(P_x, P_y, \\delta)\\\\\n",
    "return \\: best\\: of\\: (l_1, l_2),\\:(r_1, r_2),\\:(s_1, s_2)\n",
    "$$\n",
    "\n",
    "The modified ClosestSplitPair works as follows\n",
    "\n",
    "1. Find the right most point in left half of the points sorted by x coordinates. This becomes the median point. We call this point $\\bar{x}$ and since the points are already sorted by x coordinate, this operation is done in $O(1)$.\n",
    "2. We then filter points $S_y := \\{q_1, q_2\\dots q_l\\}$ which are all points by doing a linear scan of $P_y$ and with their x coordinates between $\\bar{x} - \\delta$ and $\\bar{x} + \\delta$. Thus give is a thick strip of width $2\\delta$ around the median $\\bar{x}$. This operation takes linear time.\n",
    "3. Next, we iterate through these points in $S_y$ and look no further than 7 points (or the last point in the list $S_y$, whichever is minimum) to and compute distance between them in Brute force, on each iteration we compare te distance with the best we have and continue in the loop. The upper bound of this loop iterations is $7n \\leq 7n$ and thus this is a linear time operation and the entire ClosestSplitPair is a linear time operation.\n",
    "\n",
    "Few question still are open.\n",
    "- Why does this algorithm work and guarantees to return an answer\n",
    "- Where did we get this number 7 in the step 3 of the above algorithm. We will now see a proof for the correctness of this algorithm\n",
    "\n",
    "***\n",
    "\n",
    "Proof of correctness of ClosestSplitPair.\n",
    "\n",
    "We start with the following Lemma\n",
    "\n",
    "In the ClosestSplitPair, suppose (p, q) are the closest split pair with $d(p, q) < \\delta$, then\n",
    "\n",
    "- p and q will be included in the set $S_y$\n",
    "- at most 6 points of $S_y$ have a y coordinate in between p and q, or in other words, p and q are not separated by more than 6 points between them in $S_y$\n",
    "\n",
    "The following corollary this follows, if (p, q) is a split pair such that $d(p, q) < \\delta$, then in the ClosestSplitPair algorithm, the points p and q, by above Lemma, are guaranteed to be present in $S_y$ and are no do not have more than 6 points separating them. Since the routine performs am exhaustive search and satisfies both these conditions, the ClosestSplitPair thus indeed find the correct points.\n",
    "\n",
    "We will now prove the Lemma we stated earlier. We will prove the two claims one by one. The first one is simple of the two.\n",
    "\n",
    "**p and q will be included in the set $S_y$**\n",
    "\n",
    "*Proof:*\n",
    "\n",
    "- $p = (x_1, y_1),\\:q = (x_2, y_2)$ and $\\delta $ is the minimum of the distance between points on the left of right. The value $\\bar{x}$ is the median, or the right most x value of the left half.\n",
    "- Since p and q are split with p in left half and q in right half we can conclude that $x_1 \\leq \\bar{x} < x_2$\n",
    "- $\\delta > d(p, q) = \\sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2} \\geq \\sqrt{max\\{|x_1 - x_1|, |y_1 - y_2| \\}}$, thus we can say $|x_1 - x_2|, |y_1 - y_2| < \\delta$\n",
    "- Since $x_1 \\leq \\bar{x}$ and $x_2$ is no more than $\\delta$ bigger than $x_1$, we can say $x_2 \\leq \\bar{x} + \\delta$\n",
    "- Similarly, $x_2 > \\bar{x}$ and is no more than $\\delta$ bigger than $x_2$, we can say $x_1 \\geq \\bar{x} - \\delta$\n",
    "- Thus $x_1 - \\delta \\leq \\bar{x} \\leq x_2 - \\delta$, another way of thing about it as follows. Suppose two people are tied with a rope of length delta, then the maximum they can travel around $\\bar{x}$ is, $\\bar{x} - \\delta$ and $\\bar{x} + \\delta$\n",
    "\n",
    "**In the set $S_y$, there are no more than 6 points separating the true closest points, p and q**\n",
    "\n",
    "*Proof:*\n",
    "\n",
    "To prove this point let us draw 8 boxes, each box of size $\\delta / 2$. 4 boxes will be stacked in a grid of $2\\times 2$ on each side of the $\\bar{x}$. The points p and q are split and are present on either side of $\\bar{x}$. From above proof, we know that both these points are present in $S_y$ and their x coordinates are between $\\bar{x} - \\delta$ and $\\bar{x} + \\delta$. From above proof we know that both points p and q also satisfy the condition $|x_1 - x_2|, |y_1 - y_2| < \\delta$\n",
    "\n",
    "Suppose we draw the boxes such that the lower edge of the lower box passes through the point whose y coordinate is minimum of the two points p and q, then the points p and q will appear in one of the 4 boxes on either side of $\\bar{x}$.\n",
    "\n",
    "How many points can these 8 boxes hold? Let us assume that a box(each of size $\\delta / 2 \\times \\delta/2$) can hold multiple points. The fartest these points can be placed is when they are kept on diagonally opposite corners of this box. By pythagoras theorem, the distance between these points is is $\\sqrt{2}\\cdot \\delta / 2 = \\delta / \\sqrt{2}\\: <\\: \\delta $. Since each box lies entirely on the left of right of $\\bar{x}$, and $\\delta$ is the minimum of the distance between a pair of points on left of right, we can say that no box, each of size $\\delta / 2 \\times \\delta /2$ can contain more than 1 point in it. Since these two points p and q also exist in one of the 8 boxes, we cannot have more than these remaining 6 points between these two points p and q\n",
    "\n",
    "Let us implement the FastPower function and analyse its time complexity and the finally implement the ClosestPair algorithm in python\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('2^1 =', 2)\n",
      "('3^2 =', 9)\n",
      "('4^3 =', 64)\n",
      "('5^11 =', 48828125)\n",
      "('5^0 =', 1)\n",
      "('2^6 =', 64)\n"
     ]
    }
   ],
   "source": [
    "def power_2rec(a, b):\n",
    "    if b == 0:\n",
    "        return 1\n",
    "    if b == 1:\n",
    "        return a\n",
    "    l = power_2rec(a, b//2)\n",
    "    r = power_2rec(a, b - b//2)\n",
    "    return l*r\n",
    "\n",
    "\n",
    "print('2^1 =',power_2rec(2, 1))\n",
    "print('3^2 =',power_2rec(3, 2))\n",
    "print('4^3 =',power_2rec(4, 3))\n",
    "print('5^11 =',power_2rec(5, 11))\n",
    "print('5^0 =',power_2rec(5, 0))\n",
    "print('2^6 =',power_2rec(2, 6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 215,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('2^1 =', 2)\n",
      "('3^2 =', 9)\n",
      "('4^3 =', 64)\n",
      "('5^11 =', 48828125)\n",
      "('2^6 =', 64)\n",
      "('5^0 =', 1)\n"
     ]
    }
   ],
   "source": [
    "def fastpower(a, b):\n",
    "    if b == 0:\n",
    "        return 1\n",
    "    if b == 1:\n",
    "        return a\n",
    "    else:\n",
    "        res = fastpower(a * a, b // 2)\n",
    "        return res if b % 2 == 0 else  res * a\n",
    "       \n",
    "\n",
    "print('2^1 =',fastpower(2, 1))\n",
    "print('3^2 =',fastpower(3, 2))\n",
    "print('4^3 =',fastpower(4, 3))\n",
    "print('5^11 =',fastpower(5, 11))\n",
    "print('2^6 =',fastpower(2, 6))\n",
    "print('5^0 =',fastpower(5, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On each recursive call of fastpower, the value of b halves. Thus the number of times fastpower is called is $log(b)$ times and thus the time complexity of the above function is $O(log(b))$. The log is technically of base 2 but since this is the Big-Oh notation thats not important. On the other hand the simple approach to multiply a b times will have the time complexity $O(b)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 317,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0, 0), (1, 1)]\n"
     ]
    }
   ],
   "source": [
    "def merge_sort(arr, axis = 0):  # axis = 0 for x-axis and 1 for y-axis\n",
    "    n = len(arr)\n",
    "    if n == 1:\n",
    "        return arr\n",
    "\n",
    "    if n >= 1:\n",
    "        left_tree = merge_sort(arr[ : n/2], axis = axis)   \n",
    "        right_tree = merge_sort(arr[n/2 : ], axis = axis)\n",
    "        sorted_arr = [(None, None)] * n\n",
    "        i = 0\n",
    "        j = 0\n",
    "        for k in range(n):\n",
    "            if(i < len(left_tree) and j < len(right_tree)):\n",
    "                if left_tree[i][axis] <= right_tree[j][axis]:\n",
    "                    sorted_arr[k] = left_tree[i]\n",
    "                    i+=1 \n",
    "                elif left_tree[i][axis] > right_tree[j][axis]:\n",
    "                    sorted_arr[k] = right_tree[j]\n",
    "                    j+=1 \n",
    "                \n",
    "            elif i >= len(left_tree) and j < len(right_tree):\n",
    "                sorted_arr[k] = right_tree[j]\n",
    "                j += 1\n",
    "            elif j >= len(right_tree) and i < len(left_tree):\n",
    "                sorted_arr[k] = left_tree[i]\n",
    "                i += 1\n",
    "\n",
    "    return  sorted_arr   \n",
    "\n",
    "    \n",
    "def euclidean_distance(p1, p2):\n",
    "    return ((p1[0]-p2[0])**2 + (p1[1]-p2[1])**2)**0.5\n",
    "\n",
    "\n",
    "def ClosestSplitPair(Px, Py, delta):\n",
    "    x_bar = Px[len(Px)//2]\n",
    "    Sy = []\n",
    "    best = delta\n",
    "    bestpair = None\n",
    "    bestpairx = None\n",
    "    for s in Py:\n",
    "        if  (x_bar[0] - delta) <= s[0] <= (x_bar[0] + delta):\n",
    "            Sy.append(s) \n",
    "    for i in range(0, len(Sy)-1):\n",
    "        for j in range(1, min(7, len(Sy) - i)):\n",
    "            dist_ps = euclidean_distance(Sy[i], Sy[i+j])\n",
    "            if dist_ps < delta:\n",
    "                best = dist_ps\n",
    "                bestpair = [Sy[i], Sy[i+j]]\n",
    "                bestpairx = bestpair if bestpair[0][0] <= bestpair[1][0]  else [bestpair[1], bestpair[0]] \n",
    "                \n",
    "    return bestpairx, bestpair\n",
    "    \n",
    "def ClosestPair(Px, Py):\n",
    "    n = len(Px)\n",
    "    if n == 1:\n",
    "        return None\n",
    "    if n == 2:\n",
    "        return Px , Py\n",
    "    if n == 3:\n",
    "        dist_pxy1 = euclidean_distance(Px[0], Px[1]) \n",
    "        dist_pxy2 = euclidean_distance(Px[1], Px[2]) \n",
    "        dist_pxy3 = euclidean_distance(Px[0], Px[2]) \n",
    "        if dist_pxy1 <= dist_pxy2 and dist_pxy1 <= dist_pxy3 :\n",
    "            newPx = Px[0:2][:]\n",
    "            return newPx, newPx if newPx[0][1] <= newPx[1][1]  else [newPx[1], newPx[0]]\n",
    "        elif dist_pxy1 > dist_pxy2 and dist_pxy2 <= dist_pxy3 :\n",
    "            newPx = Px[1:3][:]\n",
    "            return newPx, newPx if newPx[0][1] <= newPx[1][1]  else [newPx[1], newPx[0]]\n",
    "        else:\n",
    "            newPx = [Px[0], Px[2]]\n",
    "            return newPx, newPx if newPx[0][1] <= newPx[1][1]  else [newPx[1], newPx[0]]\n",
    "    \n",
    "    p1, p1y = ClosestPair(Px[:n//2], Py[:n//2])  # Qx, Qy\n",
    "    p2, p2y = ClosestPair(Px[n//2:], Py[n//2:])  # Rx, Ry\n",
    "    \n",
    "    closepairx = None\n",
    "    closepairy = None\n",
    "    dist_p1 = euclidean_distance(p1[0], p1[1]) \n",
    "    dist_p2 = euclidean_distance(p2[0], p2[1])\n",
    "    if dist_p1 <= dist_p2:\n",
    "        delta = dist_p1\n",
    "        closepairx = p1\n",
    "        closepairy = p1y\n",
    "    elif dist_p1 > dist_p2:\n",
    "        delta = dist_p2\n",
    "        closepairx = p2\n",
    "        closepairy = p2y\n",
    "    \n",
    "    p3, p3y = ClosestSplitPair(Px, Py, delta)\n",
    "        \n",
    "    if p3 != None:\n",
    "        dist_p3 = ((p3[0][0]-p3[1][0])**2 + (p3[0][1]-p3[1][1])**2)**0.5  # Euclidean distance\n",
    "        closepairx = closepairx if dist_p3 > delta else p3\n",
    "        closepairy = closepairy if dist_p3 > delta else p3y\n",
    "        \n",
    "    return closepairx, closepairy \n",
    "    \n",
    "    \n",
    "def ClosestPairAlgorithm(P):\n",
    "    Px = merge_sort(P, 0)\n",
    "    Py = merge_sort(P, 1)\n",
    "    return ClosestPair(Px, Py)\n",
    "\n",
    "\n",
    "print(ClosestPairAlgorithm( [(0, 0), (1, 1), (0.5, 4), (4, 0.5), (-2, 2)] )[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 318,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAADNtJREFUeJzt3U+InPd9x/HPR39MsoqNoZoGE2W1+OJLoJYZDEHFpHYIDjGmhx4cNoWGwl5KsGlLSKJDyUHXkJ4Cg5zUJROH1I4uJjUxJMYxNA6zslP/kXqo8cqqk2pMSG1loSb2p4cZ1ZbY1TwTzTOz3933C5bZ59lHM99h0JtHv3lG6yQCANSxb9EDAACmQ7gBoBjCDQDFEG4AKIZwA0AxhBsAiiHcAFAM4QaAYgg3ABRzoI07PXz4cFZWVtq4awDYldbX199M0mlybCvhXllZ0WAwaOOuAWBXsr3R9FiWSgCgGMINAMUQbgAohnADQDGEGwCKmRhu27fZfuEDX2/ZfmgewwFACf2+tLIi7ds3uu33W324iZcDJvkPSbdLku39kv5L0ulWpwKAKvp9aW1N2twcbW9sjLYlaXW1lYecdqnkHkn/maTx9YYAsKudOPF+tC/b3Bztb8m04X5A0qNb/cD2mu2B7cFwOLz+yQCggvPnp9s/A43DbfsGSfdL+petfp6kl6SbpNvpNPrUJgDUt7w83f4ZmOaM+7OSziT577aGAYByTp6Ulpau3Le0NNrfkmnC/Xlts0wCAHvW6qrU60lHj0r26LbXa+2NSUlykskH2UuSXpd0a5L/mXR8t9sN/8kUADRnez1Jt8mxjf53wCSbkv7ouqYCAMwEn5wEgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAiiHcAFAM4QaAYgg3ABRDuAGgGMINAMUQbgAohnADQDGEGwCKIdwAUAzhBoBiCDcAFEO4AaAYwg0AxRBuACimUbht32z7MdvnbJ+1/cm2BwMAbO1Aw+P+UdKTSf7C9g2SllqcCQBwDRPDbfsmSXdJ+itJSvKOpHfaHQsAsJ0mSyW3ShpK+o7t522fsn2o5bkAANtoEu4Dku6Q9K0kxyT9TtJXrj7I9prtge3BcDic8ZgAgMuahPuCpAtJnhtvP6ZRyK+QpJekm6Tb6XRmOSMA4AMmhjvJryW9bvu28a57JL3S6lQAgG01varkS5L64ytKXpX0xfZGAgBcS6NwJ3lBUrflWQAADfDJSQAohnADQDGEGwCKIdwAUAzhBoBiCDcAFEO4AaAYwg0AxRBuACiGcANAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAiiHcAFAM4QaAYhr9lnfbr0l6W9K7kn6fhN/4DgAL0ijcY3+W5M3WJgEANMJSCQAU0zTckfRj2+u219ocCABwbU2XSo4necP2H0t6yva5JM988IBx0NckaXl5ecZjAgAua3TGneSN8e1FSacl3bnFMb0k3STdTqcz2ykBAP9vYrhtH7J94+XvJX1G0kttDwYA2FqTpZKPSjpt+/Lx30vyZKtTAQC2NTHcSV6V9CdzmAUA0ACXAwJAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAiiHcAFAM4QaAYgg3ABRDuAGgGMINAMUQbgAohnADQDGEGwCKIdwAUAzhBoBiCDcAFEO4AaCYxuG2vd/287afaGWSfl9aWZH27Rvd9vutPAwWjNcZuG4Hpjj2QUlnJd008yn6fWltTdrcHG1vbIy2JWl1deYPhwXhdQZmotEZt+0jkj4n6VQrU5w48f5f5ss2N0f7sXvwOgMz0XSp5JuSvizpve0OsL1me2B7MBwOp5vi/Pnp9qMmXmdgJiaG2/Z9ki4mWb/WcUl6SbpJup1OZ7oplpen24+aeJ2BmWhyxn1c0v22X5P0fUl32/7uTKc4eVJaWrpy39LSaD92D15nYCYmhjvJV5McSbIi6QFJP0nyhZlOsboq9XrS0aOSPbrt9XjDarfhdQZmwkmaH2x/StLfJ7nvWsd1u90MBoPrHA0A9g7b60m6TY6d5nJAJXla0tN/wEwAgBnhk5MAUAzhBoBiCDcAFEO4AaAYwg0AxRBuACiGcANAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAiiHcAFAM4QaAYgg3ABRDuAGgGMINAMVMDLftD9n+he1f2n7Z9tfnMRgAYGsHGhzzv5LuTnLJ9kFJz9r+1yQ/b3k2AMAWJoY7SSRdGm8eHH+lzaEAANtrtMZte7/tFyRdlPRUkufaHQsAsJ1G4U7ybpLbJR2RdKftT1x9jO012wPbg+FwOOs5AQBjU11VkuS3kp6WdO8WP+sl6SbpdjqdGY0HALhak6tKOrZvHn//YUmflnSu7cEAAFtrclXJLZIesb1fo9D/IMkT7Y4FANhOk6tK/l3SsTnMAgBogE9OAkAxhBsAiiHcAFAM4QaAYgg3ABRDuAGgGMINAMUQbgAohnADQDGEGwCKIdwAUAzhBoBiCDcAFEO4AaAYwg0AxRBuACiGcANAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoJiJ4bb9cds/tX3W9su2H2xlkn5fWlmR9u0b3fb7rTwMAFR3oMExv5f0d0nO2L5R0rrtp5K8MrMp+n1pbU3a3Bxtb2yMtiVpdXVmDwMAu8HEM+4kv0pyZvz925LOSvrYTKc4ceL9aF+2uTnaDwC4wlRr3LZXJB2T9NwWP1uzPbA9GA6H001x/vx0+wFgD2scbtsfkfS4pIeSvHX1z5P0knSTdDudznRTLC9Ptx8A9rBG4bZ9UKNo95P8cOZTnDwpLS1duW9pabQfAHCFJleVWNLDks4m+UYrU6yuSr2edPSoZI9uez3emASALTjJtQ+w/1TSzyS9KOm98e6vJfnRdn+m2+1mMBjMbEgA2O1sryfpNjl24uWASZ6V5OueCgAwE3xyEgCKIdwAUAzhBoBiCDcAFEO4AaAYwg0AxRBuACiGcANAMYQbAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0AxhBsAiiHcAFAM4QaAYgg3ABRDuAGgmInhtv1t2xdtvzSPgQAA19bkjPufJN3b8hwAgIYmhjvJM5J+M4dZAAANsMYNAMXMLNy212wPbA+Gw+Gs7hYAcJWZhTtJL0k3SbfT6czqbgEAV2GpBACKaXI54KOS/k3SbbYv2P7r9scCAGznwKQDknx+HoMAAJphqQQAiiHcAFAM4QaAYgg3ABRDuAGgGMINAMUQbgAohnBjvvp9aWVF2rdvdNvvL3oioJyJH8ABZqbfl9bWpM3N0fbGxmhbklZXFzcXUAxn3JifEyfej/Zlm5uj/QAaI9yYn/Pnp9sPYEuEG/OzvDzdfgBbItyYn5MnpaWlK/ctLY32A2iMcGN+VlelXk86elSyR7e9Hm9MAlPiqhLM1+oqoQauE2fcAFAM4QaAYgg3ABRDuAGgGMINAMU4yezv1B5K2vgD//hhSW/OcJwKeM673157vhLPeVpHk3SaHNhKuK+H7UGS7qLnmCee8+63156vxHNuE0slAFAM4QaAYnZiuHuLHmABeM673157vhLPuTU7bo0bAHBtO/GMGwBwDTsm3La/bfui7ZcWPcs82P647Z/aPmv7ZdsPLnqmttn+kO1f2P7l+Dl/fdEzzYvt/baft/3EomeZB9uv2X7R9gu2B4ueZx5s32z7Mdvnxn+vP9naY+2UpRLbd0m6JOmfk3xi0fO0zfYtkm5Jcsb2jZLWJf15klcWPFprbFvSoSSXbB+U9KykB5P8fMGjtc7230rqSropyX2Lnqdttl+T1E2yZ67jtv2IpJ8lOWX7BklLSX7bxmPtmDPuJM9I+s2i55iXJL9Kcmb8/duSzkr62GKnaldGLo03D46/dsaZQ4tsH5H0OUmnFj0L2mH7Jkl3SXpYkpK801a0pR0U7r3M9oqkY5KeW+wk7RsvGbwg6aKkp5Ls+ucs6ZuSvizpvUUPMkeR9GPb67bXFj3MHNwqaSjpO+MlsVO2D7X1YIR7wWx/RNLjkh5K8tai52lbkneT3C7piKQ7be/qZTHb90m6mGR90bPM2fEkd0j6rKS/GS+F7mYHJN0h6VtJjkn6naSvtPVghHuBxuu8j0vqJ/nhoueZp/E/I5+WdO+CR2nbcUn3j9d8vy/pbtvfXexI7Uvyxvj2oqTTku5c7EStuyDpwgf+BfmYRiFvBeFekPEbdQ9LOpvkG4ueZx5sd2zfPP7+w5I+LencYqdqV5KvJjmSZEXSA5J+kuQLCx6rVbYPjd9w13i54DOSdvXVYkl+Lel127eNd90jqbULDXbM75y0/aikT0k6bPuCpH9I8vBip2rVcUl/KenF8ZqvJH0tyY8WOFPbbpH0iO39Gp00/CDJnrg8bo/5qKTTo3MTHZD0vSRPLnakufiSpP74ipJXJX2xrQfaMZcDAgCaYakEAIoh3ABQDOEGgGIINwAUQ7gBoBjCDQDFEG4AKIZwA0Ax/wcGEzhf9huFngAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d22cf90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(1, 2), (1, 4)]\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "points = [(3, 4), (1, 2), (1, 4), (3, 0.5), (6, 7)] \n",
    "plt.plot([x[0] for x in points], [x[1] for x in points], 'ro')\n",
    "plt.show()\n",
    "print(ClosestPairAlgorithm( points )[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 319,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAADg5JREFUeJzt3U+IJOd9xvHnmV0Zu2UJQTQxwvLM4ItyMFh/GoFREIpkjEyEyNFhEkgIzMUYCSUEx3sIPmxyM8nJ0EhOFNSScSTrIhxhQSxkQSzTI8vWn1UOMTurjezsiODoz0CEpSeH7o208uxOtbara37d3w801fXOO1W/l6KfrX2rusZJBACoY6XrAgAA0yG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAijnaxkavvPLKbGxstLFpAFhI29vbryVZbdK3leDe2NjQaDRqY9MAsJBs7zTty1QJABRDcANAMQQ3ABRDcANAMQQ3ABRzYHDbvsb2c+97vW777nkUBwAlDIfSxoa0sjJeDoet7u7A2wGT/LukayXJ9hFJ/ynp0VarAoAqhkNpa0va2xuv7+yM1yVpc7OVXU47VXKbpP9I0vh+QwBYaMeOvRfaZ+3tjdtbMm1wf0nSQ/v9wPaW7ZHt0e7u7sVXBgAVnDo1XfsMNA5u2x+RdKekf97v50kGSfpJ+qurjb61CQD1ra1N1z4D05xxf1HSs0n+q61iAKCc48elXu/ctl5v3N6SaYL7D3WeaRIAWFqbm9JgIK2vS/Z4ORi0dmFSkpzk4E52T9Irkj6d5H8O6t/v98NDpgCgOdvbSfpN+jZ6OmCSPUm/dVFVAQBmgm9OAkAxBDcAFENwA0AxBDcAFENwA0AxBDcAFENwA0AxBDcAFENwA0AxBDcAFENwA0AxBDcAFENwA0AxBDcAFENwA0AxBDcAFENwA0AxBDcAFENwA0AxBDcAFENwA0AxjYLb9hW2H7b9su0Ttj/XdmEAgP01PeP+e0mPJ/kdSZ+VdKK9kpbIcChtbEgrK+PlcNh1RQA+jDl/lo8e1MH25ZJulvQnkpTkbUlvt1rVMhgOpa0taW9vvL6zM16XpM3N7uoCMJ0OPstOcuEO9rWSBpJe0vhse1vSXUneOt/v9Pv9jEajWda5eDY2xgf4g9bXpZMn510NgA9rRp9l29tJ+k36NpkqOSrpeknfTHKdpLckfXWfnW7ZHtke7e7uNi52aZ06NV07gMOpg89yk+A+Lel0kmcm6w9rHOTnSDJI0k/SX11dnWWNi2ltbbp2AIdTB5/lA4M7yS8lvWL7mknTbRpPm+BiHD8u9XrntvV643YAdXTwWW56V8lXJA1t/0zStZL+prWKlsXmpjQYjOfB7PFyMODCJFBNB5/lAy9OfhhcnASA6cz64iQA4BAhuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIo52qST7ZOS3pD0jqRfN/1LxACA2WsU3BO/l+S11ioBADTCVAkAFNM0uCPp+7a3bW+1WRAA4MKaTpXclORV278t6QnbLyd56v0dJoG+JUlra2szLhMAcFajM+4kr06WZyQ9KunGffoMkvST9FdXV2dbJQDg/x0Y3LYvtX3Z2feSviDphbYLAwDsr8lUySckPWr7bP8HkzzealUAgPM6MLiT/FzSZ+dQCwCgAW4HBIBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKKZxcNs+Yvsnth9rpZLhUNrYkFZWxsvhsJXdoGPLdpy7HG9X+162YyzNf8xJGr0k3SPpQUmPHdT3hhtuyFQeeCDp9RLpvVevN27H4li249zleLva97Id42RmY5Y0SsM89rj/hdm+WtL9ko5LuifJHRfq3+/3MxqNmv/rsbEh7ez8Zvv6unTyZPPt4HBbtuPc5Xi72veyHWNpZmO2vZ2k36hvw+B+WNLfSrpM0l/sF9y2tyRtSdLa2toNO/sN5HxWVsb/Tv3mRqV3322+HRxuy3acuxxvV/tetmMszWzM0wT3gXPctu+QdCbJ9oX6JRkk6Sfpr66uNix1Ym1tunbUtGzHucvxdrXvZTvGUidjbnJx8iZJd9o+Kenbkm61/cBMqzh+XOr1zm3r9cbtWBzLdpy7HG9X+162Yyx1M+amk+GTKZVb1MbFyWQ8kb++ntjj5SJfzFhmy3acuxxvV/tetmOczGTMmvXFybNs36LzzHG/39QXJwFgyU0zx310mg0neVLSkx+iJgDAjPDNSQAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAo5sDgtv1R2z+2/VPbL9r++jwKAwDs72iDPv8r6dYkb9q+RNLTtv8lyY9arg0AsI8DgztJJL05Wb1k8kqbRQEAzq/RHLftI7afk3RG0hNJnmm3LADA+TQK7iTvJLlW0tWSbrT9mQ/2sb1le2R7tLu7O+s6AQATU91VkuRXkp6UdPs+Pxsk6Sfpr66uzqg8AMAHNbmrZNX2FZP3H5P0eUkvt10YAGB/Te4quUrS/baPaBz030nyWLtlAQDOp8ldJT+TdN0cagEANMA3JwGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIo5MLhtf8r2D2yfsP2i7btaqWQ4lDY2pJWV8XI4bGU3AFDd0QZ9fi3pz5M8a/sySdu2n0jy0syqGA6lrS1pb2+8vrMzXpekzc2Z7QYAFsGBZ9xJfpHk2cn7NySdkPTJmVZx7Nh7oX3W3t64HQBwjqnmuG1vSLpO0jP7/GzL9sj2aHd3d7oqTp2arh0Alljj4Lb9cUmPSLo7yesf/HmSQZJ+kv7q6up0VaytTdcOAEusUXDbvkTj0B4m+e7Mqzh+XOr1zm3r9cbtAIBzNLmrxJLuk3QiyTdaqWJzUxoMpPV1yR4vBwMuTALAPpzkwh3s35X0Q0nPS3p30vy1JN873+/0+/2MRqOZFQkAi872dpJ+k74H3g6Y5GlJvuiqAAAzwTcnAaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaCYA4Pb9rdsn7H9wjwKAgBcWJMz7n+UdHvLdQBYFMOhtLEhrayMl8Nh1xW1b85jPnpQhyRP2d5otQoAi2E4lLa2pL298frOznhdkjY3u6urTR2MmTluALNz7Nh7AXbW3t64fVF1MOaZBbftLdsj26Pd3d1ZbRZAJadOTde+CDoY88yCO8kgST9Jf3V1dVabBVDJ2tp07YuggzEzVQJgdo4fl3q9c9t6vXH7oupgzE1uB3xI0r9Jusb2adt/1lo1AGrb3JQGA2l9XbLHy8FgcS9MSp2M2UlmvtF+v5/RaDTz7QLAorK9naTfpC9TJQBQDMENAMUQ3ABQDMENAMUQ3ABQDMENAMUQ3ABQDMGN+VrGR34CM3bgY12BmVnGR34CLeCMG/OzjI/8BFpAcGN+lvGRn0ALCG7MzzI+8hNoAcGN+VnGR34CLSC4MT/L+MhPoAXcVYL52twkqIGLxBk3ABRDcANAMQQ3ABRDcANAMQQ3ABTTyh8Ltr0raedD/vqVkl6bYTkVMObFt2zjlRjztNaTrDbp2EpwXwzbo6Z/6XhRMObFt2zjlRhzm5gqAYBiCG4AKOYwBveg6wI6wJgX37KNV2LMrTl0c9wAgAs7jGfcAIALODTBbftbts/YfqHrWubB9qds/8D2Cdsv2r6r65raZvujtn9s+6eTMX+965rmxfYR2z+x/VjXtcyD7ZO2n7f9nO1R1/XMg+0rbD9s++XJ5/pzre3rsEyV2L5Z0puS/inJZ7qup222r5J0VZJnbV8maVvSHyR5qePSWmPbki5N8qbtSyQ9LemuJD/quLTW2b5HUl/S5Unu6Lqettk+KamfZGnu47Z9v6QfJrnX9kck9ZL8qo19HZoz7iRPSfrvruuYlyS/SPLs5P0bkk5I+mS3VbUrY29OVi+ZvA7HmUOLbF8t6fcl3dt1LWiH7csl3SzpPklK8nZboS0douBeZrY3JF0n6ZluK2nfZMrgOUlnJD2RZOHHLOnvJP2lpHe7LmSOIun7trdtb3VdzBx8WtKupH+YTInda/vStnZGcHfM9sclPSLp7iSvd11P25K8k+RaSVdLutH2Qk+L2b5D0pkk213XMmc3Jble0hclfXkyFbrIjkq6XtI3k1wn6S1JX21rZwR3hybzvI9IGib5btf1zNPkv5FPSrq941LadpOkOydzvt+WdKvtB7otqX1JXp0sz0h6VNKN3VbUutOSTr/vf5APaxzkrSC4OzK5UHefpBNJvtF1PfNge9X2FZP3H5P0eUkvd1tVu5L8VZKrk2xI+pKkf03yRx2X1Srbl04uuGsyXfAFSQt9t1iSX0p6xfY1k6bbJLV2o8Gh+ZuTth+SdIukK22flvTXSe7rtqpW3STpjyU9P5nzlaSvJflehzW17SpJ99s+ovFJw3eSLMXtcUvmE5IeHZ+b6KikB5M83m1Jc/EVScPJHSU/l/Snbe3o0NwOCABohqkSACiG4AaAYghuACiG4AaAYghuACiG4AaAYghuACiG4AaAYv4PKIiGqdskyRUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bf2d790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(4.5, 4), (5, 4)]\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "points = [(3, 4), (1, 2), (1, 4), (3, 0.5), (6, 7), (5, 4), (4, 4), (6, 4), (4.5, 4), (6, 6), (2, 6), (5, 1), (6, 1)] \n",
    "plt.plot([x[0] for x in points], [x[1] for x in points], 'ro')\n",
    "plt.show()\n",
    "print(ClosestPairAlgorithm( points )[0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IV. THE MASTER METHOD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem 1:** <br>\n",
    "You are given as input an unsorted array of n distinct numbers, where n is a power of 2. Give an algorithm that identifies the $2^{nd}$ largest number in the array, and that uses at most $n + log_2n - 2$ comparisons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 338,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Second-largest number : 92\n"
     ]
    }
   ],
   "source": [
    "def second_largest(arr):\n",
    "    n = len(arr)\n",
    "    #print(\"arr  {}\".format(arr))\n",
    "    if n == 2:\n",
    "        return [arr[1], arr[0]] if arr[1] > arr[0] else arr\n",
    "\n",
    "    if n > 2:\n",
    "        left_tree = second_largest(arr[ : n/2])   \n",
    "        #print(\"left_tree  {}\".format(left_tree))\n",
    "        right_tree = second_largest(arr[n/2 : ])\n",
    "        #print(\"right_tree  {}\".format(right_tree))\n",
    "        largest2_arr = [None] * 4\n",
    "        \n",
    "        if left_tree[0] >= right_tree[0]:\n",
    "            largest2_arr[0] = left_tree[0]  \n",
    "            largest2_arr[1] = right_tree[0] \n",
    "        else:\n",
    "            largest2_arr[0] = right_tree[0]\n",
    "            largest2_arr[1] = left_tree[0]\n",
    "            \n",
    "        largest2_arr[2] = left_tree[1]\n",
    "        largest2_arr[3] = right_tree[1]\n",
    "            \n",
    "    return  largest2_arr   \n",
    "\n",
    "\n",
    "def select_second_largest(array):\n",
    "    outarr = second_largest(array) \n",
    "    second_max = outarr[1]\n",
    "    if second_max < outarr[2]:\n",
    "        second_max = outarr[2]\n",
    "    if second_max < outarr[3]:\n",
    "        second_max = outarr[3]\n",
    "    print(\"Second-largest number : {}\".format(second_max))\n",
    "    \n",
    "    \n",
    "\n",
    "#select_second_largest([40, 50, 5, 10])\n",
    "#select_second_largest([ 5, 10, 40, 50])\n",
    "#select_second_largest([ 3, 20, 33, 2, 59, 65, 82, 3, 5, 10, 40, 50, 99, 45, 92, 43])\n",
    "select_second_largest([ 3, 20, 33, 99, 59, 65, 82, 3, 5, 10, 40, 50, 2, 45, 92, 43]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem 2:** <br>\n",
    "You are a given a unimodal array of n distinct elements, meaning that its entries are in increasing order up until its maximum element, after which its elements are in decreasing order. Give an algorithm to compute the maximum element that runs in O(log n) time.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 352,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9\n"
     ]
    }
   ],
   "source": [
    "def maximum_in_unimodal(array):\n",
    "    if len(array) == 1:\n",
    "        return array[0]\n",
    "    elif len(array) == 2:\n",
    "        return max(array[0], array[1])\n",
    "    \n",
    "    mid = len(array) // 2 \n",
    "    if array[mid-1] > array[mid] and array[mid] > array[mid+1]:\n",
    "        return maximum_in_unimodal(array[:mid])\n",
    "    elif array[mid-1] < array[mid] and array[mid] < array[mid+1]:\n",
    "        return maximum_in_unimodal(array[mid+1:])\n",
    "    else:\n",
    "        return array[mid]\n",
    "    \n",
    "#print(maximum_in_unimodal([8, 9, 8]))\n",
    "#print(maximum_in_unimodal([5, 6, 7, 8, 9, 8, 7, 6, 5, 4, 3, 2, 1]))\n",
    "print(maximum_in_unimodal([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 7, 6]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem 3:** <br>\n",
    "You are given a sorted (from smallest to largest) array A of n distinct integers which can be positive, negative, or zero. You want to decide whether or not there is an index i such that A[i] = i. Design the fastest algorithm that you can for solving this problem.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 384,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4\n"
     ]
    }
   ],
   "source": [
    "def index_value_equality(array, idx = 0):\n",
    "    if len(array) == 1:\n",
    "        if array[0] == idx:\n",
    "            return array[0]\n",
    "        else:\n",
    "            return None\n",
    "    elif len(array) == 2:\n",
    "        if array[0] == idx:\n",
    "            return array[0]\n",
    "        elif array[1] == idx+1:\n",
    "            return array[1]\n",
    "        else:\n",
    "            return None\n",
    "    \n",
    "    mid = len(array) // 2 \n",
    "    if array[mid] > mid:\n",
    "        return index_value_equality(array[:mid], idx)\n",
    "    elif array[mid] < mid:\n",
    "        return index_value_equality(array[mid+1:], idx = mid+1)\n",
    "    else:\n",
    "        return array[mid]\n",
    "\n",
    "#print(index_value_equality([0, 2, 70]))    \n",
    "#print(index_value_equality([-1, 0, 2]))    \n",
    "#print(index_value_equality([-50, 1, 70]))\n",
    "#print(index_value_equality([-50, 1, 70, 100]))\n",
    "print(index_value_equality([-20, -15, -10, -5, 4, 70, 80, 90, 100]))\n",
    "#print(index_value_equality([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 7, 6]))\n",
    "#print(index_value_equality([-4, -3, -2, 3, 5, 6]))\n",
    "#print(index_value_equality([-4, -3, -2, -1, 0, 5]))\n",
    "#print(index_value_equality([-4, -3, -2, -1, 0]))\n",
    "#print(index_value_equality([-1, 0, 1, 2, 4, 5, 7]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem 4:** <br>\n",
    "You are given an n by n grid of distinct numbers. A number is a local minimum if it is smaller than all of its neighbors. (A neighbor of a number is one immediately above, below, to the left, or the right. Most numbers have four neighbors; numbers on the side have three; the four corners have two.) Use the divide-and-conquer algorithm design paradigm to compute a local minimum with only O(n) comparisons between pairs of numbers. (Note: since there are $n^2$\n",
    "  numbers in the input, you cannot afford to look at all of them. Hint: Think about what types of recurrences would give you the desired upper bound.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 455,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matrix_local_minimum(matrix, setglobal):\n",
    "    global colidx\n",
    "    global matrix_main \n",
    "    if setglobal == 0:\n",
    "        colidx = matrix.shape[1] - 1\n",
    "        matrix_main = matrix\n",
    "        setglobal = 1\n",
    "    \n",
    "    m = matrix.shape[1] // 2\n",
    "    i_min = np.argmin(matrix[:,m])  # Done in O(n)\n",
    "    min_val = matrix[i_min, m]\n",
    "    #Comparison done in constant time, the recursive calls would be done log n times\n",
    "    print(\"m: {}, i_min: {}\".format(m, i_min))\n",
    "    left_val = matrix_main[i_min, m-1] if m > 0  else float('inf')\n",
    "    right_val = matrix_main[i_min, m+1] if m < colidx  else float('inf')\n",
    "    if min_val < left_val and min_val < right_val:\n",
    "        return (i_min, np.where(matrix_main == min_val)[1][0])\n",
    "    elif  min_val > left_val and m > 0 :\n",
    "        return matrix_local_minimum(matrix[:,:m], 1) \n",
    "    elif  min_val > right_val and  m < colidx:\n",
    "        return matrix_local_minimum(matrix[:,m+1:], 1)\n",
    "    else:\n",
    "        return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 456,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input is\n",
      " [[12 11 10  9]\n",
      " [13  6  5  4]\n",
      " [14  7  3  1]\n",
      " [15  8  2  0]]\n",
      "m: 2, i_min: 3\n",
      "m: 0, i_min: 3\n",
      "\n",
      "A local minima is at (3, 3) with value 0\n"
     ]
    }
   ],
   "source": [
    "#array = np.matrix([[1, 2, 5, 4], [0, -2, 9, 6], [7, -1, 3, 10], [11, 12, 13, 14]])\n",
    "#array = np.matrix([[1, 2, 5, 4], [0, 9, 9, 6], [7, 7, 7, 10], [11, 12, 13, -3]])\n",
    "#array = np.matrix([[-3, 2, 5, 4], [0, 9, 9, 6], [7, 7, 7, 10], [11, 12, 13, -3]])\n",
    "array = np.matrix([[12, 11, 10, 9], [13, 6, 5, 4], [14, 7, 3, 1], [15, 8, 2, 0]])\n",
    "print('Input is\\n {}'.format(array))\n",
    "ix = matrix_local_minimum(array, setglobal = 0)\n",
    "print('\\nA local minima is at {} with value {}'.format(ix, array[ix]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem 5:** <br>\n",
    "Download the following text file:  IntegerArray.txt\n",
    "\n",
    "This file contains all of the 100,000 integers between 1 and 100,000 (inclusive) in some order, with no integer repeated.\n",
    "\n",
    "Your task is to compute the number of inversions in the file given, where the $i^{th}$\n",
    "  row of the file indicates the $i^{th}$\n",
    "  entry of an array.\n",
    "\n",
    "Because of the large size of this array, you should implement the fast divide-and-conquer algorithm.\n",
    "\n",
    "\n",
    "[TIP: before submitting, first test the correctness of your program on some small test files or your own devising.]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 459,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/aishwaya/anaconda2/lib/python2.7/site-packages/urllib3/connectionpool.py:858: InsecureRequestWarning: Unverified HTTPS request is being made. Adding certificate verification is strongly advised. See: https://urllib3.readthedocs.io/en/latest/advanced-usage.html#ssl-warnings\n",
      "  InsecureRequestWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2407905288\n"
     ]
    }
   ],
   "source": [
    "import urllib3\n",
    "http = urllib3.PoolManager()\n",
    "\n",
    "# Test case\n",
    "r1 = http.request('GET', \"https://d3c33hcgiwev3.cloudfront.net/_bcb5c6658381416d19b01bfc1d3993b5_IntegerArray.txt?Expires=1559692800&Signature=FnkGMm2moEXaXZOXClN4Xr~7pIIkEBU1XnfHiQvX4WuTmGCpTZMowf~aMUxplVRe3ELbMWrJdjuQvuzneGkGzEZ29hiRp8I5m5x0KJks3uI456R01269xT83jK2w9Ug73gkjCGf9VtuOQrxjMxYtySd-QyiMYccV7czkJkm-~uA_&Key-Pair-Id=APKAJLTNE6QMUY6HBC5A\")\n",
    "IntegerArrayString = r1.data.split('\\r\\n')\n",
    "del IntegerArrayString[-1]\n",
    "IntegerArray = [int(i) for i in IntegerArrayString]\n",
    "#count_inversions_brute_force(IntegerArray)\n",
    "print(sort_and_count_inversions(IntegerArray , count_inversions = 0)[1])    "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}