{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mathematical Techniques"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Note that Parsl is not effective if multiple CPU cores aren't available because Parsl's ability to execute tasks in parallel is depenedent on the availability multiple cores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cores available: 4\n"
     ]
    }
   ],
   "source": [
    "import multiprocessing\n",
    "print('Cores available: {}'.format(multiprocessing.cpu_count()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing Libraries and Configuration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll be using the htex configuration for Parsl. Read more [here.]( https://github.com/Parsl/parsl/blob/master/parsl/configs/htex_local.py)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<parsl.dataflow.dflow.DataFlowKernel at 0x10c51edd0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import parsl\n",
    "import os\n",
    "from parsl.app.app import python_app, bash_app\n",
    "from parsl.providers import LocalProvider\n",
    "from parsl.channels import LocalChannel\n",
    "\n",
    "from parsl.config import Config\n",
    "from parsl.executors import HighThroughputExecutor\n",
    "\n",
    "config = Config(\n",
    "    executors=[\n",
    "        HighThroughputExecutor(\n",
    "            label=\"htex_local\",\n",
    "            cores_per_worker=1,\n",
    "            provider=LocalProvider(\n",
    "                channel=LocalChannel(),\n",
    "                init_blocks=1,\n",
    "                max_blocks=1,\n",
    "            ),\n",
    "        )\n",
    "    ],\n",
    ")\n",
    "\n",
    "parsl.load(config)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Block Method of Matrix Multiplication"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The block method of matrix multiplication between two matrices A and B of sizes sXq and rXs is as follows:\n",
    "\n",
    "![](./images/formula.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We carry out the operation above for two matrices A and B. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_a = [[2,3,4,5],[4,5,6,7],[4,5,6,7]]\n",
    "matrix_b = [[3,4],[4,5],[6,7],[8,9]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "A python app for multiplying corresponding elements\n",
    "'''\n",
    "\n",
    "@python_app\n",
    "def multiply(inputs=[]):\n",
    "    return inputs[0]*inputs[1]\n",
    "\n",
    "@python_app\n",
    "def sum_elements(inputs=[]):\n",
    "    return sum(inputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "A Function for two matrices that uses the block method to multiply them.\n",
    "'''\n",
    "\n",
    "def multiply_matrices(A,B):\n",
    "    C = [[0 for _ in range(len(B[0]))] for _ in range(len(A))] # Result matrix\n",
    "    \n",
    "    for q in range(len(A)):\n",
    "        for r in range(len(B[0])):\n",
    "            \n",
    "            s = len(A[0])\n",
    "            \n",
    "            elements = []\n",
    "            for i in range(s):\n",
    "                elements.append(multiply(inputs=[A[q][i],B[i][r]])) # Multiplying elements\n",
    "            elements = [i.result() for i in elements] # evaluation for each row\n",
    "            C[q][r] = sum_elements(inputs=elements)\n",
    "            \n",
    "    C = [[element.result() for element in row] for row in C]\n",
    "    \n",
    "    return C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[82, 96], [124, 146], [124, 146]]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''\n",
    "Example implementation of the two given matrices\n",
    "'''\n",
    "multiply_matrices(matrix_a, matrix_b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Markov Chains"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We build a simple implementation of markov chains using Parsl. A Markov Chain is a chain of samples obtained using a transition matrix. Each sample is only dependent on the sample before it. We'll run 5 different chains for different amount of steps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will use the Parsl code for the block multiplication of matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def multiply_matrices_for_markov(A,B):\n",
    "    '''\n",
    "    This function takes two matrices A and B, and uses the matrix multiplication in Parsl to compute C.\n",
    "    \n",
    "    Here, C = AxB\n",
    "    '''\n",
    "    C = [[0 for _ in range(len(B[0]))] for _ in range(len(A))] # Result matrix\n",
    "    \n",
    "    for q in range(len(A)):\n",
    "        for r in range(len(B[0])):\n",
    "            \n",
    "            s = len(A[0])\n",
    "            \n",
    "            elements = []\n",
    "            for i in range(s):\n",
    "                elements.append(multiply(inputs=[A[q][i],B[i][r]])) \n",
    "                # Multiplying elements\n",
    "                \n",
    "            elements = [i.result() for i in elements] \n",
    "            # evaluation for each row\n",
    "            \n",
    "            C[q][r] = sum_elements(inputs=elements) \n",
    "            # Implementing each element as sum of the multiplied elements\n",
    "    return C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "transition_matrix = np.array([[0.5,0.1,0.3],[0.25,0.1,0.1],[0.25,0.8,0.6]])\n",
    "\n",
    "'''\n",
    "This is the transition probability kernel that maps \n",
    "the probability for going from one state to the other\n",
    "'''\n",
    "\n",
    "initial_state = np.array([[5],[6],[1]]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrices = []\n",
    "for num_steps in [50,100,200,300,400,500]: # For different burn-in points\n",
    "    new_matrix = initial_state\n",
    "    \n",
    "    for i in range(num_steps): \n",
    "        '''\n",
    "        Evaluating the matrices sequentially\n",
    "        '''\n",
    "        new_matrix = multiply_matrices_for_markov(transition_matrix,new_matrix)\n",
    "        \n",
    "    matrices.append(new_matrix) \n",
    "    # Collecting the final matrix at the end of the markov process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results = []\n",
    "for matrix in matrices: \n",
    "    \n",
    "    matrix = [[element.result() for element in row] for row in matrix] \n",
    "    # Evaluating the final results\n",
    "    \n",
    "    final_results.append(matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[[4.048192771084336], [1.8072289156626502], [6.14457831325301]],\n",
       " [[4.048192771084336], [1.8072289156626502], [6.14457831325301]],\n",
       " [[4.048192771084336], [1.8072289156626502], [6.14457831325301]],\n",
       " [[4.048192771084336], [1.8072289156626502], [6.14457831325301]],\n",
       " [[4.048192771084336], [1.8072289156626502], [6.14457831325301]],\n",
       " [[4.048192771084336], [1.8072289156626502], [6.14457831325301]]]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "final_results # Six different final states "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Mandelbrot Fractal Set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Mandelbrot Fractal Set is a figure made by the repetition of the same pattern at different scales. Each element is a replica of the same function plotted at a smaller scale. Since the evaluation of the function and conversion to the coordinates can be done in parallel, we can use Parsl to execute them in parallel threads.\n",
    "\n",
    "Code adapted from: https://www.geeksforgeeks.org/mandelbrot-fractal-set-visualization-in-python/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run the Mandelbrot Fractal Set for a single pixel first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "@python_app\n",
    "def mandelbrot(x, y):\n",
    "    \n",
    "    from numpy import complex, array \n",
    "    import colorsys\n",
    "    \n",
    "    '''\n",
    "    Internal function for converstion to the RGB type \n",
    "    tuple that can be added to the image.\n",
    "    '''\n",
    "    \n",
    "    def rgb_conversion(i): \n",
    "        color = 255 * array(colorsys.hsv_to_rgb(i / 255.0, 1.0, 0.5)) \n",
    "        return tuple(color.astype(int)) \n",
    "    \n",
    "    c0 = complex(x, y)\n",
    "    c = 0\n",
    "    for i in range(1, 1000): \n",
    "        if abs(c) > 2:\n",
    "            return rgb_conversion(i) \n",
    "        c = c * c + c0 \n",
    "    return (0, 0, 0) \n",
    "    \n",
    "def mandelbrot_non_parallel(x,y):\n",
    "    \n",
    "    from numpy import complex, array \n",
    "    import colorsys\n",
    "    \n",
    "    '''\n",
    "    Internal function for converstion to the RGB type \n",
    "    tuple that can be added to the image\n",
    "    '''\n",
    "    def rgb_conversion(i): \n",
    "        color = 255 * array(colorsys.hsv_to_rgb(i / 255.0, 1.0, 0.5)) \n",
    "        return tuple(color.astype(int)) \n",
    "    \n",
    "    c0 = complex(x, y)\n",
    "    c = 0\n",
    "    for i in range(1, 1000): \n",
    "        if abs(c) > 2:\n",
    "            return rgb_conversion(i) \n",
    "        c = c * c + c0 \n",
    "    return (0, 0, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(127, 6, 0)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Evaluation for a single pixel.\n",
    "mandelbrot(40,10).result()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now execute the code using the Parsl function for an entire image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Parallel Execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from PIL import Image \n",
    "from numpy import complex, array \n",
    "import colorsys \n",
    "import time\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = time.time()\n",
    "# Setting up the image\n",
    "WIDTH = 256\n",
    "img = Image.new('RGB', (WIDTH, int(WIDTH / 2)))\n",
    "pixels_raw = [[0 for _ in range(img.size[1])] for _ in range(img.size[0])]\n",
    "\n",
    "for x in range(img.size[0]):\n",
    "    for y in range(img.size[1]):\n",
    "        '''\n",
    "        Looping through the image width and height to update the \n",
    "        raw pixels list with the output of the function\n",
    "        '''\n",
    "        a = (x - (0.75 * WIDTH)) / (WIDTH / 4)\n",
    "        b = (y - (WIDTH / 4)) / (WIDTH / 4)\n",
    "        pixels_raw[x][y] = mandelbrot(a, b)\n",
    "\n",
    "pixels_raw = [[i.result() for i in x] for x in pixels_raw]\n",
    "\n",
    "# Projecting the raw pixels onto the image\n",
    "pixels = img.load()\n",
    "\n",
    "for x in range(img.size[0]):\n",
    "    for y in range(img.size[1]):\n",
    "        pixels[x,y] = pixels_raw[x][y]\n",
    "\n",
    "end = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Non-Parallel Execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import time\n",
    "start2 = time.time()\n",
    "# Setting up the image\n",
    "WIDTH = 256\n",
    "img = Image.new('RGB', (WIDTH, int(WIDTH / 2)))\n",
    "pixels_raw = [[0 for _ in range(img.size[1])] for _ in range(img.size[0])]\n",
    "\n",
    "for x in range(img.size[0]):\n",
    "    for y in range(img.size[1]):\n",
    "        \n",
    "        # Looping through the image width and height to update the \n",
    "        # raw pixels list with the output of the function\n",
    "        \n",
    "        a = (x - (0.75 * WIDTH)) / (WIDTH / 4)\n",
    "        b = (y - (WIDTH / 4)) / (WIDTH / 4)\n",
    "        pixels_raw[x][y] = mandelbrot_non_parallel(a, b)\n",
    "\n",
    "# Projecting the raw pixels onto the image\n",
    "pixels = img.load()\n",
    "\n",
    "for x in range(img.size[0]):\n",
    "    for y in range(img.size[1]):\n",
    "        pixels[x,y] = pixels_raw[x][y]\n",
    "\n",
    "end2 = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print('Time taken with Parsl: {}'.format(end-start))\n",
    "print('Time taken with Parsl: {}'.format(end2-start2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the image of the projects of all the function evaluations\n",
    "plt.figure()\n",
    "plt.imshow(img)\n",
    "plt.show()"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}