From b97069372fde4dcd29f54cf1ff994460187b90a6 Mon Sep 17 00:00:00 2001 From: "Kamal Bentahar (ab3735)" Date: Tue, 26 Mar 2019 14:07:59 +0000 Subject: [PATCH] Investigating 3SAT --- Investigating 3SAT.ipynb | 1033 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 1033 insertions(+) create mode 100644 Investigating 3SAT.ipynb diff --git a/Investigating 3SAT.ipynb b/Investigating 3SAT.ipynb new file mode 100644 index 0000000..694e748 --- /dev/null +++ b/Investigating 3SAT.ipynb @@ -0,0 +1,1033 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "_Kamal Bentahar_ (individual work)\n", + "\n", + "[https://github.coventry.ac.uk/ab3735/380CT/tree/master/SAT](https://github.coventry.ac.uk/ab3735/380CT/tree/master/SAT)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "_This is a guide report for the 380CT Assignment._" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Notation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let $x_1,x_2,\\ldots,x_n$ be Boolean **variables**, and let $\\phi$ be a Boolean formula written in 3-cnf (Conjunctive Normal Form) given by\n", + "\t$$\n", + "\t\\phi = c_1\\land c_2\\land \\cdots \\land c_\\ell,\n", + "\t$$\n", + "where each **clause** $c_m$ has the form $\\alpha\\lor \\beta \\lor \\gamma$, where each of $\\alpha, \\beta, \\gamma$ is a **literal**: a variable $x_i$ or its _negation_ $\\bar x_i$.\n", + "\n", + "\n", + "The ratio $\\ell/n$ is important for experiments, and will be denoted by $\\rho$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Definition of the problem" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Given $\\phi$ as above, the versions of the 3SAT are defined as follows:\n", + "\n", + "* **Decisional 3SAT**:\n", + "\n", + " _Decide if $\\phi$ is satisfiable._\n", + "\n", + " **NP-complete**, because 3SAT $\\in$ NP and SAT $\\leq_p$ 3SAT.\n", + " \n", + " * 3SAT $\\in$ NP: once an assignment of the variables is given (a certificate) we can quickly evaluate $\\phi$ in $O(\\ell)$ time to verify it.\n", + " * SAT $\\leq_p$ 3SAT: Reduction from the SAT problem, which is known to be NP-complete by the Levin-Cook theorem (Sipser, 2013, p. 304).\n", + "\n", + "* **Computational/Search 3SAT**:\n", + "\n", + " _If $\\phi$ is satisfiable then find a satisfying assignment._\n", + "\n", + " **NP-Hard**, because we can reduce the decision version of 3SAT to it trivially: if a solution is found then return **yes**, otherwise return **no**.\n", + "\n", + "* **Optimization 3SAT (Max 3SAT)**:\n", + "\n", + " _Find an assignment that minimizes the number of non-satisfying clauses._\n", + "\n", + " **NP-Hard**, because the optimization version of (decision) NP-complete problems are automatically NP-Hard. (using the same method sketched above for **Search 3SAT**)\n", + "\n", + "The facts about the complexity classes memberships can also be found in (Garey and Johnson, 1979)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Testing methodology" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* **Exhaustive search**:\n", + " average time for instances with increasing $n$.\n", + "* **Dynamic programming**:\n", + " average time for instances with increasing $\\ell$.\n", + "* **Greedy and meta-heuristics**:\n", + " quality of approximation with increasing $\\rho$ . Quality of approximation is calculated as the ratio of satisfied clauses to $\\ell$.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Random instances sampling strategy" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Following (Hoos & Stutzler, 2005), general 3SAT instances will be generated by selecting exactly 3 different literals from\n", + "$$\\{x_1, \\bar x_1, x_2, \\bar x_2, \\ldots, x_n, \\bar x_n\\}$$\n", + "uniformly at random.\n", + "\n", + "Do not allow clauses including both $x_i$ and $\\bar x_i$ (tautological clauses).\n", + "\n", + "\n", + "For *yes* instances, a random variable assignment is fixed first, then clauses are randomly constructed making sure each is satisfiable.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Code" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First start by importing relevant libraries." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "from random import sample, choice, shuffle, randint\n", + "from itertools import product\n", + "from numpy import arange\n", + "from time import time\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Base class definition" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "class SAT():\n", + " def __init__(self):\n", + " ''' Initialize an empty 3SAT instance: phi '''\n", + " self.variables = [] # Names of the variables :1,2,... -> x1,x2,...\n", + " self.clauses = [] # Clauses, e.g. [1,2,-3] -> x1 and x2 and not x3\n", + " self.nv = 0 # Number of variabes = len(self.variables)\n", + " self.nc = 0 # Number of clauses = len(self.clauses)\n", + " self.configuration = [] # Current value assignment of the variabes" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def sat_repr(self):\n", + " ''' To print out the SAT instance in readable form '''\n", + " r = \"SAT instance with \"+str(self.nv)\n", + " r += \" variables and \"+str(self.nc)+\" clauses: \\n\"\n", + " for c in self.clauses:\n", + " r+= str(c) + \"\\n\"\n", + " return r\n", + "SAT.__repr__ = sat_repr" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def set_variable(self, variable, value):\n", + " ''' Set the given 'variable' to the giben 'value' '''\n", + " self.configuration[ abs(variable)-1 ] = value\n", + "def flip_variable(self,variable):\n", + " '''\n", + " Flip the value of 'variable'\n", + " '''\n", + " idx = abs(variable)-1\n", + " self.configuration[ idx ] = not self.configuration[ idx ]\n", + "def get_variable(self,variable):\n", + " ''' Return the value assigned to 'variable' '''\n", + " return self.configuration[ abs(variable)-1 ]\n", + "\n", + "SAT.set_variable = set_variable\n", + "SAT.flip_variable = flip_variable\n", + "SAT.get_variable = get_variable" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Evaluation" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def clause_value(self, clause):\n", + " ''' Return the value of 'clause' '''\n", + " x1 = self.get_variable( clause[0] )\n", + " x2 = self.get_variable( clause[1] )\n", + " x3 = self.get_variable( clause[2] )\n", + " # if negative then negate boolean value of the variable\n", + " if clause[0]<0: x1 = not x1\n", + " if clause[1]<0: x2 = not x2\n", + " if clause[2]<0: x3 = not x3\n", + " return (x1 or x2 or x3)\n", + "\n", + "def value(self):\n", + " ''' evaluate phi - stop as soon as a clause is False '''\n", + " for c in self.clauses:\n", + " # if any clause is False then phi is not satisfied\n", + " if self.clause_value( c ) == False:\n", + " return False\n", + " # if no clause is False then phi is satisfied\n", + " return True \n", + "\n", + "def value_full(self):\n", + " ''' evaluate phi - do not stop until we evalaute all clauses '''\n", + " val = True\n", + " for c in self.clauses:\n", + " val &= self.clause_value( c )\n", + " return val\n", + "\n", + "def ratio_unsatisfied(self):\n", + " ''' calculate ratio of satisfied clauses '''\n", + " count = 0\n", + " for c in self.clauses:\n", + " if self.clause_value( c ) == False:\n", + " count += 1\n", + " return count/self.nc\n", + "\n", + "SAT.clause_value = clause_value\n", + "SAT.value = value\n", + "SAT.value_full = value_full\n", + "SAT.ratio_unsatisfied = ratio_unsatisfied" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Methods for random constructions" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def random_clause(self):\n", + " ''' helper function to generate a clause randomly '''\n", + " c = sample(self.variables,3)\n", + " c[0] *= choice([-1,1])\n", + " c[1] *= choice([-1,1])\n", + " c[2] *= choice([-1,1])\n", + " return c\n", + "\n", + "def random_instance(self, nv, nc):\n", + " ''' build a random 3SAT instance with nv variables and nc clauses '''\n", + " self.nv = nv\n", + " self.variables = list(range(1,nv+1))\n", + " self.configuration = [False]*nv\n", + " self.nc = nc\n", + " # generate nc clauses: (... or ... or ...)\n", + " self.clauses = []\n", + " for c in range(nc):\n", + " self.clauses.append( self.random_clause() )\n", + "\n", + "def random_configuration(self):\n", + " '''\n", + " Build a random configuration (= assignment of the variables)\n", + " Return value of objective function.\n", + " '''\n", + " self.configuration = [choice([True,False]) for nv in range(self.nv)]\n", + "\n", + "def random_yes_instance(self, nv, nc):\n", + " ''' build a random 'yes' 3SAT instance with nv variables and nc clauses '''\n", + " # choose a random configuration to be satisfied\n", + " self.nv = nv\n", + " self.variables = list(range(1,nv+1))\n", + " self.random_configuration()\n", + " # choose satisfying clauses\n", + " self.nc = nc\n", + " self.clauses = []\n", + " for nc in range(self.nc):\n", + " c = self.random_clause() #sample(self.variables,3)\n", + " # alter c until it becomes satisfying\n", + " while self.clause_value(c) == False:\n", + " c[ choice([0,1,2]) ] *= -1 # by flipping terminals (negating)\n", + " self.clauses.append(c)\n", + " if self.value() != True:\n", + " print(\"something has gone wrong!\")\n", + "\n", + "SAT.random_clause = random_clause\n", + "SAT.random_instance = random_instance\n", + "SAT.random_configuration = random_configuration\n", + "SAT.random_yes_instance = random_yes_instance" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Solution methods" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exact methods -- Exhaustive search" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Pseudo-code:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "1. **for all** possible variable assignments of $x_1,x_2,\\ldots,x_n$ **do**\n", + "2. $\\quad$ **if** $\\phi(x_1,x_2,\\ldots,x_n)$ evaluates to True **then**\n", + "3. $\\qquad$ **return** True\n", + "4. $\\quad$ **end if**\n", + "5. **end for**\n", + "6. **return** False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are $2^n$ possible assignments, and each evaluation of $\\phi$ costs $O(\\ell)$. So this algorithm costs $$O(\\ell\\, 2^n).$$" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def exhaustive_search(self):\n", + " '''\n", + " Solve in the 3SAT problem instance using exhaustive search\n", + " -- stop as soon as a solution is found\n", + " '''\n", + " # iterate over all the possible Boolean variable assignments\n", + " for self.configuration in product([True,False],repeat=self.nv):\n", + " if self.value() == True:\n", + " return True\n", + " return False\n", + "\n", + "def exhaustive_search_full(self):\n", + " ''' Solve in the 3SAT problem instance using exhaustive search '''\n", + " # iterate over all the possible Boolean variable assignments\n", + " self.decision = False\n", + " for self.configuration in product([True,False],repeat=self.nv):\n", + " # if any configuration satisfies it then we get True\n", + " self.decision |= self.value_full()\n", + " return self.decision\n", + "\n", + "SAT.exhaustive_search = exhaustive_search\n", + "SAT.exhaustive_search_full = exhaustive_search_full" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Testing" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Average time for randomly generated instances with $\\rho=3$." + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "n\tExhaustive\n", + "10\t0.0016361308097839356\n", + "11\t0.002511744499206543\n", + "12\t0.004523138999938964\n", + "13\t0.009381508827209473\n", + "14\t0.017256977558135985\n", + "15\t0.02649304151535034\n", + "16\t0.054238312244415283\n", + "17\t0.07540356397628784\n", + "18\t0.1821987247467041\n", + "19\t0.4292074728012085\n", + "20\t0.6677624154090881\n" + ] + } + ], + "source": [ + "instance = SAT() # global variable... [TODO?]\n", + "pnts_n = []\n", + "pnts_t = []\n", + "def time_Exhaustive():\n", + " # test exhaustive search\n", + " print(\"n\\tExhaustive\" ) # header\n", + " max_repeats = 100\n", + " n = 10\n", + " t0 = t1 = 0\n", + " while t1-t0<60: # in seconds; if it takes too long then stop testing\n", + " t0 = time()\n", + " for repeats in range(max_repeats): # e.g. average over 100 instances\n", + " instance.random_yes_instance(n,3*n) # rho=3\n", + " decision = instance.exhaustive_search()\n", + " t1 = time()\n", + " # record average time\n", + " print( str(n)+\"\\t\"+str((t1-t0)/max_repeats) )\n", + " pnts_n.append( n )\n", + " pnts_t.append( (t1-t0)/max_repeats )\n", + " n += 1\n", + "time_Exhaustive()" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAHThJREFUeJzt3Xt0VeWd//H3lyAiXhAkiAUSULEubS2OEcfa1ltRBi/oONMfSr2iFCr60zW/VufntNOuWcyMnc7PaZCIEREvUcRBkSqKlw51qJcmWKUCQ6UqECoS8QpBrt/fH89Jcwg5yUlyznnOOfm81jorZ+/95OzvJvBhZ+9nP4+5OyIiUlx6xC5AREQyT+EuIlKEFO4iIkVI4S4iUoQU7iIiRUjhLiJShNIKdzMbY2arzWyNmd3ayvYfmNkbiddbZrbbzPpnvlwREUmHtdfP3cxKgD8Ao4F6oBa41N1Xpmh/AXCzu5+V4VpFRCRN6Zy5jwLWuPs77r4DmAuMa6P9pcAjmShOREQ6p2cabQYD65OW64FTWmtoZn2AMcDU9j50wIABPmzYsDR2LyIiTZYtW/ahu5e21y6dcO+IC4DfuPtHrW00s0nAJICysjLq6uoyvHsRkeJmZmvTaZfOZZkNwNCk5SGJda0ZTxuXZNy92t0r3L2itLTd/3hERKST0gn3WmCEmQ03s16EAF/YspGZ9QVOB57MbIkiItJR7V6WcfddZjYVWAyUALPdfYWZTU5sn5loejHwnLtvzVq1IiKSlna7QmZLRUWF65q7iEjHmNkyd69or52eUBURKUIKdxGRXKmpgWHDoEeP8LWmJmu7ynRXSBERaU1NDUyaBI2NYXnt2rAMMGFCxnenM3cRkVy47bbmYG/S2BjWZ4HCXUQkF9at69j6LlK4i4jkQllZx9Z3kcJdRCQX/uEf9l3Xpw9Mm5aV3SncRURyYWvi+c5Bg8AMysuhujorN1NBvWVERLJv926YPh1OOw2WLs3JLnXmLiKSbc88A3/8I9x4Y852qXAXEcm2ykoYPBguvjhnu1S4i4hk08qV8Pzz8P3vw3775Wy3CncRkWy6807Yf3+47rqc7lbhLiKSLZ98AvffD5ddBjmeoEjhLiKSLbNnhyEGcngjtYnCXUQkG3bvDpdkvvUtGDky57tXuIuIZMPTT8O770Y5aweFu4hIdlRWwtChMG5clN0r3EVEMm3FCnjxRbj+eugZZyAAhbuISKZNnw69e8O110YrQeEuIpJJH38MDzwA3/0uHHZYtDLSCnczG2Nmq81sjZndmqLNGWb2hpmtMLNfZ7ZMEZECce+9sG0b3HBD1DLavRhkZiXADGA0UA/UmtlCd1+Z1OZQoAoY4+7rzGxgtgoWEclbTd0fzzgDTjghainpnLmPAta4+zvuvgOYC7S8/XsZ8Li7rwNw902ZLVNEpAAsXBgmvo7U/TFZOuE+GFiftFyfWJfsGKCfmS0xs2VmdkVrH2Rmk8yszszqGhoaOlexiEi+qqwMk3BccEHsSjJ2Q7UncBJwHnAu8CMzO6ZlI3evdvcKd68ozfE4CyIiWbV8OSxZErX7Y7J0KtgADE1aHpJYl6we2OzuW4GtZvYS8DXgDxmpUkQk302fDgccABMnxq4ESO/MvRYYYWbDzawXMB5Y2KLNk8A3zKynmfUBTgFWZbZUEZE8tXkzPPQQXH459O8fuxogjTN3d99lZlOBxUAJMNvdV5jZ5MT2me6+ysyeBZYDe4BZ7v5WNgsXEckbs2bBF19E7/6YzNw9yo4rKiq8rq4uyr5FRDJm1y448kgYMSIMOZBlZrbM3Svaaxf/qr+ISCF78klYvz70b88jGn5ARKQrKith+HA477zYlexF4S4i0llvvAEvvQRTp0JJSexq9qJwFxHprOnToU8fuOaa2JXsQ+EuItIZH34INTVw5ZVw6KGxq9mHwl1EpDPuuQe2bw+XZPKQwl1EpKN27oSqKhg9Go47LnY1rVJXSBGRjlqwAOrr4a67YleSks7cRUQ6qrISjjoKxo6NXUlKCncRkY54/XVYujRca++RvxGav5WJiOSjyko48EC4+urYlbRJ4S4ikq5Nm+CRR+Cqq6Bv39jVtEnhLiKSrupq2LEjb7s/JlO4i4iko6n747nnwrHHxq6mXeoKKSKSjvnz4f33w9jtBUBn7iIi6aisDGO2jxkTu5K0KNxFRNpTWwuvvBJmWsrj7o/JCqNKEZGYpk+Hgw8Og4QVCIW7iEhbNm6EuXNDv/ZDDoldTdoU7iIibamuDj1lCqD7YzKFu4hIKjt2hMHBxo4NN1MLSFrhbmZjzGy1ma0xs1tb2X6GmX1qZm8kXj/OfKkiIjn2n/8ZLsvceGPsSjqs3X7uZlYCzABGA/VArZktdPeVLZr+t7ufn4UaRUTiqKyEL385jNteYNI5cx8FrHH3d9x9BzAXGJfdskREInvttfAqoO6PydKpeDCwPmm5PrGupa+b2XIze8bMjm/tg8xskpnVmVldQ0NDJ8oVEcmRysrQO+aKK2JX0imZ+u/odaDM3U8ApgMLWmvk7tXuXuHuFaWlpRnatYhIhv3pTzBvHlxzTejfXoDSCfcNwNCk5SGJdX/m7p+5+5bE+0XAfmY2IGNViojk0t13w+7dcP31sSvptHTCvRYYYWbDzawXMB5YmNzAzAaZmSXej0p87uZMFysiknXbt8PMmXDeeXD00bGr6bR2e8u4+y4zmwosBkqA2e6+wswmJ7bPBP4GmGJmu4BtwHh39yzWLSKSHfPmhUk5CrD7YzKLlcEVFRVeV1cXZd8iIq1yh5NPhsZGWLECwgWJvGJmy9y9or12Gs9dRKTJq6/CsmVhUo48DPaOKLzOmyIi2VJZGeZGvfzy2JV0mcJdRARgw4Yw3MC118JBB8WupssU7iIiEHrIFHj3x2QKdxGRL74IfdsvvBCGD49dTUYo3EVEHn0UGhoKvvtjMoW7iHRv7vCLX8Dxx8OZZ8auJmPUFVJEureXX4bf/S5clinw7o/JdOYuIt1bZSX06wcTJsSuJKMU7iLSfa1fD/Pnh+6PBx4Yu5qMUriLSPd1113hmvv3vx+7koxTuItI97RtG1RXw7hxMGxY7GoyTuEuIt1LTU0I8z59YPNmOO642BVlhcJdRLqPmhqYNAnWrm1ed8cdYX2RUbiLSPdx221hON9kjY1hfZFRuItI97FuXcfWFzCFu4h0H2VlHVtfwBTuItJ9TJsG++2397o+fcL6IqNwF5Hu4zvfCQ8r9e4dhhooLw/dIYvs6VTQ2DIi0p388pfwySewcCFccEHsarIqrTN3MxtjZqvNbI2Z3dpGu5PNbJeZ/U3mShQRyZAZM8L19bFjY1eSde2Gu5mVADOAvwKOAy41s316/Sfa3Q48l+kiRUS6bNUq+NWvYPJkKCmJXU3WpXPmPgpY4+7vuPsOYC4wrpV2NwDzgU0ZrE9EJDNmzoRevWDixNiV5EQ64T4YWJ+0XJ9Y92dmNhi4GLgrc6WJiGTIli0wZw787d/CwIGxq8mJTPWW+Q/gFnff01YjM5tkZnVmVtfQ0JChXYuItOPhh+Gzz4py9MdU0uktswEYmrQ8JLEuWQUw18IsJgOAsWa2y90XJDdy92qgGqCiosI7W7SISNrcw43UkSPh1FNjV5Mz6YR7LTDCzIYTQn08cFlyA3f/83ThZjYHeKplsIuIRPHyy7B8eejPXkTT6LWn3XB3911mNhVYDJQAs919hZlNTmyfmeUaRUQ6r6oK+vaFyy5rv20RSeshJndfBCxqsa7VUHf3q7pelohIBmzaBI89Fq61F9k0eu3R8AMiUrzuvRd27oQpU2JXknMKdxEpTrt3h77tZ58NX/5y7GpyTuEuIsXp6afDOO3dqPtjMoW7iBSnqioYPBguvDB2JVEo3EWk+Lz9NixeDN/7HvTsnoPfKtxFpPjMnBlC/dprY1cSjcJdRIpLYyPcdx9ccgkccUTsaqJRuItIcZk7Fz7+uNveSG2icBeR4tE0jszxx8M3vxm7mqi6550GESlOtbXw+uuhp0w3GkemNTpzF5HiMWMGHHQQfPe7sSuJTuEuIsXhww/h0Ufhiivg4INjVxOdwl1EisN998H27d3+RmoThbuIFL49e+Cuu+D008PNVFG4i0gRePZZePddnbUnUbiLSOGrqoJBg+Cii2JXkjcU7iJS2N59FxYtgkmToFev2NXkDYW7iBS2mTOhRw+47rrYleQVhbuIFK4vvgizLY0bB0OGxK4mryjcRaRwPfYYbN4M118fu5K8o3AXkcI1Y0aYQu/MM2NXknfSCnczG2Nmq81sjZnd2sr2cWa23MzeMLM6M/tG5ksVEUmybBm89lro/tjNx5FpTbsDh5lZCTADGA3UA7VmttDdVyY1exFY6O5uZicA84Bjs1GwiAgQHlrq0ycMNyD7SOfMfRSwxt3fcfcdwFxgXHIDd9/i7p5YPBBwRESy5eOP4eGHwwBhhx4au5q8lE64DwbWJy3XJ9btxcwuNrP/AZ4Grmntg8xsUuKyTV1DQ0Nn6hURgTlzYNs2mDIldiV5K2M3VN39CXc/FrgI+KcUbardvcLdK0pLSzO1axHpTprGkfn612HkyNjV5K10wn0DMDRpeUhiXavc/SXgSDMb0MXaRET29eKL8Pbb6v7YjnTCvRYYYWbDzawXMB5YmNzAzI42C7erzewvgP2BzZkuVkSEGTOgtDRMgC0ptdtbxt13mdlUYDFQAsx29xVmNjmxfSZwCXCFme0EtgH/K+kGq4hIZqxbB7/8JdxyC+y/f+xq8lpac6i6+yJgUYt1M5Pe3w7cntnSRERaqK4OX7/3vbh1FAA9oSoihWH7drjnHjj/fCgvj11N3lO4i0hhePxx2LRJE3KkSeEuIoWhqgqOPhpGj45dSUFQuItI/lu+HJYuDQ8t9VBspUN/SiKS/6qqoHdvuOqq2JUUDIW7iOS3Tz+Fhx6CSy+F/v1jV1MwFO4ikt8efBC2btUTqR2kcBeR/OUeLsmMGgUnnRS7moKS1kNMIiJRLFkCq1aFUSClQ3TmLiL5q6oqXGf/zndiV1JwFO4ikp82bIAnnoCJE+GAA2JXU3AU7iKSn+65J4zdrnFkOkXhLiL5Z+fOMEjYmDFw1FGxqylIuqEqIvlnwQJ4//1w9i6dojN3Eck/VVUwbFg4c5dOUbiLSH5ZuTJ0gZw8GUpKYldTsBTuIpJfqqrCLEsTJ8aupKAp3EUkf3z+OTzwQOjXPmBA7GoKmsJdRPJHTU0IeE3I0WUKdxHJD03jyJx4IpxySuxqCp66QopIfli6FH7/e5g1C8xiV1Pw0jpzN7MxZrbazNaY2a2tbJ9gZsvN7Pdm9rKZfS3zpYpIUauqgr59w7jt0mXthruZlQAzgL8CjgMuNbPjWjR7Fzjd3b8K/BNQnelCRaSIbdwI8+fD1VdDnz6xqykK6Zy5jwLWuPs77r4DmAuMS27g7i+7+8eJxVeBIZktU0SK2qxZYciBKVNiV1I00gn3wcD6pOX6xLpUJgLPtLbBzCaZWZ2Z1TU0NKRfpYgUr1274O67YfRoOOaY2NUUjYz2ljGzMwnhfktr29292t0r3L2itLQ0k7sWkUL11FNQX6/ujxmWTrhvAIYmLQ9JrNuLmZ0AzALGufvmzJQnIkWrpiaMH3PxxWGYgc8/j11RUUmnK2QtMMLMhhNCfTxwWXIDMysDHgcud/c/ZLxKESkuNTUwaRI0Nobl3bvDWDI9esCECXFrKxLtnrm7+y5gKrAYWAXMc/cVZjbZzCYnmv0YOAyoMrM3zKwuaxWLSOG77bbmYG/S2BjWS0aYu0fZcUVFhdfV6f8AkW6pR4/wRGpLZmH2JUnJzJa5e0V77TT8gIjkXqoOFWVlua2jiCncRSS3liyBjz7ad4iBPn1g2rQoJRUjhbuI5M4rr8D554f+7FVVUF4eQr68PMyZqpupGaOBw0QkN5YtC9PmfelL8MILcMQRoYeMZIXO3EUk+956C845B/r1gxdfDMEuWaVwF5HsWr0avv1t6N07BPvQoe1/j3SZLsuISPa8+y6cfXbo3rhkCRx1VOyKug2Fu4hkR319CPbGxhDsxx4bu6JuReEuIpn3wQch2DdvDpdiTjghdkXdjsJdRDJr8+Zwjb2+Hp57DirafZhSskDhLiKZ88knoVfM22/DokVw2mmxK+q2FO4ikhlbtsDYsWGS6wUL4KyzYlfUrSncRaTrtm2DCy6A3/4WHn00hLxEpXAXka7Zvh3++q/h17+GBx+ESy6JXZGgcBeRrti5Ey69FJ59NkxyrbFh8oaeUBWRztm9G668Ep54AiorYeLE2BVJEoW7iHTcnj1w3XXwyCNw++1www2xK5IWFO4i0jHucOONcN998I//CD/8YeyKpBUKdxFJn3sI8xkz4Ac/COEueUnhLiLp++lP4ec/h+uvD5djWs6mJHkjrXA3szFmttrM1pjZra1sP9bMXjGz7Wb2fzJfpohE97OfhXC/+upwA1XBntfa7QppZiXADGA0UA/UmtlCd1+Z1Owj4EbgoqxUKSJx3Xkn3HILjB8P99wDPfRLf75L5yc0Cljj7u+4+w5gLjAuuYG7b3L3WmBnFmoUkZjuvTf0hrnoInjgASgpiV2RpCGdcB8MrE9ark+sE5Fi9/DDocvjmDEwdy7st1/siiRNOf3dyswmmVmdmdU1NDTkctci0lGPPw5XXAFnnBHe779/7IqkA9IJ9w1A8qSHQxLrOszdq929wt0rSktLO/MRIpILixaF6+ujRsHChXDAAbErkg5KJ9xrgRFmNtzMegHjgYXZLUtEovnVr8JAYCecAM88AwcdFLsi6YR2e8u4+y4zmwosBkqA2e6+wswmJ7bPNLNBQB1wCLDHzG4CjnP3z7JYu4hk2m9+AxdeCCNGwOLF0Ldv7Iqkk9IaFdLdFwGLWqybmfR+I+FyjYgUkpoauO02WLcODj88zKRUVgbPPw+HHRa7OukCDfkr0l3V1MCkSdDYGJY3bgwPJk2dCoMGxa1NukxPIoh0R01jxDQFe/L6f//3ODVJRunMXaQ7eP99qK2Furrmrx9+2HrbdetyW5tkhcJdpNh89FEI76Ygr62FDYneyz16wFe+Em6aLlgQ2rZUVpbbeiUrFO4ihWzLFvjd75pDvLYW/vjH5u0jRsDpp8PJJ4fXiSdCnz5h21ln7X3NHcK2adNyewySFQp3kXyQ3GulrCwEbMv5SLdvh+XL9w7yVavCrEgQvq+iAq69NgT5SSfBoYem3mfT57e3XylI5u5RdlxRUeF1dXVR9i2SV1r2WoFwBv2jH8HAgc1Bvnx5mJAaoLS0+Wz85JNDqB9+eJz6JafMbJm7V7TbTuEuEtmwYbB2berthxwSwjs5zIcO1Xjq3VS64a7LMiIx7NgBr7wSHhZqK9hXr4ajj9b46dJhCneRXHAP18effz68liyBrVvD2Oj77x+up7dUXg7HHJPzUqU4KNxFsmXTJnjhheZAb+qOOGIEXHUVjB4dhtN96in1WpGMU7iLZMq2bbB0aQjy556DN98M6/v3h7PPhnPOCYFeXr7396nXimSBbqiKdNaePaEHS1OYL10KX3wRZis67bTmMD/xRE1NJxmjG6oi2bBhQ/Nlluefh6YZxY4/HqZMCWH+rW/BgQfGrVO6PYW7SLKWDxP9+Meh/3jT2fmqVaHd4YeHM/NzzoFvfxu+9KW4dYu0oHAXaXL//eHse9u2sLx2LUycGN737h3OyCdODGfnX/2q+plLXlO4S/5J51H8dG3bBh980PZr48bw9dNPW/+MgQND0Pfu3fljEskxhbvkl5aP4q9dG5ahOeC3bGk/sJten3/e+n769QuXVg4/HEaODF+nT2+9bUODgl0KjnrLSHyNjbB5cxhffMyY0D+8pV69YMiQENhbt7b+OYcd1hzYLV+DBjW/HzgwfF5LqYYBKC+H997ryhGKZIx6y0jXdebySGNjCOmmsE7nfdM17rbs2AGnnpo6vEtLQxfErpg2TQ8TSdFQuOe7TF5/Tpc7PPDAvjcXr7kGFi+GI49MHdZtBXX//uHsesCAMPDVyJHh/YABzesnTw5n5y2Vl8NDD2XneJvoYSIpImmFu5mNAX4BlACz3P1fW2y3xPaxQCNwlbu/nuFa4wRdzH2nc/3ZPTw4s2VL8+vzz/deTvVqq92uXfvWs2MHPPhg6CXSr19zKA8dGh7UaQro5LBuet+vH/RM46/b1q1xz54nTFCYS1Fo95q7mZUAfwBGA/VALXCpu69MajMWuIEQ7qcAv3D3U9r63A5fc0815nV1ddf+MbqH15494bV7d/P7pte8eXDzzXuflfbuDT/5CZx7bhhju+m1Y8fey115zZu37wTGEEJy4MDmIG6arKE9ZnDQQalfBx/c/P6f/zn1Z+zcmd0nLmP+Jy6S5zI2nruZnQr8xN3PTSz/PYC7/0tSm7uBJe7+SGJ5NXCGu7+f6nM7HO6pbnb17BlutKUK5vbWR7qhvM8x7Lffvq/161N/z8SJqYM5VXAfcED6fbN1c1EkL2XyhupgIDll6gln5+21GQzsFe5mNgmYBFDW0Ul4U83IvmtXeLikR4/mV0nJ3sup1qXb9qabWt+3Gcyf33owN7169Wp7e8+eqQO3rYCdNatjf34dpZuLIgUtpzdU3b0aqIZw5t6hby4rSx1099+fifJSu+OO1vddVgYXX5y9/cYMWN1cFClo6UzvsgEYmrQ8JLGuo226Ztq05lnbm+Qq6GLte8KEcE+hvDyc3ZeXd/0eQ0f3/9574RLWe+8p2EUKSDrhXguMMLPhZtYLGA8sbNFmIXCFBX8JfNrW9fZOiRl0sfetgBWRDkrrCdVEb5j/IHSFnO3u08xsMoC7z0x0hbwTGEPoCnm1u7d5t1RPqIqIdFxGn1B190XAohbrZia9d+D6jhYpIiLZoSnVRUSKkMJdRKQIKdxFRIqQwl1EpAhFG8/dzBqAVp4MSssA4MMMllMIdMzdg465e+jKMZe7e2l7jaKFe1eYWV06XYGKiY65e9Axdw+5OGZdlhERKUIKdxGRIlSo4V4du4AIdMzdg465e8j6MRfkNXcREWlboZ65i4hIG/I+3M1stpltMrO3ktb1N7PnzeztxNd+MWvMtBTH/G9m9j9mttzMnjCzQ2PWmGmtHXPStr8zMzezATFqy5ZUx2xmNyR+1ivM7Gex6suGFH+3R5rZq2b2hpnVmdmomDVmkpkNNbP/MrOViZ/n/06sz3qG5X24A3MIo00muxV40d1HAC8mlovJHPY95ueBr7j7CYQ5bf8+10Vl2Rz2PWbMbChwDpBiKq6CNocWx2xmZwLjgK+5+/HAzyPUlU1z2Pfn/DPgp+4+EvhxYrlY7AL+zt2PA/4SuN7MjiMHGZb34e7uLwEftVg9Dmiaful+4KKcFpVlrR2zuz/n7rsSi68SJkQpGil+zgB3AD8Eiu7mUIpjngL8q7tvT7TZlPPCsijFMTtwSOJ9X+BPOS0qi9z9fXd/PfH+c2AVYQrSrGdY3od7CocnTQayETg8ZjERXAM8E7uIbDOzccAGd38zdi05dAzwTTN7zcx+bWYnxy4oB24C/s3M1hN+Uym230oBMLNhwInAa+Qgwwo13P8sMZZ80Z3VpWJmtxF+1auJXUs2mVkf4P8Sfk3vTnoC/Qm/wv8AmJeYDKeYTQFudvehwM3AvZHryTgzOwiYD9zk7p8lb8tWhhVquH9gZkcAJL4W1a+uqZjZVcD5wAQv/j6sRwHDgTfN7D3CZajXzWxQ1Kqyrx543IPfAnsI45AUsyuBxxPvHwOK5oYqgJntRwj2GndvOs6sZ1ihhvtCwl8IEl+fjFhLTpjZGMK15wvdvTF2Pdnm7r9394HuPszdhxFC7y/cfWPk0rJtAXAmgJkdA/Si+AfV+hNweuL9WcDbEWvJqMRvXfcCq9z9/yVtyn6GuXtev4BHgPeBnYR/4BOBwwh3mN8GXgD6x64zB8e8BlgPvJF4zYxdZ7aPucX294ABsevMwc+5F/AQ8BbwOnBW7DpzcMzfAJYBbxKuR58Uu84MHu83CJdclif92x2biwzTE6oiIkWoUC/LiIhIGxTuIiJFSOEuIlKEFO4iIkVI4S4iUoQU7iIiRUjhLiJShBTuIiJF6P8DzZlm7uLy/48AAAAASUVORK5CYII=" + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "plt.plot(pnts_n, pnts_t, 'ro-')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "from scipy.optimize import curve_fit\n", + "\n", + "def func(x, a, b):\n", + " return a * x * (2 ** x) + b\n", + "\n", + "popt, pcov = curve_fit(func, pnts_n, pnts_t)" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'3.31991597899e-08 * n * 2^n + 0.0134577258739'" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "str(popt[0])+\" * n * 2^n + \"+str(popt[1])" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXJ5OEkBD2yBYgAQPKlgARFUVEi2VTcWGTulzbctFa9dfaipd6215La6utt7fScnm0XrdUqOKCioIbiiIIKPsSwh4EDGFPQpLJfH5/zCQOIcsQZnJmJp/n4zGPs2bO5yT65jvfOed7RFUxxhgTXWKcLsAYY0zwWbgbY0wUsnA3xpgoZOFujDFRyMLdGGOikIW7McZEIQt3Y4yJQhbuxhgThSzcjTEmCsU6deD27dtrWlqaU4c3xpiItGbNmsOqmlLffo6Fe1paGqtXr3bq8MYYE5FEZE8g+1m3jDHGRCELd2OMiUIW7sYYE4Uc63OvSXl5Ofn5+Zw+fdrpUkwNEhISSE1NJS4uzulSjDH1CKtwz8/PJzk5mbS0NETE6XKMH1WlsLCQ/Px80tPTnS7HGFOPsOqWOX36NO3atbNgD0MiQrt27exTlTERIqzCHbBgD2P2tzEmcgQU7iIySkS2iUieiMyoYfvPRGSt77VRRCpEpG3wyzXGmMhVUFTERU8/zY/efjvkx6o33EXEBcwGRgN9gCki0sd/H1V9QlWzVDULeAT4WFWPhKLgUBMRvve971Utu91uUlJSGDdunINVNZ7du3fzz3/+0+kyjIlKK/fvZ1thIZsPHw75sQJpuQ8B8lR1p6qWAfOAG+vYfwrwUjCKq09ODqSlQUyMd5qTc/7vmZSUxMaNGykpKQHgvffeo0uXLuf/xg3gdrsb/ZgW7saEzuf79gFwWSNkSiDh3gXY57ec71t3FhFJBEYBC86/tLrl5MC0abBnD6h6p9OmBSfgx4wZw9u+j00vvfQSU6ZMqdpWVFTE3XffzZAhQxg4cCBvvPEG4A3FYcOGMWjQIAYNGsTy5csBOHDgAFdddRVZWVn069ePZcuWAdCiRYuq93zllVe46667ALjrrruYPn06l156KT//+c9rPd6zzz7L+PHjGTlyJGlpaTz99NP86U9/YuDAgVx22WUcOeL94LRjxw5GjRrF4MGDGTZsGFu3bq06zv3338/QoUPp0aMHr7zyCgAzZsxg2bJlZGVl8dRTT53/L9MYU2XF/v0AXJqaGvqDqWqdL+BW4O9+y7cDT9ey7yTgzTreaxqwGljdrVs3rW7z5s1nratN9+6q3lg/89W9e8BvUaOkpCRdt26d3nLLLVpSUqKZmZn60Ucf6dixY1VV9ZFHHtEXXnhBVVWPHj2qGRkZeurUKS0qKtKSkhJVVc3NzdXBgwerquqTTz6pv/nNb1RV1e1264kTJ6qOU+nll1/WO++8U1VV77zzTh07dqy63e46j/d///d/2rNnTz1x4oR+88032rJlS/3b3/6mqqoPPvigPvXUU6qqes0112hubq6qqq5YsUJHjBhRdZxbb71VKyoqdNOmTdqzZ09V1TPOtSbn8jcyxnzLXVGhLX77W+VXv9KvfTnQEMBqrSe3VTWg69z3A139llN962oymTq6ZFR1LjAXIDs7WwM4dq327j239ediwIAB7N69m5deeokxY8acsW3JkiUsXLiQJ598EvBevrl37146d+7Mfffdx9q1a3G5XOTm5gJwySWXcPfdd1NeXs748ePJysqq9/gTJkzA5XLVeTyAESNGkJycTHJyMq1ateL6668HoH///qxfv55Tp06xfPlyJkyYUPXepaWlVfPjx48nJiaGPn36cOjQoYb+uowxAdhcUMCpsjK6t2pFp+TkkB8vkHBfBWSISDreUJ8M3FZ9JxFpBQwHvld9Wyh06+btiqlpfTDccMMNPPTQQyxdupTCwsKq9arKggUL6N279xn7/+pXv6JDhw6sW7cOj8dDQkICAFdddRWffPIJb7/9NnfddRc/+clPuOOOO864rLD6teNJSUn1Hm/lypU0a9asajkmJqZqOSYmBrfbjcfjoXXr1qxdu7bGc/T/eW+DwBgTKivy8wG4rDG6ZAigz11V3cB9wGJgC/AvVd0kItNFZLrfrjcBS1S1KDSlnmnWLEhMPHNdYqJ3fTDcfffd/PKXv6R///5nrP/ud7/LX/7yl6ow/OqrrwA4fvw4nTp1IiYmhhdeeIGKigoA9uzZQ4cOHfjhD3/ID37wA7788ksAOnTowJYtW/B4PLz22mu11lHb8QLRsmVL0tPTefnllwFvgK9bt67On0lOTubkyZMBH8MYE5iwC3cAVV2kqr1UtaeqzvKtm6Oqc/z2eVZVJ4eq0OqmToW5c6F7dxDxTufO9a4PhtTUVO6///6z1j/66KOUl5czYMAA+vbty6OPPgrAvffey3PPPUdmZiZbt26tan0vXbqUzMxMBg4cyPz583nggQcAePzxxxk3bhxDhw6lU6dOtdZR2/EClZOTwz/+8Q8yMzPp27dv1ReytRkwYAAul4vMzEz7QtWYIKr8MrWxwl2c+jienZ2t1R/WsWXLFi6++GJH6jGBsb+RMefu2OnTtPn974l3uTgxYwbNYhs+rJeIrFHV7Pr2C7vhB4wxJtqs8rXau8Z2pPeFsUG9N6c2YTUqpDHGRKPK/vY9n6Xi9l0IUnlvDgSvO9mftdyNMSbEKvvb3bvP7G8vLoaZM0NzTAt3Y4wJIVWtarmTf/aXqcG4N6cmFu7GGBNCeUeOcKSkhJjiJDjW6qztwbo3pzoLd2OMCaHlvsHCBqakkph45jMRgnlvTnUW7tW4XC6ysrKqXrt372b16tVV17wvXbq0alAwgNdff53Nmzef83H8Bw7zd/DgQSZPnkzPnj0ZPHgwY8aMqRrKwBgTeT719btMHtotpPfmVGdXy1TTvHnzs27XT0tLIzvbe1np0qVLadGiBUOHDgW84T5u3Dj69Olz1nudK1Xlpptu4s4772TevHkArFu3jkOHDtGrV6+A3qOioqJqXBpjjPOW+cJ9WLduXDo0dGFenbXcA7B06VLGjRvH7t27mTNnDk899RRZWVl8/PHHLFy4kJ/97GdkZWWxY8eOWofY3bVrF5dffjn9+/fnF7/4RY3H+eijj4iLi2P69G9HdcjMzGTYsGFVNVS67777ePbZZwHvPz4PP/wwgwYN4oknnmDIkCFV++3evbtqCIU1a9YwfPhwBg8ezHe/+10OHDgQ7F+VMcZPQVER2woLaR4by8A67kQPhbBtucuvfx2S99Vf/rLO7SUlJVUjN6anp58x7ktaWhrTp0+nRYsWPPTQQ4B3gLFx48Zx6623AnDttdcyZ84cMjIyWLlyJffeey8ffvghDzzwAPfccw933HEHs2fPrvHYGzduZPDgwQ06r3bt2lWNWzNv3jx27dpFeno68+fPZ9KkSZSXl/PjH/+YN954g5SUFObPn8/MmTN55plnGnQ8Y0z9PvP1t1+amkp8I3+iDttwd0pN3TKBqmuI3c8++4wFC7zPMLn99tt5+OGHz79YP5MmTaqanzhxIvPnz2fGjBnMnz+f+fPns23bNjZu3MjIkSMBb/dNXWPaGGPO3zLf0LXDQnVJTB3CNtzra2GHo/qG2PUf5rcmffv2rXoiUnWxsbF4PJ6q5bqGCZ40aRITJkzg5ptvRkTIyMhgw4YN9O3bl88//zzQ0zHGnKdPfS33Kx0Id+tzP0fVh8T1X65riN0rrrii6kvSnFoGlLjmmmsoLS1l7ty5VevWr1/PsmXL6N69O5s3b6a0tJRjx47xwQcf1Fpjz549cblcPPbYY1Ut+t69e1NQUFAV7uXl5WzatKmhvwZjTD2Kysr48sABYkQabSRIfxbu5+j666/ntddeIysri2XLljF58mSeeOIJBg4cyI4dO2odYvfPf/4zs2fPpn///uzfX/ODrESE1157jffff5+ePXvSt29fHnnkETp27EjXrl2ZOHEi/fr1Y+LEiQwcOLDOOidNmsSLL77IxIkTAYiPj+eVV17h4YcfJjMzk6ysrDMu6TTGBNfK/ftxezxkdexIS78H4zQWG/LXnBP7GxkTmP/6+GN+uXQp9w8Zwp9Hjw7a+9qQv8YY46DKm5ec6G8HC3djjAk6t8fD577BwizcfexBzeHL/jbGBGbdwYOcKiujZ5s2dEpOdqSGgMJdREaJyDYRyRORGbXsc7WIrBWRTSLycUOKSUhIoLCw0EIkDKkqhYWFJCQkOF2KMWHvE9/17U612iGA69xFxAXMBkYC+cAqEVmoqpv99mkN/BUYpap7ReSChhSTmppKfn4+BQUFDflxE2IJCQmkOnBJlzGR5qPduwG4Oi3NsRoCuYlpCJCnqjsBRGQecCPgPxTibcCrqroXQFW/aUgxcXFxpKenN+RHjTEmLFR4PFUt9xEOhnsg3TJdgH1+y/m+df56AW1EZKmIrBGRO4JVoDHGRJKvDh7keGkp6a1b0711a8fqCNbwA7HAYOBaoDnwuYisUNUzBiIXkWnANIBuDvZFGWNMqHy0axfgbKsdAmu57we6+i2n+tb5ywcWq2qRqh4GPgEyq7+Rqs5V1WxVzU5JSWlozcYYE7Yq+9tHONzFHEi4rwIyRCRdROKBycDCavu8AVwpIrEikghcCmwJbqnGGBPeyisqqh7O4XTLvd5uGVV1i8h9wGLABTyjqptEZLpv+xxV3SIi7wLrAQ/wd1XdGMrCjTEm3Kw5cIBTZWVktG1Ll5YtHa0loD53VV0ELKq2bk615SeAJ4JXmjHGRJZw6W+HMLxD1RhjItWHYdLfDhbuxhgTFKVuN5/5+tudvHmpkoW7McYEwRf791PidnNx+/Z0bNHC6XIs3I0xJhg+DKP+drBwN8aYoFiycycA1/Xs6XAlXhbuxhhzno6dPs3K/HxcImHxZSpYuBtjzHn7aNcuKlS5vGtXR56XWhMLd2OMOU9LduwA4LoePRyu5FsW7sYYc57Crb8dLNyNMea87DhyhJ1Hj9ImIYHszp2dLqeKhbsxxpyHyi6Z7/TogSsmfCI1fCoxxpgIFI5dMmDhbowxDVZeUVF189LIMPoyFSzcjTGmwb7Yv58TpaX0btfO0Ufq1cTC3RhjGmixr7893FrtYOFujDEN9vb27QCMychwuJKzWbgbY0wDfH3yJF8eOEDz2NiwGOK3Ogt3Y4xpgEW+Vvt3evSgeVycw9WczcLdGGMaoLJLZmwYdslAgOEuIqNEZJuI5InIjBq2Xy0ix0Vkre/1n8Ev1RhjwkOp2817vi9Tx/bq5XA1Nav3Adki4gJmAyOBfGCViCxU1c3Vdl2mquNCUKMxxoSVj/fsoai8nMwOHUht2dLpcmoUSMt9CJCnqjtVtQyYB9wY2rKMMSZ8vZWbC8C4MG21Q2Dh3gXY57ec71tX3VARWS8i74hI36BUZ4wxYUZVq8I9XPvbIYBumQB9CXRT1VMiMgZ4HTjrrEVkGjANoFu3bkE6tDHGNJ6thw+z69gx2icmMqRLTe3c8BBIy30/0NVvOdW3roqqnlDVU775RUCciLSv/kaqOldVs1U1OyUl5TzKNsYYZ1S22kdfeGFYjQJZXSCVrQIyRCRdROKBycBC/x1EpKOIiG9+iO99C4NdrDHGOO2NbduA8O5vhwC6ZVTVLSL3AYsBF/CMqm4Skem+7XOAW4F7RMQNlACTVVVDWLcxxjS6g6dOsXzfPpq5XIy+8EKny6lTQH3uvq6WRdXWzfGbfxp4OrilGWNMeHlj61YUGNmzJ8lh8iDs2oRvh5ExxoSZ17ZuBeCmiy5yuJL6WbgbY0wAjp0+zYe7dhEjwvVh3t8OFu7GGBOQt3NzKfd4uKp7d1KSkpwup14W7sYYE4BI6pIBC3djjKlXSXk57+TlARbuxhgTNZbs2EFxeTnZnTvTtVUrp8sJiIW7McbU41Vfl8zNEdJqBwt3Y4yp02m3m9d94X5Lnz4OVxM4C3djjKnD4rw8TpSWMrBjR3q1a+d0OQGzcDfGmDrM27QJgMn9+jlcybmxcDfGmFoUlZWx0DdQ2MS+kfWYCgt3Y4ypxdvbt1NcXs5lqamktW7tdDnnxMLdGNOk5ORAWhrExHinOTm17ztv40YAJkdYqx2C9yQmY4wJezk5MG0aFBd7l/fs8S4DTJ165r4nSktZtH07AkyIwHC3lrsxpsmYOfPbYK9UXOxdX90bW7dSWlHB8LQ0OicnN06BQWThboxpMvbuDXz9S74umUkR2GoHC3djTBPSrVtg6w+eOsXiHTuIi4nh1gi6ccmfhbsxpsmYNQsSE89cl5joXe8vZ/16PKqM7dWL9tV/IEJYuBtjmoypU2HuXOjeHUS807lzz/4y9fn16wG4MzPTgSqDI6BwF5FRIrJNRPJEZEYd+10iIm4RuTV4JRpjTPBMnQq7d4PH451WD/a1Bw+y/tAh2jVvzpiMDCdKDIp6w11EXMBsYDTQB5giImd1Qvn2+z2wJNhFGmNMY3l+3ToApvTrR7zL5XA1DRdIy30IkKeqO1W1DJgH3FjDfj8GFgDfBLE+Y4xpNOUVFeRs2ADAHRHcJQOBhXsXYJ/fcr5vXRUR6QLcBPwteKUZY0zjWrJjB98UFXFR+/Zkd+7sdDnnJVhfqP438LCqeuraSUSmichqEVldUFAQpEMbY0xwPOfrkrkzMxMRcbia8xPI8AP7ga5+y6m+df6ygXm+X0Z7YIyIuFX1df+dVHUuMBcgOztbG1q0McYEW0FREa9v3UqMCN8bMMDpcs5bIOG+CsgQkXS8oT4ZuM1/B1VNr5wXkWeBt6oHuzHGhLPn1q2j3ONhXK9epLZs6XQ5563ecFdVt4jcBywGXMAzqrpJRKb7ts8JcY3GGBNSqsrcNWsAmDZokMPVBEdAo0Kq6iJgUbV1NYa6qt51/mUZY0zj+XjPHrYfOUKX5GRGR/C17f7sDlVjTJNX2Wr//sCBxMZERyxGx1kYY0wDHS4uZsGWLQjw/SjpkgELd2NME/f8unWUVVQwOiODbq1aOV1O0Fi4G2OaLE8UfpFaycLdGNNkvb9zJ9sKC0lt2ZKxvXo5XU5QWbgbY5qs/1m5EoB7srOj5ovUStF1NsYYE6C8I0dYtH07zVwufhhlXTJg4W6MaaJmf/EFCtzWvz8pSUlOlxN0Fu7GmCbnZGkpz6xdC8CPhwxxuJrQsHA3xjQ5z69bx4nSUq7s1o2BnTo5XU5IWLgbY5oUjyp/+eILAO6P0lY7WLgbY5qYt3Jz2VZYSNeWLRl/0UVOlxMyFu7GmCbl9599BsBPLr+cuAh+Rmp9LNyNMU3GZ3v3snzfPtokJPCDKLz80Z+FuzGmyfjD8uUA/OiSS2gRH+9wNaFl4W6MaRI2FxSwcNs2EmJj+fGllzpdTshZuBtjmoQnfa32f8vK4oIovGmpOgt3Y0zU23v8OC+uX0+MCD+9/HKny2kUFu7GmKj3u2XLKPd4mNS3Lz3btnW6nEYRULiLyCgR2SYieSIyo4btN4rIehFZKyKrReTK4JdqjDHnbt/x4/zjq68Q4NGrrnK6nEZT7wOyRcQFzAZGAvnAKhFZqKqb/Xb7AFioqioiA4B/AdF7d4AxJmL87tNPKfd4mNKvHxenpDhdTqMJpOU+BMhT1Z2qWgbMA27030FVT6mq+haTAMUYYxzWVFvtEFi4dwH2+S3n+9adQURuEpGtwNvA3cEpzxhjGu7xTz+lrKKCSU2s1Q5B/EJVVV9T1YuA8cBjNe0jItN8ffKrCwoKgnVoY4w5y97jx/l7E221Q2Dhvh/o6rec6ltXI1X9BOghIu1r2DZXVbNVNTulif0raoxpXP/50UeUVVQwuV8/+jTBvAkk3FcBGSKSLiLxwGRgof8OInKhiIhvfhDQDCgMdrHGGBOIDYcO8fy6dcTGxPDYiBFOl+OIeq+WUVW3iNwHLAZcwDOquklEpvu2zwFuAe4QkXKgBJjk9wWrMcY0qv/48EMUmD54cJO5rr26esMdQFUXAYuqrZvjN/974PfBLc0YY87dsj17eCs3l6S4OH7RBPvaK9kdqsaYqKGqPPz++wA8NHQoHVq0cLgi51i4G2OixoItW/g8P5+UxMQmM4ZMbSzcjTFRoaS8nIeWLAHg11dfTXKzZg5X5CwLd2NMVPjj55+z5/hxBnTowLTBg50ux3EW7saYiJd/4gS/+/RTAP48ahSuGIs2+w0YYyLejPffp7i8nFv79OHqtDSnywkLFu7GmIj26d695GzYQDOXiydGjnS6nLBh4W6MiVhlFRX8+1tvAfDwFVeQ1rq1wxWFDwt3Y0zEenL5cjYXFJDRti2PDBvmdDlhxcLdGBORdhw5wmOffALAnHHjSIgN6Ib7JsPC3RgTcVSVexct4rTbze0DBnBNerrTJYUdC3djTMR5cf16luzYQdvmzfnjddc5XU5YsnA3xkSU/SdOcP+77wLwx+uuIyUpyeGKwpOFuzEmYqgq0956i2OnTzOuVy/uzMx0uqSwZeFujIkYz65dy6Lt22mTkMD/jhuH7xlBpgYW7saYiLD3+HEeXLwYgL+MHk3n5GSHKwpvFu7GmLBX4fFw+2uvcaK0lPEXXcRt/fs7XVLYs3A3xoS93y5bxid79tCpRQvmWndMQCzcjTFh7bO9e/nVxx8jwAs33WRXxwQooHAXkVEisk1E8kRkRg3bp4rIehHZICLLRcS+wjbGnLejJSXc9uqreFT5+RVXcG2PHk6XFDHqDXcRcQGzgdFAH2CKiPSpttsuYLiq9gceA+YGu1BjTNPiUeXO119n7/HjDOnShcdGjHC6pIgSSMt9CJCnqjtVtQyYB9zov4OqLlfVo77FFUBqcMs0xjQ1j3/6KW/m5tImIYGXbrmFOJfL6ZIiSiDh3gXY57ec71tXm+8D75xPUcaYpm3Jjh384sMPESDn5pvp0aaN0yVFnKAOoyYiI/CG+5W1bJ8GTAPo1q1bMA9tjIkSe44d47YFC1DgV8OHMzojw+mSIlIgLff9QFe/5VTfujOIyADg78CNqlpY0xup6lxVzVbV7JSUlIbUa4yJYidLS7lh3jwKS0oYk5HBo8OHO11SxAok3FcBGSKSLiLxwGRgof8OItINeBW4XVVzg1+mMSba5ORAWhrExHinL7zoYeqrr7L+0CF6tWvHizfdRIxdz95g9XbLqKpbRO4DFgMu4BlV3SQi033b5wD/CbQD/uq7ucCtqtmhK9sYE8lycmDaNCgu9i7v2QN3v/Q+7iHeL1DfmjKFNs2bO1tkhBNVdeTA2dnZunr1akeObYxxVlqaN9CrDF4D178Fnhg++rfbuTotzaHKwp+IrAmk8Wx3qBpjGt3evX4LF22BsW97598aa8EeJBbuxphGV3WxXPc9cOsCiFH4aDjdCwc5Wlc0sXA3xjS6WbMgofshmPISxFbAqsE0XzWcWbOcrix62OPCjTGNbuB1BcTvep7TFaWw+WK6bRrDb+cKU6c6XVn0sHA3xjSqbYcPc+3zz3OiopiRPXqwcObNJMRaJ0Kw2W/UGNNo8o4c4Zrnn+fgqVNck57OG5MnkxBrbcxQsN+qMaZRbCko4DsvvMDXJ09yVffuLJw8meZxcU6XFbUs3I0xIffVgQNc9+KLHC4u5qru3XlryhSS4uOdLiuqWbgbY0Jq+b59jMnJ4XhpKaMuvJAFEyeSaC32kLNwN8aEzOtbtzJlwQJOu93ccvHF/POWW4i3cdkbhYW7MSYk/rJyJQ+8+y4KfH/gQOaMG0dsjF3D0Vgs3I0xQVXh8fDz997jTytWAPDYiBHMHDYMsREeG5WFuzEmaCofaP1uXh6xMTH844YbuCMz0+mymiQLd2NMUGw9fJgbXnqJ7UeO0D4xkZcnTLBBwBxk4W6MOW//2rSJHyxcyMmyMjI7dOD1yZNJa93a6bKaNAt3Y0yDlbrd/GTxYv7qezbDxL59eeaGG+wa9jBg4W6MaZCthw8z9dVX+fLAAeJdLv503XXce8kl9sVpmLBwN8acE1Xlb6tX89CSJZS43aS3bs2/Jkwgu3Nnp0szfizcjTEByz9xgmlvvsk7eXkA3JGZyf+MGkWrhASHKzPVBXRHgYiMEpFtIpInIjNq2H6RiHwuIqUi8lDwyzTGOMmjyv+uXk2f2bN5Jy+PNgkJ/OvWW3lu/HgL9jBVb8tdRFzAbGAkkA+sEpGFqrrZb7cjwP3A+JBUaYxxzJaCAu55+20+9j3R+sbevfnr2LF0Tk52uDJTl0C6ZYYAeaq6E0BE5gE3AlXhrqrfAN+IyNiQVGmMaXSnysp47OOP+dOKFbg9Hi5ISuLp0aO5tU8f+9I0AgQS7l2AfX7L+cClDTmYiEwDpgF0q3pCrjEmnHhUyVm/nkc++ID9J08iwLRBg/jttdfSLjHR6fJMgBr1C1VVnQvMBcjOztbGPLYxpn6f7NnDTxYvZs2BAwAM7tSJv44dy5AuXRyuzJyrQMJ9P9DVbznVt84YEyW+PHCAX3z4YdVVMJ2Tk5l1zTXcPmAALhvJMSIFEu6rgAwRSccb6pOB20JalTGmUaw7eJDHPvmEBVu2ANAiPp6fDR3KTy+/3O4yjXD1hruqukXkPmAx4AKeUdVNIjLdt32OiHQEVgMtAY+IPAj0UdUTIazdGNNAK/LzmbVsGW/l5npXlMfScuslPH7DFdwzPMnZ4kxQBNTnrqqLgEXV1s3xmz+It7vGGBOmKjwe3szN5akVK/jEd1kj5bGwZjB8dgUnTibz0DvQ0gVTpzpbqzl/doeqMVHuaEkJz61bx9NffMGOo0cBaNmsGXxxCSeWXAZF37bUi4th5kwL92hg4W5MFFJVVn39NXPXrOGfGzZQ4nYDkNa6NQ9eeil3DxxIq+bNoIZr1vbubeRiTUhYuBsTRQ6dOkXOhg0889VXbCooqFr/nR49uCc7mxt69656jmm3blDZO+PPbkGJDhbuxkS4k6WlvLZ1KzkbNvD+zp141NscT0lM5PYBA/j37Gx6tWt31s/NmgXTpnm7YiolJnrXm8hn4W5MGMjJ8fZ1793rbTnPmlV3v/fx06d5MzeXlzdvZnFeHqUVFQDExsQwNiODuwcOZExGBvFwoecGAAAKYklEQVQuV63vUfn+53JcEzks3I1xWE7OmS3oPXu8y3Bm0O46epS3cnNZmJvL0t27cXs8AAhwZbduTO3fnwl9+pzTEAFTp1qYRysLd2McNnPmmV0j4F1+5JdltLt0L+/m5fFOXh65hYVV22NEGN69OxP69OGmiy+2ERrNWSzcjXFY1dUpceXQJR/S9kD6Lval5jM6x1O1X8tmzbiuZ09u6NWLMRkZNoiXqZOFuzEOOXjqFJ/v20fyzfs40XofdP4aXN+GOR7hks6d+U6PHoy+8EIuS00lro4+dGP8Wbgb0wgKior46uBBvjxwgFVff80X+/eTf8I3Okd/304KHOgIe7oT/3Uaf/lpGtPusKccmYaxcDcmiErdbnILC9lUUMD6Q4dYf+gQ6w4d+jbI/STHxzOkSxeGdu3K6e2pzHuyK/l5CXbVigkKC3djzpGq8k1REduPHGF7YSFbDx9mm2+ad+QIFXr2bZ9JcXFkdezIwI4duaRLFy7p3Jne7dsTU/lEoxHwh2mNfCImqlm4G+On8nrzPQdL6XTRcW679xjpWcfYdewYO48eZdexY+w4coSTZWU1/nyMCBlt29InJYUBHTowoEMH+l9wARe2bWvjoptGZeFumpRSt5tDRUUcPHWKAydPcvDUKfafPMn+EydYlXuSjXtOoFNOQEIpB4A/HgAOnP0+rZo1I6NdOzLatqV3u3Zc1L49vdu3p3e7djSPi2vs0zLmLBbuJmK5PR6OnT7N0ZISjpSUUFhSQmFxMYUlJRwuLuZwcTEFxcV8U1REQVERh4qKOHb6dN1vmuKblsfC8VZwrDUtKlrxyD2t6dmmDelt2tCjTRvaNW9uD4k2Yc3C3TQqt8dDcXk5RWVlFJWXc6qsjKKyMk6WlXGqrIyTpaV8tLyMNxeXcux0KS3aljHgklLadi7l+OnTHPN71dY1UheXCBckJdGxRQs6JSfTqUULOicn0yU5melTW8KJlnC8JZQ0x3vvJxQJ/Mc/g/yLMCbELNyjlEcVt8dT9SqvqKC8hqnb46GsooLyigrKfK9y37qyigpK3W7v1DdffXra7eZ0RYV36nZTUl7unfrmi8vLKXG7KfbNl/nGQKnXYO/kFLD8FJB79i4CtE5IoE3z5rRt3px2vmn7xETaJybSrnlzLkhKIiUpiZTERDq0aEHb5s2//RKzmt+Vwp6DZ6+3URJNJIq4cH83L4/VX38NeK9aqLwuoab5hk4r38Pjt76meU8N6/1f/vvU9KrweL6d91uu8NteUc/U7ZuvDPHKdTUM0x0WBEiKjycpLo6k+Hha+F5JcXEkN2vG4oXxFB2Nh9JmUFo5TSClZTNefqEZrRMSaJ2QQKuEBFo2a1ZrUDeEjZJooknEhfvCbdv42+rVTpcREWJjYoiLiSHO5aqaj/Ut+8/H+15xMTFV8we/drF2tYuKMhdUuMDtIlZiGTXSxeDMWJq5XDSLjaV5bCzNYr3LzePiaB4bS0JsLM3j4kj0LVfOJ8bF0czlqrOvOmYyNT5A4rDA8LSQ/aoAGyXRRJeAwl1ERgF/xvuA7L+r6uPVtotv+xigGLhLVb8Mcq3k5MArcy6EZgm0bCmMuBr696MqLISz5xs6Be9lbZXLq1cLb74BR44IbdsIN40XLr/Mt49vP1dMTNXPxIhUbYvxe7n81rtEcMXEVP1s5bbK+XffFX7zXzGcLhbQGPAICc1iePy3wi03e8O58j1iY75djnO5zrtFm5YGFdUe5OAGNmyEN3ef11vXyekHSNgoiSZqqGqdL7yBvgPoAcQD64A+1fYZA7yDN1MvA1bW976DBw/Wc/Hii6qJiarw7Ssx0bs+1Jw6dvfuZx6z8tW9e2iPq6oqUvOxRUJ7XCf/zsZEAmC11pOvqkogd1UMAfJUdaeqlgHzgBur7XMj8Lzv2CuA1iLS6bz/5fFT27CoM2cG8yjhdezanmXZGM+4rK2lHOoW9NSpMHcudO8OIt7p3LnWmjbmXAUS7l2AfX7L+b5157oPIjJNRFaLyOoCv+c7BsLJoHPq2E4FLHj7mquPKNtYXy5OnQq7d4PH451asBtz7hr1fmhVnauq2aqanZKSUv8P+HEy6Jw6ttMBay1oYyJXIOG+H+jqt5zqW3eu+5wXJ4POqWM7HbDWgjYmcgUS7quADBFJF5F4YDKwsNo+C4E7xOsy4Liq1jAiR8M5GXROH9sC1hhzrkRrGJ70rJ1ExgD/jffKmWdUdZaITAdQ1Tm+SyGfBkbhvRTy31S1zovRs7OzdbVdr26MMedERNaoanZ9+wV0nbuqLgIWVVs3x29egR+da5HGGGNCwwaYNsaYKGThbowxUcjC3RhjopCFuzHGRKGArpYJyYFFCoAahogKSHvgcBDLiQR2zk2DnXPTcD7n3F1V670L1LFwPx8isjqQS4GiiZ1z02Dn3DQ0xjlbt4wxxkQhC3djjIlCkRruc50uwAF2zk2DnXPTEPJzjsg+d2OMMXWL1Ja7McaYOoR9uIvIMyLyjYhs9FvXVkTeE5HtvmkbJ2sMtlrO+QkR2Soi60XkNRFp7WSNwVbTOftt+6mIqIi0d6K2UKntnEXkx76/9SYR+YNT9YVCLf9tZ4nIChFZ63uYzxAnawwmEekqIh+JyGbf3/MB3/qQZ1jYhzvwLN7RJv3NAD5Q1QzgA99yNHmWs8/5PaCfqg4AcoFHGruoEHuWs88ZEekKXAc0wjO3Gt2zVDtnERmB97GVmaraF3jSgbpC6VnO/jv/Afi1qmYB/+lbjhZu4Keq2gfv86V/JCJ9aIQMC/twV9VPgCPVVt8IPOebfw4Y36hFhVhN56yqS1TV7VtcgfeBKFGjlr8zwFPAz4Go+3KolnO+B3hcVUt9+3zT6IWFUC3nrEBL33wr4OtGLSqEVPWAqn7pmz8JbMH7CNKQZ1jYh3stOvg9DOQg0MHJYhxwN/CO00WEmojcCOxX1XVO19KIegHDRGSliHwsIpc4XVAjeBB4QkT24f2kEm2fSgEQkTRgILCSRsiwSA33Kr6x5KOuVVcbEZmJ96NejtO1hJKIJAL/gfdjelMSC7TF+xH+Z8C/fA/DiWb3AP9PVbsC/w/4h8P1BJ2ItAAWAA+q6gn/baHKsEgN90Mi0gnAN42qj661EZG7gHHAVI3+a1h7AunAOhHZjbcb6ksR6ehoVaGXD7yqXl8AHrzjkESzO4FXffMvA1HzhSqAiMThDfYcVa08z5BnWKSG+0K8/0Hgm77hYC2NQkRG4e17vkFVi52uJ9RUdYOqXqCqaaqahjf0BqnqQYdLC7XXgREAItILiCf6B9X6Ghjum78G2O5gLUHl+9T1D2CLqv7Jb1PoM0xVw/oFvAQcAMrx/g/+faAd3m+YtwPvA22drrMRzjkP2Aes9b3mOF1nqM+52vbdQHun62yEv3M88CKwEfgSuMbpOhvhnK8E1gDr8PZHD3a6ziCe75V4u1zW+/2/O6YxMszuUDXGmCgUqd0yxhhj6mDhbowxUcjC3RhjopCFuzHGRCELd2OMiUIW7sYYE4Us3I0xJgpZuBtjTBT6/14LAWjq4XSRAAAAAElFTkSuQmCC" + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = pnts_n\n", + "y = pnts_t\n", + "plt.figure()\n", + "plt.plot(x, y, 'bo', label=\"Measurement\")\n", + "xx = np.linspace(pnts_n[0],pnts_n[-1],100)\n", + "plt.plot(xx, func(xx, *popt), 'r-',\n", + " linewidth=2, color='teal', label=\"Fitted Curve\")\n", + "plt.legend(loc='best')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Discussion" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* Exhaustive search exhibits exponential running time:\n", + " linear in $\\ell$ but exponential in $n$.\n", + "\n", + "* So it is useful when $n$ is small even if $\\ell$ is large." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Approximation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Greedy" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Find the variable that appears most often and assign it accordingly to\n", + "maximize the number of satisfied clauses.\n", + "\n", + "\n", + "1. $L \\gets \\emptyset$\n", + "2. **for** $w\\in\\{x_1,\\bar x_1,\\ldots,x_n,\\bar x_n\\}$ **do**\n", + "3. $\\quad$ Count occurrences of $w$ in $\\phi$\n", + "4. $\\quad$ Append pair ($w$, count of occurrences of $w$ in $\\phi$) to $L$\n", + "5. **end for**\n", + "6. Sort $L$ with respect to the second component\n", + " Comment: We now have the literals in decreasing order of total occurrence\n", + "7. **for all** $(w,c)\\in L$ **do**\n", + "8. $\\quad$ Set $w$ to True $\\quad$ (i.e. if $w=\\bar x_i$ then set $x_i$ to False)\n", + "9. **end for**\n", + "10. **return** count of satisfied clauses" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Cost\n", + "Steps 1-5 and 7-10 are $O(n)$, and step 6 costs $O(n\\log n)$ assuming the use of an $O(n\\log n)$ sorting\n", + "algorithm.\n", + "\n", + "The total costs is the sum of these three sequential blocks, yielding a total cost of: $$(n \\log n)$$" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def greedy(self):\n", + " '''\n", + " Solve in the 3SAT problem instance using a simple greedy search.\n", + " Find the variable that appears most often and assign it accordingly\n", + " '''\n", + " literals_occurance = [[i,0] for i in range(-self.nv,self.nv+1)]\n", + " for c in self.clauses:\n", + " for v in c:\n", + " literals_occurance[ v+self.nv ][1] += 1\n", + " literals_occurance.sort( key=lambda a:a[1] )\n", + " for v in literals_occurance:\n", + " self.set_variable( v[0] , v[0]>0 )\n", + " return self.ratio_unsatisfied()\n", + "\n", + "SAT.greedy = greedy" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Greedy Randomized Adaptive Search Procedure (GRASP)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The general template pseudocode for GRASP is given by:\n", + "\n", + "1. $best\\_candidate\\gets \\emptyset$\n", + "2. **while** \"termination condition\" is not met **do**\n", + "3. $\\quad$ $greedy\\_candidate\\gets \\text{ConstructGreedyRandomizedSolution}()$\n", + "4. $\\quad$ $grasp\\_candidate\\gets \\text{LocalSearch}(greedy\\_candidate)$\n", + "5. $\\quad$ **if** $f(grasp\\_candidate) < f(best\\_candidate)$ **then**\n", + "6. $\\qquad$ $best\\_candidate\\gets grasp\\_candidate$\n", + "7. $\\quad$ **end if**\n", + "8. **end while**\n", + "9. **return** $best\\_candidate$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For our problem, we have:\n", + "\n", + "- **termination condition** is simply to repeat a fixed number of times, e.g. 100 times.\n", + "\n", + "- $f$ gives the ratio of unsatisfied clauses to $\\ell$. Objective is to minimize it.\n", + "\n", + "- **ConstructGreedyRandomizedSolution()** works like Greedy but shuffles $L$ in blocks of a given size. \\[TODO: EXPLAIN MORE\\]\n", + "\n", + "- **LocalSearch()** works by flipping one variable's assignment at a time, and seeing if $f$ improves. The best improvement is selected after trying all of them. [@greedy]" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def greedy_randomized(self, block_size=2):\n", + " '''\n", + " Solve in the 3SAT problem instance using a randomized greedy search\n", + " Find the variable that appears most often and assign it accordingly\n", + " Same as greedy() but randomizes the order of 'literals_occurance' a bit\n", + " '''\n", + " literals_occurance = [[i,0] for i in range(-self.nv,self.nv+1)]\n", + " for c in self.clauses:\n", + " for v in c:\n", + " literals_occurance[ v+self.nv ][1] += 1\n", + " literals_occurance.sort( key=lambda a:a[1] )\n", + " # randomize list block by block\n", + " for i in range(0,len(literals_occurance), block_size):\n", + " tmp = literals_occurance[i:i+block_size]\n", + " shuffle( tmp )\n", + " literals_occurance[i:i+block_size] = tmp\n", + " for v in literals_occurance:\n", + " self.set_variable( v[0] , v[0]>0 )\n", + " return self.ratio_unsatisfied()\n", + "\n", + "def GRASP(self, max_repetition=100):\n", + " ''' Solve the 3SAT problem instance using the GRASP meta-heuristic '''\n", + " best = 1 # i.e. start with the worst solution that satisfies no clauses\n", + " for i in range(max_repetition):\n", + " # random block_size from {2,3,4}\n", + " candidate = self.greedy_randomized( choice([2,3,4]) )\n", + " # local search: flip variable assignment and see if things improve\n", + " # Choose the best one from the flipped attempts\n", + " best_v = 0\n", + " for v in self.variables:\n", + " self.flip_variable(v) # try changing variable v\n", + " neighbour = self.ratio_unsatisfied()\n", + " if neighbour < candidate:\n", + " candidate = neighbour\n", + " best_v = v\n", + " self.flip_variable(v) # undo, to try flipping next variable\n", + " if candidate < best:\n", + " best = candidate\n", + " return best\n", + "\n", + "SAT.greedy_randomized = greedy_randomized\n", + "SAT.GRASP = GRASP" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Genetic Algorithm" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The general template pseudocode for Genetic Algorithms is given by:\n", + "\n", + "1. determine initial population $p$\n", + "2. **while** termination criterion not satisfied **do**\n", + "3. $\\quad$ generate set $pr$ of new candidate by recombination\n", + "4. $\\quad$ generate set $pm$ of new candidates from $p$ and $pr$ by mutation\n", + "5. $\\quad$ select new population $p$ from candidates in $p, pr, pm$\n", + "6. **end while**\n", + "7. **return** fittest candidate found" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For our problem, we have:\n", + "\n", + "- **Initial population**: generate 20 candidates by randomly assigning True/False values to the variables.\n", + "- **termination criterion** is to simply repeat 20 times.\n", + "- **Recombination** is done by cutting a pair of candidates at some point ($ab$,$cd$) then creating the possible combinations $ad$ and $cb$.\n", + "- Each candidate is **mutated** by flipping one random variable with probability 20%.\n", + "- The candidates are then sorted according to the objective function, and the fittest are kept." + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def GA(self, pop_size=20, max_generations=20):\n", + " '''\n", + " Solve the 3SAT problem instance using the Genetic Algorithm meta-heuristic\n", + " '''\n", + " def sort_key(p):\n", + " self.configuration = p\n", + " return self.ratio_unsatisfied()\n", + " population = [\n", + " [choice([True,False]) for nv in range(self.nv)]\n", + " for i in range(pop_size)\n", + " ]\n", + " for generation in range(max_generations):\n", + " # crossing\n", + " for i in range( pop_size ):\n", + " for j in range(i+1, pop_size):\n", + " cut = randint(0,self.nv-1)\n", + " population.append( population[i][:cut] + population[j][cut:] )\n", + " population.append( population[j][:cut] + population[i][cut:] )\n", + " # mutation\n", + " for i in range( len(population) ):\n", + " j = randint(0,self.nv-1) # variable to mutate\n", + " if randint(0,100)<=20: # probabilistically mutate\n", + " population[i][j] = not population[i][j]\n", + " # choose fittest\n", + " population.sort(key=sort_key)\n", + " population = population[:pop_size]\n", + " self.configuration = population[0]\n", + " return self.ratio_unsatisfied()\n", + "\n", + "SAT.GA = GA" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Testing" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "rho\tGreedy\tGRASP\tGA\n", + "1\t0.9835\t0.9995\t1.0000\n", + "2\t0.9653\t0.9932\t1.0000\n", + "3\t0.9505\t0.9788\t0.9995\n", + "4\t0.9451\t0.9736\t0.9972\n", + "5\t0.9416\t0.9678\t0.9943\n", + "6\t0.9443\t0.9665\t0.9938\n", + "7\t0.9364\t0.9609\t0.9936\n", + "8\t0.9337\t0.9575\t0.9921\n", + "9\t0.9366\t0.9589\t0.9932\n" + ] + } + ], + "source": [ + "instance = SAT() # global variable... [TODO?]\n", + "pnts_rho = []\n", + "pnts_greedy = []\n", + "pnts_grasp = []\n", + "pnts_ga = []\n", + "\n", + "def test_Approximation():\n", + " print( \"rho\\tGreedy\\tGRASP\\tGA\" )\n", + " max_repeats = 100\n", + " nv = 20\n", + " for rho in arange(1,10,1): # rho = nc/nv\n", + " nc = int(rho*nv)\n", + " q = [0]*3 # to measure 'quality'\n", + " for repeat in range(max_repeats):\n", + " instance.random_yes_instance(nv, nc)\n", + " q[0] += instance.greedy()\n", + " q[1] += instance.GRASP()\n", + " q[2] += instance.GA()\n", + " for i in range(3): q[i] = '{:1.4f}'.format(1-q[i]/max_repeats)\n", + " print( str(rho)+'\\t'+\"\\t\".join(q) )\n", + " pnts_rho.append( rho )\n", + " pnts_greedy.append( q[0] )\n", + " pnts_grasp.append( q[1] )\n", + " pnts_ga.append( q[2] )\n", + "\n", + "test_Approximation()" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4lNX1wPHvSQiQBAiLAYGsLih7kMiq/hDKorJYNzZFUQRapEitilirtuJaFW1RxA2LKagouAsFBRUkGHYREMQkBBAiKFuAQLi/P+5kn5AJmeSd5XyeZx6Zd96ZOYlw7vve5VwxxqCUUip4hDgdgFJKqeqliV8ppYKMJn6llAoymviVUirIaOJXSqkgo4lfKaWCjCZ+pZQKMpr4lVIqyGjiV0qpIFPD6QDcOeuss0xCQoLTYSillN9YtWrVL8aYaE/O9cnEn5CQQFpamtNhKKWU3xCRDE/P1a4epZQKMpr4lVIqyGjiV0qpIKOJXymlgowmfqWUCjLlJn4ReU1E9orId2W8LiLyvIhsE5H1InJRkdf6icgW12uTvBl4SSkbUkiYmkDIwyEkTE0gZUNKVX6dx3w1LqVU8PLkin8m0O80r18BnO96jAZeBBCRUGCa6/VWwFARaVWZYMuSsiGF0R+OJuNABgZDxoEMRn842vEk66txKaWCW7nz+I0xX4pIwmlOGQT8x9g9HFeISH0RaQokANuMMdsBRGSO69zvKxt0Sfcvvp+cEznFjuWcyGHcx+PY8ssWb3+dx55Pfd5tXBM+nUCt0FpEhEWc9lEjpOqWWaRsSOH+xfeTeSCTuKg4pvSawvC2w6vs+5RSvsMbmaU5sKPI8yzXMXfHO5f1ISIyGnvHQFxcXIUCyDzgft3CgeMHeOTLRyr0Wd5kcL+f8b6j+7j+nevLfX9YSBgRYRFE1ow8fSNR4/QNSMnHgh8XMGnRJI6ePApQcCcCaPJXKgj4zMpdY8wMYAZAcnJyhXaAjzscSkadvFLH4w+Hkv7USe8EeAYSpiaQ4aZRala3GQtuXEDOiRxyTuRwJPdIwZ/LfJws/vy3Y7+5Pe9M5d+JXNjoQlpFtyI8LLwyP7pSyod5I/HvBGKLPI9xHQsr47jXTVmQx+gBkFOz8FhErj3OkwZEquJry4+r1xRGfzi6WEKOCIvgyd5P0qZxG69/nzGGYyePlduIDHtvmNv37zu6j+SXkwmREC5odAHtmrQr9oitF4s49LtUSnmPNxL/B8Adrj78zsABY8xuEckGzheRRGzCHwK4zziVNPxgPHyYwf29IDMK4g7AlMUwfAPQrRv87W/Qr1+1NwD53SbV1ZcuIoSHhRMeFk4jGpV53n2L7yvzTuS5fs+xfs961u9Zz8qdK3lr41sFr0fViqJdk3a0b9K+oDFo07gNkTUjq+TnUUpVDbFjsqc5QWQ20AM4C9gDPIi9mscYM13sJeC/sTN/coCRxpg013uvBKYCocBrxpgpngSVnJxsKlSkLSUFRo+GnCJdHRERMGQILF4MGRmQnGwbgP79HbsD8BX5s41K3onMGDCjVKN08PhBvtv7Het+XmcbhL22UTicexgAQTi34bm2IWjcjvZn20YhoX4CIaLLRJSqLiKyyhiT7NG55SV+J1Q48YNN/vffD5mZEBcHU6bA8OGQmwuzZtnnP/0EHTrYBmDQoKBuACozq+eUOUXGbxkFdwbr9thGYdv+bQUD2nVq1qFt47bFuoraNm5LVO2oqvyxlApawZn4y3PihG0cpkyBbdugfXt44AH4/e8hRK9MveFI7hE2Zm8saBDyG4Xfjv1WcE5C/YSCu4P8BuG8hucRGhIK6DRTpc6UJv7TOXkSZs+GRx6BH36ANm1sA3DttRAaWjXfGcSMMWQdzCpsDFxdRVt+2UKesTOxwmuE07pxayLDIlm+YzknTp0oeH9ZXVBKqeI08XsiLw/efhv+8Q/YtAlatrQNwA03aANQDY6dPMb32d8Xuzv4Iv0LTplTpc5tFN6IFaNWcG6Dc3VWkVJl0MRfEXl58O678Pe/w8aN0KIF/PWvMHQo1PCZZQ5BIeThkDIXvQFER0TTLbYb3WK70T22Ox2bdaR2jdrVGKFSvqsiiV87t0ND7VX++vUwdy7Urg0jRtg7gJkz7diAqhZxUe5XbDet05SX+r/Eledfycbsjdy76F4uef0S6j1Wj66vduWuBXfx7vfvsvvQ7mqOWCn/pFf8JZ06BR98YO8A1qyBc86ByZPhppugZs3y36/OmKfTTPce2cs3O75h+Y7lLM9azrc7v+V43nEAEusnFtwVdIvtRtvGbQsGjpUKZNrV4w3GwEcf2QYgLQ3i4+G++2DkSG0AqtCZzOrJzctl9e7VtiHYsZxlO5bx8+GfATuttEtMF7rF2IagS0wXnVKqApImfm8yBj77DB5+GFJTITYWJk2C226DWrWcjk65YYwtgb18x3KWZS5jedZy1u9ZzylzCkFo07hNsbsCHTRWgUATf1UwBv73P9sALF8OzZvDvffCqFEQrgXNfN2h44dYuXNlQffQNzu+4cDxA0DhoHH32O50i+3mdtBY1xcoX6eJvyoZA59/bhuAr76Cs8+2DcDo0bZMhPILp8wpvs/+vqB7aPmO5WzdvxWw5bA7NutIt5hudI/rzu5Du7ln0T0elbhQyima+KvLkiV2DOCLL6BxY7j7bvjDHyBSi5b5o9MNGrvTKLwRL1z1Qrn7H4SFhFVZV5Kv3on4alyBTBN/dfvqK9sALFoEZ50Ff/kLjBsHdeo4HZmqhNy8XNbsXkOXV7tU6nNCJbTcxiGyZmSFN9RZ+ONC7lt8X8GGOuAbdyIVKQKovEcTv1OWL7cNwIIF0KgR/PnP9k7gkUdKF49TfqOsDXWa123OwpsWlrnvgdsNdk6Ws+GO632nW8jmiciwSGqE1HDk8fQ3Txerz5Qvtl4smRMzK/VzBSpv3CFp4ndaaqptAD75pPRrEREwY4Ymfz9S3Vewxhhy83LLbSCGvDukzM+4q+tdnDx1slof+bWXTqdxZGPiouKIi4ojPiq+4M/5j+iI6KCbYTVr3SzGfDSm0ndumvh9RdOm8PPPpY/Hx0N6erWHo86cL/ZZl3UnEh8VT/qd6dUejzGGPJPHOc+dw46DO0q9HlUrihta30DmgUwyD2SScSCj1HahtUJrlWoMij5i68X6/LagxhgOHj/I3iN7iz32HNlT6tjeI3vZd3Sf28+p6P9HTfy+IiTEzgIqScSuEFaqEny1L93TuIwx/Hrs14KGoOgj40AGmQcy2X1od6lur+iIaOLru+4W6pVuHBpHNi7zruFMG/DjJ4+7Tdp7j+xlb07pY7l5uW4/p2F4QxpHNi58RDTmhbQX3J4rCKce9DxPVCTxaxWyqhQXZ3f/Kik2tvQxpSqourf29HZcIkLD8IY0DG9I0tlJbj8rNy+XnQd3lm4cDmay+ZfNLNi2gCMnjhR7T63QWsRGxRY2Bq7GYdv+bUxNncqxk8cAyDiQwagPRrEpexPtmrQrM7HvObKHg8cPuo2vdo3aNIlsQuPIxjSt05T2TdoXT+yuR5PIJpwVcRZhoWGlPuPjrR+7vXMrq3aVN+gVf1VytyUkwHXXwTvvOBOTUgHkdHcN+Y9dh3Z5PFgeIiGcFXFWqavyxpGNaVKnSamEHhkWWekxCW/duekVv6/IH8DN3xIyNtY+5s6FN96Am292Nj6l/FxF7hrOff5ctw2AIGz4wwYaRzamYXjDai/q58Sdm17xV7fcXLjqKrv46+OPoU8fpyNSKij42mC4t3m9Hr+I9BORLSKyTUQmuXm9gYjME5H1IrJSRNoUeW2CiHwnIhtF5E7Pf4wAVbOm3fildWu73ePq1U5HpFRQmNJrChFhxcuqRIRFMKXXFIcick65iV9EQoFpwBVAK2CoiLQqcdpkYK0xph0wAnjO9d42wO1AJ6A90F9EzvNe+H6qXj07x79hQ3v1r1M7lapyw9sOZ8aAGcRHxSMI8VHxjs+AcoonffydgG3GmO0AIjIHGAR8X+ScVsDjAMaYzSKSICJNgJZAqjEmx/XepcA1wJPe+xH8VLNm8Omn0L079OsHy5bZ1b5KqSozvO3woEz0JXnS1dMcKLoaI8t1rKh12ISOiHQC4oEY4DvgUhFpJCIRwJWAzmXM16qV3e0rPR0GDoSjR8t9i1JKVZa39tx9HKgvImuB8cAaIM8Yswl4AlgIfAasBdyu6xaR0SKSJiJp2dnZXgrLD1x6Kbz5JnzzjZ0FlFf+snellKoMTxL/Topfpce4jhUwxhw0xow0xiRh+/ijge2u1141xnQ0xlwG/Ar84O5LjDEzjDHJxpjk6OjoM/hR/Nh118Gzz8K8eTBxovvVvkop5SWe9PF/C5wvIonYhD8EGFb0BBGpD+QYY3KBUcCXxpiDrtcaG2P2ikgctjuocjVuA9WECbBjBzz9tJ3rf/fdTkeklApQ5SZ+Y8xJEbkDWACEAq8ZYzaKyFjX69Oxg7hviIgBNgK3FfmId0WkEXACGGeMKV2vVVlPPglZWXDPPXZrx2HDyn+PUkpVkC7g8jXHj0Pfvra2/2efQc+eTkeklPIDXl/ApapRrVowfz60aAG//z2sX+90REqpAKOJ3xfVr2/n+NetC1deafv+lVLKSzTx+6rYWJv8Dx2CK66A33RoRCnlHZr4fVnbtrbb54cf4Oqrbf+/UkpVkiZ+X3f55TBzJixdass4685dSqlK0nr8/mDYMDvN8957ISYG/vlPpyNSSvkxTfz+4u67iy/wmjDB6YiUUn5KE7+/EIGpU2HnTlvWoXlzW+pBKaUqSPv4/UloqN3Ht1s3uPFG+OorpyNSSvkhTfz+Jjwc3n8fEhJg0CDYtMnpiJRSfkYTvz9q1MiWc6hVy27ismuX0xEppfxIwCT+lBR7ERwSYv+bkuJ0RFUsIcFu37h/v13de/Cg0xEppfxEQCT+lBQYPRoyMmwp+4wM+zzgk3+HDjB3LmzcaDduz811OiKllB8IiMR///2Qk1P8WE6OPR7w+vaFl1+GRYtg1CjdxEUpVa6AmM6ZmVmx4wHnllvsAq8HHrBz/KdMcToipZQPC4jEHxdnu3fcHQ8a999vF3g9+qhN/mPHOh2RUspHBURXz5QpEBFR/Fjt2kF24SsC06bBgAEwbpyd8qmUUm4EROIfPhxmzID4eJv/QkLsjMdrr3U6smpWowbMng3JyTB0KKxY4XRESikfFBCJH2zyT0+3xSs//thWNpg82emoHBAZCR99ZEs69O9vSzorpVQRAZP4i+rXD+64A5591k52CTrR0XYTl5AQ+8vYs8fpiJRSPiQgEz/AE0/AhRfaCS/79zsdjQPOO89e+e/ZA1ddBYcPOx2RUspHeJT4RaSfiGwRkW0iMsnN6w1EZJ6IrBeRlSLSpshrE0Vko4h8JyKzRaS2N3+AskRE2AVce/bAH/4QpNPbO3WCt96CNWtg8GA4edLpiJRSPqDcxC8iocA04AqgFTBURFqVOG0ysNYY0w4YATznem9z4E9AsjGmDRAKDPFe+Kd30UXw97/D228HwSresvTvD9On2/IOY8cGaQuolCrKkyv+TsA2Y8x2Y0wuMAcYVOKcVsDnAMaYzUCCiDRxvVYDCBeRGkAEUK0Vxe65By65xM5wdDfXPyjcfrtd3PXqq7YlVEoFNU8Sf3NgR5HnWa5jRa0DrgEQkU5APBBjjNkJ/BPIBHYDB4wxCysbdEWEhsKsWfZCd8QIyMurzm/3IQ8/DCNHwkMP2QZAKRW0vDW4+zhQX0TWAuOBNUCeiDTA3h0kAs2ASBG50d0HiMhoEUkTkbTs7GwvhWUlJMC//w1ffhnE29WKwEsv2do+Y8bYrRyDqpypUiqfmHL6fEWkK/CQMaav6/l9AMaYx8o4X4CfgHZAX6CfMeY212sjgC7GmD+e7juTk5NNWlpaBX+U0zMGbrjBLmhNTbWFLYPSoUPQtm3pfq+ICLsKbvhwZ+JSSlWKiKwyxiR7cq4nV/zfAueLSKKI1MQOzn5Q4gvru14DGAV8aYw5iO3i6SIiEa4GoRfgyJZRInaM86yzbG47etSJKHxA3bru+7uCppypUqrcxG+MOQncASzAJu23jTEbRWSsiORXAmsJfCciW7Czfya43psKzAVWAxtc3zfD6z+Fhxo1gpkz7W6Fk0pNSg0iO3e6Px405UyVCm7ldvU4oSq6eoqaMAGefx4WLIA+farsa3xXQoL7KU6hoXbWz623wtlnV3tYSqkz5+2unoDz+OPQqpVd1btvn9PROMBdOdOaNeGCC2x3T2ysHRD54gud969UAArKxB8ebiex/PKLneASdLmtZDnT+Hh47TW7heOWLfaWaPFi6NnT1r149tkgrXuhVGAKyq6efE8+Cffea/v9b765yr/Ovxw9avfznT4dli+3GxwMHmxX/3bubBsMpZTPqEhXT1An/rw8e1G7Zg2sWweJiVX+lf5p3Tq7BmDWLFvsLSnJNgDDhtlZQkopx2kfv4dCQ+E//7EXrzfdFMSresvTvj288ALs2mXvAIyxib9ZM1sBb/16pyNUSlVAUCd+sN3b06bBsmW2lLM6jbp17aDImjV2d69rr7X9ZO3bQ7duthUN2gUSSvmPoE/8YMc6Bw+GBx+EVaucjsYPiNh+/pkz7ZqA/MHfm2+GmBi46y7d+UspH6aJH5vHXnwRmjSBG2+0i1iVhxo2hDvvtKviPv8cfvc7u0jiggugVy87QHzihNNRKqWK0MTv0qABvPEGbN5sSzmrChKByy+3G7/s2GHXCvz4I1x/PcTFwV//qiuDlfIRmviL6NULJk60ff6ffup0NH7s7LPtTvc//mi3f0xOhkcftdOmBgyAjz/WkXSlHBTU0zndOXYMLr7YLu7asMEWdVNekJEBL78Mr7xi98OMj4fRo7U8hFJeotM5K6F2bbuqd/9+u3GVD7aL/ik+Hh55xHYDvfMOnHtuYXmIwYMLy0OkpOg+AUpVMb3iL8PTT8Nf/mI3q7r1VkdDCVxbttjSEa+/Dr/+aq/89+0rPhis+wQo5RFduesFp07ZCSrffgtr19oLVFVFjh61dwG33w65uaVfj4+H9PRqD0spf6JdPV4QEmJn+YSG2lW9J086HVEACw+3GyKXNe0zM1P73JTyIk38pxEba+f3f/ONLeWsqlhcnPvjxkCnTnaqqLbASlWaJv5yDB1qHw89ZLt9VBVyt09AeDiMHAkHDsCQIXD++XaB2OHDzsSoVADQxO+BadNsPbIbb4QjR5yOJoC52yfg5ZftXgGbN8O8edC8ud0vIC7OzgravdvpqJXyOzq466EvvrALvMaMsd0/ykHffAP//KdtCMLCbIt81112WzWlgpQO7laByy+3uWX6dLvwVDmoa1d4911bCG7UKJg9G1q3hv79YelSHQhWqhya+CvgkUegXTs7r3/vXqejUZx3nu2Hy8yEhx+GlSuhRw8dCFaqHB4lfhHpJyJbRGSbiExy83oDEZknIutFZKWItHEdv0BE1hZ5HBSRO739Q1SXWrXgzTftOKOu6vUhZ50Ff/ubLQsxfboOBCtVjnITv4iEAtOAK4BWwFARKdmZOhlYa4xpB4wAngMwxmwxxiQZY5KAjkAOMM+L8Ve7tm3hscfggw9s2RnlQ8LD7SDM5s0wf37pgeCff3Y6QqV8gidX/J2AbcaY7caYXGAOMKjEOa2AzwGMMZuBBBFpUuKcXsCPxpiMSsbsuAkT7EDvnXfC1q1OR6NKCQmBQYPg66/tRvE9e9rWOj7ejgls2uR0hEo5ypPE3xzYUeR5lutYUeuAawBEpBMQD8SUOGcIMPvMwvQtISF286maNXVVr8/r2tVuBpM/EPzf/9rZPwMG6ECwClreGtx9HKgvImuB8cAaoKDguojUBAYC75T1ASIyWkTSRCQtOzvbS2FVnZgYeOklSE21646Ujys5EJyaageCO3eGt9/W1lsFFU8S/04gtsjzGNexAsaYg8aYka6+/BFANLC9yClXAKuNMXvK+hJjzAxjTLIxJjk6OtrjH8BJN9xgp5D/4x82jyg/UHQg+KWX7EDw4ME6EKyCiieJ/1vgfBFJdF25DwE+KHqCiNR3vQYwCvjSGHOwyClDCZBunpL+/W87hnjjjZoz/Ep4uN0IZtMmHQhWQafcxG+MOQncASwANgFvG2M2ishYERnrOq0l8J2IbMFe3U/If7+IRAK9gfe8HbwviIqCWbPsLoN//rPT0agK83QgWDeIUQFESzZ4yb33wpNPwvvvw8CBTkejKmXbNnj2WbtBzNGjkJRkG4DjxwvP0Q1ilI/RjVgccPw4dOkCO3favXqblJzMqvzPL7/ACy/YweBTp0q/3qyZ3UWsTp3qj02pEjTxO2TjRujY0e7c9eGHtsCkCgAhIaef9lm3rh0jaNas7EfTpnZDZ6WqSEUSf42qDiaYtG5tu3smTLATRsaOLf89yg/ExdlZQCU1agT33AO7dhU+li2z/y3aLZSvYcPijYG7xqJJE1tx1FMpKXYwOjPTxjllinY/qXLpFb+XnToF/frZscI1a+CCC5yOSFVaSoqdAZSTU3jsdH38xtjN43ftsn1/RRuGoo/duyEvr/h7RaBx47LvHPIbi+hoW5W0InGpgKZdPQ7btcvW9DnnHDtRpCIXcMpHVcWVdV6eHUcor4HYu7d0V1ONGvZYyYYDdHP6IKWJ3wfMnQvXXw8PPAB//7vT0Si/duIE7NlT2BDkNxKPPur+fBH3g9EqoGni9xG33AJvvGG7bffu1S5Y5WUJCe7HHsDOMLjzTrjiCjs4rQKe7sDlIy65xF587dlj78ozMmyXrK79UV5R1ub0N9xg1x307w8tW9opqbpZtCpCE38VeuSR0l2zOTm2q1ipSitrc/q33oKffrKVSKOiYNw4W1Xw3nthx47yP1cFPO3qqUJlTf/WLlhVbYyxm9NPnWr3KRaB666z3UBdujgdnfIi7erxEXFxFTuulNeJQLdutvT09u0wcSJ89pndp6BLF5gzxw4eq6Ciib8KueuCBRg6tPpjUYr4eHjqKcjKsmVl9++3fxnPOQeeeMI+V0FBE38VKtkFGxtru1pfeMGWd1DKEXXq2H7/zZttbZEWLWDSJPsX9I9/tMdVQNPEX8WGD7draU6dsmt/li2zdwFXXaUl35XDQkLszJ/Fi2HdOhgyBF57zc4EuvJK+N//dGvKAKWJv5rFxcFHH0F2ti3fXHS1vVKOadcOXn21cGvK1auhTx9o08bOFDp61OkIlRdp4ndAx462zEpamr0jcLfqXilHNG5cuDXlzJlQs6ZdfBIbC3/9q10xrPyeJn6HDBxo9/qYP98WeFTKp9SqBTffbK/8lyyxqxEffdSuFr7pJli1yukIVSVo4nfQhAkwfjw884wd8FXK54jA//2fvULZutUO/s6fD8nJcOml8N57esvqhzTxO+zZZ2HAANsAfPyx09EodRrnnmsXgmVl2auVrCy49lo47zz7/MABpyNUHtLE77DQULuyPikJBg+2NfyV8mlRUXYh2LZtdjVwbCzcdZedqzxhAvz4o25O7+O0ZIOP2LXLLqTMy4PUVPtvSCm/sWqVvRt46y27Ejg0tHgXkG4QU+W8XrJBRPqJyBYR2SYik9y83kBE5onIehFZKSJtirxWX0TmishmEdkkIl09/1GCR7Nmtqvn0CE7tfrQIacjUqoCOnaEWbPsopWoqNL9/jk5cN99joSmSis38YtIKDANuAJoBQwVkVYlTpsMrDXGtANGAM8Vee054DNjzIVAe2CTNwIPRG3b2g1cvvvOVtY9edLpiJSqoGbN4OBB96/t2GFrBD3wACxdCrm51RubKuDJFX8nYJsxZrsxJheYAwwqcU4r4HMAY8xmIEFEmohIFHAZ8KrrtVxjzG9eiz4A9ekDL75o62iNH68LJ5UfKqsKYVSUnSX02GPQowc0aGBXCD/zDKxfr3/Zq5Enib85ULSId5brWFHrgGsARKQTEA/EAIlANvC6iKwRkVdEJLLSUQe422+3c/unT7f/JpTyK+6qE0ZEwLRpdhPqffvslNBbb7X7Btx1F7RvD2efbccAXn9d9w2oYt6a1fM4UF9E1gLjgTVAHlADuAh40RjTATgClBojABCR0SKSJiJp2dnZXgrLfz32mN2z9+677cQJpfyGuw1iig7sRkXBoEHwr3/ZncIyM22y793b1g269VZ713DBBXDHHbaR+C3AOwqqeRZUubN6XIOxDxlj+rqe3wdgjHmsjPMF+AloB0QAK4wxCa7XLgUmGWOuOt13BuOsHneOHoVevewUzyVLoHNnpyNSqooZY0vX/u9/sGiRHQs4csQmxIsvtnsJ9+5tp8DVquV0tN6RkmLLYhQt3HUGs6C8utm6iNQAfgB6ATuBb4FhxpiNRc6pD+QYY3JF5HbgUmPMCNdrXwGjjDFbROQhINIYc/fpvlMTf6HsbPt3/PBhWLECEhOdjkipapSba+c3L1pkG4OVK+2MoYgIuOwy2xD87nd2ZoSvbyp//Djs3GkXvhV9vPKK+yJ48fF2lpSHvJr4XR94JTAVCAVeM8ZMEZGxAMaY6a67gjcAA2wEbjPG/Op6bxLwClAT2A6MzH+tLJr4i9u82W6idPbZtqxzgwZOR6SUQw4csHcBixbZxybXJMHGje3tcX5DUN3b3OXkFCb1HTtKJ/esLHsVV1JUVNkrniu4R6vXE39108Rf2tKl9g73kkvsjJ+aNZ2OSCkfkJVlxwXyG4L8TS5atChsBC6/HOrXt8dTUuD+++24QlycHYgurzvl0CH3ibzow93uZY0a2ZWYZT2aN4e6dW2ffkZG6fc7fcVf3TTxuzdrFowYYYsmvv66vSBQSrkYA99/X9gttGRJ4fhAcjI0bWqvmo4fL3xPeDg8+KCdVVQ0kRe9ane3LqFx49Mn9ZgY+9me8MU+fido4i/bQw/ZfTL+8Q9bHl0pVYbcXDsmkN8QLF9e/ntEbJ9qWck8NtYuUvP2wPKZ3ImUCl0Tf8Ayxl7xz5pl/64MG+Z0REr5iZAQ94vERODrr21ib9oUwsKqPzYvqEjir1HVwSjvErE74WVmwsiR9gLk0kudjkopPxAX574vPS6M723pAAAZuklEQVTOzp4IIj4+/0m5U6uW3f8iMRGuvhp++MHpiJTyA2WtKJ4yxZl4HKSJ3081bGireYaEwFVXwS+/OB2RUj6uvBXFQUT7+P3cN9/Y2WodO9pZbbVrOx2RUsoJXq/Hr3xX1652oHf5crjllgqt91BKBSlN/AHg+uvhiSfs5kcPPOB0NEopX6ezegLE3XfbLVAffdTuiX3rrU5HpJTyVZr4A4SILXeekQFjxtgZar/7ndNRKaV8kXb1BJCwMHjnHWjZEq691la3VUqpkjTxB5h69ew0z8hIu6tdfs0qpZTKp4k/AMXGwocf2rn9AwbYOlVKKZVPE3+A6tgR5syB1avt+pS8PKcjUkr5Ck38AWzAAJg6Fd5/3876UUop0Fk9AW/8ePjxR3j2WTvNc9w4pyNSSjlNE38QePpp+Okn+NOfbHmS/v2djkgp5STt6gkCoaHw3/9Chw4wZIjt91dKBS9N/EEiMtLO9GnY0F7x79jhdERKKado4g8iTZvaOf6HD0P37nZ1b0iI3es5JcXp6JRS1cWjxC8i/URki4hsE5FJbl5vICLzRGS9iKwUkTZFXksXkQ0islZEtNayw9q2hT/+0V7x79hhd6LLyLB7PWvyVyo4lJv4RSQUmAZcAbQChopIqxKnTQbWGmPaASOA50q8frkxJsnTWtGqas2ZU/pYTo7d61kpFfg8ueLvBGwzxmw3xuQCc4BBJc5pBXwOYIzZDCSISBOvRqq8JjPT/fGMDNi3r3pjUUpVP08Sf3Og6FBglutYUeuAawBEpBMQD8S4XjPAIhFZJSKjy/oSERktImkikpadne1p/OoMxMWV/VpMDNx2G6xZU33xKKWql7cGdx8H6ovIWmA8sAbILxJwiTEmCdtVNE5ELnP3AcaYGcaYZGNMcnR0tJfCUu6Utef0Y4/BzTfbrqCLLoJLLrF/zs11Jk6lVNXwJPHvBGKLPI9xHStgjDlojBnpSvAjgGhgu+u1na7/7gXmYbuOlIPK2nN60iSYPh2ysuCZZ2xlz6FD7ayfhx+G3budjlwp5Q2eJP5vgfNFJFFEagJDgA+KniAi9V2vAYwCvjTGHBSRSBGp6zonEugDfOe98NWZGj4c0tPtHr3p6fZ5vgYNYOJE+OEHO/0zKQkeesh2EQ0davf3NcahwJVSlVZu4jfGnATuABYAm4C3jTEbRWSsiIx1ndYS+E5EtmC7dCa4jjcBvhaRdcBK4GNjzGfe/iFU1QgJsTX9P/nENgJ33AGffmrXAHTsCK+9BkePOh2lUqqixPjgpVtycrJJS9Mp/77o8GE73/9f/7I7fDVsCKNG2bUB8fFOR6dU8BKRVZ5OmfebxH/ixAmysrI4duyYQ1H5rtq1axMTE0NYWFi1facxsHQp/PvfMH++fT5ggL0r6NXLjh0opapPQCb+n376ibp169KoUSNEs0oBYwz79u3j0KFDJCYmOhLDjh12UHjGDLvr14UX2gZgxAioW9eRkJQKOhVJ/H5Tq+fYsWOa9N0QERo1auTonVBsrJ0iumMHvPEG1KljE3/z5rYU9JYtjoWmlHLDbxI/oEm/DL7ye6ld217lf/stpKbCoEHw0kv2DqBPH1sdVLeAVMp5fpX4nbZnzx6GDRvGOeecQ8eOHenatSvz5s3z+vfMnDmTO+64w+ufW506dYJZs2x5iEcege+/h4ED4bzz4KmnYP9+pyNUKngFbOJPSbELj7xVdtgYw9VXX81ll13G9u3bWbVqFXPmzCErK6vYeSdPnqzcFwWYJk1s8beffoJ33rEzf+65x3YDjRoF69Y5HaFSwScgE39Kii0znJHhvbLDn3/+OTVr1mTs2LEFx+Lj4xk/fjwzZ85k4MCB9OzZk169egHw1FNPcfHFF9OuXTsefPDBgve8+eabdOrUiaSkJMaMGUOeq+/j9ddfp0WLFnTq1Illy5YBFAzYnjhxAoCDBw8We+5PwsLguutgyRKb7EeMgNmz7eKwSy+Ft96CEye832ArpUrzyz1377wT1q4t+/UVK+D48eLHcnJs8bGXX3b/nqQkmDq17M/cuHEjF110UZmvr169mvXr19OwYUMWLlzI1q1bWblyJcYYBg4cyJdffkl0dDRvvfUWy5YtIywsjD/+8Y+kpKTQu3dvHnzwQVatWkVUVBSXX345HTp0oG7duvTo0YOPP/6Yq6++mjlz5nDNNddU67TNqtCune37f/xxmDkTpk2zW0JGRdn/T/ntWn6DDcVXFiulKicgr/hLJv3yjp+JcePG0b59ey6++GIAevfuTcOGDQFYuHAhCxcupEOHDlx00UVs3ryZrVu3snjxYlatWsXFF19MUlISixcvZvv27aSmptKjRw+io6OpWbMmgwcPLvieUaNG8frrrwP2rmDkyJHe+yEcVrI0xPHjhUk/X04OTJ7sTHz+QO+Q1Jnwyyv+012Zg/0HkJFR+nh8vO1qOBOtW7fm3XffLXg+bdo0fvnlF5KT7bTZyMjIgteMMdx3332MGTOm2Gf861//4uabb+axxx4rdnz+/Pllfm/37t1JT09nyZIl5OXl0aZNmzLP9Vf5pSHKapgzM2HMGDszqFcvqF+/euPzVfldmjk59rneISlPBeQVf1llh6dMOfPP7NmzJ8eOHePFF18sOJaT/y+uhL59+/Laa69x+PBhAHbu3MnevXvp1asXc+fOZe/evQDs37+fjIwMOnfuzNKlS9m3bx8nTpzgnXfeKfZ5I0aMYNiwYQF1te9OWfsEhIfb8YDrroNGjaBbN1s0bvlyCOax9MmTC5N+Pt1JTXkiIBN/WWWHK3MVJCLMnz+fpUuXkpiYSKdOnbj55pt54oknSp3bp08fhg0bRteuXWnbti3XXXcdhw4dolWrVjzyyCP06dOHdu3a0bt3b3bv3k3Tpk156KGH6Nq1K927d6dly5Ylfp7h/PrrrwwdOvTMfwA/UFaD/fLLdmewr76ySe3UKfjHP2yxuLPOgmuusWMGP/3kTNzVJTvbFsx76CG44oqyd1Ir67hS+fymZMOmTZtKJcRgMXfuXN5//31mzZpV5jmB8vtJSbHJPTPT3gFMmeK+wd6/HxYvhoULYcECu2oY7DqBPn3s4/LLoV696o3fW3Jy7C5oqamwcqV95DdsItC6tX1+5Ejp9zZoYBuJ0NDqjVk5KyBr9QRKYquo8ePH8+mnn/LJJ5/QokWLMs8L1t8P2Cm7P/xQ2Ah88YVNnDVqQNeuhQ1Bx46+mQzz8uwCt/wEv3IlbNhQuMo5NhY6d7aL4jp1sruj1a1buo8f7HjJqVNw8cX2LqhDB2d+JlX9NPEHIf39FDp+HL75xjYECxfCqlX2eMOG8LvfFTYEsbGn/5yqYIy9Oyma5NPSCq/co6IKE3ynTjaBN21a9ue5u0MKCbGzpbKzba2kv/9di+UFA038QUh/P2XLzoZFiwobgl277PELL4S+fW0j8H//B0UmZnnNr7/axJ6f5FNTYc8e+1rNmnb9SNFEf/75NnFX1m+/2cHf6dOhWTN4/nn4/e+1XHYg08QfhPT34xljbLfKggW2EVi6FI4dsyuLL7nENgJ9+0L79sUTsCdjD8eP21XJ+Ql+5UrbBZXvggtscs/vtmnXDmrVqtqfd8UKGDvWxtW/v91AJyGhar9TOUMTfxDS38+ZOXYMvv66cHxg/Xp7PDoaeve2DcGhQ3DvvcX70iMi7OyaJk0Kr+bXri1cgHb22cX75ZOTnVt/cPKkveL/299s//+DD8Kf/2wbOxU4NPEHIf39eMfu3cW7hVxLLk4rMtL2xRftsomJ8b1ulR07bJ///PnQpo3tBure3emolLcE5EYsvqCsssxLliwhKiqKpKQkLrzwQv7yl7+Ueu/VV19Nly5dih3bsmULPXr0ICkpiZYtWzLateyy6Oe1bNmShx9+uFp+PmUHUm+6yZaU3r3bTqk8nQ0b4MABO5PoiSfg2mvtoLGvJX2wcc2bB++/DwcP2q6t22/XEtnByKPELyL9RGSLiGwTkUluXm8gIvNEZL2IrBSRNiVeDxWRNSLykbcCL5eXi5iUV5b50ksvZe3ataxZs4aPPvqooMImwG+//caqVas4cOAA27dvLzj+pz/9iYkTJ7J27Vo2bdrE+PHjC17L/7y0tDTefPNNVq9eXan4VcWFhNjB17I2kY+Pt1fOvjhF9HQGDoSNG+Huu+H11+3Yw3/+Y8c/VHAoN/GLSCgwDbgCaAUMFZFWJU6bDKw1xrQDRgDPlXh9ArCp8uF6qArqMp+uLHNR4eHhJCUlsXPnzoJj7733HgMGDGDIkCHMmTOn4Pju3buJiYkpeN62bdtS3xsZGUnHjh3Ztm3bGceuKqcqSoA4rU4dePJJWL3aziS6+Wbo2RM2b3Y6MlUdPCnS1gnYZozZDiAic4BBwPdFzmkFPA5gjNksIgki0sQYs0dEYoCrgCnAn70StQN1mcsry5zv119/ZevWrVx22WUFx2bPns3f/vY3mjRpwrXXXstkV7nJiRMn0rNnT7p160afPn0YOXIk9UuMAO7bt48VK1bwwAMPlPvdqmrkz97xZEWxv2nXzg5uv/KKHcBu187+d/JkWyNJBSZPunqaAzuKPM9yHStqHXANgIh0AuKB/EvZqcA9wKlKRVoR1VCXuWRZ5q+++or27dvTvHlz+vbty9lnnw3YcYGtW7dyySWX0KJFC8LCwvjuu+8AGDlyJJs2beL6669nyZIldOnSheOuGL/66is6dOhAnz59mDRpEq1bt/Za7Krihg+H9HQ7KyY9PTCSfr6QEHtDvGULDB5st8ps29YObqsAZYw57QO4DnilyPObgH+XOKce8DqwFpgFfAskAf2BF1zn9AA+Os33jAbSgLS4uDhT0vfff1/qWJni442xnTzFH/Hxnn9GCYsWLTKXXXZZsWPZ2dkmPj7efPHFF+aqq64yxhizfft207hxY7NmzRpjjDHPP/+8qVevnomPjzfx8fGmQYMGZvLkyW6/o3Xr1iYtLa3Y53mqQr8fpU5j8WJjWrSw/2SGDDFm1y6nI1KeANJMOfk8/+HJFf9OoOji9hjXsaKNx0FjzEhjTBK2jz8a2A50BwaKSDowB+gpIm+W0QDNMMYkG2OSo6OjPQjrNKqgU9bTssyJiYlMmjSpoGrn7Nmz+eyzz0hPTyc9Pb1gUBjgs88+K9hG8eeff2bfvn00b17yZkqp6tWzp13P8PDDdhbQhRfCCy8U1g5S/s+TxP8tcL6IJIpITWAI8EHRE0Skvus1gFHAl67G4D5jTIwxJsH1vs+NMTd6MX73qqAuc0XKMo8dO5Yvv/yS9PR0MjIyik3jTExMJCoqitTUVBYuXEibNm1o3749ffv25amnniroIlLKSbVq2QVfGzbYdQnjxtl9EMqb3qr8g0cLuETkSmxffSjwmjFmioiMBTDGTBeRrsAbgAE2ArcZY34t8Rk9gL8YY/qX9326gKvi9PejqooxMGeOFn7zdV5fwGWM+cQY08IYc64xZorr2HRjzHTXn79xvX6BMeaakknfdc4ST5K+Usq3iMDQoXaq55gx8Nxz0LIlvPeezv33V7pyVynlkfr1bV//8uV257Nrr7WLwdLTnY5MVZQmfqVUhXTpYktNP/20LVXRurVdDJZfoE75Pk38SqkKq1HDVvjctMlWML33XrszWJFKJaoCvFxhplya+JVSZ6yswm8zZlRvIvNnVVBhplya+JVSlVa08Nurr9pB4OpMZP5s8uTiez2AfX7//VX3nZr4K6Csssz57rzzTpo3b86pU9VXnUIpX5Ff+M3dUpScHNso6CIwyM215cSeeQauu87Wf3KnrOPeELCJP2VDCglTEwh5OISEqQmkbKjassynTp1i3rx5xMbGsnTpUm/8CEr5pZ9/dn989267mXyPHnDPPTB3rk1ugT4ldO9e2xV2771w6aVQrx507Qp33WWro5YsMpAvLq7qYvKkOqffSdmQwugPR5Nzwt4/ZRzIYPSHdpOT4W3PbPVueWWZlyxZQuvWrRk8eDCzZ8/m8ssvr+RPoZR/iouz3TslNWoEw4bZbSqfe85e+YLdvjJ/57LOne02lQ0aVG/M3nLqlN3TefnywsfWrfa1sDDo2LFwFXTXrtCsWWEff8mtPauy7LdfJv47P7uTtT+XXZZ5RdYKjucVr8SZcyKH296/jZdXuS/LnHR2ElP7nXlZ5tmzZzN06FAGDRrE5MmTOXHiBGG6qakKQlOmuE9kzz1XWDXl+HFbDyh/v+KVK+HDDwvPb9Gi+FaWSUlVvzH9mTh0yMa+fLmd0bRihd2RDey+zd26wahR9r/JyVC7dunPcKLst18m/vKUTPrlHT8T48aN4+uvv6ZmzZosW7aMTz75hGeeeYa6devSuXNnFixYQP/+ulBZBR9PElmtWnaf4osvtlfAYBNmWppNpKmpdu/jN10lHcPCbPIv2hi0aGFnDVWX/IHq/CS/fLltvE6dsqubW7e2Za27d7eJ/txzPd+Cc/jw6i317ZeJ/3RX5gAJUxPIOFD6XjM+Kp4ltyw5o+9s3bo17777bsHzadOm8csvv5CcnMyCBQv47bffCnbQysnJITw8XBO/ClpnksiioqBXL/sAm2h37ix+V/DGGzBtWuH5JTe5b9rUez9Dbq7tgy/abbN7t32tTh27kO2vf7VJvnNnu7LZX/hl4i/PlF5TivXxA0SERTClV+XKMk+ePJkXX3yRP/zhD0BhWebZs2fzyiuvMHToUACOHDlCYmIiOTk5RJQ1cqOUOi0RiImxj2uuscfy8mzNoKKNwZNPwsmT9vWYmMKxgk6dbJ960WJyKSll34ns3QvffFOY5L/9tnDvpsREW666Wzf7aNPGLmLzVx5V56xu3qjOmbIhhfsX30/mgUziouKY0mvKGQ/s5tu9ezcTJ04kNTWV6OhoIiMjueWWW5g4cSLp6enUq1ev4NxrrrmGwYMHM3jw4Ep9p6e0OqcKVkeP2p1YU1MLG4Mff7SviUCrVrYRAJg9G44dK3xvzZq2kfj559KDsPlJPn8Q1tdVpDpnwCb+YKO/H6UK7dtnr9jzG4LUVPjlF/fnhoTAgAE2yXfvbpO+u0FYX1eRxO/HNytKKeVeo0bQr599gB0vCA11v2bAGJg/v3rjc1rALuBSSql8ImUviKrKhVK+ShO/UiooVMFW3H7LrxK/L45H+AL9vShVvirYittv+U0ff+3atdm3bx+NGjVCPF0VEQSMMezbt4/a/jgapVQ1q+6FUr7KbxJ/TEwMWVlZZGdnOx2Kz6lduzYxMTFOh6GU8hN+k/jDwsJITEx0OgyllPJ7ftXHr5RSqvI08SulVJDRxK+UUkHGJ0s2iEg24GYrB4+cBZSxONtRGlfFaFwVo3FVTCDGFW+MifbkRJ9M/JUhImme1quoThpXxWhcFaNxVUywx6VdPUopFWQ08SulVJAJxMQ/w+kAyqBxVYzGVTEaV8UEdVwB18evlFLq9ALxil8ppdRpBEziF5HXRGSviHzndCz5RCRWRL4Qke9FZKOITHA6JgARqS0iK0VknSuuh52OqSgRCRWRNSLykdOxFCUi6SKyQUTWikha+e+oHiJSX0TmishmEdkkIl19IKYLXL+n/MdBEbnT6bgARGSi6+/9dyIyW0R8osKhiExwxbSxqn9XAdPVIyKXAYeB/xhj2jgdD4CINAWaGmNWi0hdYBVwtTHme4fjEiDSGHNYRMKAr4EJxpgVTsaVT0T+DCQD9Ywx/Z2OJ5+IpAPJxhifmv8tIm8AXxljXhGRmkCEMeY3p+PKJyKhwE6gszHmTNfneCuW5ti/762MMUdF5G3gE2PMTIfjagPMAToBucBnwFhjzLaq+L6AueI3xnwJ7Hc6jqKMMbuNMatdfz4EbAKaOxsVGOuw62mY6+ETVwAiEgNcBbzidCz+QESigMuAVwGMMbm+lPRdegE/Op30i6gBhItIDSAC2OVwPAAtgVRjTI4x5iSwFLimqr4sYBK/rxORBKADkOpsJJarO2UtsBf4nzHGJ+ICpgL3AKecDsQNAywSkVUiMtrpYFwSgWzgdVf32CsiEul0UCUMAWY7HQSAMWYn8E8gE9gNHDDGLHQ2KgC+Ay4VkUYiEgFcCcRW1Zdp4q8GIlIHeBe40xhz0Ol4AIwxecaYJCAG6OS61XSUiPQH9hpjVjkdSxkucf3OrgDGuboXnVYDuAh40RjTATgCTHI2pEKurqeBwDtOxwIgIg2AQdgGsxkQKSI3OhsVGGM2AU8AC7HdPGuBvKr6Pk38VczVh/4ukGKMec/peEpydQt8AfRzOhagOzDQ1Zc+B+gpIm86G1Ih19Uixpi9wDxsf6zTsoCsIndsc7ENga+4AlhtjNnjdCAuvwN+MsZkG2NOAO8B3RyOCQBjzKvGmI7GmMuAX4Efquq7NPFXIdcg6qvAJmPMM07Hk09EokWkvuvP4UBvYLOzUYEx5j5jTIwxJgHbPfC5McbxqzEAEYl0DdDj6krpg709d5Qx5mdgh4hc4DrUC3B08kAJQ/GRbh6XTKCLiES4/n32wo69OU5EGrv+G4ft3/9vVX2X3+zAVR4RmQ30AM4SkSzgQWPMq85GRXfgJmCDqz8dYLIx5hMHYwJoCrzhmm0RArxtjPGpqZM+qAkwz7Xfcw3gv8aYz5wNqcB4IMXVrbIdGOlwPEBBA9kbGON0LPmMMakiMhdYDZwE1uA7q3jfFZFGwAlgXFUO0gfMdE6llFKe0a4epZQKMpr4lVIqyGjiV0qpIKOJXymlgowmfqWUCjKa+JVSKsho4ldKqSCjiV8ppYLM/wNoti1I+Z8LbgAAAABJRU5ErkJggg==" + }, + "execution_count": 35, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "plt.figure()\n", + "p1 = plt.plot(pnts_rho, pnts_greedy, 'bo-', label=\"Greedy\")\n", + "p2 = plt.plot(pnts_rho, pnts_grasp, 'ro-', label=\"GRASP\")\n", + "p3 = plt.plot(pnts_rho, pnts_ga, 'go-', label=\"GA\")\n", + "plt.legend(loc='best')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discussion" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For $1\\leq \\rho\\leq 9$ we notice that:\n", + "\n", + "- Greedy gets about 94-98% of the clauses satisfied.\n", + "- GRASP improves this to about 97-100%.\n", + "- GA performs better than GRASP at about 99-100%.\n", + "- For all of them, the quality decreases as $\\rho$ increases." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Special cases" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* $n=1$.\n", + "\n", + " Each clause will fall into one of 3 categories as follows:\n", + "\n", + " 1. Only $x$ appears. The clause reduces to: $x\\lor x\\lor x = x$\n", + " 2. Only $\\bar x$ appears. The clause reduces to: $\\bar x\\lor \\bar x\\lor \\bar x = \\bar x$.\n", + " 3. Both $x$ and $\\bar x$ appear. The clause becomes a tautology: $x\\lor \\bar x \\lor \\ldots = T$, and the clause can be ignored.\n", + "\n", + " So $\\phi$ simplifies to a conjunction of $x$'s and $\\bar x$'s, whose satisfiability is easy to establish:\n", + " * If both $x$ and $\\bar x$ appear then return **no**.\n", + " * Otherwise, return **yes**.\n", + "\n", + "* $n=2$.\n", + "\n", + " We get 2-SAT which is known to be in **P**, (Krom, 1967).\n", + "\n", + "* $\\ell=1$.\n", + "\n", + " Always satisfiable.\n", + " \n", + " [TODO: true for $\\ell\\leq n$?]\n", + "\n", + "* Nested satisfiability, (Knuth, 1990).\n", + "\n", + " If the clauses have a \"hierarchical structure\", as defined in (Knuth, 1990), then it is solvable in linear time." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Conclusion" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* If instance is a special case then can be solved in polynomial time.\n", + "* Exhaustive search useful when $n$ is small.\n", + "* Dynamic programming useful when $\\ell$ is small.\n", + "* Otherwise, use GRASP, GA, or other metaheuristics for approximate solutions." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Reflection" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Reflection ...." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# References" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* Garey, S. and Johnson, D. (1979)\n", + "**Computers and Intractability: A Guide to the Theory of NP-Completeness.**\n", + "Freeman.\n", + "\n", + "* Hoos, H. and Stutzler, T. (2005)\n", + "**Stochastic Local Search: Foundations and Applications.**\n", + "Morgan Kaufmann.\n", + "\n", + "* Knuth, D.E. (1990).\n", + "**Nested satisfiability.**\n", + "Acta Informatica, 28(1), 1-6.\n", + "\n", + "* Koutsoupias, E., & Papadimitriou, C. H. (1992).\n", + "**On the greedy algorithm for satisfiability.**\n", + "Information Processing Letters, 43(1), 53-55.\n", + "\n", + "* Krom, M.R. (1967).\n", + "**The Decision Problem for a Class of First-Order Formulas in Which all Disjunctions are Binary**\n", + "Zeitschrift für Mathematische Logik und Grundlagen der Mathematik, 13: 15–20, [doi:10.1002/malq.19670130104](doi:10.1002/malq.19670130104).\n", + "\n", + "* Sipser, M. (2013).\n", + "**Introduction to the theory of computation**\n", + "(3rd international ed.). Cengage Learning." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (Ubuntu Linux)", + "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.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}