Permalink
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
TSP-Guidance/Investigating 3SAT.ipynb
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
1033 lines (1033 sloc)
74.6 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "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 | |
| } |