Z3 Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction \n",
    "\n",
    "## Hi, I'm Phil!\n",
    "\n",
    "I like math, computers, physics, and long walks on the beach. You can check out what I'm up to at https://www.philipzucker.com/ or give me a twoot at [@SandMouth](https://twitter.com/SandMouth). The video corresponding to this notebook is available here https://youtu.be/56IIrBZy9Rc\n",
    "\n",
    "I work at Draper Labs where we use Z3 as a backend for tools like the Comparative Binary Analysis (CBAT) project https://github.com/draperlaboratory/cbat_tools which formally verifies binary program transformations.\n",
    "\n",
    "## Formal Methods \n",
    "\n",
    "What are formal methods? Well, let's google it:\n",
    "\n",
    "\"In computer science, specifically software engineering and hardware engineering, formal methods are a particular kind of mathematically based techniques for the specification, development and verification of software and hardware systems.[1] The use of formal methods for software and hardware design is motivated by the expectation that, as in other engineering disciplines, performing appropriate mathematical analysis can contribute to the reliability and robustness of a design.[2]\" - https://en.wikipedia.org/wiki/Formal_methods\n",
    "\n",
    "Broadly speaking, Formal Methods are the application of mathematical thinking to computational systems.\n",
    "\n",
    "More narrowly, Formal Methods make the connection between specifications and implementations rigid and reliable.\n",
    "\n",
    "A specification is a precise, clear, and succinct statement of what a program or system should do. It often has a logical or mathematical flavor.\n",
    "\n",
    "An implementation is also precise, but often significantly more verbose. It needs much more detail to actually describe how the system achieves it's goals. It often has a programming flavor. \n",
    "\n",
    "The two should be equivalent in some sense. The process of showing this is verification. This reduces the possibility of bugs in the implementation to bugs in the specification and in this sense a program can be proven correct.  \n",
    "\n",
    "How do we prove that other than by waving our hands and hoping? Well, it's pretty hard. One astonishingly powerful and versatile tool we'll explore today is the Z3 theorem prover.\n",
    "\n",
    "## What is Z3?\n",
    "\n",
    "Well, let's google it.\n",
    "\n",
    "\"Z3 is a state-of-the art theorem prover from Microsoft Research. It can be used to check the satisfiability of logical formulas over one or more theories. Z3 offers a compelling match for software analysis and verification tools, since several common software constructs map directly into supported theories.\"\n",
    "\n",
    "### What the heck does that mean?\n",
    "\n",
    "There are a number of different interesting things you can do with Z3. It has some simplification capabilities that overlap with computer algebra systems like Mathematica, some capabilities that overlap with constraint satisfaction solvers like Minizinc, and a capability to manipulate and express logical statements. The sum total of these capabilities enables z3 to effectively and automatically find solutions to a wide variety of problems or (importantly!) effectively and automatically prove that such solutions _do not exist_.\n",
    "\n",
    "Z3 was designed with the intention that it is going to be used for hardware and software verification tasks.\n",
    "Many verification tools extensively use Z3 or equivalent solvers as the heavy lifting. They are often interfaces that translate the task at hand into a bunch of queries to Z3. But Z3 is also useful on it's own.\n",
    "\n",
    "Z3 was open sourced in the beginning of 2015. The source code is licensed under MIT License and hosted on GitHub.\n",
    "\n",
    "Z3 has official bindings for C, C++, Python, .NET, OCaml, Java and often bindings available in other languages as well. It also has a command line interface that takes in a standardized logical formula input format called [smtlib](http://smtlib.cs.uiowa.edu/). The python bindings are full featured, and python is a familiar language to many, so that is our choice today.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Walking through a Z3 problem\n",
    "\n",
    "To get started, let's do some simple examples\n",
    "\n",
    "First of you need Z3 installed. Z3 is available on pip, the python package installer. \n",
    "Oh, not _that_ Z3. Amusingly, `pip install z3` actually installs a completely unrelated tool.\n",
    "The command you actually need is\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install z3-solver"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a simple problem by first importing the library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from z3 import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then declare our variables, state a set of constraints we wish to hold, and then call the convenience function `solve` to get a solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[y = 0, x = 7]\n"
     ]
    }
   ],
   "source": [
    "x = Int('x')\n",
    "y = Int('y')\n",
    "solve(x > 2, y < 10, x + 2*y == 7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can ask the same question using the longer form syntax of building a solver object and adding constraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sat\n",
      "[y = 0, x = 7]\n"
     ]
    }
   ],
   "source": [
    "x = Int('x')\n",
    "y = Int('y')\n",
    "\n",
    "s = Solver()\n",
    "s.add(x > 2)\n",
    "s.add(y < 10)\n",
    "s.add(x + 2 *y == 7)\n",
    "\n",
    "print(s.check())\n",
    "print(s.model())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use z3 to prove properties and theorems. We'll get more into what this means later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "p &rArr; p"
      ],
      "text/plain": [
       "Implies(p, p)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = Bool(\"p\")\n",
    "my_true_thm = Implies(p, p)\n",
    "my_true_thm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "proved\n"
     ]
    }
   ],
   "source": [
    "prove(my_true_thm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the property is not true, z3 can supply a counterexample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "q &rArr; p"
      ],
      "text/plain": [
       "Implies(q, p)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = Bool(\"q\")\n",
    "my_false_thm = Implies(q, p)\n",
    "my_false_thm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "counterexample\n",
      "[p = False, q = True]\n"
     ]
    }
   ],
   "source": [
    "prove(my_false_thm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Z3 Variables\n",
    "\n",
    "You can create z3 variables of different sorts. Here is an incomplete list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Bool(\"x\")\n",
    "x = Const(\"x\" , BoolSort())\n",
    "p, q, r = Bools(\"p q r\") # convenience function for multiple definitions\n",
    "x = Real(\"x\")\n",
    "y = Int(\"x\") \n",
    "v = BitVec(\"n\", 32) # 32 bit bitvector\n",
    "f = FP(\"f\", Float64()) #Floating point values\n",
    "a = Array(\"a\", IntSort(), BoolSort()) # arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Z3 constants\n",
    "\n",
    "Often values are cast in Z3 expressions for you automatically, but sometimes they aren't. You can explicitly cast constants when this is the case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(BoolVal(True))\n",
    "print(RealVal(2.5))\n",
    "print(IntVal(7))\n",
    "print(Q(1,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constraints / operations\n",
    "\n",
    "Operations and constraints in Z3 are built mostly by overloading the ordinary python operations. There is a good chance you can guess what the right operator is. All operations are also available as capitalized functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Common numerical operations\n",
    "print(x + x) # addition\n",
    "print(x * x) # multiplication\n",
    "print(x ** 4) # fourth power of x \n",
    "\n",
    "# Common boolean operations\n",
    "print(And(p,q))\n",
    "print(And(p,q,r)) # some can accept more than one argument for convenience.\n",
    "print(Or(p,q))\n",
    "print(Implies(p,q))\n",
    "print(Xor(p,q))\n",
    "print(Not(p))\n",
    "\n",
    "# Constraints \n",
    "print(x == x) #equal\n",
    "print(x != x) # not equal\n",
    "print(x <= RealVal(3)) # inequality\n",
    "print(Or( x < 3,  x == 3, x > 3  )) # Constraints are operations into booleans\n",
    "\n",
    "# And many more! Look in the API docs for others https://z3prover.github.io/api/html/namespacez3py.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've shown some of the basics, you're ready to try something simple on your own\n",
    "\n",
    "### Exercise 1\n",
    "Get a root of the polynomial `x**3 + 3*x**2 + 4*x + 2 == 0` using z3. Can you use z3 to show this is the only solution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's take a look\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "xs = np.linspace(-3,3,200)\n",
    "plt.plot(xs,xs**3 + 3*xs**2 + 4*xs + 2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FILL IN\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### * General purpose python exploration tricks\n",
    "\n",
    "I find it extremely useful to explore libraries interatively in python. It is often faster and more honest than looking up the documenation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "p = Real(\"p\") # Define the z3 expression of Real constant value 3.\n",
    "type(p) # What is the type of this thing\n",
    "help(p) # What is the documentation. z3 has quite lovely documentation actually.\n",
    "dir(p) # Does this thing have any suspicious looking functions?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### * Getting Values out of Z3 models\n",
    "\n",
    "Z3 has a number of helper functions like `solve` above, but a more powerful interface is through the Solver class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Int(\"x\")\n",
    "y = Real(\"y\")\n",
    "z = Real(\"z\")\n",
    "p = Bool(\"p\")\n",
    "a = BitVec(\"a\",32)\n",
    "\n",
    "s = Solver()\n",
    "s.add(x == 2)\n",
    "s.add(y == 1/3)\n",
    "s.add(p == True)\n",
    "s.add(a == 0xDEADBEEF)\n",
    "print(s.check())\n",
    "m = s.model()\n",
    "\n",
    "print(m[x].as_long())\n",
    "print(m[y].numerator_as_long() / m[y].denominator_as_long())\n",
    "print(m[a].as_long())\n",
    "print(bool(m[p]))\n",
    "\n",
    "m.eval(2*x + 1) # you can also evaluate z3 expressions in the context of the model\n",
    "      \n",
    "s = Solver()     \n",
    "s.add(z ** 2 == 2)\n",
    "print(s.check())\n",
    "m = s.model()\n",
    "print(m[z].approx())\n",
    "print(m[z].approx().numerator_as_long() / m[z].approx().denominator_as_long())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### * Conversion Functions\n",
    "\n",
    "It is useful to sometimes convert between sorts. There are conversion functions for that for many different sorts. Here is a simple example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Int(\"x\")\n",
    "y = Real(\"y\")\n",
    "ToReal(x)\n",
    "ToInt(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### * Simplifying expressions with z3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Z3 can be used to simplify expressions similar to how one might use sympy or mathematica."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Int('x')\n",
    "y = Int('y')\n",
    "print(simplify(x + y + 2*x + 3))\n",
    "print(simplify(x < y + x + 2))\n",
    "print(simplify(And(x + 1 >= 3, x**2 + x**2 + y**2 + 2 >= 5))) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### * Fresh Variable Names\n",
    "It can be useful to generate variables with fresh names, especially when building problems programmatically. Z3 has built in functionality for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(FreshReal(prefix=\"x\"))\n",
    "print(FreshBool(prefix=\"b\"))\n",
    "print(FreshConst(IntSort() , prefix=\"f\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Useful Documentation\n",
    "\n",
    "Now that we have a simple example, let's look at where we can go to find out more\n",
    "\n",
    "- The API docs: https://z3prover.github.io/api/html/namespacez3py.html\n",
    "- https://rise4fun.com/z3/tutorial has an interactive tutorial in z3. It uses the smt-lib2 input language, not python\n",
    "- https://ericpony.github.io/z3py-tutorial/guide-examples.htm A tutorial for the python bindings upon which this tutorial is heavily based\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Z3 To Find Stuff\n",
    "\n",
    "Z3 is from one perspective just executing a flexible and fast search similar to constraint satisfaction solvers, SAT solvers, and computer algebra systems. You can use Z3 to solve puzzles like Sudoku and Nonograms. \n",
    "It isn't just fun and games, there are many practical usages of such things (it's a whole industry).\n",
    " - Graph Coloring (Useful for register allocation for example)\n",
    " - Scheduling Problems\n",
    " - Routing problems\n",
    " - etc.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8 Queens\n",
    "\n",
    "The eight queens puzzle is the problem of placing eight chess queens on an 8x8 chessboard so that no two queens attack each other. Thus, a solution requires that no two queens share the same row, column, or diagonal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Q_5 = 8,\n",
      " Q_8 = 6,\n",
      " Q_3 = 2,\n",
      " Q_2 = 7,\n",
      " Q_6 = 1,\n",
      " Q_4 = 4,\n",
      " Q_7 = 3,\n",
      " Q_1 = 5]\n"
     ]
    }
   ],
   "source": [
    "# We know each queen must be in a different row.\n",
    "# So, we represent each queen by a single integer: the column position\n",
    "Q = [ Int(f\"Q_{row + 1}\") for row in range(8) ]\n",
    "\n",
    "# Each queen is in a column {1, ... 8 }\n",
    "val_c = [ And(1 <= Q[row], Q[row] <= 8) for row in range(8) ]\n",
    "\n",
    "# At most one queen per column\n",
    "col_c = [ Distinct(Q) ]\n",
    "\n",
    "# Diagonal constraint\n",
    "diag_c = [ If(i == j,\n",
    "              True,\n",
    "              And(Q[i] - Q[j] != i - j, Q[i] - Q[j] != j - i))\n",
    "           for i in range(8) for j in range(i) ]\n",
    "\n",
    "solve(val_c + col_c + diag_c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f87f2b8e6d0>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAECCAYAAADXWsr9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAALjUlEQVR4nO3d34tc9RnH8c+n68Y1MRLQVKwbTIU2IFKTENJKQNrEH7FKetOLBBQqLelFK4YWRHtT/AfEXhRBolYwRjQaKNIagz8QoY1N4qrRjUFDxG3UjYrEHzTR+PRiTsJ2Xbtns+d7dnae9wuGzM5OzvPsLp/5njlzZh5HhAD0tm/NdAMAyiPoQAIEHUiAoAMJEHQgAYIOJNAVQbe91vYbtt+0fVvhWvfZHrW9r2SdMfUW2X7W9rDt12zfUrjegO0Xbb9c1bujZL2qZp/tl2w/UbpWVe+Q7VdtD9neXbjWAtvbbO+v/oaXF6y1pPqZTl6O2t7UyMYjYkYvkvokvSXpYklzJL0s6ZKC9a6QtFzSvpZ+vgskLa+uz5d0oPDPZ0lnV9f7Je2S9KPCP+PvJD0k6YmWfqeHJJ3XUq0HJP2quj5H0oKW6vZJek/SRU1srxtW9JWS3oyIgxFxXNLDkn5WqlhEPC/po1Lbn6DeuxGxt7r+iaRhSRcWrBcR8Wn1ZX91KXZWlO1BSddJ2lyqxkyxfY46C8O9khQRxyPi45bKr5H0VkS83cTGuiHoF0p6Z8zXIyoYhJlke7GkZeqssiXr9NkekjQqaWdElKx3l6RbJX1VsMZ4Iekp23tsbyxY52JJRyTdXz012Wx7XsF6Y62XtLWpjXVD0D3BbT13Xq7tsyU9JmlTRBwtWSsiTkTEUkmDklbavrREHdvXSxqNiD0ltv9/rIqI5ZKulfQb21cUqnOGOk/z7o6IZZI+k1T0GJIk2Z4jaZ2kR5vaZjcEfUTSojFfD0o6PEO9FGG7X52Qb4mIx9uqW+1mPidpbaESqySts31Inadcq20/WKjWKRFxuPp3VNJ2dZ7+lTAiaWTMHtE2dYJf2rWS9kbE+01tsBuC/i9J37P93eqRbL2kv85wT42xbXWe4w1HxJ0t1Ftoe0F1/SxJV0raX6JWRNweEYMRsVidv9szEXFDiVon2Z5ne/7J65KullTkFZSIeE/SO7aXVDetkfR6iVrjbFCDu+1SZ9dkRkXEl7Z/K2mHOkca74uI10rVs71V0o8lnWd7RNIfI+LeUvXUWfVulPRq9bxZkv4QEX8rVO8CSQ/Y7lPngfyRiGjlZa+WnC9pe+fxU2dIeiginixY72ZJW6pF6KCkmwrWku25kq6S9OtGt1sdygfQw7ph1x1AYQQdSICgAwkQdCABgg4k0FVBL3w644zVoh71ZrpeVwVdUpu/zFb/cNSj3kzW67agAyigyAkzc3xmDGjqb/L5QsfUrzMb72ema0233vd/8PmU/8+RD09o4bl9p1XvwCtzp/x/ZtPvs5fr/Uef6Xgc+9obxYqcAjugefqh15TYdEo7dgxNfqcGXfOdpa3WQ3N2xdMT3s6uO5AAQQcSIOhAAgQdSICgAwkQdCABgg4kQNCBBGoFvc2RSQCaN2nQqw8Z/LM6H0F7iaQNti8p3RiA5tRZ0VsdmQSgeXWCnmZkEtCr6ryppdbIpOqN8hslaUBTf/cTgHLqrOi1RiZFxD0RsSIiVrT5dj4Ak6sT9J4emQRkMOmue9sjkwA0r9YHT1RzwkrNCgNQGGfGAQkQdCABgg4kQNCBBAg6kABBBxIg6EACBB1IoMiklrbtONzbk0yYnILpYkUHEiDoQAIEHUiAoAMJEHQgAYIOJEDQgQQIOpAAQQcSIOhAAnVGMt1ne9T2vjYaAtC8Oiv6XyStLdwHgIImDXpEPC/poxZ6AVAIz9GBBBp7myqz14Du1diKzuw1oHux6w4kUOflta2S/iFpie0R278s3xaAJtUZsrihjUYAlMOuO5AAQQcSIOhAAgQdSICgAwkQdCABgg4kQNCBBHpi9hqzydDN2pwNuPKazye8nRUdSICgAwkQdCABgg4kQNCBBAg6kABBBxIg6EACBB1IgKADCdT5cMhFtp+1PWz7Ndu3tNEYgObUOdf9S0m/j4i9tudL2mN7Z0S8Xrg3AA2pM3vt3YjYW13/RNKwpAtLNwagOVN6jm57saRlknaVaAZAGbXfpmr7bEmPSdoUEUcn+D6z14AuVWtFt92vTsi3RMTjE92H2WtA96pz1N2S7pU0HBF3lm8JQNPqrOirJN0oabXtoery08J9AWhQndlrL0hyC70AKIQz44AECDqQAEEHEiDoQAIEHUiAoAMJEHQgAYIOJNATs9cwu7U5m0xqf1Zfm/UOxIcT3s6KDiRA0IEECDqQAEEHEiDoQAIEHUiAoAMJEHQgAYIOJEDQgQTqfArsgO0Xbb9czV67o43GADSnzrnuxyStjohPq893f8H23yPin4V7A9CQOp8CG5I+rb7sry5RsikAzao7qaXP9pCkUUk7I4LZa8AsUivoEXEiIpZKGpS00val4+9je6Pt3bZ3f6FjTfcJYBqmdNQ9Ij6W9JyktRN8j9lrQJeqc9R9oe0F1fWzJF0paX/pxgA0p85R9wskPWC7T50Hhkci4omybQFoUp2j7q9IWtZCLwAK4cw4IAGCDiRA0IEECDqQAEEHEiDoQAIEHUiAoAMJuPMu1GatuGwgXtyxqPHtfpO2Z2kB3WpXPK2j8ZHH386KDiRA0IEECDqQAEEHEiDoQAIEHUiAoAMJEHQgAYIOJEDQgQRqB70a4vCSbT4YEphlprKi3yJpuFQjAMqpO5JpUNJ1kjaXbQdACXVX9Lsk3Srpq4K9ACikzqSW6yWNRsSeSe53avbakQ9PNNYggOmrs6KvkrTO9iFJD0tabfvB8XcaO3tt4bl9DbcJYDomDXpE3B4RgxGxWNJ6Sc9ExA3FOwPQGF5HBxKoM2TxlIh4Tp2xyQBmEVZ0IAGCDiRA0IEECDqQAEEHEiDoQAIEHUiAoAMJTOmEmboOvDK3p+eh7Tg81Gq9Xv5doh2s6EACBB1IgKADCRB0IAGCDiRA0IEECDqQAEEHEiDoQAIEHUig1imw1Uc9fyLphKQvI2JFyaYANGsq57r/JCI+KNYJgGLYdQcSqBv0kPSU7T22N5ZsCEDz6u66r4qIw7a/LWmn7f0R8fzYO1QPABslaUBzG24TwHTUWtEj4nD176ik7ZJWTnCfU7PX+nVms10CmJY601Tn2Z5/8rqkqyXtK90YgObU2XU/X9J22yfv/1BEPFm0KwCNmjToEXFQ0mUt9AKgEF5eAxIg6EACBB1IgKADCRB0IAGCDiRA0IEECDqQQJHZa72u12ehMVuu97CiAwkQdCABgg4kQNCBBAg6kABBBxIg6EACBB1IgKADCRB0IIFaQbe9wPY22/ttD9u+vHRjAJpT91z3P0l6MiJ+bnuOxIQGYDaZNOi2z5F0haRfSFJEHJd0vGxbAJpUZ9f9YklHJN1v+yXbm6tBDv/D9kbbu23v/kLHGm8UwOmrE/QzJC2XdHdELJP0maTbxt+JkUxA96oT9BFJIxGxq/p6mzrBBzBLTBr0iHhP0ju2l1Q3rZH0etGuADSq7lH3myVtqY64H5R0U7mWADStVtAjYkjSisK9ACiEM+OABAg6kABBBxIg6EACBB1IgKADCRB0IAGCDiTA7DV8Ta/PQss4W44VHUiAoAMJEHQgAYIOJEDQgQQIOpAAQQcSIOhAAgQdSGDSoNteYntozOWo7U1tNAegGZOeAhsRb0haKkm2+yT9W9L2wn0BaNBUd93XSHorIt4u0QyAMqYa9PWStpZoBEA5tYNefab7OkmPfsP3mb0GdKmprOjXStobEe9P9E1mrwHdaypB3yB224FZqVbQbc+VdJWkx8u2A6CEuiOZPpd0buFeABTCmXFAAgQdSICgAwkQdCABgg4kQNCBBAg6kABBBxIg6EACjojmN2ofkXQ671k/T9IHDbfTDbWoR7226l0UEQvH31gk6KfL9u6IWNFrtahHvZmux647kABBBxLotqDf06O1qEe9Ga3XVc/RAZTRbSs6gAIIOpAAQQcSIOhAAgQdSOC/XyaSWkdXkugAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "s = Solver()\n",
    "s.add(val_c + col_c + diag_c)\n",
    "s.check()\n",
    "m = s.model()\n",
    "board = np.zeros((8,8))\n",
    "\n",
    "for i in range(8):\n",
    "    board[i, m[Q[i]].as_long() - 1] = 1\n",
    "\n",
    "plt.matshow(board)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Send More Money\n",
    "\n",
    "There is an assignment of unique digits to letters such that the equation `send + more = money` holds. Find an assignment of the digits to make this true.\n",
    "\n",
    "Here we show the longer form interface to the Solver as compared to the `solve` convenience method used above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "[m = 1, s = 1]"
      ],
      "text/plain": [
       "[m = 1, s = 1]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "digits = Ints('s e n d m o r y')\n",
    "s,e,n,d,m,o,r,y = digits\n",
    "send = Sum([10**(3-i) * d for i,d in enumerate([s,e,n,d])]) # convert digits of number to number\n",
    "more = Sum([10**(3-i) * d for i,d in enumerate([m,o,r,e])])\n",
    "money = Sum([10**(4-i) * d for i,d in enumerate([m,o,n,e,y])])\n",
    "\n",
    "solver = Solver()\n",
    "solver.add([s > 0, m > 0]) # first digits are nonzero\n",
    "\n",
    "# Constrain all digits constrained to be 0 through 9\n",
    "# FILL IN\n",
    "\n",
    "# Constrain all digits to be unique (Hint: Look at the 8-Queens constraints. Anything useful?)\n",
    "# FILL IN\n",
    "\n",
    "# Constrain send + more == money\n",
    "# FILL IN\n",
    "\n",
    "solver.check()\n",
    "solver.model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## * Robot Arm Inverse Kinematics\n",
    "\n",
    "The inverse kinematics problem is to find the angles of the joints of a robot arm in order to place the hand in some position. \n",
    "You can describe the joint positions using polynomial constraints. Pairs of joints are constrained to be a distance apart equal the the length of that piece of the robot arm. You can also enter inequality constraints (like staying out of the floor) with varying amounts of success."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos = np.array([  [ Real(f\"x_{(i,j)}\")  for j in range(2)] for i in range(3)])\n",
    "s = Solver()\n",
    "s.add( np.sum(pos[0]**2) == 1 )\n",
    "s.add( np.sum(pos[1]**2) == 1 )\n",
    "s.add( np.sum((pos[2] - pos[0])**2) == 1 )\n",
    "s.add( np.sum((pos[2] - pos[1])**2) == 1 )\n",
    "s.add( pos[2][0] == 0 )\n",
    "s.add( pos[2][1] == 1 )\n",
    "s.check()\n",
    "s.model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## * Discrete Tomography\n",
    "\n",
    "Tomographic reconstruction is an interesting class of problems. One example is x-ray tomography. The data you get out is the integrated extinction of an x-ray along a line going thorugh a patient. From the data of many such lines, one can attempt to back out an image of the x-ray absorption through out the patient. \n",
    "A simplified schematic version of this can be done in a discrete setting. We have a binary image and can collect the sums of horizontal, verticals, and diagonals. Then we ask z3 to find an assignment that satisfies these sum constraints. There may or may not be a unique solution.\n",
    "The problem is related to Picross/Nonograms puzzles.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "bin_img = np.array([\n",
    "[0,0,0,0,0,0,0],\n",
    "[0,0,1,0,1,0,0],\n",
    "[0,0,1,0,1,0,0],\n",
    "[0,0,0,0,0,0,0],\n",
    "[0,1,0,0,0,1,0],\n",
    "[0,0,1,1,1,0,0],\n",
    "[0,0,0,0,0,0,0]])\n",
    "bin_img = np.hstack([bin_img,bin_img.T])\n",
    "bin_img = np.vstack([bin_img,bin_img])\n",
    "plt.imshow(bin_img)\n",
    "plt.title(\"Original Image\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rows, cols = bin_img.shape\n",
    "bvars = [[ If(Bool(f\"x_({row},{col})\"), IntVal(1),IntVal(0)) for col in range(cols)] for row in range(rows)]\n",
    "\n",
    "s = SolverFor(\"QF_FD\")\n",
    "# for each row, column sum constraint\n",
    "s.add([Sum(bvars[row]) == IntVal(int(np.sum(bin_img[row,:]))) for row in range(rows)]) \n",
    "# for each column, row sum constraint\n",
    "s.add([Sum([bvars[row][col] for row in range(rows)]  ) == IntVal(int(np.sum(bin_img[:,col]))) for col in range(cols)])\n",
    "\n",
    "# diagonal sum constraints\n",
    "for diag in range(rows+cols): # Challenge mode: Can you do this more elegantly and clearly? \n",
    "    s.add(Sum([ bvars[row][col] for row in range(rows) for col in range(cols) if col + row == diag]) == IntVal(int(sum([ bin_img[row,col] for row in range(rows) for col in range(cols) if col + row == diag]))))\n",
    "    s.add(Sum([ bvars[row][col] for row in range(rows) for col in range(cols) if col - row == diag - rows +1]) == IntVal(int(sum([ bin_img[row,col] for row in range(rows) for col in range(cols) if col - row == diag - rows + 1]))))\n",
    "for i in range(50):  # loop through all solutions\n",
    "    res = s.check()\n",
    "    if res == sat:\n",
    "        m = s.model()\n",
    "        sol = [[m.eval(s).as_long() for s in row] for row in bvars]\n",
    "        plt.figure()\n",
    "        plt.imshow(sol)\n",
    "        plt.title(\"Reconstructed Image\")\n",
    "        s.add(Not(And([  And([v == m.eval(v) for v in row]) for row in bvars])))\n",
    "    else:\n",
    "        print(\"No more models\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interlude: A Rough Sketch of How Z3 works\n",
    "\n",
    "It is very useful to have at least a simplistic mental model of how a tool works, so that you can predict the places where it might run into difficulty, or formulate your problem in a way amenable to the solver.\n",
    "\n",
    "Z3 (and other related SMT solvers) are built around the core concept of a SAT solver, so let's discuss that first.\n",
    "\n",
    "## SAT Solving\n",
    "\n",
    "A propositional logical formula is the combination of variables using the boolean connectives And, Or, Implies, and Not. Here is a propositional formula over variables p, q, r, and w in mathematical notation.\n",
    "\n",
    "$$ (p \\lor q) \\land (q \\lor \\neg w) \\land (w \\lor r \\lor p) \\land (r \\lor \\neg w) $$\n",
    "\n",
    "If we assign values to the variables, we can evaluate the formula to determine if the whole thing evaulates to true or false. In this sense, formulas are like boolean circuits. Here for example, is a python function that evaluates a boolean formula.\n",
    "\n",
    "```python\n",
    "def foo(p,q,r,w):\n",
    "    return (p or q) and (q or not w) and (w or r or p) and (r or not w)\n",
    "```\n",
    "\n",
    "One useful question we can ask about a propositional formula is if it is satisfiable. Satisfiable means there exists an assignment of the boolean variables to true and false values that makes the entire formula evaluate to true. For a circuit, this is deductively running it backwards to find an input, something that they are physically designed not to do.\n",
    "\n",
    "Specialized solvers for this question are called SAT solvers, one of which is at the core of z3. Here is a Z3 query that just sticks to the SAT fragment of its capabilities.\n",
    "\n",
    "```python\n",
    "p,q,r,w = Bools(\"p q r w\")\n",
    "solve( And( [Or(p,q), Or(q, Not(w)), Or(w, r, p), Or(r, Not(w))  ] ))\n",
    "```\n",
    "\n",
    "Given a satisfying solution, it is quite easy and efficient to check if it is indeed a satisfying solution by merely plugging the solution. It is this property that puts the SAT problem into complexity class NP.\n",
    "\n",
    "In principle, one can answer the SAT question by a brute force search (a giant nested set of `for` loops), enumerating every possible combination of true/false the variables can take on.\n",
    "\n",
    "However, for every new variable we add the size of the space we are searching is multiplied by 2. Thus brute force becomes an unacceptable approach for a number of variables well below 100.\n",
    "\n",
    "Instead SAT solvers search the space more carefully, trying to deduce and propagate information in a manner somewhat similar to how a human solves a Sudoku puzzle.\n",
    "\n",
    "A revolution in the power of SAT solvers occurred in the 90s/00s where new techniques were found, some of the most prominent of which are backjumping, the two watched literal technique, conflict directed clause learning, and random restarts. An interesting resource to learn more https://www.youtube.com/watch?v=II2RhzwYszQ&ab_channel=SimonsInstitute\n",
    "\n",
    "There are surprising number of useful and interesting questions one can encode into this form to leverage the power of SAT solvers. One very natural application is hardware verification of logical circuits. It is often very useful to even _try_ to twist a problem into a SAT solvable form, because off the shelf solvers have become so good. SAT solvers have even been used to prove previously unproven mathematical conjectures.\n",
    "\n",
    "## SMT solving\n",
    "\n",
    "There are two layers the SMT solver is working at. There may be complex boolean logical structure amenable to a pure SAT solver, but then there is also domain specific structure about the integers, reals, uninterpeted formulas, etc.\n",
    "\n",
    "We can abstract all complex domain specific facts like the truth of $x + y \\ge 7$ as an opaque boolean variables $p$. We then use the technques of SAT solvers to find a satisfying assignment of these variables. But there are implicit constraints required by the actual interpretation of the boolean variables that may not have been apparent to just the SAT search. So we take the truth assignment and query a theory specific solver to see if that particular assignment is possible. If it is, we are done. If it is not, there is a reason it is not. We add this reason to our boolean formula and iterate.\n",
    "\n",
    "This is where the intimidating name Satisfiability Modulo Theories comes from. The solver is solving SAT problems, but only accepting solutions that also make sense in underlying domain specific theories (\"modulo-ing\" those solutions out).\n",
    "\n",
    "In actual SMT solvers, there is no reason to not use domain specific information much earlier than after a full boolean assignment. The only objective is to do things that are logically correct and as fast as possible. The SAT solving and the domain specific querying may be very tightly interwoven.\n",
    "\n",
    "It is helpful to know these things so that you can have some intuition for what Z3 might like or not like. Z3 will likes conjunction (ANDs) formulas. Every disjunction (Or) tends to add branching behavior where Z3 has to pick which one it wants to be true. Z3 has built in procedures and understanding of\n",
    "\n",
    "- propositional logic via SAT solving\n",
    "- linear equations $3x+ 4y == 7$ \n",
    "- linear inequalities  $ 3x + 4y <= 7$\n",
    "- Bitvectors, which can also be reduced to SAT by bit-blasting as a last resort.\n",
    "- Arrays\n",
    "- Algebraic Data Types\n",
    "- Uninterpeted Functions\n",
    "- polynomial equalities $x^2 + y^2 == 2$\n",
    "\n",
    "These pieces give a very powerful and surprisingly flexible base to work from. \n",
    "\n",
    "## What Z3 Can't Do\n",
    "\n",
    "Z3 seems incredibly powerful, and it is. However, you can express far more things to Z3 than it can solve. I've been cherry picking problems I know it handles well.\n",
    "\n",
    "As an example, factoring integers is at the heart of RSA cryptography. While Z3 can factor integers, it can't magically factor even moderately big numbers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[y = 7, x = 3]\n"
     ]
    }
   ],
   "source": [
    "x,y = Ints(\"x y\")\n",
    "pubkey = 3\t* 7\n",
    "solve(x * y == pubkey, x > 1, y > 1) # easy peasy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "failed to solve\n"
     ]
    }
   ],
   "source": [
    "x,y = Ints(\"x y\")\n",
    "pubkey = 1000000993\t* 1000001011\n",
    "solve(x * y == pubkey, x > 1, y > 1) # nope"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nonlinear equations in general are going to be tough. Not impossible, but tough. Z3 gives up immediately if you try to find the solution to a seemingly solvable problem involving an exponential. Z3 also can't really understand sines, cosines, and logarithms. The main exception is polynomial equality constraints, for which z3 has an intrinsic understanding, however these routines can be computationally expensive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Real('x')\n",
    "s = Solver()\n",
    "s.add(2**x == 3)\n",
    "print(s.check())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Proof = Exhaustively Not Finding Counterexamples\n",
    "\n",
    "So far we've been using z3 to find solutions to a set of constraints/formulas. Z3 is so effective at this process that it can be exhaustive and know when there are no solutions to be found. This is useful for theorem proving. If we state a theorem like `Implies(And(p,q),p)` that must be true for any value of `p` and `q`, we can prove its truth by exhaustively search for a counterexample and failng to find one. A counterexample is an assignment of the variables for which the theorem evaluates to false, or equivalently for which the negation `Not(th)` evaluates to true.\n",
    "\n",
    "`solver.check()` may return \n",
    "   - sat - \"I have found a solution. You may ask for it with `solver.model()`\"\n",
    "   - unknown - This means \"I give up\". There may be a solution out there somewhere, there may not.\n",
    "   - unsat - \"I KNOW there is no solution. I've tried EVERYTHING\".\n",
    "\n",
    "Here is a simplified (not even by that much) version of the Z3 convenience function `prove` in which we perform the negation and intepret `unsat` as a proof."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prove2(f): # a simplified version of the z3py function \"prove\"\n",
    "    s = Solver()\n",
    "    s.add(Not(f))\n",
    "    if s.check() == unsat:\n",
    "        print(\"proved\")\n",
    "    else:\n",
    "        print(\"failed to prove\")     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Propositional Logic / Boolean Logic\n",
    "\n",
    "Propositional formula are logical formula with boolean variables and logical connectives like and ($\\land$) , or ($\\lor$),  negation ( $ \\neg $ ), and implication $\\implies$. Conceptually, z3 checks all possible assingments of true and false to the variables. If the given formula always evaluates to true under any assignment, z3 states `proved`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A refresher on basic boolean operations\n",
    "p, q = Bools('p q')\n",
    "print(And(p,q))\n",
    "print(Or(p,q))\n",
    "print(Xor(p,q))\n",
    "print(Not(p))\n",
    "print(Implies(p,q))\n",
    "print(p == q)\n",
    "\n",
    "# a reminder from earlier\n",
    "my_true_thm = Implies(And(p,q), p)\n",
    "prove(my_true_thm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "Prove:\n",
    "- De Morgan's Law `p & q == ~ (~p | ~q)`\n",
    "- `p -> q ==  ~ p | q`\n",
    "- Peirce's Law `((p -> q) -> p) -> p` is always true\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FILL IN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proving properties of arithmetic\n",
    "\n",
    "Z3 is capable of understanding many of the properties that they taught you in arithmetic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise : Properties of Ints, Reals, and Floats\n",
    "Prove or find counterexamples that the mathematical reals and 16 bit floats\n",
    "- Addition commutes `x + y == y + x`\n",
    "- Addition is associative `(x + y) + z == x + (y + z)`\n",
    "- `x` is either less than, equal to, or greater than 0\n",
    "- `x*x` is always non negative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solution\n",
    "x = Real(\"x\")\n",
    "y = Real(\"y\")\n",
    "z = Real(\"z\")\n",
    "\n",
    "# Prove Properties here\n",
    "# FILL IN\n",
    "\n",
    "x = FP(\"x\", Float16())\n",
    "y = FP(\"y\", Float16())\n",
    "z = FP(\"z\", Float16())\n",
    "\n",
    "# FILL IN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "            \n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## * Properties of Vectors\n",
    "\n",
    "I find a sometimes useful trick is to embed python objects into numpy arrays. When you do this, you get numpy array slicing and broadcasting, vectorized addition, matrix multiplication and scalar multiplication for free. Here for example, this makes for succinct expression of some properties of linear algebra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "x,y = Reals(\"x y\")\n",
    "q = np.array([x,y]) \n",
    "\n",
    "np.dot(q, q) # dot product\n",
    "np.array([x]) * q # scalar product\n",
    "q * x # scalar product\n",
    "\n",
    "np.RealVal = np.vectorize(RealVal) # simple mapping of functions over a vector\n",
    "\n",
    "import operator as op\n",
    "def vec_eq(x,y):\n",
    "    return And(np.vectorize(op.eq)(x,y).tolist())\n",
    "def NPArray(n, prefix=None, dtype=RealSort()): # a vector of fresh variables.\n",
    "    return np.array( [FreshConst(dtype, prefix=prefix) for i in range(n)] )\n",
    "\n",
    "\n",
    "u = NPArray(2)  \n",
    "v = NPArray(2)\n",
    "w = NPArray(2)\n",
    "\n",
    "prove(vec_eq(   v + u       ,  u + v            )) #commutativity\n",
    "prove(vec_eq(   v + (u + w) , (v + u) + w       )) # associativity\n",
    "prove(vec_eq(  (v + w) * x  , (w * x) + (v * x) ))# scalar distrubitivity\n",
    "prove(np.dot(v, u)**2 <= np.dot(v,v) * np.dot(u,u)) #cauchy schwartz inequality\n",
    "\n",
    "\n",
    "# proving properties of a linear operator defined by a matrix\n",
    "\n",
    "def Z(x):\n",
    "    z = np.RealVal( np.arange(4).reshape(2,2) )\n",
    "    return z @ x\n",
    "\n",
    "def linear(Z):\n",
    "    x = NPArray(2)\n",
    "    y = NPArray(2)\n",
    "    return vec_eq(Z(x + y), Z(x) + Z(y))\n",
    "\n",
    "\n",
    "p = linear(Z)\n",
    "prove(p)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## * Geometric Theorem Proving\n",
    "\n",
    "You can convert geometric theorems into systems of polynomial equations.\n",
    "\n",
    "Consider the following. You have a circle with a diameter. Pick an arbitrary point on the circle and connect it with lines to the points on the diameter. Is the angle formed by these lines a perpendicular angle?\n",
    "\n",
    "Two useful definitions:\n",
    "\n",
    "- A point can be constrained to a circle by `x**2 + y**2 == R**2`. \n",
    "- Perpendicularity can be defined as the dot product of the vectors going along the two lines is 0. \n",
    "\n",
    "![my_circle](images/circle_diagram.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = Reals(\"x y\") # point on circle\n",
    "# assume wlog horizontal diameter y = 0. Hits circle at x = 1, x = -1\n",
    "d1 = (x - 1, y) # vector along line 1\n",
    "d2 = (x + 1, y) # vector along line 2\n",
    "perp = d1[0]*d2[0]+d1[1]*d2[1] == 0 #perpendicularity constraint\n",
    "prove(Implies(x**2 + y**2 == 1, perp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## * Axiomatic Groups\n",
    "\n",
    "A [Group](https://en.wikipedia.org/wiki/Group_(mathematics\\)) is a mathematical abstraction of the properties of symmettries. It is a set with a binary operation that is associative and an inverse operation. Common examples include integers, matrices, and polynomials.\n",
    "Z3 has sufficient expressive power to directly express the axiomatic formulation of a group. Z3's solving capabilities become less reliable in the presence of quantifiers so your mileage may vary.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<b>unsat</b>"
      ],
      "text/plain": [
       "unsat"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G = DeclareSort(\"G\")\n",
    "e,a,b,c = Consts(\"e a b c\", G)\n",
    "times = Function('*', G,G, G)\n",
    "inv = Function('~', G,G)\n",
    "\n",
    "axioms = [\n",
    "  ForAll([a], times(inv(a),a) == e ),\n",
    "  ForAll([a], times(a,inv(a)) == e ),\n",
    "  ForAll([a], times(e,a) == a ),\n",
    "  ForAll([a], times(a,e) == a ),\n",
    "  ForAll([a,b,c] , times(a,times(b,c)) == times(times(a,b),c))\n",
    "]\n",
    "\n",
    "# A simple theorem\n",
    "s = Solver()\n",
    "s.add(axioms)\n",
    "theorem = inv(times(a,b)) == times(inv(b), inv(a))\n",
    "s.add(Not(theorem))\n",
    "s.check()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## * An Inductive Fact About a Sum\n",
    "\n",
    "Z3 can crush finite structures, but isn't particularly designed for proving facts about the infinite. Questions like these need some kind of inductive reaosning.\n",
    "Here is an example trying to prove a closed form for the sum $\\sum_{i=0}^n i = \\frac{n(n+1)}{2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = Int(\"n\")\n",
    "Sumn = Function(\"sumn\", IntSort(), IntSort())\n",
    "s = Solver()\n",
    "s.add( Sumn(0) == 0)\n",
    "s.add(ForAll([n], Sumn(n+1) == n + 1 + Sumn(n)))\n",
    "s.add(Not(ForAll([n], Implies(n >= 0, 2 * Sumn(n) == n * (n + 1)))))\n",
    "s.check()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By maually instantiating an induction principle, we can get this to go through forever, however this process is not at all fool proof. Ultimately, it appears that complicated questions like these are best treated in interactive proof assistants like Coq."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = Int(\"n\")\n",
    "Sumn = Function(\"sumn\", IntSort(), IntSort())\n",
    "s = Solver()\n",
    "# defining sum(n) = 1 + 2 + 3 + ... n\n",
    "s.add( Sumn(0) == 0)\n",
    "s.add(ForAll([n], Sumn(n+1) == n + 1 + Sumn(n)))\n",
    "def prop_sumn(n):\n",
    "    return 2 * Sumn(n) == n * (n + 1)\n",
    "def induction(p):\n",
    "    n = Int(\"n\")\n",
    "    return Implies(And(  p(0),    ForAll([n],Implies(And(n >= 0, p(n)),p(n+1)))),\n",
    "                     # --------------------------------------------------------\n",
    "                       ForAll( [n] , Implies(n>=0, p(n))))\n",
    "s.add(induction(prop_sumn))\n",
    "s.add( Not(ForAll([n], Implies( n >= 0, prop_sumn(n)))))\n",
    "if s.check() == unsat:\n",
    "    print(\"proved\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Proving Programs\n",
    "\n",
    "So far we've been dealing with things that seem like puzzles or pure math. Now we have enough tools in our toolbox to attack problems that feel more directly related to programming.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bit Tricks\n",
    "\n",
    "Some low level hacks are very popular with C programmers. Here is a list of cool bit tricks\n",
    "http://graphics.stanford.edu/~seander/bithacks.html\n",
    "\n",
    "### Power of two\n",
    "This hack is frequently used in C programs (Z3 included) to test whether a machine integer is a power of two. We can use Z3 to prove it really works. To understand this hack, it helps to look at the boolean form of these numbers. It seems that x - 1 and x share no 1s if and only if x is a power of 2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0010000000000000\n",
      "0001111111111111\n"
     ]
    }
   ],
   "source": [
    "print(format(2**13, '016b'))\n",
    "print(format(2**13-1, '016b'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0000000000101010\n",
      "0000000000101001\n"
     ]
    }
   ],
   "source": [
    "# A non power of two\n",
    "print(format(42, '016b'))\n",
    "print(format(42-1, '016b'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True]\n",
      "[]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def is_pow2_slow(x):\n",
    "    return x in [2 ** i for i in range(32)]\n",
    "def is_pow2_fast(x):\n",
    "    return not x & (x - 1)\n",
    "\n",
    "# Tests\n",
    "print([is_pow2_fast(2**i) for i in range(32)])\n",
    "print([ i for i in np.random.randint(2**32, size=100000) if is_pow2_fast(i) != is_pow2_slow(i)  ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "counterexample\n",
      "[x = 0]\n"
     ]
    }
   ],
   "source": [
    "x      = BitVec('x', 32)\n",
    "fast   = x & (x - 1) == 0\n",
    "powers = [ 2**i for i in range(32) ]\n",
    "slow   = Or([ x == p for p in powers ])\n",
    "prove(fast == slow)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Really? But my spiel sounded so good."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "print(is_pow2_slow(0))\n",
    "print(is_pow2_fast(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes, this is a counter example. Some funny behavior happens subtracting at 0 when looking at the bit representation. What about it we exclude that case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "proved\n"
     ]
    }
   ],
   "source": [
    "def is_pow2_fast(x):\n",
    "    return x != 0 and not x & (x - 1)\n",
    "\n",
    "fast = And(x != 0, x & (x - 1) == 0)\n",
    "prove(fast == slow)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Exactly One Negative\n",
    "The following bit trick gives a fast way to determine if exactly one of 2 given integers is negative. The reason it works is because in the two's comnplement encoding of negative numbers the most significant bit is a [sign bit](https://en.wikipedia.org/wiki/Sign_bit)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(-2, 0), (-2, 1), (-1, 0), (-1, 1), (0, -2), (0, -1), (1, -2), (1, -1)]\n",
      "[(-2, 0), (-2, 1), (-1, 0), (-1, 1), (0, -2), (0, -1), (1, -2), (1, -1)]\n"
     ]
    }
   ],
   "source": [
    "def opposite_signs_fast(x,y):\n",
    "    return x ^ y < 0\n",
    "def opposite_signs_slow(x,y):\n",
    "    if x < 0 and y >= 0:\n",
    "        return True\n",
    "    if x >= 0 and y < 0:\n",
    "        return True\n",
    "    return False\n",
    "\n",
    "print( [ (i,j) for i in range(-2,2) for j in range(-2,2) if opposite_signs_fast(i,j) ])\n",
    "print( [ (i,j) for i in range(-2,2) for j in range(-2,2) if opposite_signs_slow(i,j) ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = BitVec('x', 32)\n",
    "y = BitVec('y', 32)\n",
    "\n",
    "# FILL IN\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Symbolic Execution and Bounded Model Checking\n",
    "It is fairly straightforward with a few tricks to directly translate a program into a Z3 expression. Depending on the community, goals, and differences in search strategy this technique goes by the name of symbolic execution or bounded model checking.\n",
    "\n",
    "Here are some real world projects that are interesting to check out\n",
    "\n",
    "- https://klee.github.io/ Klee Symbolic Execution Engine\n",
    "- http://www.esbmc.org/ ESMBC \n",
    "- https://apalache.informal.systems/ Apalache - Bounded Model Checking for TLA+\n",
    "\n",
    "### Assignment and Equality\n",
    "\n",
    "Consider this simple program"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def foo(x):\n",
    "    y = 2*x + 4\n",
    "    z = y*y\n",
    "    assert(z >= 0)\n",
    "    return z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can translate this to Z3 with very few changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "x,y,z = Ints(\"x y z\")\n",
    "prog = [\n",
    "  y == 2 * x + 4,\n",
    "  z == y * y\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can ask Z3 if this always returns a positive number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = Solver()\n",
    "s.add(prog)\n",
    "s.add( Not(z >= 0) )\n",
    "s.check()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider this program"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def foo(x):\n",
    "    x = 2*x + 4\n",
    "    x = x*x\n",
    "    assert(x >= 0)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is not an accurate translation of this program. Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Int(\"x\")\n",
    "prog = [\n",
    "  x == 2 * x + 4,\n",
    "  x == x * x\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The issue is subtle differences between incompatible notions of equality. This is not just a Z3 issue.\n",
    "Consider the difference between these two statements.\n",
    "\n",
    "```\n",
    "x = 3\n",
    "x == 3\n",
    "```\n",
    "The first assigns `x` to the number 3. The second checks if `x` holds the number 3. Both could be called \"equals\" but are very different. You could imagine a beginner programmer being confused on this point.\n",
    "\n",
    "Z3 is static. It has no intrinsic notion of time or mutation. Its expressions are like those taught to you in algebra class rather than in your programming class.\n",
    "\n",
    "When we're modelling a programming language with mutation in Z3, we need to translate the mutational `=` to Z3's notion of `==`. One way to do so is to make every new assignment go to a fresh variable.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def foo(x):\n",
    "    x1 = 2*x + 4\n",
    "    x2 = x1*x1\n",
    "    return x2\n",
    "\n",
    "x,x1,x2 = Ints(\"x x1 x2\")\n",
    "prog = [\n",
    "  x1 == 2 * x + 4,\n",
    "  x2 == x1 * x1\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This transformation is commonly used in compiler's where it is known as Single Static Assignment https://en.wikipedia.org/wiki/Static_single_assignment_form . The sameness of `x` at different points in your code is an abstraction. Your language may not be guarantee that variables named `x` really are the same location in memory at different points in your function language. They may be flipping things in and out of registers and memory as the compiler sees fit.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Branches\n",
    "Describing branching which is a dynamic process to the static Z3 also takes a little ingenuity. Now the value a variable takes on depends on which branch control came from"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def foo(x):\n",
    "    if x <= 0:\n",
    "        y = x\n",
    "    else:\n",
    "        y = -x\n",
    "    z = y*y\n",
    "    assert(z >= 0)\n",
    "    return z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x,y,z = Ints(\"x y z\")\n",
    "prog = [\n",
    "    If(x <= 0, y == x, y == -x),\n",
    "    z == y * y\n",
    "]\n",
    "prove( Implies(And(prog), z >= 0) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise: Two Sort\n",
    "\n",
    "Consider the following program which sorts x and y. This is the simplest possible sorting algorithm that deserves the name.\n",
    "Translate this to z3 and prove that the result of this program `(a,b)` has the property `a <= b`. Does this property specify the function is a sorting algorithm?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sort(x,y):\n",
    "    if y < x:\n",
    "        temp = x\n",
    "        x = y\n",
    "        y = temp\n",
    "    else:\n",
    "        x = x\n",
    "        y = y\n",
    "    assert(x <= y)\n",
    "    return (x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FILL IN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## * Loops\n",
    "\n",
    "The simplest and most effective way to deal with looping behavior is to simply unroll the loops. This is simple when the number of executions are known. The following two programs are equivalent.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def foo(x):\n",
    "    for i in range(3):\n",
    "        x = x*x\n",
    "    return x\n",
    "\n",
    "def foo(x):\n",
    "    x = x*x\n",
    "    x = x*x\n",
    "    x = x*x\n",
    "    return x\n",
    "\n",
    "x,x1,x2,x3 = Ints(\"x x1 x2 x3\")\n",
    "prog = [\n",
    "    x1 == x*x,\n",
    "    x2 == x1*x1,\n",
    "    x3 == x2*x2,\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the bound on the loop isn't known (which is very common) a strategy is to just unroll the loop a bounded amount. Counterexamples returned can be confirmed as counterexamples. Sometimes false counterexamples may be returned, in which case one may increase the number of loop executions. This is the strategy from which Bounded Model Checking receives it's name."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### * Verifying the Babylonian Square Root Method\n",
    "\n",
    "The babylonian method is an iterative procedure for computing square roots. You take an initial guess and then compute the next iteration using the formula\n",
    "$ x_{n+1} = \\frac{1}{2}(\\frac{r}{x_n} + x_n)$\n",
    "Intuitively, if your current $x$ is low for the square root of $r$, $\\frac{r}{x}$ is high. If $x$ is high, $\\frac{r}{x}$ is low. So if you take the average of the two, you're hopefully getting closer to the real answer.\n",
    "We can also see that if you are precisely at the square root $x = \\sqrt{r}$, (r/x + x)/2 is a fixed point of the iteration. $\\frac{1}{2}(\\frac{r}{x} + x) = \\frac{1}{2}(\\frac{r}{\\sqrt{r}} + \\sqrt{r}) = \\sqrt{r}$\n",
    "\n",
    "How good is the method after 7 iterations though? What precision can we rely on? Well, we can ask z3 to prove a bound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "proved\n"
     ]
    }
   ],
   "source": [
    "def babylonian(x):\n",
    "    res = 1\n",
    "    for i in range(7):\n",
    "        res = (x / res + res) / 2\n",
    "    return res    \n",
    "\n",
    "def bab_iter(x, y):\n",
    "    return (x / y + y) / 2\n",
    "\n",
    "x, y = Reals(\"x y\")\n",
    "prove(Implies(And(y**2 == x, y >= 0, 0 <= x, x <= 10), babylonian(x) - y <= 0.01))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAXFklEQVR4nO3df5Dcd13H8efrdvcugIMp0yuGJDURDjSDWOKRRhBFfmgSmR6gaDJCamQ8gw0Co0KQmarj6DD8EIhmchMhSsZKRCxyg6ehU3/wj8FcIYSmbexNLORIaI9xDIyZ9Hq5t3/sd6/bzd7uJ7m9214+r8fMzu1+v5/P7uczydzrPp/P97sfRQRmZpafnm43wMzMusMBYGaWKQeAmVmmHABmZplyAJiZZarc7QZcjRtvvDHWrVvX7WaYmS0r991333cior/x+LIKgHXr1jE+Pt7tZpiZLSuSvtHseNIUkKQtkk5LmpC0t8l5SdpXnD8paWO7upJukXRM0glJ45I2XUvHzMzs2rQNAEklYD+wFdgA7JC0oaHYVmCgeAwDBxLqfhD4w4i4BbizeG1mZkskZQSwCZiIiDMRMQ0cAYYaygwBh6PqGLBS0qo2dQN4dvH8+4FzC+yLmZldhZQ1gNXA2brXk8CtCWVWt6n7LuCopA9TDaKXN/twScNURxXcfPPNCc01M7MUKSMANTnW+AVC85VpVfftwLsjYi3wbuCTzT48Ig5GxGBEDPb3X7GIbWZm1yglACaBtXWv13DldM18ZVrVvR24u3j+d1Sni8zMbImkBMBxYEDSekm9wHZgtKHMKLCzuBpoM3AhIs63qXsO+Oni+auBhxfYFzMzuwpt1wAiYkbSHuAoUAIORcQpSbuL8yPAGLANmAAuArta1S3e+teBj0sqA5co5vkXw70PPsrpR7/Hb77qBYv1EWZmy46W034Ag4ODcS03gt35+fsZ/do5Ttz5s4vQKjOzpzdJ90XEYOPxLL4LqK/cw+NPzHa7GWZmTyuZBECJx2cus5xGO2Zmiy2LAFhR6WE2YGbWAWBmVpNFAPSVSwA8PuNpIDOzmjwCoFLt5uNPXO5yS8zMnj7yCIByEQAeAZiZzckkADwFZGbWKJMAqI0APAVkZlaTRwDMrQF4BGBmVpNHAHgKyMzsCpkEgKeAzMwaZRIAxQjAU0BmZnPyCICKLwM1M2uURwB4CsjM7AqZBIAXgc3MGmUSAP4qCDOzRnkEQLEGcMkjADOzOUkBIGmLpNOSJiTtbXJekvYV509K2tiurqS/lXSieDwi6URnunSl3pJvBDMza9R2T2BJJWA/8DpgEjguaTQiHqgrthUYKB63AgeAW1vVjYhfrvuMjwAXOtSnK5RLPZR75EVgM7M6KSOATcBERJyJiGngCDDUUGYIOBxVx4CVklal1JUk4JeATy+wLy31lXu8CGxmViclAFYDZ+teTxbHUsqk1H0l8GhEPNzswyUNSxqXND41NZXQ3Ob6KiWPAMzM6qQEgJoca9xbcb4yKXV30OKv/4g4GBGDETHY39/fsqGteGN4M7OnarsGQPWv9rV1r9cA5xLL9LaqK6kMvAn48fQmXxtPAZmZPVXKCOA4MCBpvaReYDsw2lBmFNhZXA20GbgQEecT6r4WeCgiJhfckzb6yp4CMjOr13YEEBEzkvYAR4EScCgiTknaXZwfAcaAbcAEcBHY1apu3dtvZ5EXf2v6Kh4BmJnVS5kCIiLGqP6Srz82Uvc8gDtS69ad+9XUhi7UinLJawBmZnWyuBMYaiMATwGZmdXkEwBeBDYze4qMAqDkADAzq5NRAHgKyMysXj4BUPGNYGZm9fIJAE8BmZk9RUYB4CkgM7N6mQXALNVbFszMLJ8AqJSIgCcuOwDMzCCnAKjtC+xpIDMzIMsA8EKwmRlkFQAlwAFgZlaTTwBUahvDewrIzAxyCoBiCuiSbwYzMwOyCoDaFJBHAGZmkFUAeBHYzKxeUgBI2iLptKQJSXubnJekfcX5k5I2ptSV9I7i3ClJH1x4d+Y3twbgADAzAxJ2BJNUAvYDr6O6+ftxSaMR8UBdsa3AQPG4FTgA3NqqrqSfAYaAl0TE45Ju6mTHGs1NAXkR2MwMSBsBbAImIuJMREwDR6j+4q43BByOqmPASkmr2tR9O/CBiHgcICIe60B/5uUpIDOzp0oJgNXA2brXk8WxlDKt6r4QeKWkL0v6d0kva/bhkoYljUsan5qaSmhuc74PwMzsqVICQE2ONX6hznxlWtUtAzcAm4HfBT4j6YryEXEwIgYjYrC/vz+huc09uQbgKSAzM0hYA6D6V/vautdrgHOJZXpb1J0E7o7q13P+p6RZ4Ebg2v/Mb2FuCsj3AZiZAWkjgOPAgKT1knqB7cBoQ5lRYGdxNdBm4EJEnG9T9x+AVwNIeiHVsPjOgns0jxUVTwGZmdVrOwKIiBlJe4CjQAk4FBGnJO0uzo8AY8A2YAK4COxqVbd460PAIUn3A9PA7bGIX9bfW/IUkJlZvZQpICJijOov+fpjI3XPA7gjtW5xfBp4y9U0diF6ekRvqccjADOzQjZ3AkOxK5jXAMzMgNwCoOJ9gc3MavIKgHLJU0BmZoXMAsBrAGZmNVkFQG+5x98FZGZWyCoA+iqeAjIzq8krAMpeBDYzq8kwADwCMDOD7AKg5PsAzMwKeQWA7wMwM5uTVwB4CsjMbE5mAeCrgMzMajILgB4u+T4AMzMgtwCoeArIzKwmrwAol5iemWURtx0wM1s2MguA2qYwHgWYmSUFgKQtkk5LmpC0t8l5SdpXnD8paWO7upL+QNK3JJ0oHts606X5OQDMzJ7UNgAklYD9wFZgA7BD0oaGYluBgeIxDBxIrPvRiLileFyxa1in9c3tC+yFYDOzlBHAJmAiIs4U2zgeAYYaygwBh6PqGLBS0qrEuktmbgTgu4HNzJICYDVwtu71ZHEspUy7unuKKaNDkm5o9uGShiWNSxqfmppKaO78PAVkZvaklABQk2ONl9HMV6ZV3QPA84FbgPPAR5p9eEQcjIjBiBjs7+9PaO78VngKyMxsTjmhzCSwtu71GuBcYpne+epGxKO1g5L+AvhCcquvkUcAZmZPShkBHAcGJK2X1AtsB0YbyowCO4urgTYDFyLifKu6xRpBzRuB+xfYl7b6ysUIwGsAZmbtRwARMSNpD3AUKAGHIuKUpN3F+RFgDNgGTAAXgV2t6hZv/UFJt1CdEnoE+I1OdqyZvkptBOApIDOzlCkgiks0xxqOjdQ9D+CO1LrF8bdeVUs7wFNAZmZPyuxO4NoisAPAzCyzAKjdB+ApIDOzvAKg4ikgM7OavALAU0BmZnMyCwBfBWRmVpNnAPg+ADOzvAJAEr3eGN7MDMgsAKA6CvAUkJlZlgFQ8gjAzIwsA6DHawBmZuQYAJUeLnkKyMwswwAolzwCMDMjywDwIrCZGWQbAB4BmJnlFwAVXwVkZgY5BkC5x98GamZGpgEw7RGAmVlaAEjaIum0pAlJe5ucl6R9xfmTkjZeRd3fkRSSblxYV9L4RjAzs6q2ASCpBOwHtgIbgB2SNjQU2woMFI9h4EBKXUlrgdcB31xwTxKtqPgqIDMzSBsBbAImIuJMREwDR4ChhjJDwOGoOgaslLQqoe5HgfdQ3Rh+Sfg+ADOzqpQAWA2crXs9WRxLKTNvXUm3Ad+KiK+1+nBJw5LGJY1PTU0lNLe1voovAzUzg7QAUJNjjX+xz1em6XFJzwTeD9zZ7sMj4mBEDEbEYH9/f9vGttNX7mH68iyzs0s26DAze1pKCYBJYG3d6zXAucQy8x1/PrAe+JqkR4rjX5H0A1fT+GtR2xZy+rJHAWaWt5QAOA4MSFovqRfYDow2lBkFdhZXA20GLkTE+fnqRsTXI+KmiFgXEeuoBsXGiPh2pzo2H+8KZmZWVW5XICJmJO0BjgIl4FBEnJK0uzg/AowB24AJ4CKwq1XdRelJor5K/b7AlW42xcysq9oGAEBEjFH9JV9/bKTueQB3pNZtUmZdSjs6oTYF5IVgM8tdlncCA74XwMyyl20AXPIagJllLr8AqHgKyMwMcgwATwGZmQFZB4BHAGaWtwwDoJgC8hqAmWUuvwCoeArIzAxyDABPAZmZAVkGQG0KyCMAM8tbfgFQ8QjAzAxyDABPAZmZARkGQG+p9m2gngIys7xlFwCS6Ct7VzAzs+wCAHAAmJmRawBUSr4PwMyyl2UArKj0+E5gM8teUgBI2iLptKQJSXubnJekfcX5k5I2tqsr6Y+KsickfVHS8zrTpfb6yiVPAZlZ9toGgKQSsB/YCmwAdkja0FBsKzBQPIaBAwl1PxQRL4mIW4AvAHcuvDtpqmsAngIys7yljAA2ARMRcSYipoEjwFBDmSHgcFQdA1ZKWtWqbkR8t67+s4BYYF+SeRHYzCwtAFYDZ+teTxbHUsq0rCvpjyWdBX6FeUYAkoYljUsan5qaSmhue33lktcAzCx7KQGgJsca/1qfr0zLuhHx/ohYC9wF7Gn24RFxMCIGI2Kwv78/obnt9VU8BWRmlhIAk8DautdrgHOJZVLqAvwN8AsJbekITwGZmaUFwHFgQNJ6Sb3AdmC0ocwosLO4GmgzcCEizreqK2mgrv5twEML7EsyXwVkZgbldgUiYkbSHuAoUAIORcQpSbuL8yPAGLANmAAuArta1S3e+gOSXgTMAt8Adne0Zy30lXv8XUBmlr22AQAQEWNUf8nXHxupex7AHal1i+NLNuXTqLoG4BGAmeUtyzuBPQVkZpZtAPgqIDOzTAOgxBOXg8uzS3bvmZnZ006eAVBsCzntaSAzy1ieATC3LaSngcwsX5kGQAnwvsBmlrdMA6Da7Uu+F8DMMpZnAFRqU0AeAZhZvvIMgNoUkL8R1MwylmkAeBHYzCzzAPAIwMzylWcAVGpXAXkEYGb5yjMAaiMArwGYWcayDIAVFd8HYGaWZQB4EdjMLPsA8AjAzPKVFACStkg6LWlC0t4m5yVpX3H+pKSN7epK+pCkh4ryn5O0sjNdam9uEdhrAGaWsbYBIKkE7Ae2AhuAHZI2NBTbCgwUj2HgQELde4AXR8RLgP8C3rfg3iTyFJCZWdoIYBMwERFnImIaOAIMNZQZAg5H1TFgpaRVrepGxBcjYqaofwxY04H+JCn3iB55CsjM8pYSAKuBs3WvJ4tjKWVS6gL8GvBPzT5c0rCkcUnjU1NTCc1tT5K3hTSz7KUEgJoca9xKa74ybetKej8wA9zV7MMj4mBEDEbEYH9/f0Jz0/RVenjc3wZqZhkrJ5SZBNbWvV4DnEss09uqrqTbgdcDr4mIJd2fsbovsEcAZpavlBHAcWBA0npJvcB2YLShzCiws7gaaDNwISLOt6oraQvwXuC2iLjYof4k8xSQmeWu7QggImYk7QGOAiXgUESckrS7OD8CjAHbgAngIrCrVd3irf8c6APukQRwLCJ2d7JzrVRHAJ4CMrN8pUwBERFjVH/J1x8bqXsewB2pdYvjL7iqlnZYdQ3AIwAzy1eWdwKDp4DMzDIOAE8BmVneMg8AjwDMLF8ZB0DJawBmlrV8A6DiKSAzy1u+AVDu4ZJHAGaWsYwDoOQRgJllLeMA8CKwmeUt3wCoOADMLG/5BkC5xOXZYOayQ8DM8pRxAHhfYDPLmwPAAWBmmco2AFbUNob3lUBmlqlsA6CvUowAfC+AmWUq3wAo10YADgAzy1PGAVBbA/AUkJnlKSkAJG2RdFrShKS9Tc5L0r7i/ElJG9vVlfRmSackzUoa7Ex30nkEYGa5axsAkkrAfmArsAHYIWlDQ7GtwEDxGAYOJNS9H3gT8KWFd+PqeQ3AzHKXMgLYBExExJmImAaOAEMNZYaAw1F1DFgpaVWruhHxYESc7lhPrpKngMwsdykBsBo4W/d6sjiWUialbld4CsjMcpcSAGpyLBLLpNRt/eHSsKRxSeNTU1NXU7UljwDMLHcpATAJrK17vQY4l1gmpW5LEXEwIgYjYrC/v/9qqrbkNQAzy11KABwHBiStl9QLbAdGG8qMAjuLq4E2Axci4nxi3a7wFJCZ5a7crkBEzEjaAxwFSsChiDglaXdxfgQYA7YBE8BFYFerugCS3gj8GdAP/KOkExHxc53u4Hw8BWRmuWsbAAARMUb1l3z9sZG65wHckVq3OP454HNX09hOmgsATwGZWaayvRO4XOqh1CNPAZlZtrINAKhtC+kpIDPLkwPAIwAzy1TmAVDi0hMeAZhZnvIOAG8Mb2YZyzsAyj2+CsjMspV5AJS8CGxm2co8ADwFZGb5yjsAvAZgZhnLOwA8BWRmGcs6AFZUvAhsZvnKOgCqIwAHgJnlKfMA8FdBmFm+HAAeAZhZpvIOgErJawBmlq28A6CYAqpuZ2BmlpfsA2A2YGbWAWBm+UkKAElbJJ2WNCFpb5PzkrSvOH9S0sZ2dSU9R9I9kh4uft7QmS6l877AZpaztgEgqQTsB7YCG4AdkjY0FNsKDBSPYeBAQt29wL0RMQDcW7xeUn2V2raQvhLIzPKTsifwJmAiIs4ASDoCDAEP1JUZAg4XewMfk7RS0ipgXYu6Q8CrivqfAv4NeO8C+3NVavsC/+LIf1Du0VJ+tJnZVfmTN/0oL1v3nI6+Z0oArAbO1r2eBG5NKLO6Td3nRsR5gIg4L+mmZh8uaZjqqIKbb745obnpfnKgnzfc8jymL3sKyMye3p5RKXX8PVMCoNmfxo2rpvOVSanbUkQcBA4CDA4OdnS1dvXKZ/Cx7S/t5FuamS0bKYvAk8DautdrgHOJZVrVfbSYJqL4+Vh6s83MbKFSAuA4MCBpvaReYDsw2lBmFNhZXA20GbhQTO+0qjsK3F48vx34/AL7YmZmV6HtFFBEzEjaAxwFSsChiDglaXdxfgQYA7YBE8BFYFerusVbfwD4jKS3Ad8E3tzRnpmZWUtaTnfBDg4Oxvj4eLebYWa2rEi6LyIGG49nfSewmVnOHABmZplyAJiZZcoBYGaWqWW1CCxpCvjGNVa/EfhOB5uzHLjPeXCf87CQPv9gRPQ3HlxWAbAQksabrYJfz9znPLjPeViMPnsKyMwsUw4AM7NM5RQAB7vdgC5wn/PgPueh433OZg3AzMyeKqcRgJmZ1XEAmJllKosAaLep/fVG0lpJ/yrpQUmnJL2z221aCpJKkr4q6QvdbstSKLZe/aykh4p/65/odpsWm6R3F/+n75f0aUkrut2mTpN0SNJjku6vO/YcSfdIerj4eUMnPuu6D4DETe2vNzPAb0fEjwCbgTsy6DPAO4EHu92IJfRx4J8j4oeBH+M677uk1cBvAYMR8WKqXzG/vbutWhR/BWxpOLYXuDciBoB7i9cLdt0HAHWb2kfENFDbmP66FRHnI+IrxfPvUf3FsLq7rVpcktYAPw98otttWQqSng38FPBJgIiYjoj/7W6rlkQZeIakMvBMrtydcNmLiC8B/9NweAj4VPH8U8AbOvFZOQTAfBvWZ0HSOuClwJe725JF9zHgPcBstxuyRH4ImAL+spj2+oSkZ3W7UYspIr4FfJjqBlLnqe48+MXutmrJPLfYZZHi502deNMcAmDBG9MvV5K+D/h74F0R8d1ut2exSHo98FhE3NfttiyhMrAROBARLwX+jw5NCzxdFfPeQ8B64HnAsyS9pbutWt5yCICUTe2vO5IqVH/53xURd3e7PYvsFcBtkh6hOsX3akl/3d0mLbpJYDIiaiO7z1INhOvZa4H/joipiHgCuBt4eZfbtFQelbQKoPj5WCfeNIcASNnU/roiSVTnhh+MiD/tdnsWW0S8LyLWRMQ6qv++/xIR1/VfhhHxbeCspBcVh14DPNDFJi2FbwKbJT2z+D/+Gq7zhe86o8DtxfPbgc934k3bbgq/3LXZmP569QrgrcDXJZ0ojv1eRIx1sU3Wee8A7ir+sDkD7OpyexZVRHxZ0meBr1C90u2rXIdfCSHp08CrgBslTQK/D3wA+Iykt1ENwjd35LP8VRBmZnnKYQrIzMyacACYmWXKAWBmlikHgJlZphwAZmaZcgCYmWXKAWBmlqn/B3+KAgL1vjasAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This is also pretty convincing\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "xs = np.linspace(0,10)\n",
    "plt.plot( xs, babylonian(xs)- np.sqrt(xs))\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## * Weakest Precondition\n",
    "\n",
    "Another useful and popular technique for translating programs to logical queries is the method of [predicate transformers](https://en.wikipedia.org/wiki/Predicate_transformer_semantics) of which Weakest Precondition is the most commonly encountered. In this method, the program statements are interpreted as acting upon logical formulas of the states of the computer, for example replacing a variable in the formula can be used to interpret assignment. A logical formula can be thought of as representing the set of states of the computer for which the formula is true. We can determine whether sets of states in the program remain in the \"safe\" set of states of the specification via a Z3 query.\n",
    "\n",
    "For a more in depth example of using Z3Py to perform a Weakest Precondition analysis in the style of this tutorial, check out this blog post https://www.philipzucker.com/weakest-precondition-z3py/\n",
    "\n",
    "For systems that implement variations of this technique over real languages see\n",
    "\n",
    "- Dafny https://rise4fun.com/Dafny/tutorial\n",
    "- ADA/Spark https://www.adacore.com/about-spark\n",
    "- Why3 http://why3.lri.fr/\n",
    "- Frama-C https://frama-c.com/\n",
    "- JML https://www.openjml.org/\n",
    "- Viper https://www.pm.inf.ethz.ch/research/viper.html\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sorting Networks\n",
    "\n",
    "[Sorting networks](https://en.wikipedia.org/wiki/Sorting_network) are a circuit flavored take on sorting. Although you can build circuits for any size input, any paritcular circuit works for a fixed sized input. They are like an unrolling of the loops or recursion of regular style sorting algorithms.\n",
    "\n",
    "\n",
    "Here's an interesting thing. We can go to Wikipedia and get a little python snippet for the comparison order of a [Batcher even-odd mergesort](https://en.wikipedia.org/w/index.php?title=Batcher_odd%E2%80%93even_mergesort&oldid=969926478) Kind of a confusing algorithm. Why does it even work? Is it even right?\n",
    "\n",
    "<img src=\"images/batcher.png\" alt=\"Drawing\" style=\"width: 500px;\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, 1),\n",
       " (2, 3),\n",
       " (0, 2),\n",
       " (1, 3),\n",
       " (1, 2),\n",
       " (4, 5),\n",
       " (6, 7),\n",
       " (4, 6),\n",
       " (5, 7),\n",
       " (5, 6),\n",
       " (0, 4),\n",
       " (2, 6),\n",
       " (2, 4),\n",
       " (1, 5),\n",
       " (3, 7),\n",
       " (3, 5),\n",
       " (1, 2),\n",
       " (3, 4),\n",
       " (5, 6)]"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# from https://en.wikipedia.org/w/index.php?title=Batcher_odd%E2%80%93even_mergesort&oldid=969926478\n",
    "def oddeven_merge(lo: int, hi: int, r: int):\n",
    "    step = r * 2\n",
    "    if step < hi - lo:\n",
    "        yield from oddeven_merge(lo, hi, step)\n",
    "        yield from oddeven_merge(lo + r, hi, step)\n",
    "        yield from [(i, i + r) for i in range(lo + r, hi - r, step)]\n",
    "    else:\n",
    "        yield (lo, lo + r)\n",
    "\n",
    "def oddeven_merge_sort_range(lo: int, hi: int):\n",
    "    \"\"\" sort the part of x with indices between lo and hi.\n",
    "\n",
    "    Note: endpoints (lo and hi) are included.\n",
    "    \"\"\"\n",
    "    if (hi - lo) >= 1:\n",
    "        # if there is more than one element, split the input\n",
    "        # down the middle and first sort the first and second\n",
    "        # half, followed by merging them.\n",
    "        mid = lo + ((hi - lo) // 2)\n",
    "        yield from oddeven_merge_sort_range(lo, mid)\n",
    "        yield from oddeven_merge_sort_range(mid + 1, hi)\n",
    "        yield from oddeven_merge(lo, hi, 1)\n",
    "\n",
    "def oddeven_merge_sort(length: int):\n",
    "    \"\"\" \"length\" is the length of the list to be sorted.\n",
    "    Returns a list of pairs of indices starting with 0 \"\"\"\n",
    "    yield from oddeven_merge_sort_range(0, length - 1)\n",
    "\n",
    "def compare_and_swap(x, a, b) -> None:\n",
    "    if x[a] > x[b]:\n",
    "        x[a], x[b] = x[b], x[a]\n",
    "        \n",
    "list(oddeven_merge_sort(8)) # compare and swap sequence for 8 items"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well we can confirm this relatively straightforwardly using z3 by replacing the implementation of compare_and_swap with it's z3 equivalent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_and_swap_z3(x,y):\n",
    "    x1, y1 = FreshInt(), FreshInt()\n",
    "    c = If(x <= y, And(x1 == x, y1 == y) , And(x1 == y, y1 == x) )\n",
    "    return x1, y1, c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's somewhat interesting to note that the output of odd_even_merge is a sequence of instructions, we can think of this as intepreting a very small 1 instruction (compare_and_swap) programming language."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "proved\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# predicates of interest\n",
    "def sorted_list(x): # list is sorted\n",
    "    return And([x <= y for x,y in zip(x , x[1:])])\n",
    "\n",
    "def verify_network(pairs_to_compare, N):\n",
    "    s = Solver()\n",
    "\n",
    "    a = [Int(f\"x_{i}\") for i in range(N)] #build initial array in z3 variables\n",
    "\n",
    "    #a_orig = a.copy() # keep around copy for correctness predicate\n",
    "    for i,j in pairs_to_compare:\n",
    "       x = a[i]\n",
    "       y = a[j]\n",
    "       x1, y1, c = compare_and_swap_z3(x,y) \n",
    "       a[i] = x1\n",
    "       a[j] = y1\n",
    "       s.add(c)\n",
    "\n",
    "    s.add(Not(sorted_list(a)))\n",
    "    if s.check() == unsat:\n",
    "        print(\"proved\")\n",
    "        return True\n",
    "    else:\n",
    "        return False\n",
    "N = 8\n",
    "verify_network(list(oddeven_merge_sort(N)), N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a reasonable sounding question: Can we prove that the Batcher even-odd sort is correct for all possible sizes N? \n",
    "Actually, this is very difficult to express in Z3 because here we run into a usage of induction, something Z3 is rather weak at.\n",
    "In fact such a question is best verified by a tool like Coq, which Cody will talk about tomorrow."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Synthesis of Optimal Sorting Networks\n",
    "An optimal sorting network is one that uses the minimal number of compare and swap operations. The batcher even-odd merge sort is good, but perhaps not optimal. What is? We we can find it by enumerating all possible combinations of compare and swaps and the finding one that passes verification as a sorting algorithm. This is an example of program synethesis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def all_swaps(m): # all pairs of integers from 0 to m-1\n",
    "    return [ [(i, j)] for i in range(m) for j in range(i+1, m) ]\n",
    "    \n",
    "# All list of length n on swaps on m wires \n",
    "def all_networks(m, n): \n",
    "   if n == 0:\n",
    "     return []\n",
    "   elif n == 1:\n",
    "     return all_swaps(m)\n",
    "   else:\n",
    "     return [ c + swap for c in all_networks(m,n-1) for swap in all_swaps(m)]\n",
    "\n",
    "\n",
    "def synthesize(N):\n",
    "    for n in range(N**2): # we can definitely do it in N*2 gates.\n",
    "       print(f\"trying network size: {n}\")\n",
    "       for pairs_to_compare in all_networks(N,n):   \n",
    "           if verify_network(pairs_to_compare, N):\n",
    "               return pairs_to_compare\n",
    "\n",
    "\n",
    "synthesize(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The space we are naively looking through is very very large. Consider what else might you want to do to speed things up. Early pruning of bad networks with testing, avoiding ever looking at obviously bad networks, maybe the completely optimal sorting network isn't necessary for you. Maybe a randomized search might be faster. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verifying a Neural Network\n",
    "\n",
    "Neural networks are all the rage these days. One way of looking at neural networks is that they are a flexibly parametrized class of functions. What makes them useful is: \n",
    "\n",
    "1. They can be used at insane scale due to their simplicity and excellent available implementations.\n",
    "2. There are intuitive ways to input abstract structure and symmetry expected of a problem, for example translation symmetry, or a hierarchy of small scale pattern recognition combining into large scale strutures\n",
    "3. Inspirational analogies can be drawn from nature.\n",
    "\n",
    "Neural networks made out of just ReLUs (rectified linear units, $ReLU(x) = max(0,x)$ ) and linear layers are particularly amenable to formal analysis. Regarding the weights as fixed, a neural network is a piecewise linear function and piecewise linear functions are among the functions most amenable to rigorous anaysis.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAfUElEQVR4nO3deXxU5dn/8c9lZJF9C/uqhE3LEiIIVYtLFXB7bKuy2Wq1FASX2qq0rk99+rTqz12EUmtrC4KiqNTiXrdqVZIAgbCGPaxhX0O2+/dHhj7TEJKT5EzOLN/365UXM3OW+XLmcOXmnpnrmHMOERGJfacEHUBERPyhgi4iEidU0EVE4oQKuohInFBBFxGJEyroIiJxQgVdpBxm9pCZzQw6h0hVqKCLiMQJFXRJSGZ2atAZRPymgi4Jw8w2mNk9ZpYFHDazzmb2upnlmdl6M7vtJNsNM7PccvZ1ca0EF/FIBV0SzWjgMqAF8AawBOgAXATcYWaXBphNpEZU0CXRPOOc2wycBSQ7537tnCtwzq0D/gCMCjaeSPVpHlESzebQn12A9ma2L2xZEvB57UcS8YcKuiSa4+1FNwPrnXMpHrY5DDQ4fsfMkoDkCGQTqRFNuUii+gY4EHqT9DQzSzKzs8zs7HLWXQ3UN7PLzKwOcB9Qr1bTiniggi4JyTlXDFwB9AfWA7uAF4Cm5ay7H7gltHwLpSP23LLriQTNdIELEZH4oBG6iEicUEEXEYkTKugiInFCBV1EJE4E9jn0Vq1aua5duwb19CIiMSkjI2OXc67c70EEVtC7du1Kenp6UE8vIhKTzGzjyZZpykVEJE6ooIuIxAkVdBGROKGCLiISJ1TQRUTiRKUF3cxeNLOdZrbsJMvNzJ4xsxwzyzKzVP9jiohIZbyM0P8MDK9g+QggJfQzHphW81giIlJVlRZ059xnwJ4KVrkK+Isr9RXQzMza+RVQRCSePP3hGrK37o/Ivv34YlEH/u+yXlDaJ7oDsK3simY2ntJRPJ07d/bhqUVEYsfrGbk8+eFqCotLOLP9Ca33a8yPN0WtnMfKbbLunJvhnEtzzqUlJ+sKXiKSOFZuP8C9by5lyOktueNiL1c+rDo/Cnou0Cnsfkdgqw/7FRGJCwfyC5k4M5Mm9evwzOgBnJoUmQ8Y+rHX+cAPQ592OQfY75w7YbpFRCQROee4e24Wm/Yc4bkxqSQ3jtzlaCudQzez2cAwoJWZ5QIPAnVCQacDC4CRQA5wBLgxUmFFRGLNH/+5nnezt3PvyN4M6tYios9VaUF3zo2uZLkDJvmWSEQkTizcsIffvrOS4We25ebzukX8+fRNURGRCMg7eIxJszLp1Pw0Hr2mL2blfX7EX4H1QxcRiVdFxSXcNnsRB/ILeenHg2hSv06tPK8KuoiIzx7/YDX/Wreb/3dNP3q3a1Jrz6spFxERH32wfAfTPlnL6EGd+cHAjrX63CroIiI+2bT7CHe+upizOjThwSv61Przq6CLiPggv7CYibMyOMWMaWMHUr9OUq1n0By6iIgPHpqfTfbWA7x4QxqdWjQIJING6CIiNTQ3fTNzFm5m0gVncGGvNoHlUEEXEamB5VsPcN+byxh6Rkvu/G7PQLOooIuIVNOB/EJumZVBswalTbeSTon8l4cqojl0EZFqcM7xi1eXkLv3KHPGn0OrRpFruuWVRugiItXwh8/X8f7yHUwZ0Yu0rpFtuuWVCrqISBV9vW43j7y7ipHfastN50a+6ZZXKugiIlWw80A+k2cvokuLBjzy/dppuuWV5tBFRDwqKi5h8uxFHMwv5K83DaJxLTXd8koFXUTEo8feX8U36/fw5HX96NW29ppueaUpFxERD97P3s7vP13H2MGduXpA7Tbd8koFXUSkEht3H+bnc5fQt2NTHgig6ZZXKugiIhXILyxmwsxMTjFj6phU6p1a+023vNIcuohIBR54axkrth3gTzecHVjTLa80QhcROYlXF27m1fRcbr2wOxf0ah10nEqpoIuIlCN7637uf2sZ53ZvxR0X9wg6jicq6CIiZew/WsjEmZk0b1CXp0f1D7zplleaQxcRCeOc4xdzl7B131Fe+ekQWkZB0y2vNEIXEQnz+8/W8cHyHfxqZG8GdmkedJwqUUEXEQn519rdPPruSi7r244bv9016DhVpoIuIkJp061bZy+ia6uGUdd0yyvNoYtIwissLmHyy4s4fKyIl38ymEb1YrM0xmZqEREfPfbeKr7ZsIenrutPjzaNg45TbZpyEZGE9u6y7cz4bB3Xn9OF/xrQIeg4NaKCLiIJa/2uw9w1dwn9OjXjvst7Bx2nxjwVdDMbbmarzCzHzKaUs7ypmf3NzJaYWbaZ3eh/VBER/xwtKGbizAySkoypYwZEddMtryot6GaWBEwFRgB9gNFmVrZ/5CRguXOuHzAMeNzM6vqcVUTEF8457n9rGat2HOSp6/rTsXl0N93yyssIfRCQ45xb55wrAOYAV5VZxwGNrfRzPo2APUCRr0lFRHzyysLNvJaRy60XpjCsZ/Q33fLKS0HvAGwOu58beizcc0BvYCuwFLjdOVdSdkdmNt7M0s0sPS8vr5qRRUSqb9mW/TwwP5vzUlpx+0UpQcfxlZeCXt6n612Z+5cCi4H2QH/gOTM74YJ7zrkZzrk051xacnJylcOKiNTE/iOFTJyVQcuGdXl61ICYabrllZeCngt0CrvfkdKReLgbgXmuVA6wHujlT0QRkZorKXH8fO5itu/PZ+rYVFo0jL+3+bwU9IVAipl1C73ROQqYX2adTcBFAGbWBugJrPMzqIhITUz7dC0frtjJvSN7k9o5tppueVXpN0Wdc0VmNhl4D0gCXnTOZZvZhNDy6cDDwJ/NbCmlUzT3OOd2RTC3iIhnX67dxePvr+KKfu350dCuQceJGE9f/XfOLQAWlHlsetjtrcAl/kYTEam57fvzuW32Irq1asjvvvetmGy65ZV6uYhI3CptupXJkYJiZv/kHBrGaNMtr+L7byciCe2Rd1aSvnEvz4weQEoMN93ySr1cRCQuvbN0Gy/8cz0/GtKFK/u1DzpOrVBBF5G4sy7vEHe9lkX/Ts2497KynUrilwq6iMSVowXF3DIrkzpJxtSxqdQ9NXHKnObQRSRuOOe4982lrNpxkJduHESHZqcFHalWJc6vLhGJe7O/2cy8zC3cflEK5/dIvPYiKugiEheycvfx0Pxszu+RzG0XxlfTLa9U0EUk5u07UsDEmZm0alSXp67rzylx1nTLK82hi0hMKylx/OyVxew8mM/cCUPjsumWVxqhi0hMe/6THD5elcf9l/ehf6dmQccJlAq6iMSsL3J28cQHq7myX3uuP6dL0HECp4IuIjHpeNOtM5Ib8ds4b7rllQq6iMScwuISJr2cSX5hMdPGDYz7plte6SiISMz57YKVZGzcy3NjBtC9daOg40QNjdBFJKb8PWsbL36xnhuGduXyvonRdMsrFXQRiRlr8w5x92tLSO3cjF+N7B10nKijgi4iMeFIQRETZ2ZQr05SwjXd8kpz6CIS9Zxz3PvGMtbsPMRffjyIdk0Tq+mWV/oVJyJRb9bXm3hj0RZ+dnEPzktJvKZbXqmgi0hUW7J5H7/+23KG9Uxm8gXdg44T1VTQRSRq7T1cwC2zMkluXI8nr03cplteaQ5dRKJSSYnjZ68uJu/gMeZOGELzBG665ZVG6CISlZ77OIdPVuVx/xV96JfgTbe8UkEXkajz+Zo8nvxwNVcP6MC4wZ2DjhMzVNBFJKps3XeU2+csJqV1I35z9VlqulUFKugiEjUKikqbbhUUlTBt3EAa1NXbfFWhoyUiUeN/F6xg0aZ9PD82lTOS1XSrqjRCF5Go8LclW/nzlxv48be7MfJb7YKOE5NU0EUkcDk7DzHl9SwGdmnOL0f2CjpOzPJU0M1suJmtMrMcM5tyknWGmdliM8s2s0/9jSki8erwsdKmW/XrJDF1TCp1kjTOrK5K59DNLAmYCnwXyAUWmtl859zysHWaAc8Dw51zm8ysdaQCi0j8cM7xy3lLWZt3iL/eNJi2TesHHSmmeflVOAjIcc6tc84VAHOAq8qsMwaY55zbBOCc2+lvTBGJR3/9aiPzl2zlzu/24NvdWwUdJ+Z5KegdgM1h93NDj4XrATQ3s0/MLMPMfljejsxsvJmlm1l6Xl5e9RKLSFxYtGkvD7+9nAt7teaWYWq65QcvBb28T/W7MvdPBQYClwGXAvebWY8TNnJuhnMuzTmXlpysFpgiiWrP4QImzcqkTZP6PHFtPzXd8omXz6HnAp3C7ncEtpazzi7n3GHgsJl9BvQDVvuSUkTiRnGJ445XFrPrUAGvTRxCswZquuUXLyP0hUCKmXUzs7rAKGB+mXXeAs4zs1PNrAEwGFjhb1QRiQfP/mMNn63O48Er+9C3o5pu+anSEbpzrsjMJgPvAUnAi865bDObEFo+3Tm3wszeBbKAEuAF59yySAYXkdjz6eo8nv5oDd9L7cCYQWq65Tdzrux0eO1IS0tz6enpgTy3iNS+LfuOcvkzn9OmSX3euOXbnFY3KehIMcnMMpxzaeUt0yf4RSTiCopKmDQrk8Jix/NjU1XMI0TNuUQk4n7z9+Us3ryP6eNSOV1NtyJGI3QRiaj5S7by0r82cvO53Rh+lppuRZIKuohEzJodB5nyehZnd23OPSPUdCvSVNBFJCIOHStiwswMGtRN4jk13aoVmkMXEd8555jyehbrdx1m5s2DadNETbdqg35liojvXvpyA29nbePnl/Rk6BlqulVbVNBFxFeZm/bymwUruKhXayZ+54yg4yQUFXQR8c3uQ8eYNCuTtk3r88S1/dV0q5ZpDl1EfHG86dbuwwXMmziUpg3qBB0p4WiELiK+ePqjNXy+Zhe/vvJMzurQNOg4CUkFXURq7JNVO3n2H2v4wcCOXHd2p8o3kIhQQReRGsnde4Q7XllMzzaNefiqszDTvHlQVNBFpNqOFRUzaVYmxcWO6eMGqulWwPSmqIhU2/+8vYIlufuZPm4gXVs1DDpOwtMIXUSq5c1FW/jrVxsZf/7pDD+rbdBxBBV0EamG1TsO8st5SxnUtQV3X9oz6DgSooIuIlVyvOlWw3qn8tyYAZyqpltRQ6+EiHjmnOOe17LYsOswz44eQGs13YoqKugi4tmfvtjA35du465LezHkjJZBx5EyVNBFxJOMjXv43wUruLh3GyZ85/Sg40g5VNBFpFK7Dh1j0qxFtG92Go9f209fHopS+hy6iFSouMRx+5xF7D1SwLxbhtL0NDXdilYq6CJSoac+XM0XObt59Pt9ObO9mm5FM025iMhJfbxyJ8/+I4dr0zpyrZpuRT0VdBEp1+Y9pU23+rRrwq+vOivoOOKBCrqInCC/sJhbZmVS4hzTxqVSv46absUCzaGLyAl+/fZylm7Zz4zrB9KlpZpuxQqN0EXkP8zLzOXlrzfx0++cziVnqulWLFFBF5F/W7n9AL96YymDu7XgrkvUdCvWqKCLCAAH8wuZODOTxvXr8KyabsUkT6+YmQ03s1VmlmNmUypY72wzKzazH/gXUUQizTnH3a9lsWnPEZ4bPYDWjdV0KxZVWtDNLAmYCowA+gCjzazPSdZ7BHjP75AiEll//Od63lm2nbsv7cng09V0K1Z5GaEPAnKcc+uccwXAHOCqcta7FXgd2OljPhGJsPQNe/jdOyu5pE8bxp+vpluxzEtB7wBsDrufG3rs38ysA3A1ML2iHZnZeDNLN7P0vLy8qmYVEZ/tOnSMSS9n0qH5aTx2jZpuxTovBb28V9iVuf8UcI9zrriiHTnnZjjn0pxzacnJyV4zikgEFJc4bpu9iH1HCpk2dqCabsUBL18sygXCmzh0BLaWWScNmBP67d4KGGlmRc65N31JKSK+e+KDVXy5djeP/aAvfdo3CTqO+MBLQV8IpJhZN2ALMAoYE76Cc67b8dtm9mfgbRVzkej10YodTP14LaPO7sQ1aWq6FS8qLejOuSIzm0zpp1eSgBedc9lmNiG0vMJ5cxGJLpt2H+FnryzmzPZNeOjKM4OOIz7y1MvFObcAWFDmsXILuXPuhprHEpFIyC8s5paXMwCYNnagmm7FGTXnEkkg//23bJZtOcALP0yjc8sGQccRn+m7vSIJ4rWMXGZ/s5mJw87g4j5tgo4jEaCCLpIAVmw7wL1vLGXI6S35+Xd7BB1HIkQFXSTOHcgvZOLMDJqeVodnRqvpVjzTHLpIHHPOcffcLDbvPcrsn5xDcuN6QUeSCNKvapE49sLn63k3eztThvdiULcWQceRCFNBF4lT36zfw+/eXcnwM9ty83ndKt9AYp4Kukgc2nkwn8kvZ9K5RQMevaavmm4lCBV0kThTVFzCbbMXcSC/kGnjUmlSX023EoXeFBWJM49/sJqv1u3h8Wv60autmm4lEo3QReLIB8t3MO2TtYwe1JnvD+wYdBypZSroInFi4+7D3PnqYs7q0IQHrzjhKpGSAFTQReJAfmExE2dmcoqZmm4lMM2hi8SBB9/KZvm2A7x4QxqdWqjpVqLSCF0kxr2avplX0jcz6YIzuLCXmm4lMhV0kRiWvXU/97+5jKFntOTO7/YMOo4ETAVdJEbtP1rILbMyadagtOlW0in68lCi0xy6SAxyznHX3CVs2XuUOePPoVUjNd0SjdBFYtKMz9bx/vIdTBnRi7SuarolpVTQRWLM1+t28+h7qxj5rbbcdK6absn/UUEXiSE7D+QzefYiurRowCPfV9Mt+U+aQxeJEUXFJUyevYhD+UXMvGkwjdV0S8pQQReJEY+9t4pv1u/hyev60bNt46DjSBTSlItIDHgvezu//2wdYwd35uoBarol5VNBF4lyG3Yd5hevLqFvx6Y8oKZbUgEVdJEoll9YzMRZmZxyijF1TCr1TlXTLTk5zaGLRLH731zGim0H+NMNZ6vpllRKI3SRKPXKwk3Mzcjl1gu7c0Gv1kHHkRiggi4ShZZt2c/9b2VzbvdW3HFxj6DjSIxQQReJMsebbrVoUJenR/VX0y3xzFNBN7PhZrbKzHLMbEo5y8eaWVbo50sz6+d/VJH4V1Li+PmrS9i67yhTx6bSUk23pAoqLehmlgRMBUYAfYDRZlb2s1Prge845/oCDwMz/A4qkgh+/9k6Plyxg3sv683ALs2DjiMxxssIfRCQ45xb55wrAOYAV4Wv4Jz70jm3N3T3K0DffBCpon+t3c1j763ksr7tuGFo16DjSAzyUtA7AJvD7ueGHjuZm4B3yltgZuPNLN3M0vPy8rynFIlzOw7kc+vsTLq1aqimW1JtXgp6eWeWK3dFswsoLej3lLfcOTfDOZfmnEtLTk72nlIkjhUWlzD55UwOHytm2riBNKqnr4dI9Xg5c3KBTmH3OwJby65kZn2BF4ARzrnd/sQTiX+PvruShRv28vSo/vRoo6ZbUn1eRugLgRQz62ZmdYFRwPzwFcysMzAPuN45t9r/mCLx6d1l2/jD5+u5/pwuXNW/oplMkcpVOkJ3zhWZ2WTgPSAJeNE5l21mE0LLpwMPAC2B50Nzf0XOubTIxRaJfet3HeauuVn069SM+y7vHXQciQOeJuuccwuABWUemx52+2bgZn+jicSvowXFTJyZQVKSMXXMADXdEl/o3ReRWuac4743l7Fqx0H+dMPZdGyuplviD331X6SWzVm4mdczc7n1whSG9VTTLfGPCrpILVq2ZT8Pzs/mvJRW3H5RStBxJM6ooIvUkv1HCpkwM4OWDevy9KgBarolvtMcukgtKClx3PnqYnYcyOeVnw6hRcO6QUeSOKQRukgtmPbpWj5auZP7LutDamc13ZLIUEEXibAvcnbx+PuruKJfe344pEvQcSSOqaCLRND2/fncNnsRpyc34nff+5aabklEaQ5dJEKON906WljMK+NSaaimWxJhOsNEIuR376wkfeNenhk9gO6t1XRLIk9TLiIRsGDpNv74z/X8aEgXruzXPug4kiBU0EV8ti7vEHe/lkX/Ts2497KyV2sUiRwVdBEfHSkoYuLMTOokGVPHplL3VP0Tk9qjOXQRnzjnuO+NZazeeZCXbhxEh2anBR1JEoyGDyI+efmbTcxbtIXbL0rh/B66xKLUPhV0ER9k5e7jv+cv5/weydx2oZpuSTBU0EVqaN+RAibOzKRVo7o8dV1/TlHTLQmI5tBFaqCkxHHHK4vZeTCfuROGqumWBEojdJEamPpxDp+syuOBy/vQv1OzoONIglNBF6mmf67ZxRMfruaq/u0Zd46abknwVNBFqmHb/qPcNmcR3ZMb8Vs13ZIooYIuUkUFRSVMmpXJscJipo0bSIO6eitKooPORJEq+u07K8jctI/nxgyge+tGQccR+TeN0EWq4O2srfzpiw3cMLQrl/dV0y2JLiroIh7l7DzEPa9lkdq5Gb8a2TvoOCInUEEX8eBIQRG3zMqgXp0kNd2SqKU5dJFKOOf41bylrNl5iL/8eBDtmqrplkQnDTNEKjHz6028uXgrP7u4B+elqOmWRC8VdJEKLN68j4f/tpxhPZOZfEH3oOOIVEgFXeQk9h4uYNKsTJIb1+PJa9V0S6Kf5tBFynG86VbewWO8NnEIzdV0S2KApxG6mQ03s1VmlmNmU8pZbmb2TGh5lpml+h9VpPY8+48cPl2dxwNX9KFvRzXdkthQaUE3syRgKjAC6AOMNrOyV74dAaSEfsYD03zOKVJrPludx1MfrebqAR0YO7hz0HFEPPMy5TIIyHHOrQMwsznAVcDysHWuAv7inHPAV2bWzMzaOee2+R3409V5/M/byytfUaSaNu45QkrrRvzm6rPUdEtiipeC3gHYHHY/FxjsYZ0OwH8UdDMbT+kIns6dqzfyaVTvVFLaqH+GRE6/Ts24/aIUNd2SmOPljC1viOKqsQ7OuRnADIC0tLQTlnsxsEtzBnYZWJ1NRUTimpc3RXOBTmH3OwJbq7GOiIhEkJeCvhBIMbNuZlYXGAXML7POfOCHoU+7nAPsj8T8uYiInFylUy7OuSIzmwy8ByQBLzrnss1sQmj5dGABMBLIAY4AN0YusoiIlMfTuz7OuQWUFu3wx6aH3XbAJH+jiYhIVeir/yIicUIFXUQkTqigi4jECRV0EZE4YaXvZwbwxGZ5wMZqbt4K2OVjHL9Eay6I3mzKVTXKVTXxmKuLc67cK60EVtBrwszSnXNpQecoK1pzQfRmU66qUa6qSbRcmnIREYkTKugiInEiVgv6jKADnES05oLozaZcVaNcVZNQuWJyDl1ERE4UqyN0EREpQwVdRCRORG1BN7NrzCzbzErM7KQf7znZBazNrIWZfWBma0J/NvcpV6X7NbOeZrY47OeAmd0RWvaQmW0JWzaytnKF1ttgZktDz51e1e0jkcvMOpnZx2a2IvSa3x62zNfjVZMLnle2bYRzjQ3lyTKzL82sX9iycl/TWso1zMz2h70+D3jdNsK57grLtMzMis2sRWhZJI/Xi2a208yWnWR5ZM8v51xU/gC9gZ7AJ0DaSdZJAtYCpwN1gSVAn9CyR4EpodtTgEd8ylWl/YYybqf0ywAADwG/iMDx8pQL2AC0qunfy89cQDsgNXS7MbA67HX07XhVdL6ErTMSeIfSq3CdA3ztddsI5xoKNA/dHnE8V0WvaS3lGga8XZ1tI5mrzPpXAP+I9PEK7ft8IBVYdpLlET2/onaE7pxb4ZxbVclq/76AtXOuADh+AWtCf74Uuv0S8F8+Ravqfi8C1jrnqvutWK9q+vcN7Hg557Y55zJDtw8CKyi9Jq3fKjpfwvP+xZX6CmhmZu08bhuxXM65L51ze0N3v6L0qmCRVpO/c6DHq4zRwGyfnrtCzrnPgD0VrBLR8ytqC7pHJ7s4NUAbF7pqUujP1j49Z1X3O4oTT6bJof9uvejX1EYVcjngfTPLsNKLdld1+0jlAsDMugIDgK/DHvbreFV0vlS2jpdtI5kr3E2UjvKOO9lrWlu5hpjZEjN7x8zOrOK2kcyFmTUAhgOvhz0cqePlRUTPr0Ava25mHwJty1l0r3PuLS+7KOexGn8Os6JcVdxPXeBK4JdhD08DHqY058PA48CPazHXt51zW82sNfCBma0MjSqqzcfj1YjSf3h3OOcOhB6u9vEq7ynKeczrBc8jcq5V8pwnrmh2AaUF/dywh31/TauQK5PS6cRDofc33gRSPG4byVzHXQF84ZwLHzVH6nh5EdHzK9CC7py7uIa7qOji1DvMrJ1zblvovzQ7/chlZlXZ7wgg0zm3I2zf/75tZn8A3q7NXM65raE/d5rZG5T+V+8zAj5eZlaH0mI+yzk3L2zf1T5e5ajJBc/retg2krkws77AC8AI59zu449X8JpGPFfYL16ccwvM7Hkza+Vl20jmCnPC/5AjeLy8iOj5FetTLhVdwHo+8KPQ7R8BXkb8XlRlvyfM3YWK2nFXA+W+Gx6JXGbW0MwaH78NXBL2/IEdLzMz4I/ACufcE2WW+Xm8anLBcy/bRiyXmXUG5gHXO+dWhz1e0WtaG7nahl4/zGwQpTVlt5dtI5krlKcp8B3CzrkIHy8vInt+ReKdXj9+KP3HmwscA3YA74Uebw8sCFtvJKWfilhL6VTN8cdbAh8Ba0J/tvApV7n7LSdXA0pP7KZltv8rsBTICr1g7WorF6XvoC8J/WRHy/GidPrAhY7J4tDPyEgcr/LOF2ACMCF024CpoeVLCfuE1cnONZ+OU2W5XgD2hh2f9Mpe01rKNTn0vEsofbN2aDQcr9D9G4A5ZbaL9PGaDWwDCimtXzfV5vmlr/6LiMSJWJ9yERGREBV0EZE4oYIuIhInVNBFROKECrqISJxQQRcRiRMq6CIiceL/A/x5OzVG6JyvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "xs = np.linspace(-1,1,1000)\n",
    "plt.plot(xs, np.maximum(0,xs))\n",
    "plt.title(\"ReLU\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most machine learning tasks don't have a mathematically precise specification. What is the mathematically precise definition of a picture of a horse? We could try to come up with something, but it would be rather suspect. \n",
    "\n",
    "The surprising power of neural networks has often been in areas where a mathematically precise rule based specification has not worked. Tasks that _do_ have a spec are questionable areas for machine learning techniques, because how can you know that the network meets the spec? Also, one would suspect such a problem would have exploitable structure such that you might be better off with a less general purpose approach. \n",
    "\n",
    "However, there a a couple areas where one does have reasonable formal questions one might want to ask about a neural network:\n",
    "- Robustness around training examples and valdiation examples. Adversarial examples.\n",
    "- Control Systems - There are many control systems which we do have a reasonable formal spec one could write such as walking robots. These systems are so high dimensional that it is difficult to derive a good controller from the spec, and hence reinforcement learning may be a reasonable option\n",
    "- Neural networks as computational accelerators. There are problems which we know how to solve, but are very slow. It may be useful to presolve a large number of examples offline and train a neural network to give good estimates. Neural networks can be evaluated very quickly and easily thanks to modern frameworks.\n",
    "- Dead neuron pruning. Can we show that some neuron never fires over the domain of interest?\n",
    "  \n",
    "\n",
    "Here we'll train a network to fit a simple 1-D function and then verify it's accuracy. This is of course a complete toy. How often do you see 1-D input space neural networks? Not often I bet.\n",
    "\n",
    "But it's nice for a couple reasons:\n",
    "1. Because we can visualize it. \n",
    "2. It drives home the point about neural nets being a space of piecewise linear function approximators, and how similar training is to curve fitting.\n",
    "3. It's simple enough that z3 can crush it. There is a big question if z3 or other formal methods can scale to modern neural nets. Modern top of the line neural nets are insanely huge. There are special purpose SMT solvers being built for this task. There are other approaches being pursued such as the use of mixed integer programming solvers. So this is an area of research. Come see tomorrow's talk for more!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3ic5ZXw/+9R75LVpZFluci9yEa40UMzJmCbkgC7BEiAkIRskmuzuyTZ3Tfv5rfvJtmSTU8gIYGEFooLYELvYGNZ7l2WbVWr9y7N/ftjZshESFabmWfK+VzXXJp5yjxHo9GcubsYY1BKKRW6wqwOQCmllLU0ESilVIjTRKCUUiFOE4FSSoU4TQRKKRXiIqwOYDLS09NNQUGB1WEopVRA2b17d6MxJmP49oBMBAUFBZSUlFgdhlJKBRQROTPSdq0aUkqpEKeJQCmlQpwmAqWUCnGaCJRSKsRpIlBKqRDnkUQgIg+LSL2IHBxlv4jIT0SkTET2i8gKt33rROSYc98DnohHKaXU+HmqRPB7YN059l8DFDpv9wK/BBCRcODnzv0LgVtFZKGHYlJKKTUOHhlHYIx5R0QKznHIBuBR45jzeoeIpIhIDlAAlBljygFE5EnnsYc9Eddwbxyt4+jZDr586RxvPL0KQT39QxyubedkQyfNXf30DdhJiIkgIzGahTmJzEpPICxMrA5TBYH69l4e/fAMN6ywMSsjwaPP7asBZTag0u1xlXPbSNtXjfQEInIvjtIE+fn5kwri3RONPF1SpYlATUlP/xAv7K9h274adp5qpn/QPuqxybGRXLkwi41FNi6Yk4aIJgU1OScbuvjZm2WsmZ0WsIlgpHe/Ocf2T2405kHgQYDi4uJJraaTlRRDZ98gnX2DJEQH5KBqZaGuvkEefu8UD71bTnvvIAVpcdy+egarZ6VRmJlAZlI00RHhdPYOUtvew4GqNj482cTLB8/yzO4q5mcncv+n5nDtkhxNCGrC6jt6AchKivb4c/vq07AKmO72OA+oAaJG2e4Vrhewvr2XBA9nVBXcXjpQy79sPURjZx9XLMjk7otmsWpm6ogf6MlxkSTHRTI/O4mbi6fTOzDE8/tq+PU75dz/+B4eLTjDv21cxPzsJAt+ExWoGjr6AMhIjPH4c/uq++g24HPO3kOrgTZjTC2wCygUkZkiEgXc4jzWK7KcL2Bde5+3LqGCTFvPAF95rJQvPVZKVlI0z35pLb+543xWzxp/NU9MZDg3F0/n5a9fzH/csISyhk6u/+n7PPzeKXSpWDVe9R19REeEkRTj+e/vHnlGEXkCuBRIF5Eq4P8AkQDGmF8B24H1QBnQDdzl3DcoIvcDLwPhwMPGmEOeiGkkma4SgbOIpdS5HK/r4N5HS6hu7eGbV83li5fMJjJ88t+dwsOEW1fmc9XCLP7xmf382wuH+eBkEz++pYh4rapUY6hr7yUzKdor1Yqe6jV06xj7DfCVUfZtx5EovC4zyVEiqNcSgRrDW8fq+fJjpcRHR/DEPaspLkj12HOnJUTzmzuK+d37p/n/XjzMzb/6kIfvPJ/sZM8X+VXwqG/vI9ML1UIQYiOLE6MjiI0Mp65dSwRqdC/ur+WeR0uYmR7P8/df6NEk4CIifP7Cmfz2zvM509TFpl+8z+nGLo9fRwWP+o5eMhM931AMIZYIRISspGjqOrREoEa2ZU81X32ilKLpKTxx72qvf0u/bF4mf7pvDb0DQ9z60A5NBmpU9R19ZCVpicAjMhNjtESgRvTG0Tr+/ul9rJqZxqOfX0VSTKRPrrsoN5nH71n9cTKobu3xyXVV4OgdGKKjd5AMLRF4RmZS9MfdsJRy2XW6mS/9sZRFuUk8dEcxsVHhPr3+gpwkHrt7NZ29g9z1u49o6xnw6fWVf3O1a2rVkIdkJTlKBNptT7mcaeri7kdKsE2L5Xd3nm/ZYMOFuUn8+vbzONXYxRf/UHLOEcsqtNQ5ezpmatWQZ2QlRdPdP0Rn36DVoSg/0Nk3yD2PliACv79zJWkJ3vnGNV5r56TznzctY0d5M997wStTbqkApCUCD3M1tuigMmW3G/7+T3spq+/kZ7euID8tzuqQANi43MYXL57FH3ac4dndVVaHo/yAa+yTJgIPcTW26KAy9dC75bx8qI5vr1/AhYXpVofzV/7h6nmsnpXKtzcf4FBNm9XhKIudbe8lMlxIjY/yyvOHXCLI0kFlCthf1cp/vnyMqxdl8YULZ1odzidEhIfxs9tWMC0uiq8+vofufq3KDGWuwWTemqwwZBOBdiENXV19g3ztyb2kJ0TzgxuX+u1MoOkJ0fzPZ5dxqqmLf3/xiNXhKAudbev16piWkEsECdERxEeFaxtBCPu35w9zuqmLH322iJQ47xS1PWXt7HTuvnAmj+2s4I2jdVaHoyxS195Ltpd6DEEIJgJwdMGq0zaCkPTO8QaeKqnkixfPZs3sNKvDGZdvXj2P+dmJ/OMz+2np6rc6HGWBuvZer40qhlBNBInRNGiJIOR09Q3y7c0HmJUez9evKLQ6nHGLjgjnR58torV7gH/frlVEoaajd4Cu/iGvLEjjEpKJIEtLBCHpv145RlVLDz+4aSkxkb4dOTxVC3KSuPfiWTyzu4oPyhqtDkf5kKs9U9sIPCwrKVpHF4eY0ooWfv/BaT63Zgbne2E2UV/4u8sLKUiL41ubD9A7MGR1OMpHzrY5ai+0asjDMhNj6B2w096rXfJCwZDd8M+bD5KdFMM/rptvdTiTFhMZzv/btIQzTd38+PUTVoejfOTjEoG/JwIRWScix0SkTEQeGGH/P4jIXuftoIgMiUiqc99pETng3FfiiXjG4lqprEGrh0LCEx9VcLi2nX++dqFl8wh5yto56dy4Io/fvFvOKZ2yOiScbXctWu/HiUBEwoGfA9cAC4FbRWSh+zHGmP80xhQZY4qAbwFvG2Oa3Q65zLm/eKrxjIdOMxE6Wrr6+a9XjrFmVhrrl2RbHY5H/NM184iOCOffX9S5iEJBXXsvSTERXp0R1xMlgpVAmTGm3BjTDzwJbDjH8bcCT3jgupOmg8pCx/+8epyO3kG+e/0ivx04NlGZiTHc/6k5vHaknneON1gdjvIybw8mA88kAhtQ6fa4yrntE0QkDlgHPOu22QCviMhuEbl3tIuIyL0iUiIiJQ0NU3vzuyZu0hJBcDtS285jO89w++oZzMtOtDocj7rrggJmpMXxvRcOMzCk01UHszovrkzm4olEMNLXrNG641wHvD+sWugCY8wKHFVLXxGRi0c60RjzoDGm2BhTnJGRMaWA46MjSIiO0BJBkPvBn4+SGBPJN66Ya3UoHhcdEc531i/gRH0nj+04Y3U4yovq2rw7mAw8kwiqgOluj/OAmlGOvYVh1ULGmBrnz3pgM46qJq/LTIrWGUiD2AcnG3nrWANfuWw2yXG+WXLS165cmMWaWWn89I0yXV8jSA3ZDQ2dfV7tMQSeSQS7gEIRmSkiUTg+7LcNP0hEkoFLgK1u2+JFJNF1H7gKOOiBmMaUkxxDbZsmgmBkjOEHLx0lNzmGz60psDocrxER/nHdPJq6+nn4vVNWh6O8oLGzjyG7Icvf2wiMMYPA/cDLwBHgT8aYQyJyn4jc53boJuAVY4x7n7cs4D0R2Qd8BLxojPnzVGMaj+ykWM5qIghKLx08y76qNr5x5dyAG0E8Ucvzp3H1oiwefKecZp2HKOi4qq+zvLQgjYtHOlUbY7YD24dt+9Wwx78Hfj9sWzmwzBMxTFRuSgz1HX0MDtmJCA/JcXVBaWDIzn++fIy5WQncsCLP6nB84ptXzePVw+/wy7fK+M61C8c+QQUM15fVQOg1FJCyk2MYshsaO/VbVDB5ZncVpxq7+Ier5xMeFhzdRcdSmJXIDSvyeOTDM9S09lgdjvIgX4wqhhBOBDnODFvTpv84wWJgyM7P3yxjWV4yVyzItDocn/r6FYVg4Kdv6NQTweRsey/hYUJagnerhkI2EWQnxQJoO0EQ2VxaTVVLD393eWHQDB4br7xpcXz2/Ok8s7uKai0VBI269j4yEqK9XroN2UTgKhFoz6HgMDhk52dvlrHElsyn5odWacDlvktnA/Crt05aHInylLr2Xq/3GIIQTgQpcZFER4RxVquGgsKWvTVUNHeHZGnAxZYSy03n5fFUSaUOlgwSZ9t6yfbigjQuIZsIRITclFgtEQSBQWfbwMKcpJBrGxjuS5fMYchu+PXb5VaHojzAkQi0ROBV2Uk6qCwYPL+/hlONXSFdGnDJT4tjY5GNx3aeoaFD59IKZB29A3T0DZKTEuv1a4V0IshJjtHG4gBntxt+8eZJ5mcnctXCLKvD8QtfuWw2A0N2fvOulgoCmetLaq4mAu/KTo6hrr2XIbsuWRmo3jpez4n6Tr54ySzCQmTcwFhmZSSwfkkOj+2soL13wOpw1CS5xoTkamOxd+UkxzBoNzR1ahE6UP367XJyk2P49NJcq0PxK1+8eDadfYM8sbPC6lDUJLlKBFo15GXZyY4XWNsJAtO+ylZ2nmrm8xfOJFKnCfkrS/KSWTMrjd+9f5r+QV2vIBDVtvYQJt6fZwhCPBH8ZSyBdiENRA++U05iTAS3rMy3OhS/dO8lszjb3svz+0abFV75s5q2XjITY3wyF5omArREEIgqmrp56WAtf7NqRsAvSO8tl87NYG5WAg+9W44x2g4WaGpae8hJ8X77AIR4IkiNjyIqPEx7DgWg37xXTniYcNcFBVaH4rdEhHsumsXRsx28c6LR6nDUBNW29ZKb7P32AQjxRCAiZOsCNQGnuaufP5VUsrHI5vUl/ALdhiIbWUnRPPiOTjsRSIwxjhKBD3oMQYgnAnB0IdUSQWB54qMKegfs3HPxLKtD8XtREWHcuXYm75c1caS23epw1Di1dA/QN2j3SY8h8FAiEJF1InJMRMpE5IER9l8qIm0istd5+9fxnuttuckxOhV1ABkcsvPHHWe4cE46c7MSrQ4nINy6cjoxkWE8+uFpq0NR4+QaQ2ALlDYCEQkHfg5cAywEbhWRkZZJetcYU+S8/dsEz/Wa7ORY6tp7seugsoDwyuE6att6uWNtgdWhBIyUuCg2FtnYvKea1m5diCkQfDyGIIDaCFYCZcaYcmNMP/AksMEH53pETnIMA0OGxi4dVBYIfv/BafKmxYbsVNOTdcfaAnoH7Dy1q9LqUNQ4uLq0B1KvIRvg/u6qcm4bbo2I7BORl0Rk0QTPRUTuFZESESlpaGjwQNgOrnk8alq1ncDfHalt56NTzXxuzYyQWYbSUxbkJLFqZip/2HFGp1QJADWtvUSGC+nx3h9MBp5JBCP9Rw5/p5UCM4wxy4CfAlsmcK5jozEPGmOKjTHFGRkZkw52OJszEVS3aDuBv3vkg9PERIbxmeLpVocSkO5cW0BVSw+vH6mzOhQ1htq2HrKTY3w2f5YnEkEV4P6fmQf81VBGY0y7MabTeX87ECki6eM519tsH5cINBH4s9bufrbsrWbTchspcVFWhxOQrlyYRW5yDI98eNrqUNQYalt7fdY+AJ5JBLuAQhGZKSJRwC3ANvcDRCRbnBPFi8hK53WbxnOutyXFRpAQHaHrvPq5p3ZV0jtg10biKYgID+Nv18zg/bImTtR1WB2OOoeath6fzDrqMuVEYIwZBO4HXgaOAH8yxhwSkftE5D7nYTcBB0VkH/AT4BbjMOK5U41pIkQEW0osVVo15LeG7IY/7DjDqpmpzM9OsjqcgHbL+flERYTx+w9OWx2KGsWQ3VDX3uuzMQQAHpmkxVnds33Ytl+53f8Z8LPxnutrtmmxWiLwY28craeqpYfvrF9gdSgBLzU+iuuW5rJ1bw3fXr+AeJ2nye80dvYxMGQCq0QQDGwpsdpG4Mee+KiCzMRortAVyDzitlX5dPYNsk1nJfVLrs+iQGsjCHi5KbG09QzQ2TdodShqmOrWHt46Vs9niqfrmgMesiI/hfnZiTyui9b4JV8uUemi/1k4qoZAu5D6o6d2VWKAz56vXUY9RUS4bVU+B6rbOFDVZnU4apiqlm7gL59LvqCJALexBK3dFkei3A0O2fnTrkouLsxgemqc1eEElY3LbcRGhvP4R2esDkUNU9XSQ1JMBMmxkT67piYCIM9VItDRxX7lzWMNnG3v5bZVugKZpyXFRHLdshy27q2hQxe49yvVLT3Ypvn2i48mAiAjIZrIcNGqIT/jaiTWeYW847ZVM+juH2LrXm009idVLT0ffzn1FU0EQFiYkJOsXUj9iTYSe9+yvGQW5CTx+M4KXcrSTxhjqGrp1kRgFVtKLNUt2kbgL7SR2PtcjcaHa9vZp43GfqG1e4Cu/iHytGrIGrZpsToDqZ/QRmLf2ViUS1xUOI/v1EZjf+Ca4UBLBBbJTYmlrqOX/kG71aGEPG0k9p3EmEg+vTSHF/bX0qXjaCzn6jqqicAieSmxGIOuX+wHntpVQYY2EvvMZ4qn090/xPYDtVaHEvI+LhGkaNWQJVyDN6p0LIGl6jt6efNYAzessGkjsY+cN2Mas9LjebqkyupQQl51aw+J0REkxfp2Dij9T3Oy6UplfmHLnmqG7Iabz9NGYl8REW4qzuOj082UN3RaHU5Iq2rpxjYtFues/T6jicDJtTaojiWwjjGGp0uqWJ6fwpzMBKvDCSk3rsgjTOCZ3VoqsJJjDIHvO0hoInCKjggnKymaSu1Capl9VW2cqO/U0oAFspJiuHReJs+WVumaxhZxjCHw/WAy8FAiEJF1InJMRMpE5IER9v+NiOx33j4QkWVu+06LyAER2SsiJZ6IZ7LyU+OoaNZEYJVndlcSExnGp5flWB1KSLr5vDzq2vt450SD1aGEJNcMyAGZCEQkHPg5cA2wELhVRBYOO+wUcIkxZinwPeDBYfsvM8YUGWOKpxrPVExPjaNKE4ElegeG2La3hnWLskmK8d1kW+ovLl+QRWp8FE+XVFodSkj6yxiCwKwaWgmUGWPKjTH9wJPABvcDjDEfGGNanA934Fik3u/kp8ZR295L3+CQ1aGEnFcO19HeO8jNxVotZJWoiDA2Ftl49XAdzV39VocTcqwaQwCeSQQ2wP0rRJVz22i+ALzk9tgAr4jIbhG5d7STROReESkRkZKGBu8UXadPi8MYbTC2wtMlldhSYlkzK83qUELaZ87PY2DIsHVvtdWhhByrRhWDZxLBSP2cRmxtEpHLcCSCf3LbfIExZgWOqqWviMjFI51rjHnQGFNsjCnOyMiYaswjyk9zFMm0ncC3alp7eK+skRvPyyMszLfd5tRfm5+dxNK8ZMdcTzoRnU9VtfSQEO3bdQhcPJEIqgD38nwe8Il5bUVkKfAbYIMxpsm13RhT4/xZD2zGUdVkiXznvDaVWiLwqedKqzDG0ViprHfzeXkcPdvBoZp2q0MJKa4eQ74eQwCeSQS7gEIRmSkiUcAtwDb3A0QkH3gOuN0Yc9xte7yIJLruA1cBBz0Q06RkJEQTHRFGpZYIfMYYwzO7q1g9K1UnmPMT1y3LJSo8jOdKtXrIlyqbuy1pKAYPJAJjzCBwP/AycAT4kzHmkIjcJyL3OQ/7VyAN+MWwbqJZwHsisg/4CHjRGPPnqcY0WWFhwvTUOCqaNBH4yq7TLZxu6taxA34kJS6KT83PZNu+agaHdBJGXzDGUNHczYw0axKBRya0MMZsB7YP2/Yrt/t3A3ePcF45sGz4ditNnxarbQQ+9OzuKuKjwrlmSbbVoSg3N6yw8edDZ3n3RCOX6eR/XtfQ0UfPwJBliUBHFg+TnxpHZXO3NpT5QO+AY8bLqxdnExfl20m21LldOi+TaXGRPFuqU074whnnl898i6pHNREMMz01jo6+Qdp6dEFvb3v9SD0dfYPcsFwbif1NVEQY1y3L5dXDdbTr4vZed8ZZHT0jLd6S62siGMaVkbV6yPs276kmKymaNbN17IA/umFFHn2Ddl7SdQq8rqKpizD5yyzIvqaJYBgdS+AbzV39vHWsng1FNsJ17IBfWpaXzKz0eO095ANnmrvJTYklKsKaj2RNBMNMn6aJwBde2F/DoN2wafm5BqErK4kIN6ywsfNUs3ap9rIzTdb1GAJNBJ8QHx1BWnwUlc06qMybNu+pZn52IgtykqwORZ3DhiJHotYpJ7zrTFMX+anWtA+AJoIRTXf2HFLecaqxiz0VrVoaCADTU+NYNTOV50qrtSedl7T3DtDSPaAlAn+j6xJ415Y91Yj85dum8m83rLBR3tjFvqo2q0MJSq4BrDMsHFmviWAE+alxVLf2MKCjKj3OGMOWvdWsnZ1GdnKM1eGocbhmSQ7REWE8p2MKvMLVdTRfSwT+pSA9niG7+XhaWOU5pRUtnGnqZpOOHQgYSTGRXLkwi+f31dA/qF+OPO1Mcxdg3RgC0EQwopnpjsx8urHL4kiCz+Y91cREhrFusU4pEUhuXJFHS/cAbx2rtzqUoFPR1E1afBQJ0daNrtdEMIICZ2Yu10TgUf2Ddl7YX8tVC7MtfdOribuoMJ30hCgdU+AFZ5q6La0WAk0EI0qNjyIxJkJLBB721rF6WrsH2LRCG4kDTUR4GJ9emssbx+p1+hUPq2jutrShGDQRjEhEmJkez+kmTQSetHlPNekJUVw0J93qUNQkbFxuo3/QzssHz1odStDoGxyipq2HfAvbB0ATwahmpsdzSksEHtPWPcDrR+q5blkuEeH6tgtEy/KSKUiLY4sOLvOYqpYejLG26yhoIhhVQVo81a099A4MWR1KUNh+sJb+IbvONBrARIQNRTY+LG/ibFuv1eEEBVf1c0F6EJQIRGSdiBwTkTIReWCE/SIiP3Hu3y8iK8Z7rlVmpsdjDDrC2EM2l1YzOyOexTadUiKQbVxuwxh4ft8nliVXk+CqdZidEeCJQETCgZ8D1wALgVtFZOGww64BCp23e4FfTuBcS7gytFYPTV1lczcfnW7mhhV5lizMrTxnZno8y/KStXrIQ042dDEtLpKUuChL4/BEiWAlUGaMKTfG9ANPAhuGHbMBeNQ47ABSRCRnnOdaYqaz8UYbjKfONWHZhqJciyNRnrChyMahmnbK6jusDiXgnWrsZFZGgtVheCQR2IBKt8dVzm3jOWY85wIgIveKSImIlDQ0NEw56LEkx0WSGh+lJYIpckwpUcPKglTyplnbIKY849PLcggT2LJHq4em6lRjFzMtbh8AzySCkcr6w6cpHO2Y8Zzr2GjMg8aYYmNMcUZGxgRDnJyCtDhNBFPk+ObYyYblWhoIFpmJMVwwJ52t+3RG0qno7Bukrr0vaBJBFTDd7XEeMPyrwmjHjOdcyxSkx3O6URuLp2Lr3moiw4X1i3OsDkV50IYiG5XNPZRWtFodSsA67ScNxeCZRLALKBSRmSISBdwCbBt2zDbgc87eQ6uBNmNM7TjPtczMtHjOtvfS3T9odSgBachu2LavhkvmZjAt3trGMOVZVy/KIjoiTBesmYKTDZ0AzEwPgjYCY8wgcD/wMnAE+JMx5pCI3Cci9zkP2w6UA2XAQ8CXz3XuVGPyFFfPIS0VTM7OU03UtffpugNBKDEmkisWZvHC/lqdrn2STjV2IYKlC9K4eGTmL2PMdhwf9u7bfuV23wBfGe+5/sJVd3e6qYuFudr/faK27qkhPiqcKxZkWR2K8oKNRTZe3F/LeycauWx+ptXhBJxTjV3YUmKJiQy3OhQdWXwurkRQ7izCqfHrHRhi+8Farl6cTWyU9W905XmXzM0gJS5SxxRMUnmDf/QYAk0E5xQfHUFucgxl9ZoIJuqtY/V09A6yUauFglZURBjrl+TwyqE6uvq0HW0ijDGcauxith+MIQBNBGOak5XICU0EE7ZlTw3pCdGsnZ1mdSjKizYW2egZGOLVw3VWhxJQGjr66Owb1BJBoJiTkcDJhk7sdu0vPV5tPQO8cbSe65bl6EyjQa54xjRsKbFaPTRBrkWvZvlB11HQRDCmOZkJ9A7YqW7V9YvH68/OmUa1Wij4hYUJ1xfl8u6JRpo6+6wOJ2D8peuoJoKAUJjlqMMr0wbjcduyp4aCtDiW5iVbHYrygY1FNobshhcP1FodSsA4UddJfFQ4tpRYq0MBNBGMaY6zMeekthOMy9m2XnacamJDkU1nGg0R87ITmZ+dyJY9Wj00XifqO5iTmeA3/yOaCMYwLT6KtPgoTtRpIhiP5/fVYIxj3noVOjYU2SitaKWiSQdfjseJuk7mZCZaHcbHNBGMw+zMBK0aGqcte6tZlpfsN3Wfyjeud04xrlNOjK2te4D6jj7mZvlH11HQRDAuhZkJlNV36kyLYyir7+BQTbtOKRGCbCmxrJyZypa9OiPpWMoaHOs4FGoiCCxzMhNo6xmgsbPf6lD82pY9NYSJY756FXo2Ftk42dDFoZp2q0Pxa65q5kKtGgosczIdmfuErsg0KmMMW/dVc8GcdDITY6wOR1lg/ZJsIsNFG43HcKK+k5jIML/pMQSaCMbFlbm159DoSitaqGzu0bEDISwlLopL52WybV8NQzoAc1TH6xw9hsLC/KPHEGgiGJespGgSoiN0zqFz2LKnhuiIMK5enG11KMpCG4ts1Hf0sbO8yepQ/FZZfadfVQuBJoJxERHtOXQOA0N2XjxQyxULs0iI9sjM5ipAXb4gk4ToCJ1yYhQdvQPUtvV+XN3sLzQRjFNhZgLHzmoiGMl7Jxpp7urXaiFFTGQ46xZn89KBs/QODFkdjt9x1SoUBlMiEJFUEXlVRE44f04b4ZjpIvKmiBwRkUMi8jW3fd8VkWoR2eu8rZ9KPN40PzuRxs4+GnU+lU/YsrealLhILpmbYXUoyg9sLLLR0TfIG0frrQ7F77hmMi7MCq6qoQeA140xhcDrzsfDDQJ/b4xZAKwGviIiC932/8gYU+S8+eVKZQALchwrlB07qz2H3HX1DfLKoTrWL8khKkILmArWzE4jKyma50q1emi4Y2c7iI4IIz/V+uUp3U31P3cD8Ijz/iPAxuEHGGNqjTGlzvsdONYmDrg6hHnZjgx+VBPBX3n1cB09A0NaLaQ+Fh4mbCiy8daxepq7dOyNu6Nn25mXnUi4H/UYgqkngixjTC04PvCBcy5cKiIFwHJgp1glqhYAABuESURBVNvm+0Vkv4g8PFLVktu594pIiYiUNDQ0TDHsiUtPiCY9IZqjtTpYxt2WvdXYUmIpnjHqn06FoE3LbQzaDS/ur7E6FL9hjOFIbQfzs/2rWgjGkQhE5DUROTjCbcNELiQiCcCzwNeNMa5P018Cs4EioBb479HON8Y8aIwpNsYUZ2RYUxe9ICdRSwRuGjv7ePdEI9cX5fpVn2hlvQU5SczPTuQ5HVz2sYaOPpq7+j+uZvYnYyYCY8wVxpjFI9y2AnUikgPg/Dli65CIROJIAo8ZY55ze+46Y8yQMcYOPASs9MQv5S3zsxM5Xtehg2WcXtxfy5DdaLWQGtGm5Tb2VLRy2rkaV6g74vwSOT87ABPBGLYBdzjv3wFsHX6AOCbc/i1wxBjzP8P2uU9Kswk4OMV4vGpedhJ9g3ZON+kbGxzVQvOzEz9uP1HKnWNNCtispQIAjjirlRfk+N//y1QTwfeBK0XkBHCl8zEikisirh5AFwC3A58aoZvoD0XkgIjsBy4DvjHFeLzKVbd3tFarh041drGnopVNuu6AGkV2cgxrZ6fpjKROR2vbyUmOISUuyupQPmFKw0CNMU3A5SNsrwHWO++/B4xYgWyMuX0q1/e1OZkJhIcJR8+2c+3S0J5hc3NpFWGiC9Coc9u0PI9vPr2P0opWzgvxDgVHajv8sn0AdGTxhMREhjMzPT7kG4ztdsNzexwzjWYl6UyjanTrFmcTExnG5j1VVodiqb7BIU42dPpljyHQRDBh87MTOXo2tLuQlpxpoaqlhxtWaGlAnVtCdARXLczmhf219A/arQ7HMifruxi0Gy0RBIsFOUlUNvfQ1jNgdSiW2byniriocK5epDONqrFtWm6jtXuAt46F7pQT/txQDJoIJmyxLRmAQzVtFkdijd6BIV7YX8u6xdnERelMo2psFxWmk54QFdIzkh6pbScqIoyCNP9cy1sTwQQtcSaCg9WhmQheO1JHR+8gNyzPszoUFSAiwsO4blkurx2pD9mS9MGaNhbkJBER7p8fuf4ZlR9LjY/ClhLLwerQbCfYXFpNdlIMa2anWR2KCiCbltvoH7Sz/UCt1aH4nN1uOFTdzhKbf7YPgCaCSVmUmxSSJYLGzj7eOt7AhuW5fjdplvJvS2zJzM6ID8nBZaebuujoG2SpLcXqUEaliWASltiSKW/soqM3tIq5zzvXotVqITVRIsKm5TY+OtVMZXO31eH41AHnl0ZX+6I/0kQwCYvzXA3GoVU99FxpNYtyk3RKCTUpG5xzUm3bF1ozkh6sbiM6IozCLP9alcydJoJJWJwbeg3GJ+o6OFDdxg0rtDSgJmd6ahwrC1J5dndVSE05sb/K0VAc6acNxaCJYFIyEqPJTooJqUTw3J5qwsOE65flWh2KCmA3FedR3thFaUWL1aH4hN1uOFTTztI8/60WAk0Ek7bYlvxx3V+wG7Ibtuyp5uLCdDISo60ORwWwa5fkEBcVztMloTHlxKmmLjr7Bv26fQA0EUzaYlsS5Y2OP3Kwe7+skdq2Xm48T6uF1NTER0ewfkkOL+yvpbs/+P93XLUGSzQRBKdl01MwBg5UBX+p4KmSSlLiIrlyYZbVoagg8Jni6XT2DfLSgbNWh+J1+6ucDcWZ/ttQDJoIJm35dEef4GCv62zp6ufVQ3VsLLIRHRFudTgqCJxfMI2CtDie3l1pdShet7+qlUW5/jui2GVK0YlIqoi8KiInnD9HnHBcRE47F6DZKyIlEz3fH6XERTErI549Fa1Wh+JVW/dW0z9k5zPF060ORQUJEeGm8/LYUd5MRVPwjinoH7Szv6qN5fn+/7E21TT1APC6MaYQeN35eDSXGWOKjDHFkzzf7yyfPo09FS1B3RXuTyVVLLYlsTDXf4fHq8Bzw4o8ROCZ0uBtND5S207foD0gFuSZaiLYADzivP8IsNHH51tqeX4KTV39VDb3WB2KVxysbuNwbTuf1dKA8rDclFguKszg2d1V2O3B+UVq9xlHtfGKECgRZBljagGcPzNHOc4Ar4jIbhG5dxLnIyL3ikiJiJQ0NDRMMWzPcP2Bg7Wd4OmSSqIiwrh+mS5Aozzv5vPyqG7t4YOTTVaH4hWlFS3YUmLJTvb/VfzGTAQi8pqIHBzhtmEC17nAGLMCuAb4iohcPNFAjTEPGmOKjTHFGRkZEz3dK+ZmJRAXFc6eIEwEvQNDbNlbw7pF2STHRVodjgpCVy7MIikmImgbjUvPtLA8338nmnM35soixpgrRtsnInUikmOMqRWRHGDEJYici9ljjKkXkc3ASuAdYFzn+6uI8DCW5aVQGoQNxq8crqOtZ0AbiZXXxESGs3G5jad2VdLWM0BybPB84aht66GmrZd7AqB9AKZeNbQNuMN5/w5g6/ADRCReRBJd94GrgIPjPd/fLc9P4UhtOz39Q1aH4lFPl1RiS4llra47oLzo5vOm0zdoZ1uQrV5Wesbx5TAQ2gdg6ong+8CVInICuNL5GBHJFZHtzmOygPdEZB/wEfCiMebP5zo/kKzIn8ag3bC/KnhKBVUt3bxX1sjNxXmE6boDyouW5CWz2JbEYzsrgqr3XWlFC9ERYX67WP1wU1p01hjTBFw+wvYaYL3zfjmwbCLnB5LigmmIwM5TzayaFRzfnp/a5aizvUmnlFA+cNvKGXx78wH2VLYGzDfosZRWtLAsL4WoCP8eSOYSGFH6sZS4KOZnJ7HzVHD0fBgYsvPkrko+NS+TvGlxVoejQsD1RbkkREfw2I4Kq0PxiK6+QQ5UtXFeQeAkNU0EHrBqZiq7z7TQP2i3OpQpe+1wHQ0dffzN6nyrQ1EhIiE6gg1Fubywv4a27sBf9W/3mRYG7YY1AVRDoInAA1bPSqN3wB4U7QR/3HkGW0osl8wddUiHUh5326p8+gbtPBsEI413lDcRESYBMaLYRROBB6ycmQo42gkCWXlDJ++XNXHryum6OL3yqUW5yRRNT+HxjwK/0XhHeRNL85KJj55SE6xPaSLwgNT4KOZnJ7KjPLDbCZ74qIKIMNGxA8oSt63Kp6y+k12nA3eAZlffIPur2lgTYN2uNRF4yKqZqZScbmFgKDDbCXoHhnh6dxVXLcoiM8n/h8Sr4HPd0lwSYyJ4bOcZq0OZNFf7wOoAah8ATQQes3pWGj0DQ+wP0IVqXjpYS2v3AH+7aobVoagQFRsVzo0r8njpwFmau/qtDmdSArF9ADQReMyqWWmIOJZ1DER/3FHBrPT4gCvSquBy26p8+ofsPLkrMLuSvl/WyLLpKcRFBU77AGgi8JjU+CiW2pJ557h/zIw6EQeq2th9poXbVuUjoo3EyjpzsxK5YE4af/jwTMBVszZ39bO/uo2LC/1jUsyJ0ETgQZfMzaC0oiXg+kL/7v1TxEeF85nztZFYWe+utTOpbevl5UOBtabxuycaMAYumaeJIKRdMi8Du4H3Aqh6qL69l+f313Bz8XSSYoJn9kcVuD41P5MZaXE8/N4pq0OZkLePN5ASF8kSW7LVoUyYJgIPWpaXQlJMBG8fD5zZtP+44wyDdsOdawusDkUpAMLChDvXFlBa0cq+ysAYpGm3G9453siFc9IDcgyOJgIPiggP46LCDN4+3hAQg2J6B4Z4bGcFl8/PpCA93upwlPrYTeflkRAdwe/eD4xSwZGz7TR29nHJ3MCrFgJNBB53ydwM6tr7OFbXYXUoY9q2r4amrn7uumCm1aEo9VcSYyK5uTiPF/bXUtfea3U4Y3rb2UnkYk0ECv7yRnjzqH/3HjLG8PB7p5iXlaiLzyi/dOfaAoaM4Q8f+v8As1cP17HYlkRWgA7G1ETgYdnJMSy2JfHKYf/u8fD28QaOnu3gCxfO1C6jyi/NSIvnygVZ/GHHGTr7Bq0OZ1T17b3sqWjl6oXZVocyaVNKBCKSKiKvisgJ589PDKcTkXkistft1i4iX3fu+66IVLvtWz+VePzFukXZ7Klo9esi7S/fOklOcgwbl9usDkWpUX3p0tm09QzwxE7/HWD26pE6AK5aFKKJAHgAeN0YUwi87nz8V4wxx4wxRcaYIuA8oBvY7HbIj1z7jTHbh58fiNYtdrwhXvHTftC7z7Sw81Qzd180K2BWUFKhaXn+NNbOTuM375XTN+if64K/fKiOGWlxzM1KsDqUSZvqp8AG4BHn/UeAjWMcfzlw0hjj/5V+UzAnM5HZGfH82U8TwS/fOklKXCS36AAyFQC+fOkc6tr7eK7U/xa4b+8d4MOTjVy1MCugq1inmgiyjDG1AM6fY61mcgvwxLBt94vIfhF5eKSqJRcRuVdESkSkpKHBvxtiwVEq2FHeTIufTZ51vK6D147UcceagoCaL12FrgvmpLE0L5lfv32SIbt/dct+40g9A0MmoKuFYByJQEReE5GDI9w2TORCIhIFXA887bb5l8BsoAioBf57tPONMQ8aY4qNMcUZGf7fRWvdohyG7Mbvhsn/8q2TxEaG6wAyFTBEhC9fOpvTTd1sP1BrdTh/ZeveanKTYzgvP7BmGx1uzERgjLnCGLN4hNtWoE5EcgCcP881pPYaoNQYU+f23HXGmCFjjB14CFg5tV/Hfyy2JTErPZ7Ne/ynOFtW38nWvdXcvmYG0+KjrA5HqXG7amE2szPi+dkbZX5TKmjq7OOdE41cV5RLWACOJnY31aqhbcAdzvt3AFvPceytDKsWciURp03AwSnG4zdEhE3Lbew81UxVS7fV4QDwv68dJyYynC9ePMvqUJSakLAw4etXzOVYXQcv7K+xOhwAXthfy5DdsCkIet5NNRF8H7hSRE4AVzofIyK5IvJxDyARiXPuf27Y+T8UkQMish+4DPjGFOPxK66umVv3Wv/GPXq2nRf213LXBQWkJURbHY5SE3btkhwW5CTxo1eP+8UU1Vv2VjM/O5H52UlWhzJlU0oExpgmY8zlxphC589m5/YaY8x6t+O6jTFpxpi2YeffboxZYoxZaoy53tXwHCymp8axcmYqz5ZWWT730I9ePU5idAT3XKSlARWYwsKEv79yLqebunlmd5WlsZxu7GJPRSsbigK/NAA6stjrblxho7yhi9IK62ZR3FvZysuH6vjCRTNJidO2ARW4Ll+QyfL8FH7y+gl6B6wbV/D4RxWEhwk3rNBEoMbh00tzSYyO4A8fnrbk+sYYvvfCYTISo7lbSwMqwIkI/3DVPGrbei2bg6h3YIinSyq5amFWwM4tNJwmAi+Lj47gxvPy2H7gLA0dfT6//vYDZ9l9poVvXjWXBB03oILA2jnpXDovg5+8foLGTiv+p2pp6R7gb1fP8Pm1vUUTgQ/cvmYG/UN2nvLxgty9A0P8x0tHmJ+dyE3n6ShiFTz+5dML6RkY4r9ePubza/9xxxlmpccH1ay9mgh8YHZGAhcVpvPYzgr6B33X2+G3752iqqWHf752YUCumqTUaGZnJHDn2gKeKqnkYHXb2Cd4SGlFC6UVrfzt6hkBPaXEcJoIfOQLFzoW5N68xze9HSqauvnpGye4amEWFxam++SaSvnS311RSGpcFP9n2yHsPhpk9os3HfN0fTbI5unSROAjl8zNYFleMj97s8zrfaCNMXxnywEiwsL4vxsWefVaSlklKSaSb61fwO4zLfxxp/cbjo+ddczTdefa4JunSxOBj4gIX/1UIZXNPV4fYLZtXw3vnmjkm1fNJSc51qvXUspKN66wcfHcDH7w0lGvj+D/8evHiYsK5441BV69jhU0EfjQ5QsyWZSbxI9fP+61PtC1bT3869ZDFE1P4fYgfMMq5U5E+H+bFgPwrecOeG3gZmlFC9sPnOWei2YF5Txdmgh8SET41jULqGzu4eH3T3n8+Yfshm88tZeBITv/+9kibSBWISFvWhwPrF/Auyca+e17nv+/Msbw/e1HSU+I5p4gnadLE4GPXViYzpULs/jZG2WcbfPsUpa/evskO8qb+e71iyhIj/focyvlz/52VT5XL8ri+y8dZU9Fi0ef+7nSaj463cw3riwM2rE4mggs8M/XLmDIbvj2Zs8VZd86Vs9/v3KMTy/N4ebz8jzynEoFChHhhzcuIzs5hvsf30OzhxaEauzs43svHmZFfgq3np/vkef0R5oILDAjLZ4HrpnPG0freXJX5ZSfr6y+k68+vod52Un84MalQdW/WanxSo6L5Oe3raCxs4+7H9k15XY4Ywzfeu4AXX2D/ODGpQG/5sC5aCKwyB1rCrhgThr/9/lDHKia/ICYyuZuPvfbnURHhvHQ584Lum5tSk3Esukp/O9ni9hT2crXn9w7pa7av36nnFcP1/GtaxZQmJXowSj9jyYCi4SFCT++ZTlp8dHc/egualp7JvwcVS3d3PrQDrr6h3jk8yvJmxbnhUiVCizXLMnhX65dyJ8PneX+x0vpG5x4yeDF/bX88M9HuXZpDnddUOD5IP2MJgILpSdE89s7i+nuG+Izv/6QM01d4z63tKKFjT//gPaeAf74hVUsyk32YqRKBZbPXziT7163kJcP1XHX73bRNIHJ6f588Cxfe3IPK/Kn8Z83hUZV65QSgYjcLCKHRMQuIsXnOG6diBwTkTIRecBte6qIvCoiJ5w/A3sF6EmYn53E4/espqtvkBt+8QFvHjvXss8wMGTnp6+f4LO//pC4qHCe+/JaluRpElBquDsvmMn/fGYZJWdauPYn7/HG0bpzHj8wZOfHr53gS4/tZrEtmYfvOp+4qNCoapWp9FoRkQWAHfg18E1jTMkIx4QDx3EsVVkF7AJuNcYcFpEfAs3GmO87E8Q0Y8w/jXXd4uJiU1LyiUsFtJMNnXzlsVKOnu3gyoVZ3HVBASsLUokId+Tqho4+XjpYy4PvlFPV0sN1y3L5t+sXBeXgFqU86WB1G197cg8nG7q4cE46t6+ZwcWFGcRGhQPQ3jvAa4fr+NXbJzle18nGoly+f+NSYiLDLY7c80RktzHmE1/ap5QI3J78LUZPBGuA7xpjrnY+/haAMeY/ROQYcKkxpta5kP1bxph5Y10vGBMBOKaNfuidch58p5yOvkGiIsLISoqmd8D+8VoGy6an8PXLC7lsfqbF0SoVOPoH7TzywWl+8145de19RIQJWUkxGGM4296L3cCczAT+8ep5XLUo2+pwvWa0ROCLco8NcO8jWQWsct7Pcq1T7EwGo366ici9wL0A+fnB2Z83JjKcr15eyN0XzeKNo/XsrWyhoaOPqIgwZmUkcHFhBgtyEkOizlIpT4qKCOOei2dx1wUFvFfWyK7TzdS29iIi5E2L5YI56RTPmBbUXUTPZcxEICKvASOlyO8YY7aO4xojvbITLoYYYx4EHgRHiWCi5weS2Khwrl2aw7VLc6wORamgEhEexqXzMrl0npao3Y2ZCIwxV0zxGlWA++TdeYBr+s06Eclxqxo6d0upUkopj/NF99FdQKGIzBSRKOAWYJtz3zbgDuf9O4DxlDCUUkp50FS7j24SkSpgDfCiiLzs3J4rItsBjDGDwP3Ay8AR4E/GmEPOp/g+cKWInMDRq+j7U4lHKaXUxHmk15CvBWuvIaWU8qbReg3pyGKllApxmgiUUirEaSJQSqkQp4lAKaVCXEA2FotIA3BmkqenA40eDMdTNK6J0bgmRuOaGH+NC6YW2wxjTMbwjQGZCKZCREpGajW3msY1MRrXxGhcE+OvcYF3YtOqIaWUCnGaCJRSKsSFYiJ40OoARqFxTYzGNTEa18T4a1zghdhCro1AKaXUXwvFEoFSSik3mgiUUirEBWUiEJGbReSQiNhFZNRuViKyTkSOiUiZc81k1/ZUEXlVRE44f07zUFxjPq+IzBORvW63dhH5unPfd0Wk2m3fel/F5TzutIgccF67ZKLneyMuEZkuIm+KyBHn3/xrbvs8+nqN9n5x2y8i8hPn/v0ismK853o5rr9xxrNfRD4QkWVu+0b8m/oorktFpM3t7/Ov4z3Xy3H9g1tMB0VkSERSnfu88nqJyMMiUi8iB0fZ7933ljEm6G7AAmAe8BZQPMox4cBJYBYQBewDFjr3/RB4wHn/AeAHHoprQs/rjPEsjkEgAN/FsTa0p1+vccUFnAbSp/p7eTIuIAdY4byfCBx3+zt67PU61/vF7Zj1wEs4VuVbDewc77lejmstMM15/xpXXOf6m/oorkuBFyZzrjfjGnb8dcAbPni9LgZWAAdH2e/V91ZQlgiMMUeMMcfGOGwlUGaMKTfG9ANPAhuc+zYAjzjvPwJs9FBoE33ey4GTxpjJjqIer6n+vpa9XsaYWmNMqfN+B441L2weur67c71f3ON91DjsAFLEsfLeeM71WlzGmA+MMS3OhztwrBLobVP5nS19vYa5FXjCQ9celTHmHaD5HId49b0VlIlgnGxApdvjKv7yAZJljKkFxwcN4KkFTif6vLfwyTfh/c6i4cOeqoKZQFwGeEVEdovIvZM431txASAiBcByYKfbZk+9Xud6v4x1zHjO9WZc7r6A45uly2h/U1/FtUZE9onISyKyaILnejMuRCQOWAc867bZW6/XWLz63hpzzWJ/JSKvAdkj7PqOMWY8S17KCNum3Jf2XHFN8HmigOuBb7lt/iXwPRxxfg/4b+DzPozrAmNMjYhkAq+KyFHnN5lJ8+DrlYDjH/brxph25+ZJv14jXWKEbcPfL6Md45X32hjX/OSBIpfhSAQXum32+N90AnGV4qj27HS232wBCsd5rjfjcrkOeN8Y4/5N3Vuv11i8+t4K2ERgjLliik9RBUx3e5wH1Djv14lIjjGm1ln8qvdEXCIykee9Big1xtS5PffH90XkIeAFX8ZljKlx/qwXkc04iqXvYPHrJSKROJLAY8aY59yee9Kv1wjO9X4Z65iocZzrzbgQkaXAb4BrjDFNru3n+Jt6PS63hI0xZruI/EJE0sdzrjfjcvOJErkXX6+xePW9FcpVQ7uAQhGZ6fz2fQuwzblvG3CH8/4dwHhKGOMxkef9RN2k88PQZRMwYg8Db8QlIvEikui6D1zldn3LXi8REeC3wBFjzP8M2+fJ1+tc7xf3eD/n7OGxGmhzVmmN51yvxSUi+cBzwO3GmONu28/1N/VFXNnOvx8ishLH51HTeM71ZlzOeJKBS3B7z3n59RqLd99bnm799ocbjn/6KqAPqANedm7PBba7HbceRy+TkziqlFzb04DXgRPOn6keimvE5x0hrjgc/xDJw87/A3AA2O/8Y+f4Ki4cvRL2OW+H/OX1wlHNYZyvyV7nbb03Xq+R3i/AfcB9zvsC/Ny5/wBuPdZGe6956HUaK67fAC1ur0/JWH9TH8V1v/O6+3A0Yq/1h9fL+fhO4Mlh53nt9cLxpa8WGMDx2fUFX763dIoJpZQKcaFcNaSUUgpNBEopFfI0ESilVIjTRKCUUiFOE4FSSoU4TQRKKRXiNBEopVSI+/8B2MnvCRwroFAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def poly(x):\n",
    "    return 8 * x ** 4 - 8 * x**2 + 1\n",
    "xs = np.linspace(-1,1,1000)\n",
    "ys = poly(xs)\n",
    "plt.plot(xs, ys)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'tensorflow'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-14-99e444df5507>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mtensorflow\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mtf\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mtensorflow\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mkeras\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;31m# Helper libraries\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'tensorflow'"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "\n",
    "# Helper libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "model = keras.Sequential([\n",
    "    keras.layers.Dense(20, activation='relu', input_shape=[1]),\n",
    "    keras.layers.Dense(25, activation='relu'),\n",
    "    keras.layers.Dense(1)\n",
    "])\n",
    "\n",
    "optimizer = tf.keras.optimizers.Adam()\n",
    "\n",
    "model.compile(loss='mse',\n",
    "              optimizer=optimizer,\n",
    "              metrics=['mae', 'mse'])\n",
    "\n",
    "model.fit(xs, ys, epochs=100, verbose=0)\n",
    "plt.plot(xs,model.predict(xs))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w1, b1, w2, b2, w3, b3 = model.get_weights() # unpack weights from model\n",
    "\n",
    "def Relu(x):\n",
    "    return np.vectorize(lambda y: If(y >= 0 , y, RealVal(0)))(x)\n",
    "\n",
    "# reinterpret network symbolically using z3 variables.\n",
    "def net(x):\n",
    "    x1 = w1.T @ x + b1\n",
    "    y1 = Relu(x1)\n",
    "    x2 = w2.T @ y1 + b2\n",
    "    y2 = Relu(x2)\n",
    "    x3 = w3.T @ y2 + b3\n",
    "    return x3\n",
    "\n",
    "x = np.array([Real('x')])\n",
    "y_true = poly(x)\n",
    "y_pred = net(x)\n",
    "s = Solver()\n",
    "s.add(-1 <= x[0], x[0] <= 1) # constraint input to [-1,1] interval\n",
    "\n",
    "def Abs(x):\n",
    "    return If(x <= 0, -x, x)\n",
    "s.add(Abs( y_pred[0] - y_true[0] ) >= 0.5) # Is there an input with absolute error greater than 0.5?\n",
    "\n",
    "res = s.check()\n",
    "print(res)\n",
    "if res == sat:\n",
    "    m = s.model()\n",
    "    print(\"Bad x value:\", m[x[0]])\n",
    "    x_bad = m[x[0]].numerator_as_long() / m[x[0]].denominator_as_long() \n",
    "    print(\"Error of prediction: \", abs(model.predict(np.array([x_bad])) - cheb(x_bad)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more on neural verification, I would check out\n",
    "- https://sat-smt.in/programme.html  A 3 Part tutorial at SAT+SMT school by Guy Katz\n",
    "- https://github.com/sisl/NeuralVerification.jl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Takeaways\n",
    "\n",
    "- Z3 is fun.\n",
    "- Z3 is a powerful solver.\n",
    "- It can find satisfy constraints and solve puzzles\n",
    "- Exhaustively failing to find counterexamples = Proof\n",
    "- Mathematics\n",
    "- Applications of Z3 to Programming\n",
    "\n",
    "Z3 is not a panacea. It has some weaknesses in particular in regards to abstraction and induction. These are the strengths of interactive theorems provers like Coq, which Cody will talk about tomorrow. Nevertheless, working to Z3's strengths and making small compromises you can do astounding things with minimal effort.\n",
    "\n",
    "# Thanks for Coming!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "- rise4fun https://rise4fun.com/z3/tutorialcontent/guide\n",
    "- https://ericpony.github.io/z3py-tutorial/guide-examples.htm\n",
    "- Programming Z3 - https://theory.stanford.edu/~nikolaj/programmingz3.html\n",
    "- Nikolaj Bjorner's tutorial https://youtu.be/nGwyNmsxX6I\n",
    "- Hakank's examples http://www.hakank.org/z3/\n",
    "- Yurichev's book \"SMT by Example\" https://yurichev.com/writings/SAT_SMT_by_example.pdf\n",
    "- http://hackage.haskell.org/package/sbv\n",
    "- https://www.youtube.com/watch?v=ruNFcH-KibY Tikhon Jelvis - Analyzing Programs with Z3\n",
    "- https://www.youtube.com/watch?v=rvPWDgJc0O4&ab_channel=ACMSIGPLAN - Nadia Polykarpova on Z3\n",
    "- SAT SMT school https://sat-smt.in/\n",
    "- Emina Torlak's course https://courses.cs.washington.edu/courses/cse507/19au/calendar.html\n",
    "- Lindsey Kuper - SMT Solving and Solver-Aided Systems http://composition.al/CSE290Q-2019-09/\n",
    "- http://www.sc-square.org/CSA/school/lectures.html"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}