PowersetGenerationPost_Part3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Power sets in Python: an object-oriented approach (part 3)\n",
    "\n",
    "\n",
    "## Contents\n",
    "<style>\n",
    "  ol\n",
    "  {\n",
    "    list-style-type: none\n",
    "  }\n",
    "</style>\n",
    "<ol>\n",
    "    <li><div style=\"text-align: left\"><a href=\"#Overview-main\"><cn1>Overview</cn1></a></div></li>\n",
    "    <li><div style=\"text-align: left\"><a href=\"#Preliminaries-main\"><cn1>Preliminaries</cn1></a></li>\n",
    "    <li><div style=\"text-align: left\"><a href=\"#Binary-Mask-power-set-generator-main\"><cn1>Binary Mask power set generator</cn1></a></div>\n",
    "    <ol>\n",
    "        <li><div style=\"text-align: left\"><a href=\"#Overview-binary\"><cn2>Overview</cn2></a></div></li>\n",
    "        <li><div style=\"text-align: left\"><a href=\"#Flow-of-control-binary\"><cn2>Flow of control</cn2></a></div>\n",
    "        <ol>\n",
    "            <li><div style=\"text-align: left\"><a href=\"#Aside:-initializing-the-mask-binary\"><cn3>Aside: initializing the mask<cn3></a></div></li>\n",
    "        </ol>\n",
    "        <li><div style=\"text-align: left\"><a href=\"#Testing-the-binary-mask-routine-binary\"><cn2>Testing the binary mask routine</cn2></a></div></li>\n",
    "    </ol>\n",
    "    </li>\n",
    "<li><div style=\"text-align: left\"><a href=\"#Recursive-power-set-generator-main\"><cn1>Recursive power set generator</cn1></a></div>\n",
    "    <ol>\n",
    "    <li><div style=\"text-align: left\"><a href=\"#Overview-recursive\"><cn2>Overview</cn2></a></div></li>\n",
    "    <li><div style=\"text-align: left\"><a href=\"#Flow-of-control-recursive\"><cn2>Flow of control</cn2></a></div></li>\n",
    "    <li><div style=\"text-align: left\"><a href=\"#Testing-the-recursive-routine-recursive\"><cn2>Testing the recursive routine</cn2></a></div></li>\n",
    "    </ol>\n",
    "    </li>\n",
    "    <li><div style=\"text-align: left\"><a href=\"#Quick-and-dirty-efficiency-comparison-main\"><cn1>Quick-and-dirty efficiency comparison</cn1></div></li>\n",
    "    <li><div style=\"text-align: left\"><a href=\"#Summary-main\"><cn1>Summary</cn1></div></li>\n",
    "    <li><div style=\"text-align: left\"><a href=\"#Sources-(part-3)-main\"><cn1>Sources (part 3)</cn1></div></li>\n",
    "    <li><div style=\"text-align: left\"><a href=\"#The-complete-ModSet()-class-main\"><cn1>The complete ModSet() class</cn1></a></div></li>\n",
    "</ol>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Overview-main\"><h2>Overview</h2></div>\n",
    "\n",
    "Last time, in <a href=\"https://nullexit.org/blog/index.php/2021/01/27/power-sets-in-python-an-object-oriented-approach-part-2/\">part two</a>, we built and tested our set container class, `ModSet()`, with a variety of methods that will prove useful when we construct the power set.  Methods included copying, nesting/unnesting, filtering to keep only unique members as well as set-union and set-subtraction operations. You will see that most of these will be used, either implicitly or explicitly, in both Binary Mask and Recursive methods of generating power sets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Preliminaries-main\"><h2>Preliminaries</h2></div>\n",
    "\n",
    "When we are ready to build our power set, we'll measure the time it takes each of the two routines to complete, so we can compare their performance. As a preliminary, let's import the `time` module so it stands at ready."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course we will need an example source set to operate on. Let's define it now and be done with it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clumsyPrez = ModSet(set([ 6.1, 2021, \"fight like hell\", \"we love you\",('impeachment','#2')]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Binary-Mask-power-set-generator-main\"><h2>Binary Mask power set generator</h2></div>\n",
    "\n",
    "In this iterative-based approach, we will populate the power set by applying binary masks to extract out member subsets, one-by-one, from the source set in `ModSet().val`. If you recall, a verbal description of power set generation using a binary mask we discussed way back in part one; you can find it <a href=\"https://nullexit.org/blog/index.php/2021/01/27/power-sets-in-python-an-object-oriented-approach-part-1/#The-Binary-Mask-approach\">here</a>. Below, we first provide an overview of our binary mask power set generator routine, followed by the routine itself, with heavy commenting, concluded with a detailed walk-through after that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Overview-binary\"><h3>Overview</h3></div>\n",
    "\n",
    "The method, named `.powerSet_bin()`, contains 19 lines and calls `.pushDown()` and `.union()` sibling methods.  Though we could use NumPy here to make mask generation and member extraction more compact, we did not want to confound the comparison (to our recursive approach to follow) by incorporating functionality of an external module.  Already, the code has been somewhat compressed by use of generator expressions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Flow-of-control-binary\"><h3>Flow of control</h3></div>\n",
    "\n",
    "The routine below can be divided into two main blocks: variable initialization, and the run-time loop that builds each binary mask for extracting the corresponding subset on each iteration. The loop populates the power set member-by-member and terminates when all $2^n$ members have been added. Because indexing is crucial in this approach, you'll notice heavy use of `range()` and `list()` data-types.  As you'll see later, our recursive approach requires no such crutch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class ModSet():\n",
    "    .\n",
    "    .\n",
    "    .\n",
    "    # Generate powerset via direct, binary mask approach\n",
    "    def powerSet_bin(self):\n",
    "        \n",
    "        ## Initialize local variables ##\n",
    "        S = list(self.val)         # convert to list for indexing\n",
    "        setSize = len(self.val)    # count number of members in source set\n",
    "        psetSize = pow(2, setSize) # calculate the number of elements in the power set\n",
    "        lastIndex = setSize - 1    # index value of last member\n",
    "        setIndices = range(0, setSize)  # make indices list for source set\n",
    "        psetIndices = range(0, psetSize) # make indices list for power set to be built\n",
    "        bMasks = [[False for s in setIndices] for p in psetIndices] # Initialize binary mask       \n",
    "        pSet = ModSet(set())  # initialize power set as empty ModSet() instance\n",
    "        pSet.pushDown(1)  # and nest it down one level for later joining\n",
    "        \n",
    "        ## Populate powerset with each subset, one at a time ##\n",
    "        for i in psetIndices: # loop through each member of power set\n",
    "            \n",
    "            ## Generate binary mask for subset i of power set ##\n",
    "            val = i  # assign current pSet index as current \"value\" of mask \n",
    "            for j in setIndices: # loop through each bit-digit of mask\n",
    "                \n",
    "                if (val >= pow(2, lastIndex - j)):  # if mask value >= value of current bit,\n",
    "                    bMasks[i][lastIndex - j] = True # then set corresp. mask bit to \"true\"\n",
    "                    val -= pow(2, lastIndex - j)    # subtract value of current bit from\n",
    "                                                    # mask value to determine next bit-digit\n",
    "            ## Form subset i of power set ##\n",
    "            # Use generator expression for compactness\n",
    "            dummySet = ModSet(set([S[k] for k in setIndices if bMasks[i][k] == True]))\n",
    "            dummySet.pushDown(1)  # nest ModSet instance down one level for union join\n",
    "            pSet.union(dummySet)  # include new subset in power set\n",
    "            \n",
    "        return pSet, bMasks  # return complete power set and binary masks as output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Aside:-initializing-the-mask-binary\"><h4>Aside: initializing the mask</h4></div>\n",
    "\n",
    "To have each and every element alterable in our list of binary masks, `bMasks`, we need to individually assign each and every element within the nested list. If you attempt to initialize `bMasks` like this, `bMasks = [[False]*len(setIndices)]*len(psetIndices)]`, (as *we* did on our first attempt!), you will find that you cannot change individual elements of `bMasks` later on. That is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[False, True, False],\n",
       " [False, True, False],\n",
       " [False, True, False],\n",
       " [False, True, False],\n",
       " [False, True, False]]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bMasks =[[False]*3]*5\n",
    "bMasks[3][1] = True\n",
    "bMasks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "changes to `True` the second column of *all* rows, rather than just the fourth, since all point to the same instance of `True`.  However, if you initialize `bMasks` using a list completion, the nested list assembles with a hash to a unique instance for each and every element therein"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[False, False, False],\n",
       " [False, False, False],\n",
       " [False, False, False],\n",
       " [False, True, False],\n",
       " [False, False, False]]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bMasks = [[False for j in range(0,3)] for i in range(0, 5)]\n",
    "bMasks[3][1] = True\n",
    "bMasks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so that we can make the individual bit flips necessary for our routine. Our apologies for the digression; we felt the need to address this \"rookie\" Python mistake for those who may be unaccustomed to using `list` objects in this, possibly unusual, manner."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Testing-the-binary-mask-routine-binary\"><h3>Testing the binary mask routine</h3></div>\n",
    "\n",
    "With the binary mask routine complete, we're ready to build our power set.  Below we've sandwiched the call within two `time.time()` reads, so we can measure its runtime duration. Let's examine `bMasks` output first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[False, False, False, False, False],\n",
       " [True, False, False, False, False],\n",
       " [False, True, False, False, False],\n",
       " [True, True, False, False, False],\n",
       " [False, False, True, False, False],\n",
       " [True, False, True, False, False],\n",
       " [False, True, True, False, False],\n",
       " [True, True, True, False, False],\n",
       " [False, False, False, True, False],\n",
       " [True, False, False, True, False],\n",
       " [False, True, False, True, False],\n",
       " [True, True, False, True, False],\n",
       " [False, False, True, True, False],\n",
       " [True, False, True, True, False],\n",
       " [False, True, True, True, False],\n",
       " [True, True, True, True, False],\n",
       " [False, False, False, False, True],\n",
       " [True, False, False, False, True],\n",
       " [False, True, False, False, True],\n",
       " [True, True, False, False, True],\n",
       " [False, False, True, False, True],\n",
       " [True, False, True, False, True],\n",
       " [False, True, True, False, True],\n",
       " [True, True, True, False, True],\n",
       " [False, False, False, True, True],\n",
       " [True, False, False, True, True],\n",
       " [False, True, False, True, True],\n",
       " [True, True, False, True, True],\n",
       " [False, False, True, True, True],\n",
       " [True, False, True, True, True],\n",
       " [False, True, True, True, True],\n",
       " [True, True, True, True, True]]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tStart = time.time() # clock start timestamp\n",
    "pSet_bin, bMasks = clumsyPrez.powerSet_bin() # call our binary mask power set generator\n",
    "duration = time.time() - tStart # calc run duration by subtracting tStart from current time.\n",
    "bMasks # show list of binary masks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, you can see that masks progress from all `False` to all `True` in a logical pattern. This ordering would be a great feature if we cared about how members are ordered in our power set. But sets, strictly defined, are not distinguished by the ordering of their members <a href=\"#ref1\">(Van Dalen, Doets, De Swart; 2014)</a>.  That is, set $\\{a, b, c\\}$ is equivalent to set $\\{c, a, b\\}$ and $\\{b, a, c\\}$, and so on.\n",
    "\n",
    "So, after taking a look at our power set in `pSet_bin`, we see that the nice ordering was all for naught."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Binary Mask approach took 0.008102 seconds to complete.\n",
      "The power set contains 32 subset elements\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{set(),\n",
       " {'fight like hell', 'we love you', ('impeachment', '#2'), 2021},\n",
       " {'fight like hell', 'we love you', ('impeachment', '#2'), 6.1},\n",
       " {'fight like hell', 'we love you', ('impeachment', '#2')},\n",
       " {'fight like hell', 'we love you', 2021, 6.1},\n",
       " {'fight like hell', 'we love you', 2021},\n",
       " {'fight like hell', 'we love you', 6.1},\n",
       " {'fight like hell', 'we love you'},\n",
       " {'fight like hell', ('impeachment', '#2'), 2021, 6.1},\n",
       " {'fight like hell', ('impeachment', '#2'), 2021},\n",
       " {'fight like hell', ('impeachment', '#2'), 6.1},\n",
       " {'fight like hell', ('impeachment', '#2')},\n",
       " {'fight like hell', 2021, 6.1, 'we love you', ('impeachment', '#2')},\n",
       " {'fight like hell', 2021, 6.1},\n",
       " {'fight like hell', 2021},\n",
       " {'fight like hell', 6.1},\n",
       " {'fight like hell'},\n",
       " {'we love you', ('impeachment', '#2'), 2021, 6.1},\n",
       " {'we love you', ('impeachment', '#2'), 2021},\n",
       " {'we love you', ('impeachment', '#2'), 6.1},\n",
       " {'we love you', ('impeachment', '#2')},\n",
       " {'we love you', 2021, 6.1},\n",
       " {'we love you', 2021},\n",
       " {'we love you', 6.1},\n",
       " {'we love you'},\n",
       " {('impeachment', '#2'), 2021, 6.1},\n",
       " {('impeachment', '#2'), 2021},\n",
       " {('impeachment', '#2'), 6.1},\n",
       " {('impeachment', '#2')},\n",
       " {2021, 6.1},\n",
       " {2021},\n",
       " {6.1}}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Report duration of binary mask power set generation\n",
    "print('The Binary Mask approach took %0.6f seconds to complete.'%(duration))\n",
    "print('The power set contains %i subset elements'%(len(pSet_bin.val)))\n",
    "pSet_bin.val # show power set output; remember, this is a set!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the power set output string above, it looks like we have many sub-`set` instances as members of one big super-`set`.  But, if you recall from all the way back in part one, this is *not* the case; our power set `pSet_bin` is actually a `ModSet()` object, whose `.val` attribute is a `set`, each element of which is another `ModSet()` instance that contains one unique sub-`set` of the original source set (as its `.val` attribute). We defined the `ModSet.__repr__()` method to return the code representation of its `.val` attribute; that's the reason why the return string from calling `pSet_bin.val` appears in this set-resembling format. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Recursive-power-set-generator-main\"><h2>Recursive power set generator</h2></div>\n",
    "In part one, we went over a recursive algorithm that we used to design our recursive power set generator below. You can find it <a href=\"https://nullexit.org/blog/index.php/2021/01/27/power-sets-in-python-an-object-oriented-approach-part-1/#A-Recursive-approach\">here</a> if you would like a refresher."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Overview-recursive\"><h3>Overview</h3></div>\n",
    "\n",
    "Our recursion-based method for generating power sets contains just 11 lines of code but requires five sibling methods, as well as itself, to run. Though `.powerSet_rec()` has fewer lines than `.powerSet_bin()`, it relies more heavily upon behavior of sibling methods in the `ModSet()` class, i.e. `.__copy__()`, `.pushDown()`, `.diffFunc()`, `.union()` and `.removeDuplicates()` to perform the respective obligatory processing tasks of duplication, nesting, extraction, joining and filtering. Unlike the Binary Mask approach above, indexing is not needed here, and hence, not employed. Through recursion, we can generate the identical power set using thoughtfully-placed set-subtraction, -union and -filtering operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class ModSet():\n",
    "    .\n",
    "    .\n",
    "    .\n",
    "    # Generate power set recursively.\n",
    "    def powerSet_rec(self):\n",
    "        \n",
    "        pSet = self.__copy__()      # Preserve self instance; its copy, pSet\n",
    "                                    # will be altered\n",
    "        pSet.pushDown(1)            # Nest pSet for later joining.\n",
    "              \n",
    "        if len(self.val) > 0:       # Recursion termination condition\n",
    "            for elSet in self.val:  # Iterate through members of set self.val\n",
    "                # Generate subset that remains after removing current\n",
    "                # element, elSet, from set self.val\n",
    "                dummySet = self.diffFunc(ModSet(set([elSet]))) \n",
    "                # To current powerset, append the powerset of the\n",
    "                # subset from previous step\n",
    "                pSet.union(dummySet.powerSet_rec())  # Self-call power set method,\n",
    "                                                     # union join power set of \n",
    "                                                     # dummySet with pSet\n",
    "            return pSet             # Return power set at current \n",
    "                                    # level of recursion\n",
    "        else:\n",
    "            dummySet = ModSet(set())  # Generate instance of ModSet of empty set\n",
    "            dummySet.pushDown(1)      # Nest empty set down one level so it can\n",
    "            return dummySet           # be union-joined in the recursion level\n",
    "                                      # above (that called this current run)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Flow-of-control-recursive\"><h3>Flow of control</h3></div>\n",
    "\n",
    "First a duplicate of the calling instance is made to serve as the (local) power set within our routine. This instance is promptly pushed down one level so it can be joined by union later, if necessary. Next, the number of elements in the calling instance are counted.  If empty, an empty `ModSet()` instance is returned on exit of the routine. This case path is the exit condition for the recursive flow; eventually all calling instances will dwindle down to empty as members are stripped from them in the alternative, non-empty, case that we'll describe next.  \n",
    "\n",
    "If the calling instance is not empty, `.powerSet_rec()` iterates over the elements therein, each time subtracting out the current element, `elSet`, and calling itself again to build the power set from elements that remain. Output from this recursive call is then joined by union with `pSet`.  Remember that every time `pSet.union()` runs, it calls the filtering method `removeDuplicates()` to retain only unique members in `pSet` as it is assembled.  \n",
    "\n",
    "When the loop completes, and all subsets have been included at the present level, the local instance of `pSet` is returned so that it can be union-joined with the `pSet` instance inside the calling-function one level above. The method continues in this fashion until `pSet` at the top-most level of recursion is fully populated, at which point output is returned in response to the initial method call.\n",
    "\n",
    "\\[Phew..!\\]\n",
    "\n",
    "Let's give it a whirl and gauge its run-time duration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Testing-the-recursive-routine-recursive\"><h3>Testing the recursive routine</h3></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Recursive approach took 0.006932 seconds to complete.\n",
      "The power set contains 32 subset elements\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{set(),\n",
       " {'fight like hell', 'we love you', ('impeachment', '#2'), 2021},\n",
       " {'fight like hell', 'we love you', ('impeachment', '#2'), 6.1},\n",
       " {'fight like hell', 'we love you', ('impeachment', '#2')},\n",
       " {'fight like hell', 'we love you', 2021, 6.1},\n",
       " {'fight like hell', 'we love you', 2021},\n",
       " {'fight like hell', 'we love you', 6.1},\n",
       " {'fight like hell', 'we love you'},\n",
       " {'fight like hell', ('impeachment', '#2'), 2021, 6.1},\n",
       " {'fight like hell', ('impeachment', '#2'), 2021},\n",
       " {'fight like hell', ('impeachment', '#2'), 6.1},\n",
       " {'fight like hell', ('impeachment', '#2')},\n",
       " {'fight like hell', 2021, 6.1, 'we love you', ('impeachment', '#2')},\n",
       " {'fight like hell', 2021, 6.1},\n",
       " {'fight like hell', 2021},\n",
       " {'fight like hell', 6.1},\n",
       " {'fight like hell'},\n",
       " {'we love you', ('impeachment', '#2'), 2021, 6.1},\n",
       " {'we love you', ('impeachment', '#2'), 2021},\n",
       " {'we love you', ('impeachment', '#2'), 6.1},\n",
       " {'we love you', ('impeachment', '#2')},\n",
       " {'we love you', 2021, 6.1},\n",
       " {'we love you', 2021},\n",
       " {'we love you', 6.1},\n",
       " {'we love you'},\n",
       " {('impeachment', '#2'), 2021, 6.1},\n",
       " {('impeachment', '#2'), 2021},\n",
       " {('impeachment', '#2'), 6.1},\n",
       " {('impeachment', '#2')},\n",
       " {2021, 6.1},\n",
       " {2021},\n",
       " {6.1}}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tStart = time.time()\n",
    "pSet_rec = clumsyPrez.powerSet_rec()\n",
    "duration = time.time() - tStart\n",
    "print('The Recursive approach took %0.6f seconds to complete.'%(duration))\n",
    "print('The power set contains %i subset elements'%(len(pSet_rec.val)))\n",
    "pSet_rec.val\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As promised, we see the same power set as we saw from the Binary Mask approach. But unlike that approach, there is no mask output for us to examine this time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Quick-and-dirty-efficiency-comparison-main\"><h2>Quick-and-dirty efficiency comparison</h2></div>\n",
    "\n",
    "For a given programming objective, recursive algorithms tend to be less efficient than their iterative equivalents when used in imperative-based languages, like Python, where iteration is preferred <a href=\"#ref2\">(Recursion (computer science), Wikipedia)</a>.  Below we hobbled together a quick script to measure, and statistically compare, processing times of the two methods of power set generation of `clumsyPrez`. We instruct the two power set generators to each run 500 times in a randomly alternating sequence to destroy potential processing-related biases (e.g. sequential, timing) that could otherwise arise by running the two in separate blocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# routine to repeatedly collect run-time durations\n",
    "# of functions genFunc1 and genFunc2\n",
    "def genDurationsDist(genFunc1, genFunc2, nReps):\n",
    "    durations = np.zeros(2*nReps)  # initialize output array\n",
    "    genFuncs = [genFunc1, genFunc2] # make list of the two functions to run\n",
    "    # initialize an array to contain the run-sequence of the two functions\n",
    "    funcSeq = np.concatenate([np.zeros(nReps), np.ones(nReps)]).astype(int)\n",
    "    np.random.shuffle(funcSeq)  # shuffle the ordering of the sequence\n",
    "    for i in np.arange(0, funcSeq.shape[0]): # iterate over random sequence\n",
    "        tStart = time.time() # start timestamp\n",
    "        genFuncs[funcSeq[i]]() # run one of the two functions based on funcSeq\n",
    "        durations[i] = time.time() - tStart # end timestamp, calc duration\n",
    "    # separate and return the two sets of durations\n",
    "    return durations[funcSeq == 0], durations[funcSeq == 1] \n",
    "\n",
    "durations_bin, durations_rec = genDurationsDist(clumsyPrez.powerSet_bin, \n",
    "                                                clumsyPrez.powerSet_rec, 500)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# cell for reloading nonparametric_stats module\n",
    "import importlib\n",
    "importlib.reload(nps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAFNCAYAAACUvLFdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfZxWdZ3/8ddbYAUVUWB0FVAw0QQVsJFkrU0jk7wDTVcqFUvTytzql5lkrWzpZuVm2rIWeQOmiSzeUWlCeLd5R6CICJqsYo4gICqCigp8fn+c7+DFODcXMGeumTnv5+Mxjznne+4+17muud5z7hURmJmZFcE2lS7AzMyspTj0zMysMBx6ZmZWGA49MzMrDIeemZkVhkPPzMwKw6FnjZL0K0k/qNCyQ9LeLbi8uySNaanlNUXSoZKelbRG0qhK12PWHjj02hhJiyW9nb4IX5Y0UdIOzTTv0yX9pbQtIr4SET9qjvm3JpLGSbqhtC0iPhMRkypVUz1+CPxXROwQEbeXDpC0raRrJL0gabWkxyV9pkJ1VlT6G7i40nU0h/o+l9a8HHpt07ERsQMwGBgCjK1wPa2KpI6VrqGZ7Ak81cCwjsCLwCeAbsAPgCmS+rZIZVupHb1HZWuJ11zE9brZIsI/begHWAx8qqT/p8AfS/rvA84s6T8d+EtJfwBfAZ4FXgPGAwL2A9YC64E1wOtp/InAxan7MKAGOB9YDiwFRgFHAX8DXgW+V7KsbYALgP8DVgJTgO6NvLbvpHkuAb6Uat17M17XOel1PZ/ariALhjeAOcDHU/sI4F3gvfRan6i7jFT794EX0mu9HuiWhvVNyxsD/B14BbiwpJahwOy03GXAzxt5zV8GFqV1Nw3YPbX/H7ABeDvVuG0Zn415wGcbGHY68CDwS2AV8DQwvGT47mn5r6Z6vpzaO6caeqb+7wPrgB1T/8XAL1L3tsBlaZ0sA34FdKnz2fku8DLw23pq3Bu4P9X3CnBzybAPAzNSfc8A/5Laz0rv47tpPf2+gdf/6TTdKuC/03JKP09fAhaS/U3cDezZ1N/MZky7tZ/Let+bNGwcMBW4Ic3vTDbj81fEn4oX4J/NfMNKQg/oDTwJXFEy/D6aDoc/ADsBewArgBH1jZvaJrJp6K0D/g3oRPaFvQL4HdAVGEgWnHul8b8JPJLq3Bb4NXBTA69rRPoD3R/YPs1zc0NvBtCd979oTwF6kG0VfZvsy7ZzGjYOuKFODRuXkb7IFgF7ATsAt5K+qHk/9H4DdAEGAe8A+6XhDwOnpu4dgEMaeM2fJPtyPyitn18CD9T3Xpfxudg1rfsPNzD89PTefSu9dyeTBUD3NPx+sjDoTLYHYQUpFIEHSGEKTCcL5M+UDDs+df+C7Mu5e/o8/B74cZ3Pzk/Sa+1ST403AReS/cPRGfhYat+eLCS+mN7Lg9J6G1j3M9rAa+9JFgAnpOm/QRYste/1qPRe75eGfx94qMy/mXKm3drPZWPvzbj0Wkal9daFMj9/Rf2peAH+2cw3LPsiXAOsTn9QM4GdSobfR9Ph8LGS/inABfWNm9o2fqGkL663gQ6pv2ua30dLxp8DjErdC9l0a2K39AfasZ7XdS1waUn/Pmx+6H2yiXX3GjAoddf35bJxGWm9fq1k2L61tfN+6PUuGT4LGJ26HwD+nbR11Eg91wA/LenfIS2jb8l73WTokYXYn4FfNzLO6WRb0KpT86lAH7It/K4lw34MTEzdPwKuTK/9ZbLQuJSSrUCyvQVvAh8qmccw3t+6OYxsK6ZzIzVeD0woXa+p/WTgf+u0/Rq4qO5ntIH5ngY8XNIvshCtfa/vAs4oGb4N8BZpi43G/2bKmXaLP5dlvDfjKPlHaXM+f0X98TG9tmlURHQl+yL5MNmXzuZ4uaT7LbIv23KtjIj1qfvt9HtZyfC3S+a3J3CbpNclvU4WguvJtkrq2p3si6jWC5tRU63S6ZH0bUkLJa1Ky+9G+etq9zo1vED2pV9ae0Pr8Qyy0H5a0l8lHVPOMiJiDdlu4F5l1oikbYDfkgXK15sY/aVI34rJC6mG3YFXI2J1nWG1ddxP9lk7iGzPwgyyY4mHAIsi4hWgCtgOmFPyfv8ptddaERFrG6nvfLJAmiXpKUlfSu17Ah+tnW+a9xeAf2zi9dba5LOV1kFNyfA9gStK5v1qqqP0fWjovS5n2q35XDb13nxg/pT/+SskH/RswyLifkkTyY6j1J7S/ibZl0+tcr8YIPuvtDm9CHwpIh4sY9ylZP/V1tqjzvByXtfG+iV9nOz40XDgqYjYIOk1si+kTcZtwBKyL7TSetaRBXzvxiaMiGeBz6VAOgGYKqlHRLzZ2DIkbU+22+ulJmqrHV9kW4u7AkdFxHtNTNJLkkqCbw+y3ZFLgO6SupZ8ue5RUsdDZFu6xwP3R8QCSXsAR5MFImS7G98m2+XYUP2NrvOIeJlslzmSPgb8WdIDZJ+j+yPiiC2ZL9lna+N7ltZb6Xv4InBJRNzYxHzqU860W/O5bOq9+cA0m/H5KyRv6bV9vwCOkDQ49c8FTpC0XbrG7YzNmNcyoLekf2im2n4FXCJpTwBJVZJGNjDuFOB0SQMkbQdcVGf45r6urmQhtQLoKOnfgB1Lhi8D+qYvhvrcBHxLUr90Sch/kJ1Ysa6J5SLpFElVEbEBeD01r69n1N8BX5Q0WNK2aRmPRsTippaRXEV2LOnYiHi7qZGBXYB/ldRJ0klp2jsj4kWyYPuxpM6SDiRbvzcCRMRbZLutz+H9kHsIOLu2P73W3wCXS9olrYdeko4s87Ug6SRJtWH0GtmX+Xqy42n7SDo11d5J0sGS9kvjLiM79tqQPwIHSBqVzm48h03/afoVMFbSwFRHt7R+yrG5027W57Kp96Y+m/H5KySHXhsXESvIjoXUXkB+OdmurmXAJBr546jHPWSnyL8s6ZVmKO8Ksi2J6ZJWk53U8tH6RoyIu8gC/B6yEwPuqTPK5r6uu8mOt/yNbHfQWjbdDfQ/6fdKSY/VM/21ZLsNHwCeT9Of28Qya40AnpK0hmwdjK5vt15EzCR7324h2xr5EDC6nAWkfyTOJjux4eV03eYaSV9oZLJHgf5kW2WXACdGxMo07HNkxyqXALeRHS+bUTLt/WTHDmeV9HclWz+1vkv23j0i6Q2y44z7lvN6koOBR9N6mwZ8IyKeT1s4nyZbN0vIdjXWnhAD2dbugLSL8fa6M027X08iO9N5JTCA7OzGd9Lw29L8Jqe65wNlXfO4BdNuyeeyqfemrrI+f0WlTXfxm1l7JOl0shM3PlbpWiotbUXVAF+IiHsrXY+1LG/pmVm7J+lISTul3cjfIzuG9kiFy7IKcOiZWREMI7u+8BXgWLIzoMs5DmrtjHdvmplZYeS2pZfONJol6Yl0zc2/p/aJkp6XNDf9DE7tknSlpEWS5kk6KK/azMysmPK8Tu8dsjsRrJHUCfiLpLvSsO9ExNQ643+G7Myy/mRn+F1FA2f6mZmZbYncQi9dALsm9XZKP43tSx0JXJ+meyQddN4tIpY2NEHPnj2jb9++zVWymZm1A3PmzHklIqrqG5brHVkkdSC7qHVvYHxEPCrpq2QXLP8b2f0NL4iId8huq1N6vUpNamsw9Pr27cvs2bNzq9/MzNoeSQ3exjDXszcjYn1EDCa75c9QSfuTPfvtw2QXonYnu6AV3r8NzyazqNsg6SxJsyXNXrFiRU6Vm5lZe9QilyxExOtkd7AfERFLI/MOcB3Zs58g27Irvfdib7I7ENSd14SIqI6I6qqqerdezczM6pXn2ZtVknZK3V2AT5Hd9Xu31CaymyTPT5NMA05LZ3EeAqxq7HiemZnZ5srzmN5uwKR0XG8bYEpE/EHSPZKqyHZnziV7IjHAnWRP4F5E9uiOL+ZYm5lZi3vvvfeoqalh7VrfCrM5dO7cmd69e9OpU6eyp8nz7M15wJB62j/ZwPhBdvdzM7N2qaamhq5du9K3b1+ynV22pSKClStXUlNTQ79+/cqezrchMzNrIWvXrqVHjx4OvGYgiR49emz2VrNDz8ysBTnwms+WrEuHnplZgXTo0IHBgwczaNAgDjroIB566CEAlixZwoknnrjF8x03bhySWLRo0ca2yy+/HElbdD316aefztSpdW/ctfVyvTjdzMwadvmMvzXr/L51xD5NjtOlSxfmzp0LwN13383YsWO5//772X333bc6ZA444AAmT57M97//fQCmTp3KgAEDtmqezc1bemZmBfXGG2+w8847A7B48WL2339/ACZOnMgJJ5zAiBEj6N+/P+eff35Z8xs1ahR33HEHAM899xzdunWj9Hrqr371q1RXVzNw4EAuuuiije0XXHABAwYM4MADD+S88877wHx/8IMfcPrpp7Nhw4Ytfq21vKVnZlYgb7/9NoMHD2bt2rUsXbqUe+65p97x5s6dy+OPP862227Lvvvuy7nnnkufPn3qHbfWjjvuSJ8+fZg/fz533HEHJ598Mtddd93G4Zdccgndu3dn/fr1DB8+nHnz5tG7d29uu+02nn76aSTx+uuvbzLP888/n1WrVnHdddc1y/FQhx6N72IoZ3eBmVlbUbp78+GHH+a0005j/vz5Hxhv+PDhdOvWDYABAwbwwgsvNBl6AKNHj2by5MncfffdzJw5c5PQmzJlChMmTGDdunUsXbqUBQsWMGDAADp37syZZ57J0UcfzTHHHLNx/B/96Ed89KMfZcKECVv7sjfy7k0zs4IaNmwYr7zyCvXdx3jbbbfd2N2hQwfWrVtX1jyPPfZYfvvb37LHHnuw4447bmx//vnnueyyy5g5cybz5s3j6KOPZu3atXTs2JFZs2bx2c9+lttvv50RI0ZsnObggw9mzpw5vPrqq1vxKjflLT0zs4J6+umnWb9+PT169OCtt95qlnl26dKFn/zkJ+yzz6Z7yd544w223357unXrxrJly7jrrrs47LDDWLNmDW+99RZHHXUUhxxyCHvvvffGaUaMGMGRRx7J0UcfzfTp0+natetW1+fQMzMrkNpjepDd1WTSpEl06NCh7OnPPPNMvvKVr1BdXd3gOKNHj/5A26BBgxgyZAgDBw5kr7324tBDDwVg9erVjBw5krVr1xIRXH755ZtMd9JJJ7F69WqOO+447rzzTrp06VJ2rfVRdvevtqm6ujqa43l6PqZnZi1h4cKF7LfffpUuo12pb51KmhMR9aayj+mZmVlhOPTMzKwwHHpmZlYYDj0zMysMh56ZmRWGQ8/MzArDoWdmViC1jxbaf//9OfbYYz9wr8uWcuaZZ7JgwYIWX64vTjczq5R7f9y88zt8bJOjlN57c8yYMYwfP54LL7yweetI1q9f3+CF71dffXUuy2yKt/TMzApq2LBhvPTSSxv7f/azn3HwwQdz4IEHbvLon+uvv54DDzyQQYMGceqppwIffMjrDjvsAMB9993H4Ycfzuc//3kOOOAA3nzzTY4++mgGDRrE/vvvz8033wzAYYcdxuzZs7nqqqs2eXTRxIkTOffccwG44YYbGDp0KIMHD+bss89m/fr1W/2avaVnZlZA69evZ+bMmZxxxhkATJ8+nWeffZZZs2YRERx33HE88MAD9OjRg0suuYQHH3yQnj17lnXz51mzZjF//nz69evHLbfcwu67784f//hHAFatWrXJuCeeeCLDhg3jpz/9KQA333wzF154IQsXLuTmm2/mwQcfpFOnTnzta1/jxhtv5LTTTtuq1+3QMzMrkNp7by5evJiPfOQjHHHEEUAWetOnT2fIkCEArFmzhmeffZYnnniCE088kZ49ewLQvXv3JpcxdOhQ+vXrB2RPUz/vvPP47ne/yzHHHMPHP/7xTcatqqpir7324pFHHqF///4888wzHHrooYwfP545c+Zw8MEHb6x7l1122erX792bZmYFUntM74UXXuDdd99l/PjxQHbz6bFjxzJ37lzmzp3LokWLOOOMM4iIeh/e2rFjx41PMo8I3n333Y3Dtt9++43d++yzD3PmzOGAAw5g7Nix/PCHP/zAvE4++WSmTJnCLbfcwvHHH48kIoIxY8ZsrOeZZ55h3LhxW/36HXpmZgXUrVs3rrzySi677DLee+89jjzySK699lrWrFkDwEsvvcTy5csZPnw4U6ZMYeXKlQAbd2/27duXOXPmAHDHHXfw3nvv1bucJUuWsN1223HKKadw3nnn8dhjj31gnBNOOIHbb7+dm266iZNPPhnIHmI7depUli9fvnG5L7zwwla/bu/eNDMrqCFDhjBo0CAmT57MqaeeysKFCxk2bBiQnZhyww03MHDgQC688EI+8YlP0KFDB4YMGcLEiRP58pe/zMiRIxk6dCjDhw/fZOuu1JNPPsl3vvMdttlmGzp16sRVV131gXF23nlnBgwYwIIFCxg6dCiQPa394osv5tOf/jQbNmygU6dOjB8/nj333HOrXrMfLYQfLWRmLcOPFmp+frSQmZlZAxx6ZmZWGA49MzMrDIeemVkLasvnUbQ2W7Iucws9SZ0lzZL0hKSnJP17au8n6VFJz0q6WdI/pPZtU/+iNLxvXrWZmVVC586dWblypYOvGUQEK1eupHPnzps1XZ6XLLwDfDIi1kjqBPxF0l3A/wMuj4jJkn4FnAFclX6/FhF7SxoN/AQ4Ocf6zMxaVO/evampqWHFihWVLqVd6Ny5M717996saXILvcj+lVmTejulnwA+CXw+tU8CxpGF3sjUDTAV+C9JCv9LZGbtRKdOnTbenssqI9djepI6SJoLLAdmAP8HvB4R69IoNUCv1N0LeBEgDV8F9MizPjMzK5ZcQy8i1kfEYKA3MBSo76rM2i25D97c7f1hG0k6S9JsSbO9i8DMzDZHi5y9GRGvA/cBhwA7SardrdobWJK6a4A+AGl4N+ADz7CIiAkRUR0R1VVVVXmXbmZm7UieZ29WSdopdXcBPgUsBO4FTkyjjQHuSN3TUj9p+D0+nmdmZs0pz7M3dwMmSepAFq5TIuIPkhYAkyVdDDwOXJPGvwb4raRFZFt4o3OszczMCijPszfnAUPqaX+O7Phe3fa1wEl51WNmZuY7spiZWWE49MzMrDAcemZmVhgOPTMzKwyHnpmZFYZDz8zMCsOhZ2ZmheHQMzOzwnDomZlZYTj0zMysMBx6ZmZWGA49MzMrDIeemZkVhkPPzMwKw6FnZmaF4dAzM7PCcOiZmVlhOPTMzKwwHHpmZlYYDj0zMysMh56ZmRWGQ8/MzArDoWdmZoXh0DMzs8Jw6JmZWWE49MzMrDAcemZmVhgOPTMzKwyHnpmZFUZuoSepj6R7JS2U9JSkb6T2cZJekjQ3/RxVMs1YSYskPSPpyLxqMzOzYuqY47zXAd+OiMckdQXmSJqRhl0eEZeVjixpADAaGAjsDvxZ0j4RsT7HGs3MrEBy29KLiKUR8VjqXg0sBHo1MslIYHJEvBMRzwOLgKF51WdmZsXTIsf0JPUFhgCPpqavS5on6VpJO6e2XsCLJZPV0HhImpmZbZbcQ0/SDsAtwDcj4g3gKuBDwGBgKfCftaPWM3nUM7+zJM2WNHvFihU5VW1mZu1RrqEnqRNZ4N0YEbcCRMSyiFgfERuA3/D+LswaoE/J5L2BJXXnGRETIqI6IqqrqqryLN/MzNqZPM/eFHANsDAifl7SvlvJaMcD81P3NGC0pG0l9QP6A7Pyqs/MzIonz7M3DwVOBZ6UNDe1fQ/4nKTBZLsuFwNnA0TEU5KmAAvIzvw8x2dumplZc8ot9CLiL9R/nO7ORqa5BLgkr5rMzKzYfEcWMzMrDIeemZkVhkPPzMwKw6FnZmaF4dAzM7PCcOiZmVlhOPTMzKwwHHpmZlYYDj0zMysMh56ZmRWGQ8/MzArDoWdmZoXh0DMzs8Jw6JmZWWE49MzMrDDyfIis2ea798cNDzt8bMvVYWbtkrf0zMysMBx6ZmZWGA49MzMrDIeemZkVhkPPzMwKw6FnZmaF4dAzM7PCcOiZmVlhOPTMzKwwHHpmZlYYDj0zMysMh56ZmRWGQ8/MzArDoWdmZoWRW+hJ6iPpXkkLJT0l6RupvbukGZKeTb93Tu2SdKWkRZLmSToor9rMzKyY8tzSWwd8OyL2Aw4BzpE0ALgAmBkR/YGZqR/gM0D/9HMWcFWOtZmZWQHlFnoRsTQiHkvdq4GFQC9gJDApjTYJGJW6RwLXR+YRYCdJu+VVn5mZFU+LHNOT1BcYAjwK7BoRSyELRmCXNFov4MWSyWpSm5mZWbPIPfQk7QDcAnwzIt5obNR62qKe+Z0labak2StWrGiuMs3MrAByDT1JncgC78aIuDU1L6vdbZl+L0/tNUCfksl7A0vqzjMiJkREdURUV1VV5Ve8mZm1O3mevSngGmBhRPy8ZNA0YEzqHgPcUdJ+WjqL8xBgVe1uUDMzs+bQMcd5HwqcCjwpaW5q+x5wKTBF0hnA34GT0rA7gaOARcBbwBdzrM3MzAqorNBL19hdB6wGriY7KeWCiJje0DQR8RfqP04HMLye8QM4p5x6zMzMtkS5uze/lE5C+TRQRbYVdmluVZmZmeWg3NCr3WI7CrguIp6g4a04MzOzVqnc0JsjaTpZ6N0tqSuwIb+yzMzMml+5J7KcAQwGnouItyT1wCeamJlZG1Pult6MiHgsIl4HiIiVwOX5lWVmZtb8Gt3Sk9QZ2A7omZ6GUHscb0dg95xrMzMza1ZN7d48G/gmWcDN4f3QewMYn2NdZmZmza7R0IuIK4ArJJ0bEb9soZrMzMxyUdaJLBHxS0n/BPQtnSYirs+pLjMzs2ZX7h1Zfgt8CJgLrE/NATj0zMyszSj3koVqYEC6VZiZmVmbVO4lC/OBf8yzEDMzs7yVu6XXE1ggaRbwTm1jRByXS1VmZmY5KDf0xuVZhJmZWUso9+zN+/MuxMzMLG/lnr25muxsTYB/ADoBb0bEjnkVZmZm1tzK3dLrWtovaRQwNJeKzMzMclLu2ZubiIjbgU82cy1mZma5Knf35gklvduQXbfna/bMzKxNKffszWNLutcBi4GRzV6NmZlZjso9pucHxpqZWZtX1jE9Sb0l3SZpuaRlkm6R1Dvv4szMzJpTuSeyXAdMI3uuXi/g96nNzMyszSg39Koi4rqIWJd+JgJVOdZlZmbW7MoNvVcknSKpQ/o5BViZZ2FmZmbNrdzQ+xLwL8DLwFLgRMAnt5iZWZtS7iULPwLGRMRrAJK6A5eRhaGZmVmbUO6W3oG1gQcQEa8CQ/IpyczMLB/lht42knau7UlbeuVuJZqZmbUK5YbefwIPSfqRpB8CDwE/bWwCSdem6/rml7SNk/SSpLnp56iSYWMlLZL0jKQjt+TFmJmZNabcO7JcL2k22U2mBZwQEQuamGwi8F/A9XXaL4+Iy0obJA0ARgMDya4F/LOkfSJifTn1mZmZlaPsXZQp5JoKutLxH5DUt8zRRwKTI+Id4HlJi8geXfRwucszMzNryhY9WmgrfV3SvLT7s/Y4YS/gxZJxalKbmZlZs2np0LsK+BAwmOx6v/9M7apn3HofXSTpLEmzJc1esWJFPlWamVm71KKhFxHLImJ9RGwAfsP7T1+vAfqUjNobWNLAPCZERHVEVFdV+U5oZmZWvhYNPUm7lfQeD9Se2TkNGC1pW0n9gP7ArJaszczM2r/crrWTdBNwGNBTUg1wEXCYpMFkuy4XA2cDRMRTkqaQnSizDjjHZ26amVlzyy30IuJz9TRf08j4lwCX5FWPmZlZJc7eNDMzqwiHnpmZFYZDz8zMCsOhZ2ZmheHQMzOzwnDomZlZYTj0zMysMBx6ZmZWGA49MzMrDIeemZkVhkPPzMwKw6FnZmaF4dAzM7PCcOiZmVlhOPTMzKwwHHpmZlYYDj0zMysMh56ZmRWGQ8/MzArDoWdmZoXh0DMzs8Jw6JmZWWE49MzMrDAcemZmVhgOPTMzKwyHnpmZFYZDz8zMCsOhZ2ZmheHQMzOzwnDomZlZYeQWepKulbRc0vyStu6SZkh6Nv3eObVL0pWSFkmaJ+mgvOoyM7PiynNLbyIwok7bBcDMiOgPzEz9AJ8B+qefs4CrcqzLzMwKKrfQi4gHgFfrNI8EJqXuScCokvbrI/MIsJOk3fKqzczMiqmlj+ntGhFLAdLvXVJ7L+DFkvFqUpuZmVmzaS0nsqietqh3ROksSbMlzV6xYkXOZZmZWXvS0qG3rHa3Zfq9PLXXAH1KxusNLKlvBhExISKqI6K6qqoq12LNzKx9aenQmwaMSd1jgDtK2k9LZ3EeAqyq3Q1qZmbWXDrmNWNJNwGHAT0l1QAXAZcCUySdAfwdOCmNfidwFLAIeAv4Yl51mZlZceUWehHxuQYGDa9n3ADOyasWMzMzaD0nspiZmeXOoWdmZoXh0DMzs8Jw6JmZWWE49MzMrDAcemZmVhgOPTMzKwyHnpmZFYZDz8zMCsOhZ2ZmheHQMzOzwnDomZlZYeR2w2mzFnfvjxsedvjYlqvDzFotb+mZmVlhOPTMzKwwHHpmZlYYDj0zMysMh56ZmRWGz960tqOxszPNzMrgLT0zMysMh56ZmRWGQ8/MzArDoWdmZoXh0DMzs8Jw6JmZWWE49MzMrDAcemZmVhgOPTMzKwyHnpmZFUZFbkMmaTGwGlgPrIuIakndgZuBvsBi4F8i4rVK1GdmZu1TJbf0Do+IwRFRnfovAGZGRH9gZuo3MzNrNq1p9+ZIYFLqngSMqmAtZmbWDlUq9AKYLmmOpLNS264RsRQg/d6lQrWZmVk7ValHCx0aEUsk7QLMkPR0uROmkDwLYI899sirvo0un/G3Bod964h9cl++mZk1n4ps6UXEkvR7OXAbMBRYJmk3gPR7eQPTToiI6oiorqqqaqmSzcysHWjx0JO0vaSutd3Ap4H5wDRgTBptDHBHS9dmZmbtWyV2b+4K3Capdvm/i4g/SforMEXSGcDfgZMqUJuZmbVjLR56EfEcMKie9pXA8Jaux8zMiqM1XbJgZmaWK4eemZkVhkPPzMwKw6FnZmaF4dAzM7PCqNQdWdoF363FzKxt8ZaemZkVhkPPzMwKw7s3rRju/XHDww4f23J1mFlFeUvPzMwKw6FnZmaF4dLjy0gAAAgUSURBVNAzM7PCcOiZmVlhOPTMzKwwHHpmZlYYDj0zMysMh56ZmRWGQ8/MzArDd2SxLdZubrjtu7WYFYa39MzMrDAcemZmVhgOPTMzKwyHnpmZFYZDz8zMCsNnb5o1xmd2mrUrDr2ctJvT+c3M2hGHXgU4EM3MKsOhZ5YH7xY1a5Uceq1MY1uB4C3BVqWxYLP2wf+8tDsOvXbEu03NzBrX6kJP0gjgCqADcHVEXFrhkgqtqS1PM7O2pFWFnqQOwHjgCKAG+KukaRGxoLKVtR4OoXZgS3eL5rU7zbvwtozXW5vUqkIPGAosiojnACRNBkYCDr0c5RGkje5q7XhLsy+vEJoKyzy+aLf0i72lp2tKazr+2tr+6SmY1hZ6vYAXS/prgI9WqJZ2pTVtIT783Motmm7YXj22eL5NTdvc82zp6YBcAqPReg7folm2rgDKSxFe45aq8BayIiL3hZRL0knAkRFxZuo/FRgaEeeWjHMWcFbq3Rd4JodSegKv5DDfovN6bX5ep83P67T5tfQ63TMiquob0Nq29GqAPiX9vYElpSNExARgQp5FSJodEdV5LqOIvF6bn9dp8/M6bX6taZ22thtO/xXoL6mfpH8ARgPTKlyTmZm1E61qSy8i1kn6OnA32SUL10bEUxUuy8zM2olWFXoAEXEncGeFy8h192mBeb02P6/T5ud12vxazTptVSeymJmZ5am1HdMzMzPLjUOvhKQ+ku6VtFDSU5K+Uema2gtJHSQ9LukPla6lPZC0k6Spkp5On9dhla6prZP0rfR3P1/STZI6V7qmtkjStZKWS5pf0tZd0gxJz6bfO1eqPofeptYB346I/YBDgHMkDahwTe3FN4CFlS6iHbkC+FNEfBgYhNftVpHUC/hXoDoi9ic7kW50ZatqsyYCI+q0XQDMjIj+wMzUXxEOvRIRsTQiHkvdq8m+SHpVtqq2T1Jv4Gjg6krX0h5I2hH4Z+AagIh4NyJer2xV7UJHoIukjsB21LlG2MoTEQ8Ar9ZpHglMSt2TgFEtWlQJh14DJPUFhgCPVraSduEXwPnAhkoX0k7sBawArku7jK+WtH2li2rLIuIl4DLg78BSYFVETK9sVe3KrhGxFLKNC2CXShXi0KuHpB2AW4BvRsQbla6nLZN0DLA8IuZUupZ2pCNwEHBVRAwB3qSCu4vag3SMaSTQD9gd2F7SKZWtyvLg0KtDUieywLsxIm6tdD3twKHAcZIWA5OBT0q6obIltXk1QE1E1O6FmEoWgrblPgU8HxErIuI94FbgnypcU3uyTNJuAOn38koV4tArIUlkx0kWRsTPK11PexARYyOid0T0JTsx4J6I8H/QWyEiXgZelLRvahqOH7+1tf4OHCJpu/Q9MByfHNScpgFjUvcY4I5KFdLq7shSYYcCpwJPSpqb2r6X7hJj1pqcC9yY7lH7HPDFCtfTpkXEo5KmAo+RncX9OK3oLiJtiaSbgMOAnpJqgIuAS4Epks4g+wfjpIrV5zuymJlZUXj3ppmZFYZDz8zMCsOhZ2ZmheHQMzOzwnDomZlZYTj0zMysMBx6ZptJ0jhJ5zXTvHaS9LWS/t3T9WLNIj1+aK/mmt8WLL9K0p8qtXyzuhx6ZjlLd+1vyE7AxtCLiCURcWIzLXcg0CEinmuO+W2JiFgBLJV0aKVqMCvl0DMrg6QLJT0j6c/AviXt90mqTt090z1GkXS6pP+R9HtguqQdJM2U9JikJyWNTLO4FPiQpLmSfiapb+3DNyV1lnRdGv9xSYeXzPtWSX9KD+X8aQNlf4GS2z1JWiPpJ5LmSPqzpKGp/uckHZfGGShpVqpnnqT+qf2UkvZfS+qQ2kek1/SEpJkN1HF7qsWs4nwbMrMmSPoI2X1Dh5D9zTwGlPPUiGHAgRHxatraOz4i3pDUE3hE0jSypyPsHxGD07L6lkx/DkBEHCDpw2ThuU8aNjjV8w7wjKRfRsSLdZZ/KHBTSf/2wH0R8V1JtwEXA0cAA8iecTYN+ApwRUTU3uKsg6T9gJOBQyPiPUn/DXxB0l3Ab4B/jojnJXVvYD3MTssyqziHnlnTPg7cFhFvAaSwKseMiKh9mKaA/5D0z2TPFewF7NrE9B8DfgkQEU9LegGoDb2ZEbEq1bMA2BOoG3q7kT13r9a7QO3xtSeBd1KIPQn0Te0PAxemB//eGhHPShoOfAT4a3YvZrqQ3SX/EOCBiHg+1Vj3waG1lpM9rses4hx6ZuVp6Ca163j/MEHnOsPeLOn+AlAFfCQFzeJ6xq9LjQx7p6R7PfX/Lb9dZxnvxfs3291QO4+I2FB73DEififpUbIn3d8t6cxUx6SIGLtJcdku0XJu3ts51WJWcT6mZ9a0B4DjJXWR1BU4tmTYYrKtIIDGTkDpRvYw3ffSsbk9U/tqoGsjy/0CQNqtuQfwzGbUvRDYezPGJ53p+VxEXEm2u/NAYCZwoqRd0jjdJe1JtlX4CUn9atsbmO0+wPzNqcMsLw49syZExGPAzcBcsgcM/2/J4MuAr0p6COjZyGxuBKolzSYLsqfTvFcCD0qaL+lndab5b7Jjak+m5Z8eEe9Qvj+SPeJlc5wMzE+P1vowcH1ELAC+T3ZMcR4wA9gtnZl5FnCrpCdSjUiqlnR1yTwPT7WYVZwfLWTWTknqAtxLdgLK+grW8QAwMiJeq1QNZrUcembtmKQjgYUR8fcKLb+KLHRvr8Tyzepy6JmZWWH4mJ6ZmRWGQ8/MzArDoWdmZoXh0DMzs8Jw6JmZWWH8f3G6f3cFvzBhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAEKCAYAAAA/9Q3BAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZwddZnv8c/T6SaQhOydBRKIiIIQyEKzCIEACZg0gjh3ZsIsiogEiKOOV2euzp1FvVevs7xm9N6ZQXIRCCQsShKU1aAg4BWQBBBCosZKWGL2DULkpJc894+qkxRNL6e7T506Vef7fr361X22qqeTU/2c5/f7VT3m7oiIiEj1q0s7ABERESmNkraIiEhGKGmLiIhkhJK2iIhIRihpi4iIZISStoiISEYoaSfMzL5jZn+XdhylMLNXzGx22nGISGnM7CEzuzLtOKRylLT7KUp0b5vZW2a228weMLOJxcfd/Tp3/x9l2peb2VYzq4/dV29m28ysZk64N7MPm9kvzGyfme00s8VmdnQPrzEz+yszWxf9f71mZt8ws8MqFbdkS4dje4uZ3WpmQ9KOK87d57r7orTjKCcd391T0i6PS919CDAe2Ar8nwT3tQeYG7vdDOxOcH9Vxcz+ELgD+DYwGjgZaAGeNLPh3bz0fwPzgY8DRxL+G84G7ko0YMm64rE9FZgGfLmSO49/QK8FOr5L4O766scX8AowO3a7GfhN7PatwP+Mfj4f2Ah8AdgGbAau6sW+HPhb4Pux++4B/nv4X3nwvquAtcBeYD1wbeyx0cD9hMl/F/AkUNfxdwFOBDYAV3QRy8XAr4E3gP8EHgc+FXv8k1EMu4EfAcd2+D2uA9ZFj/8HYCX8/ga8Cvx1h/vrgNXAP3TxuvcB7cAZHe6fCOwHZqb9PtJX9X11cmz/E/BA7PZA4F+A1wg/rH8HOCL2+EeAF4A3gQCY08V2vwIsjn6eFB0fV0fbfQI4HFgM7IyO22eBsdHzfwp8KoplDzA5tt1G4G1gTHT7w1E8e4CfA6d287vr+K7SL1XaZWRmg4B5wNPdPG0cMAw4mvDA/A8zG9GL3dwLnGdmw6NPnucCP+jwnG2EB+hQwgT+b2Y2PXrsC4QfHBqBscDfEB5k8d9jOrAC+Iy7v+uTqpmNJvyw8GVgFOHBfXbs8cuj7f5BtJ8ngTs7bObDwOnAFOCPgQ+V8LufABwDfD9+p7sfAJYS/qHpzCxgo7v/osPrXif8v+rqdSIAmNkEwurtt7G7/xF4P2EVfjzhMf330fPPAG4D/goYDpxHmKxLNRP4AOFxcSXh34yJhMfbdYTJ+CB33w8sA/4kdvcfA4+7+7bomL4ZuDbaxo3AD81sYCe/q47vKqakXR73mtkewk/UFwH/3M1zW4GvuXuruz8IvEX4Zi1VAbiP8MPBFcAPo/sOcvcH3D3w0OOECfjc2P7HE34ybnX3Jz36WBo5N9rmle5+fxcxNAMvu/syd28jHJraEnv8WuB/ufva6PFvAFPN7NjYc77p7nvc/TXgMcI/fD0ZHX3f3Mljmwn/gHT1us5e09PrRO41s73A64Qfhv8BwjlU4Brg8+6+y933Er7Pr4hedzVws7s/4u4H3P137v6rXuz3K+6+z93fJjxmRwHHu3u7u69y9zc7ec0dvDNp/2l0H1GsN7r7M9E2FhFWoWd1sh0d31VMSbs8Lnf34YRDVH8BPG5m47p47s7ojV70e6C3i1tuI5y7+Xj08zuY2Vwze9rMdkUfJpo5dED8M2G1sMLM1pvZlzq8/Drg5+7+WDf7P4rwjxgQjcuH1XvRscC3zWxPtP9dhENf8cUk8T8Cpf4b7Ii+j+/ksfHAdgAzezlaPPSWmZ0bva6z17zjdSKduNzdjySc2jqRQ8dRIzAIWBV7nz/MoQQxkXBIvK9ej/18O+EQ9F1mtsnM/snMGjp5zaPAEWZ2ZpRApwLLo8eOBb5QjDWKdyLhsdyRju8qpqRdRtEn2GWE8yszEtzVk4RvxrHAz+IPRMNdSwnn2sZGHyYeJDyocPe97v4Fdz8OuBT4r2Y2K7aJ64BjzOzfutn/ZmBCbJ8Wv014wF/r7sNjX0e4+8/7+PsW/Zrwj8cfxe80szrgvxDOu+HuJ7v7kOjrScI/ZhOjIcv46yYSVhqP9zMuybloxOpWwuMKwkTxNnBy7D0+zMNFaxAeA+/tYnP7CBN+UWcf8A+OfkUjYl9195MIh6k/TPiBvWOMB4DvEVbbfwrcH40AFOP5eodjcpC7dxzWBh3fVU1Ju4yi0w4+AowgXKSRiOiT76XAZR2GtgEOI6z4twNtZjaX2JxOdDrF8dGB+CbhB4z22Ov3AnMI582/2UUIDwCnmNnl0erWT/POPzzfAb5sZidH+xxmZn/UyXZ6Jfpdvwj8rZn9qZkdEY1o3ERYAXW6at/dfxPFtMTMzjKzAVFsSwkX5Py4v7FJTfgWcJGZTY0S5P8lXC8yBsDMjjaz4tztd4GrzGyWmdVFj50YPfYCcIWZNZhZE/CH3e3UzC4ws1PMbADhMdvKO4/ZuDsIp87+jEND40SxXhdV4WZmg83sEjM7spNt6PiuYkra5XGfmb1FeEB9nXA++OXebsTM/sbMHirlue7+cmf7iD5Zf5bwE/duwk/cP4w95X2Eb+K3gKeA/3T3n3bYxh7Cufm5Zvauc8zdfQfhp+F/IlzRehKwknCODHdfTrhI5y4ze5Nw5efcjtvpjJmdG/1bdvV73w18DPh8tO/NhAteZrp7V/NaEE5b3ES4Cvf3UUyvEg5/HiglNqlt7r6dcDqqeLGk/0Y41fR09D7/MdH6lGhR1FXAvxGuwH6ccFiZ6PXvJTw+v8o7k2tnxhEuDHuTsBh4nPB93FmMzxBW8kcBD8XuX0k4r/3v0X5/C3yii23o+K5i9u5CTaR3ouGrjcCf9TAXnsS+LyZcuTrL3V/oxeu+BlwOnBd9SBGRTuj4ri6qtKVPzOxDFp52NpDw9A+j+1PdEuHuKwgrhs5WwXb3ur8HFvb2dSK1QMd39VKlLX1iZl8BPkM4h74G+Gw0NCciGafju3opaYuIiGSEhsdFREQyQklbREQkI6qqg8zo0aN90qRJaYchkhmrVq3a4e5VfZlGHdcivdPdcV1VSXvSpEmsXLky7TBEMsPMXk07hp7ouBbpne6Oaw2Pi4iIZISStoiISEYoaYuIiGSEkraIiEhGKGmLiIhkhJK2iIhIRihpi1SZAwdy1UlQRCLlOLaVtEWqSKFQ4JZbbuGVV15JOxQRKaMgCLjpppv6nbiVtEWqRKFQYPHixYwbN45jjz027XBEpEyCIGDZsmXMmTOHurr+pV0lbZEq8eijjzJ+/Hiam5sxs7TDEZEy2L9/P/fddx/z5s3jmGOO6ff2quoyppKceTc+BcDd134w5Uiko0KhQHt7OxdddBH19fVK2FJ5t1wSfr/qgXTjyJldu3YxYsQIFixYwGGHHVaWbarSFklRoVBgyZIlPPfcczQ0NFRVwjazm81sm5mtjt030sweMbN10fcRacYoUq2CIOC73/0ub7zxRtkSNqjSzr1ihf3Mhl3vuK2KO33FhD1u3DhmzJiRdjiduRX4d+C22H1fAn7i7t80sy9Ft/9bCrFJORQr7Fd/9s7bqrj7pTiHPW/ePIYPH17WbavSFknJxo0bOeqoo6p2DtvdnwB2dbj7I8Ci6OdFwOUVDUokA15++eWyzWF3pEo754oVtSrs6lEoFPjtb3/L5MmTOf7449MOp7fGuvtmAHffbGZjOnuSmc0H5gOJ/OGSMilW1Kqwy2L9+vUMHz6cyy67LLF9qNIWqaDikPhrr72Gu6cdTmLcfaG7N7l7U2NjY9rhiCQuCAKWLl3Kvn37Et2PKu0aoQo7ffE57Llz51blkHgJtprZ+KjKHg9sSzsgKQNV2P0Sn8OeOHFiovtSpV0j5t341MEhckmHmXHyySdX7Rx2iX4IXBn9fCXwgxRjkXK55ZJDQ+TSawMGDEhsDrsjJW2RhBUKBe69917cnbPOOiszCdvM7gSeAk4ws41mdjXwTeAiM1sHXBTdFqlJQRDw5JNPMmnSpIqt3dDweM7plK90xYfEBw4cmHY4veLuf9LFQ7MqGogkR6d89Vl8SLySVGmLJMTdufPOOxk3blzWh8RFJGbTpk0HE3alz46walrB2tTU5CtXrkw7jFxShV1Z7e3tDBgwgM2bNzNu3LjEEraZrXL3pkQ2XiY6rjNAFXbJ2tvbAdixYwdjx45NZB/dHdeqtEXKrFAocOutt7JhwwbGjx+vClskJ4rtNc0ssYTdEyVtkTKKz2FPmjQp7XBEpEyKc9hz587td3vN/lDSFimjxx57THPYIjnT0tLC/fffn8ocdkdaPZ5zWj1eGYVCgba2NmbPnq32mpIdWj3eo507dzJy5EgWLFhAQ0ND2uGo0hbpr2purykifRcEATfffDN79uypioQNqrRzTw1DkhWfwz733HPTDkekd9QwpEvx87BHjKietvFK2iL9sGnTJo4++mg+9KEPqcIWyZE1a9ZUxRx2R0raIn1QKBRYt24dp5xyCscdd1za4YhImaxfv55hw4Zx6aWXph1KpxKd0zaz4WZ2j5n9yszWmpnGZiXzikPir7/+eq7ba4rUmkq11+yPpCvtbwMPu/sfmtlhwKCE9ycdaPV4eeWkvaZISKvHD4rPYVfbkHhcYknbzIYC5wGfAHD3FqAlqf2JVEJdXR2TJ0/mjDPOUMIWyZH6+vqqT9iQbKV9HLAduMXMpgCrgM+5+zvGHcxsPjAfqPp/rCzS6vHyKBQKPPjgg8ydO5czzzwz7XBEykOrxwmCgI0bNzJz5sy0QylJknPa9cB04AZ3nwbsA77U8UnuvtDdm9y9qbGxMcFwRPqmOCQ+cOBADj/88LTDEZEyKQ6Jv+c970k7lJIlmbQ3Ahvd/Zno9j2ESVwkM9ydu+66S5cmFcmZzZs3Z2IOu6PEhsfdfYuZvW5mJ7j7r4FZwJqk9idSbsX2mnPmzGHs2LFK2CI50d7ezpgxY7jyyisZM2ZM2uH0StKrxz8DLIlWjq8Hrkp4f9KBVo/3TXFI/IILLtB52JJfNbh6PAgCHnnkEebPn5+5hA0JJ213fwHotJG3SLWKn9aVpbkuEele/LSuNNtr9oeuiJZzWj3ee48//rjmsKU21NDq8ZaWFh544IHMzWF3pKQtEikUCrS2tjJr1iwGDBighC2SEzt27GDUqFFcf/31VdOtq6+yOT4gXZp341MHq+q4u6/9oKrsbhSHxJ9//nn1wxbJkSAIuOWWW9i9e3fmEzYoaYuovWYvmdnnzGy1mb1sZn+ZdjwiXYnPYY8cOTLtcMpCw+M5oVXifbd582a11yyRmU0GrgHOILws8cNm9oC7r0s3MumTnK8e/9WvfpX5OeyOlLSlZhUKBX79618zZcoUrRIv3QeAp9399wBm9jjwUeCfUo1KJGb9+vUMHTqUSy65JO1Qyk5JOye0Srx3ikPi48eP59RTT1WFXbrVwNfNbBTwNtAMrOz4JPUUyIgcrh6PD4nnkea0pebs379f7TX7yN3XAv8IPAI8DPwSaOvkeeopIBWXlfaa/aFKO2dUYfesrq6OU045hdNPP10Juw/c/bvAdwHM7BuEfQYky3JQYQMcdthhuU7YoEo7d7o65UvCIfFly5bR2tqqftj9YGZjou/HAH8A3JluRFKyWy45NBSeI0EQ8NhjjzFx4sRcJ2xQ0pYaEW+vecQRR6QdTtYtNbM1wH3Ap919d9oBSe0qDom/973vTTuUitDweE7olK+uqb1mebm7TmbPmpye2pXV9pr9oaQtuVZsrzl37lzGjBmjhC2SE+3t7YwdOzaT7TX7Q0k7J3TK17sVh8RnzpzJ8ccfn3Y4IunJ2aldQRCwYsUKrr322ppK2KA5bcmp+KVJa2WuS6QWFOewL7nkksy21+wPVdo5EK+uVWGHnnjiCc1hS23qrprOeIWdl/aa/aGkLblSKBRoaWnhwgsvVHtNkRzZvn07o0ePZsGCBdTX127qqt3fPAe0YvydikPixx9/PDNnzkw7HJHKyukKcTg0JH711VfnpltXX9XehIDkUqFQYPHixYwbN47zzjsv7XBEpEzy2F6zP1RpZ5hWjB+yZcsWJkyYoPaaUrtytkK86De/+U1Nz2F3pKQtmVYoFFi7di3Tpk1j0qRJaYcjImWyfv16hgwZwty5c9MOpapoeFwyqziHvWnTJtw97XBE0tHxeuJXPZD5KjsIApYuXUqhUEg7lKqjpC2ZFG+vqdO6RPKjFtpr9oeGxzOslleP19XVceqpp9LU1KSELbUpp6vFBw4cqITdDVXakimFQoGlS5fS2tqqftgiORIEAY8++igTJkxQwu6GKu0Mq7XV4/FLk6q9ptS8HK0Wjw+JS/dUaUsmuDt333235rBFcmbLli2aw+4FVdoZVisVdltbG/X19cydO5fGxkYlbKlNXVXUGa6w29raGDNmDJ/4xCdobGxMO5xMUKUtVa1QKLBo0SLWrVunftgiORIEAQsXLsTdlbB7QZV2BtXKqvH4HLb6YUvNyuEq8fgc9oABA9IOJ1MSTdpm9gqwF2gH2ty9Kcn9Sb48+eSTmsMWyZnW1lYefPBBzWH3kSV5JakoaTe5+45Snt/U1OQrV65MLJ68yXOFvX//fgYPHqz2mj0ws1XV/mFYx3WZ5KDC3rZtG42NjbS3t9d0e82edHdca05bqkpxSPz555+nvr5eCVskJ4IgYNGiRezatUsJux+S/pdzYIWZOXCjuy/s+AQzmw/MBzRUUuPic9jqhy01qbtqOsMVdnwOe9SoUWmHk2lJV9rnuPt0YC7waTN7V6Njd1/o7k3u3qQVhLVt27ZtTJgwQXPYGWBmnzezl81stZndaWaHpx2TVK9169ZpDrtMEq203X1T9H2bmS0HzgCeSHKftSBvq8cLhQJr1qxh+vTpOqgzwMyOBj4LnOTub5vZ94ArgFtTDSzLcrhCHML2moMHD2bOnDlph5IbiVXaZjbYzI4s/gxcDKxOan+STcUh8c2bN6u9ZrbUA0eYWT0wCNiUcjxSZYrtNffv3592KLmSZKU9FlgeDXPWA3e4+8MJ7q9m5OWa42qvmU3u/jsz+xfgNeBtYIW7r4g/R2tVeilH1xGHsMLWpUmTkVjSdvf1wJSkti/ZN2DAAKZMmcJpp52mhJ0hZjYC+AjwHmAP8H0z+3N3X1x8TrTodCGEp3ylEqik5vDDD1fCTojW3WdYVivsQqHAfffdR3NzM01NVX2KsXRuNrDB3bcDmNky4Gxgcbevkp5lvMIOgoANGzYwe/bstEPJLZ2nnWHzbnzq4BB5VhTnsAcNGsSgQYPSDkf65jXgLDMbZOEQySxgbcoxZd8tlxwaHs+g4mld73//+9MOJdeUtKVi1F4zH9z9GeAe4DngJcK/I++6BoPUjq1bt2oOu0ISvYxpb+lyh6XpeMrXme8ZCVT3cHlbWxsDBgxgx44djB49Wgm7THQZ0xzoeLrXsTPC7xkZKm9ra6Ouro5du3YxevTotMPJBV3GVFIVb6+pftgi+REEATfeeCPuroRdIVqIlkFZOuUrfmnS973vfWmHI1J9Mnq6l9prpkOVtiTqZz/7meawRXKmtbWVhx56SHPYKVClnUFZqbALhQIXXHABdXV1Stgi3clQlb1161bGjBnD9ddfrwo7Baq0pezi7TXVD1skP4Ig4LbbbmPnzp1K2ClRpZ0hWWgUEp/DPv/889MOR6S6ZahRSHwOW4vO0qOkLWW1fft2JkyYwMUXX6wKWyRHgiDQHHYVUNLOkGpeNV4oFFi9ejWnnXYaEydOTDsckWzIwMrxIAgYNGgQF198cdqhCJrTljIoDolv3bo17VBEpIyKQ+Ktra1phyIRJW3pF7XXFMmn+By2hsSrR0nD42Z2OHA1cDJwePF+d/9kQnFJRtTX1zN16lSmT5+uhC2SI4MGDVLCrkKlzmnfDvwK+BDwNeDPUFefiqum1ePF9ppz587ltNNOq/j+RXKhClePB0FAEASaw65SpQ6PH+/ufwfsc/dFwCXAKcmFJdUs3l5z8ODBaYcjImVSHBI/8cQT0w5FulBqpV1chbDHzCYDW4BJiUQkXaqG1ePuzve+9z3NYYuUQxWtHt+2bZvmsDOg1KS90MxGAH8L/BAYAvxdYlFJVWptbaW+vp7m5mZGjRqlhC2SE62trTQ2NnLVVVfpwilVrtSk/RN33w08ARwHYGbvSSwq6VTac9hLlixhxowZnHDCCRXfv0hupVxlB0HAQw89xPXXX6+EnQGlzmkv7eS+e8oZiFSv+KVJ3//+96cdjoiUSXEO+7LLLtO1xDOi20rbzE4kPM1rmJn9QeyhocRO/ZJkpb1q/Oc//7nmsEXKLeWV421tbTz88MOaw86YnobHTwA+DAwHLo3dvxe4JqmgpDoUCgXefvttZs6cqfaaIjmydetWGhsbue6661RhZ0y3SdvdfwD8wMw+6O5PVSgm6SCNVePFIfHjjjuOCy64IPH9idSclFaOF4fEtegsm0pdiPa8mX0aXRGtJqi9pkg+qb1m9umKaBlSqTnsnTt3MnHiRC666CINicu7mNkJwN2xu44D/t7dv5VSSNlWwVXjGzZs0Bx2xpWatI939z8ys4+4+yIzuwP4UZKByTtVYmi8UCjw0ksv0dTUxNFHH53YfiTb3P3XwFQAMxsA/A5YnmpQWVaB4fEgCDjiiCOYPXt2YvuQyij1lK+OV0Qbhq6IlivFIfFt27alHYpkyywgcPdX0w5EOlccEm9ra0s7FCkDXRGtylXidC+115R+uAK4M+0gMqkCp3ytX79elybNmZ7O0/6vsZtXRd//I/quThE5UV9fz7Rp05g2bZoStpTMzA4DLgO+3Mlj84H5gJJFigYPHqyEnTM9VdpHRt9PAE4nrLIhPGf7iVJ2EM15rQR+5+4f7kuQtSzJ070KhQI/+MEPaG5uZvr06WXbrtSMucBz7r614wPuvhBYCNDU1OSVDiwTEjzlKwgC1q1bx5w5c8q2TakO3c5pu/tX3f2rwGhgurt/wd2/AJwGTChxH59DK82rTnEOe8iQIQwZMiTtcCSb/gQNjVed4hz2SSedlHYokoBSF6IdA7TEbrdQwkI0M5tA2Hv7pl5HJolxd77//e9rDlv6zMwGARcBy9KOJdPKXGWrvWb+9eY87V+Y2XLAgY8Ci0p43beAv+bQMPu7aO6rsuLtNUeOHKmELX3i7r8HRqUdhxxSbK/5yU9+klGj9F+TVyUlbXf/upk9BJwb3XWVuz/f3WvM7MPANndfZWbnd7NtzX11o5yrx4tD4ueccw4nnnhi+YIUkd4p88rxIAh48MEHuf7665Wwc67USht3fw54rhfbPge4zMyaCS99OtTMFrv7n/cyRimD+KVJ1Q9bJD/ilyatry/5T7pkVGL/w+7+ZaJTQaJK+4tK2L1XrtXjTz31lOawRapFmVaOt7W18aMf/Uhz2DVEH8tyrlAo8Pvf/56ZM2diZkrYIjmxZcsWxowZw7XXXqv2mjWk1NXj/eLuP9U52l2bd+NTByvprtx97Qd7XWUXh8RfeOEF9cMWqTb9qLKDIOD2229n586dStg1piJJWyovPoetftgi+RGfw25sbEw7HKkwDY+nKMnriu/atUvtNUWqUT9Xjr/66quaw65hSto5UygUePHFFzn99NM56qij0g5HRMokCAIGDhzIhRdemHYokiIl7RSV+7ri8SFxEalSfVg5Hh8Sl9qmOe2cUHtNkXxSe02JU6VdId1V0+WYw25oaGD69OlMnTpVCVukmvVyDnvIkCFK2HKQKu2MKxQK3H333ezbt0/9sEVyJAgCHnroIcaMGaOELQep0k5YkivE43PYaq8pUuV6sWpcc9jSFVXaGaX2miL5tH37ds1hS5fMvXoaazU1NfnKlSvTDiMR5aywi+01d+/ezYgRI5Swa5iZrXL3prTj6E6ej+s+6abCbmlpoaGhgV27dqlbVw3r7rhWpZ0xhUKB2267jbVr16oftkiOBEHAd77zHdrb25WwpUtK2hkSn8P+wAc+kHY4IlImxTnsyy+/XO01pVtK2hny9NNPaw5bJGfa29tZsWKF5rClJPpIl7ByrB4vFArs27eP8847T+01RbKqk9Xjm1sGMfaau7n22mupq1MNJT3Tu6TKFYfEf/nLX6q9pkiOBG8PZfG2E9ixY4cStpRMlXbC+nN9cbXXFMmR2DXHg7eHsmzfdOZdOY8xY8akG5dkipJ2Fdu9e7faa4rk0Gv7j9QctvSJknaF9LbCfuGFFzjzzDMZP358glGJ9I2ZDQduAiYDDnzS3Z9KN6rqFwQBh83+DhdMnJh2KJJRmkipgHk3PnVweLwnxSHxnTt3JhyVSL98G3jY3U8EpgBrU46n6gVBwLI7bsHv/2LaoUiGqdKuImqvKVlgZkOB84BPALh7C9CSZkzV7mB7zdHrOObwt9IORzJMSTtBvT3dq6GhgdNOO40pU6YoYUs1Ow7YDtxiZlOAVcDn3H1fumFVr6E//ivmHVnPMVseDu/oZXtOkSINj1eBQqHAXXfdxVtvvaV+2JIF9cB04AZ3nwbsA74Uf4KZzTezlWa2cvv27WnEWBWCIOCBBx5gdEOBYwaqwpb+U6WdoFJO94qf1nXkkUdWND6RPtoIbHT3Z6Lb99Ahabv7QmAhhA1DKhtedXhHe81jDp3uBajClj5TpZ2ye+65R3PYkinuvgV43cxOiO6aBaxJMaSqs2PHDrXXlESo0k5YV1V2sQXfJZdcwvDhw5WwJWs+Aywxs8OA9cBVKcdTNVpaWhg1ahRXX301I0eODO9UhS1loko7BYVCgdtvv501a9aoH7Zkkru/4O5N7n6qu1/u7rvTjqkaBEHADTfcQFtb26GELVJGqrQT0tXK8UVXTjs4h33SSSelFp+IlFd8DruhoSG8s5MmIYAqbukzVdoV9otf/EJz2CI5097eziOPPKI5bEmcKu2EdFw5vujKaezdu5dRo0apvaZIjmzevJkxY8Ywf/78d3frukqrxqW8VGlXQN2BVpYsWcKLL76o9poiORIEAYsXL2bnzp1qr3+VlTAAAA9MSURBVCkVkVilbWaHA08AA6P93OPu/5DU/qpNZ3PYF154YcpRiUi5xOewu22vqSpbyijJ4fH9wIXu/paZNQA/M7OH3P3pBPdZdfbs2cMxxxzD7NmzVWGL5MjGjRs1hy0Vl1jSdncHitfta4i+cn9lpGKF/dyGbby/fgefu9fBhnCRErZILgRBQH19PTNnzuz+iVo5LglIdBLGzAaY2QvANuCR2GUP48/J3TWK6w60cvHAdRxp+9MORUTKqDgkrlEzSUuiq8fdvR2YambDgeVmNtndV3d4Tq6uUXz7Vadx++238+KQkfiRJ3L3dWenHZKIlMHB9pqlDolr5bgkoCLLHd19D/BTYE4l9pem+vp6mpqa2HLkiaBP4yK5MWzYMM1hS+qSXD3eCLS6+x4zOwKYDfxjUvtLW6FQYPny5TQ3N/ONp38PZl32zRaR7AiCgDVr1nDppZcyatSo0l+oClsSkGSlPR54zMxeBJ4lnNO+P8H9pabYXnPo0KEMHTo07XBEpEyKc9hTpkxJOxQRINnV4y8C05LafjVZunQp48aN47bXR3Dbwqffdb1xVdwi2bNz586+tdfUqnFJkC7h0w8tLS24O83NzTQ3N2sOWyQnWlpaGDlyJJ/61Kc0hy1VRdce76PikPiZZ57J5MmTgXdfb1wVtkj2BEHA/fffz4IFCxgxYkTvN6BV45IgVdp9UEzY48aN4+STT047HBEpk+Ic9kc/+tFD7TVFqogq7T549tlnu2yvqSpbJJva29v58Y9/3PfTujpW1qqwJQFK2r1QKBR48803mTFjBoCuiiSSE5s2bWLs2LFcc8016tYlVU1Ju0TFIfFjjz2W2bNnv+vxYoWtleMi2VIcEv/4xz/O2LFje78BrRaXCtJHyhLE57BnzZqVdjgiUibx9pp9StgiFaZKuwRvvPEGxx57LLNmzepySFwrx0WyZ9OmTf2/NKlWi0sFKWl3o1AosGrVKs4+++xOh8RFJJuCIKCuro5zzz037VBEekVJuwvxIfHeUIUttcDMXgH2Au1Am7s3pRtR6eJD4mWlClsqQEm7Ey0tLQcTdmendXWkIXGpURe4+460g+iNXrfX7I6GwyUFStqdaGho4IwzzmDy5Mk6rUskR4YPH84VV1zBxIkT0w5FpE+UtGMKhQLLli2jubmZU045pcfn6zQvqWEOrDAzB25094VpB9SdIAh4+eWXueyyyxg5cmT/NqZTvCRFOuUrUpzDHjZsGMOGDUs7HJFqd467TwfmAp82s/PiD5rZfDNbaWYrt2/fnk6EkeIc9tSpU1ONQ6QcVGlHli1bVvIcdpFO85Ja5e6bou/bzGw5cAbwROzxhcBCgKamJk8lSPrRXrM7OsVLUlTzSbulpYWGhgaam5sZNmyY5rBFemBmg4E6d98b/Xwx8LWUw3qX/fv3M3LkSK655hqGDx+edjgiZVHTSbs4JH7GGWeUNIcdF6+uVWFLjRkLLI8+4NYDd7j7w+mG9E5BEHDfffexYMGC/ifsripqVdiSgppN2vHzsIv9sEWkZ+6+HpiSdhxdiZ+Hfdhhh6UdjkhZ1WzSXrVqVa/nsEErxkWq2YEDB/jJT35S3vOwtUpcqkjNJe1CocAbb7zB2WefDai9pkhe/O53v2Ps2LF86lOfUntNya2aSto9tdcshVaMi1SffrfX7IxWiUsVqpmPo2qvKZJPaq8ptaRmKu0333yzx/aanemqolaFLVIdNm/eXPocdl+qZlXYUkVyn7QLhQLPPvssM2bMUHtNkRwJggCAGTNmpByJSOXkOmn3tb0maJW4SDXrdXtNrQSXnMjtnHZv22uKSDZs2LCh/JcmFcmI3FbaDQ0NnHnmmZx88sl9SthaJS5SnfrUXlMrwSUncldpF4fE9+zZo37YIjkSBAH33nsvI0aMUD9sqVm5qrTjc9i9ud5wd9W0KmyR9HU5h93bylkVtmRcrirt5cuXaw5bJGd27dqlOWyRSGKVtplNBG4DxgEHgIXu/u0k9tXS0kJ9fT3Nzc0MHTq05IStFeIi1a1QKDBy5Ejmz5/PsGHDDj2g1eBSo5KstNuAL7j7B4CzgE+b2Unl3kmhUOD2229n9erV6octkiNBEHDDDTfQ0tLyzoQtUsMSq7TdfTOwOfp5r5mtBY4G1pRrH/E57N72wwatEBepVj2219RqcKlRFZnTNrNJwDTgmU4em29mK81s5fbt23u13eeee05z2CI5c+DAAR599FHNYYt0wtw92R2YDQEeB77u7su6e25TU5OvXLmyx20WCgX27NlzsDmAErbUKjNb5e5NacfRnVKPazjUXrOurk7tNaVmdXdcJ3pUmFkDsBRY0lPCLlVxSPyll17CzJSwRXIiCALuuOMOduzYoYQt0oXEjgwLs+l3gbXu/q/l2GZ8DlvNP0TyIz6H3ZdeASK1IsmPs+cAHwMuNLMXoq/m/mxw7969TJo0SXPYIjmzdetWzWGLlCDJ1eM/A8qaWRsbG5k1a1Y5NykiVeDss89OOwSRTNDEkYiISEYoaYtIr5nZADN73szuTzsWkVqipC0iffE5YG3aQYjUGiVtEekVM5sAXALclHYsIrVGSVtEeutbwF8TNgISkQpS0haRkpnZh4Ft7r6qh+f1+fLEItI1JW0R6Y1zgMvM7BXgLsLrMCzu+CR3X+juTe7e1NjYWOkYRXJLSVtESubuX3b3Ce4+CbgCeNTd/zzlsERqRmIXV+mLVatW7TCzV3v5stHAjiTiKQPF1jeKrXTHph1ATzoc19X271cp+r1rS39/7y6P68S7fCXNzFZWa5cjxdY3ii2/avXfT793bUny99bwuIiISEYoaYuIiGREHpL2wrQD6IZi6xvFll+1+u+n37u2JPZ7Z35OW0REpFbkodIWERGpCZlN2mY20cweM7O1ZvaymX0u7ZiKzOxwM/uFmf0yiu2raccUV80dmszsFTN7ycxeMLOVacdTZGbDzeweM/tV9J77YNoxZVE1v/eSUovvHTP7fPS3b7WZ3Wlmh6cdUxLM7GYz22Zmq2P3jTSzR8xsXfR9RDn3mdmkDbQBX3D3DwBnAZ82s5NSjqloP3Chu08BpgJzzOyslGOKq/YOTRe4+9QqO1Xk28DD7n4iMIXq/verZtX+3ktCTb13zOxo4LNAk7tPBgYQXognj24F5nS470vAT9z9fcBPottlk9mk7e6b3f256Oe9hAfC0elGFfLQW9HNhuirKhYPqENT75nZUOA84LsA7t7i7nvSjSp7avG9V8PvnXrgCDOrBwYBm1KOJxHu/gSwq8PdHwEWRT8vAi4v5z4zm7TjzGwSMA14Jt1IDomGAV8AtgGPuHu1xFbtHZocWGFmq8xsftrBRI4DtgO3REO7N5nZ4LSDyqBqf+8loebeO+7+O+BfgNeAzcAb7r4i3agqaqy7b4awuATGlHPjmU/aZjYEWAr8pbu/mXY8Re7e7u5TgQnAGWY2Oe2YSu3QlLJz3H06MJdwyuO8tAMirBqmAze4+zRgH2Ue8sq7jLz3klBz751oDvcjwHuAo4DBZqbr05dJppO2mTUQJuwl7r4s7Xg6Ew2F/ZR3z3ukoaQOTWly903R923AcuCMdCMCYCOwMTZacg/hH2IpXdW/9xJSi++d2cAGd9/u7q3AMuDslGOqpK1mNh4g+r6tnBvPbNI2MyOcJ1rr7v+adjxxZtZoZsOjn48gfBP/Kt2oqr9Dk5kNNrMjiz8DFwOru39V8tx9C/C6mZ0Q3TULWJNiSJlT7e+9pNToe+c14CwzGxT9nZ5FzhffdfBD4Mro5yuBH5Rz41XV5auXzgE+BrwUzR0D/I27P5hiTEXjgUVmNoDwg9H33L1mTnHph7HA8vA4px64w90fTjekgz4DLDGzw4D1wFUpxyPZUVPvHXd/xszuAZ4jPMvneXJ6ZTQzuxM4HxhtZhuBfwC+CXzPzK4m/ADzR2Xdp66IJiIikg2ZHR4XERGpNUraIiIiGaGkLSIikhFK2iIiIhmhpC0iIpIRStoVFHX7WRC7fX6lux1F+zw7dvs6M/t4H7f1Vs/Pqjwz+6mZVVOzEck5HduVoWNbSbvShgMLenxWP0UX6e/K+cSuTuTu33H325KOqVQ9xC5SrXRs90DHdnkoaVfWN4H3Rr2i/zm6b0is1+6S6ApCmNlpZvZ41DjjR7HL4k01s6fN7EUzW17s1Rp9Av2GmT0OfC66KttSM3s2+jonaqxyHfD5KIZzzewrZvbFaBvHm9mPLewD/pyZvdfMhpjZT6LbL5nZR7r7Bc1skoU9g/+vhf10V0RXhetN7Lea2Q0W9ktfb2YzLexbu9bMbo3t6wYzW2lV2LNcao6ObR3bleHu+qrQFzAJWB27fT7wBmFTkTrgKWAGYSvPnwON0fPmATdHP78IzIx+/hrwrejnnwL/Gdv2HcCM6OdjCC/3CvAV4Iux5x28Tdgl7aPRz4cTttSrB4ZG940Gfsuhi/K81cXv2AZMjW5/D/jzXsZ+K+H1qY2w8cCbwCnRv9Gq2LZHRt8HRNs4Nba9prT/v/VVO186tnVsV+pLwxXp+4W7bwSw8HKsk4A9wGTgkejD+QBgs5kNA4a7++PRaxcB349t6+7Yz7OBk6LXAwy16LrenYkeO9rdlwO4eyG6vwH4hoXdtg4Q9iwfC2zp5nfa4O7FS8uuAib1MnaA+9zdzewlYKu7vxTF8zLhv9ELwB9b2L6znvDSsScR/vEQqQY6tt8dO+jY7hcl7fTtj/3cTvh/YsDL7v7B+BOjg6M7+2I/1wEfdPe3O2yjq9d29cCfAY3Aae7eamGXpsN7iKPj73RED8+Hd8Ye38aBDts7ANSb2XuALwKnu/vuaGitp7hEKknHdkjHdhlpTruy9gJdfiKO+TXQaGYfhPATsZmd7O5vALvN7NzoeR8DHu9iGyuAvyjeMLOp3cXgYS/yjWZ2efT8gWY2CBhG2Ae51cwuAI4tIf536WXspRhK+MfgDTMbS9h/WyQtOrZ1bFeEknYFuftO4P+Z2erYYpXOntcC/CHwj2b2S8LhouKq0CuBfzazF4GphPNHnfks0BQtDFlDuEgF4D7go8XFKh1e8zHgs9G2fw6MA5ZE21lJ+Mm8Py1GS429R+7+S8LuQS8DNwP/rx9xifSLjm0d25WiLl8iIiIZoUpbREQkI5S0RUREMkJJW0REJCOUtEVERDJCSVtERCQjlLRFREQyQklbREQkI5S0RUREMuL/A2DKY2sqLhg3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import nonparametric_stats as nps\n",
    "plt.rcParams['figure.figsize'] = [7, 5]\n",
    "durs = (10**3)*np.array([durations_bin, durations_rec])\n",
    "\n",
    "# Obtain color cycle that matplotlib uses\n",
    "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
    "colors = prop_cycle.by_key()['color']\n",
    "\n",
    "# plot histograms of duration distributions\n",
    "labels=['Bin. Mask', 'Recursive']\n",
    "ax = nps.histPlotter(50, *durs, labels=labels, colors=colors)\n",
    "ax.set_title('Runtime durations of 2 power set generators')\n",
    "ax.set_xlabel('duration (msec.)')\n",
    "ax.set_ylabel('counts')\n",
    "ax.legend()\n",
    "\n",
    "# Examine Q-Q plots to compare distributions to\n",
    "# their corresponding theoretical normal dists\n",
    "fig, axs = plt.subplots(nrows=1, ncols=2)\n",
    "fig.tight_layout()\n",
    "\n",
    "for i in range(0,2):\n",
    "    nps.qqPlotter_normal(durs[i], 30, axes=axs[i], color=colors[i])\n",
    "    axs[i].set_title(labels[i] + ' gen. Q-Q')\n",
    "    axs[i].set_xlabel('theoretical normal')\n",
    "    if i == 0:\n",
    "        axs[i].set_ylabel('data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two distributions certainly do not *appear* to be normally distributed. Below, we import and run the Mann-Whitney non-parametric comparison test to assess if the difference between the two reaches statistical significance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Median Bin. Mask duration: 1.97 msec., Median Recursive duration: 4.46 msec.\n",
      "Mann-Whitney U statistic: 2260.50, p-value: 2.02e-159\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import mannwhitneyu as mwu\n",
    "med_bin = (10**3)*np.median(durations_bin) # compute median duration of Bin. Mask ps durations\n",
    "med_rec = (10**3)*np.median(durations_rec) # compute median duration of Recursive ps durations\n",
    "\n",
    "print('Median Bin. Mask duration: %0.2f msec., Median Recursive duration: %0.2f msec.'\\\n",
    "      %(med_bin, med_rec))\n",
    "# Run non-parametric test to determine if differences between \n",
    "# distributions are statistically significant\n",
    "uStat, pVal = mwu(durations_bin, durations_rec)\n",
    "print('Mann-Whitney U statistic: %0.2f, p-value: %0.2e'%(uStat, pVal))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "variables": {
     "print('%0.2f'%med_bin)": "1.97",
     "print('%0.2f'%med_rec)": "4.46",
     "print('%0.2f'%uStat)": "2260.50"
    }
   },
   "source": [
    "Measuring the time elapsed for each of the two approaches, we find that, on average, our recursive algorithm does indeed take about twice the time to complete on our machine as the binary-mask method ({{print('%0.2f'%med_rec)}} ms versus {{print('%0.2f'%med_bin)}} ms, U = {{print('%0.2f'%uStat)}}, p < 0.01). Of course, ours is certainly not a very controlled comparison; so please take the result with \"a grain of salt\" if you will."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Summary-main\"><h2>Summary</h2></div>\n",
    "\n",
    "In <a href=\"https://nullexit.org/blog/index.php/2021/01/27/power-sets-in-python-an-object-oriented-approach-part-1/\">part one</a> we discussed power sets and their relation to the binary theorem. In addition, we mapped out a recursive procedure for generating the power set. A brief investigation revealed that Python's `set()` class could not accommodate our need to include subsets within a set as in the case of power sets.  To solve the problem, we introduced a custom container class with a set-valued attribute, whose instances are \"hashable\" and thus, eligible for membership in `set` objects.\n",
    "\n",
    "Starting with the container strategy from part one, in <a href=\"https://nullexit.org/blog/index.php/2021/01/27/power-sets-in-python-an-object-oriented-approach-part-2/\">part two</a> we expanded our custom class, `ModSet()`, to include methods that perform copying, nesting, uniqueness-filtering, set-union and set-difference operations as well as others. We verified that all `ModSet()` methods functioned as desired.  Many of the methods added would prove necessary for both methods of power set generation to follow.\n",
    "\n",
    "In <a href=\"https://nullexit.org/blog/index.php/2021/01/29/power-sets-in-python-an-object-oriented-approach-part-3/\">part three</a>, we realized our goal of generating power sets. First, we detailed program flow of both the Binary Mask and Recursive approaches of power set generation. The mask routine used list indexing and took more lines of code than recursive method, but called fewer class methods.  Our Recursive approach was shorter in terms of lines of code because it incorporated more sibling methods of the `ModSet()` class. Explaining its program flow was \"interesting\" to the least; with effort on the part of the reader, it *might* strengthen understanding of recursive algorithms. Without effort, it has potential to confuse.  Finally, we compared relative run times of the two approaches. It turned out that, as expected, the Binary Mask routine took less time to run than its recursive counterpart."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"Sources-(part-3)-main\"><h2>Sources (part 3)</h2></div>\n",
    "\n",
    "<a name=\"ref1\">1.</a> D. Van Dalen; H. C. Doets; H. De Swart (9 May 2014). <a href=\"https://books.google.com/books?id=PfbiBQAAQBAJ&pg=PA1#v=onepage&q&f=false\">Sets: Naïve, Axiomatic and Applied: A Basic Compendium with Exercises for Use in Set Theory for Non Logicians, Working and Teaching Mathematicians and Students</a>. Elsevier Science. ISBN 978-1-4831-5039-0.\n",
    "\n",
    "<a name=\"ref2\">2.</a> <a href=\"https://en.wikipedia.org/wiki/Recursion_(computer_science)#Recursion_versus_iteration\">Recursion (computer science), Recursion versus Iteration</a>,  Wikipedia."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"The-complete-ModSet()-class-main\"><h2>The complete ModSet() class</h2></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can download `ModSet()`'s complete definition from github: <a href=\"https://github.com/rweaton/NullExitBlogPosts/blob/main/PowersetGenerationPost/modset.py\">modset.py</a>"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\" /></a><br /><span xmlns:dct=\"http://purl.org/dc/terms/\" property=\"dct:title\">The ModSet() class</span> by <a xmlns:cc=\"http://creativecommons.org/ns#\" href=\"nullexit.org\" property=\"cc:attributionName\" rel=\"cc:attributionURL\">nullexit.org</a> is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">Creative Commons Attribution 4.0 International License</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "class ModSet():\n",
    "    \n",
    "    # Instance constructor\n",
    "    def __init__(self, setElement):\n",
    "        self.val = setElement  # set arg saved in .val attribute\n",
    "    \n",
    "    # String eval representation of self instance \n",
    "    def __repr__(self):\n",
    "        return self.val.__repr__() # Return string eval represent.\n",
    "                                   # of string object in .val\n",
    "    \n",
    "    # Method to make a copy of self instance\n",
    "    def __copy__(self):\n",
    "        return ModSet(self.val)\n",
    "    \n",
    "    # Modify .val to contain the set of itself, of itself, ...\n",
    "    # nesting .val \"depth\" number of levels.\n",
    "    def pushDown(self, depth):\n",
    "        while depth > 0:\n",
    "            self.val = set([ModSet(self.val)])\n",
    "            depth -= 1     \n",
    "    \n",
    "    # Remove one nesting level from set in self.val.\n",
    "    # If un-nested, ignore.\n",
    "    def pullUpOneLevel(self):\n",
    "        listSet = list(self.val)\n",
    "        if len(listSet) == 1:\n",
    "            self.val = listSet[0].val\n",
    "        else:\n",
    "            pass\n",
    "    \n",
    "    # Remove \"height\" nesting levels from set in\n",
    "    # self.val by repeatedly calling above method\n",
    "    def pullUp(self, height):\n",
    "        while height > 0:\n",
    "            self.pullUpOneLevel()\n",
    "            height -= 1\n",
    "    \n",
    "    # Within a single set, multiple ModSets with \n",
    "    # equivalent .val attributes can exist.  \n",
    "    # This can occur because each instance of ModSet\n",
    "    # allocated its own location in memory. For sets\n",
    "    # to have unique members in terms of ModSet.val\n",
    "    # values, we define the following filtering method.\n",
    "    def removeDuplicates(self):\n",
    "        uniqueSet = set()  # initialize as empty the unique set to be built\n",
    "        for s in self.val:  # s is a member of the ModSet().val set\n",
    "            inUniqueSet = False  # initialize match detection flag as false\n",
    "            sTest = s  # default conditional testing value for s\n",
    "            for us in uniqueSet: # us is a member of the uniqueSet set\n",
    "                usTest = us  # default conditional testing value for us\n",
    "                if isinstance(us, ModSet): # if member us is a ModSet\n",
    "                    usTest = us.val        # change testing value to its attribute\n",
    "                if isinstance(s, ModSet):  # if member s is a ModSet\n",
    "                    sTest = s.val          # change testing value to its attribute\n",
    "                if usTest == sTest:  # compare us and s testing values on this run\n",
    "                    inUniqueSet = True  # if match, set existence flag to true\n",
    "            if not inUniqueSet:    # only add member s to uniqueSet if\n",
    "                uniqueSet.add(s)   # match is NOT detected\n",
    "        self.val = uniqueSet    # set .val to the uniqueSet from above\n",
    "    \n",
    "    # union join multiple items, enforce that ModSet members be unique\n",
    "    def union(self, *modSets):\n",
    "        for modSet in modSets:\n",
    "            self.val = self.val.union(modSet.val)  # this is union method from set class\n",
    "            self.removeDuplicates()  # removes duplicate-valued instances of ModSet members\n",
    "    \n",
    "    # set intersect multiple items\n",
    "    def intersection(self, *modSets):\n",
    "        for modSet in modSets:\n",
    "            self.val = self.val.intersection(modSet.val)\n",
    "    \n",
    "    # set difference multiple items. note: arg order matters here! \n",
    "    def difference(self, *modSets):\n",
    "        for modSet in modSets:\n",
    "            self.val = self.val.difference(modSet.val)\n",
    "    \n",
    "    # version of above method that returns a new instance,\n",
    "    # only takes two arguments.\n",
    "    def diffFunc(modSet1, modSet2):\n",
    "        return ModSet(set.difference(modSet1.val, modSet2.val))\n",
    "    \n",
    "    # Generate powerset via direct, binary mask approach\n",
    "    def powerSet_bin(self):\n",
    "        \n",
    "        ## Initialize local variables ##\n",
    "        S = list(self.val)         # convert to list for indexing\n",
    "        setSize = len(self.val)    # count number of members in source set\n",
    "        psetSize = pow(2, setSize) # calculate the number of elements in the power set\n",
    "        lastIndex = setSize - 1    # index value of last member\n",
    "        setIndices = range(0, setSize)  # make indices list for source set\n",
    "        psetIndices = range(0, psetSize) # make indices list for power set to be built\n",
    "        bMasks = [[False for s in setIndices] for p in psetIndices] # Initialize binary mask       \n",
    "        pSet = ModSet(set())  # initialize power set as empty ModSet() instance\n",
    "        pSet.pushDown(1)  # and nest it down one level for later joining\n",
    "        \n",
    "        ## Populate powerset with each subset, one at a time ##\n",
    "        for i in psetIndices: # loop through each member of power set\n",
    "            \n",
    "            ## Generate binary mask for subset i of power set ##\n",
    "            val = i  # assign current pSet index as current \"value\" of mask \n",
    "            for j in setIndices: # loop through each bit-digit of mask\n",
    "                \n",
    "                if (val >= pow(2, lastIndex - j)):  # if mask value >= value of current bit,\n",
    "                    bMasks[i][lastIndex - j] = True # then set corresp. mask bit to \"true\"\n",
    "                    val -= pow(2, lastIndex - j)    # subtract value of current bit from\n",
    "                                                    # mask value to determine next bit-digit\n",
    "            ## Form subset i of power set ##\n",
    "            # Use generator expression for compactness\n",
    "            dummySet = ModSet(set([S[k] for k in setIndices if bMasks[i][k] == True]))\n",
    "            dummySet.pushDown(1)  # nest ModSet instance down one level for union join\n",
    "            pSet.union(dummySet)  # include new subset in power set\n",
    "            \n",
    "        return pSet, bMasks  # return complete power set and binary masks as output\n",
    "    \n",
    "    # Generate powerset recursively.\n",
    "    def powerSet_rec(self):\n",
    "        \n",
    "        pSet = self.__copy__()      # Preserve self instance; its copy, pSet\n",
    "                                    # will be altered\n",
    "        pSet.pushDown(1)            # Nest pSet for later joining.\n",
    "              \n",
    "        if len(self.val) > 0:       # Recursion termination condition\n",
    "            for elSet in self.val:  # Iterate through members of set self.val\n",
    "                # Generate subset that remains after removing current\n",
    "                # element, elSet, from set self.val\n",
    "                dummySet = self.diffFunc(ModSet(set([elSet]))) \n",
    "                # To current powerset, append the powerset of the\n",
    "                # subset from previous step\n",
    "                pSet.union(dummySet.powerSet_rec())  # Self-call powerset method,\n",
    "                                                     # union join powerset of \n",
    "                                                     # dummySet with pSet\n",
    "            return pSet             # Return powerset at current \n",
    "                                    # level of recursion\n",
    "        else:\n",
    "            dummySet = ModSet(set())  # Generate instance of ModSet of empty set\n",
    "            dummySet.pushDown(1)      # Nest empty set down one level so it can\n",
    "            return dummySet           # be union joined in the recursion level\n",
    "                                      # above that called this run-through."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}