notebooks/boolean-satisfiability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differentiable Boolean Satisfiability\n",
    "\n",
    "Differentiable implementation of the [Boolean Satisfiability Problem](https://en.wikipedia.org/wiki/Boolean_satisfiability_problem), encoded in [Conjunctive Normal Form](https://en.wikipedia.org/wiki/Conjunctive_normal_form).\n",
    "\n",
    "$True$ is encoded as floating point $1.0$ and $False$ as floating point $0.0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.executing_eagerly()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[PhysicalDevice(name='/physical_device:GPU:0', device_type='GPU')]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.config.list_physical_devices('GPU')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://github.com/tensorflow/tensorflow/issues/29781#issuecomment-542496509\n",
    "def compute_hessian(f, *x):\n",
    "    grad_grads = []\n",
    "    with tf.GradientTape(persistent=True) as hess_tape:\n",
    "        with tf.GradientTape() as grad_tape:\n",
    "            y = f(*x)\n",
    "        grad = grad_tape.gradient(y, x)\n",
    "        for g in grad:\n",
    "            grad_grads += list(hess_tape.gradient(g, x))\n",
    "    hessian = tf.stack(grad_grads)\n",
    "    return hessian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differentiable AND gate\n",
    "\n",
    "Using floating point product as an approximation for $AND$ gate.\n",
    "\n",
    "$$0.0 * 1.0 = 0.0$$\n",
    "$$1.0 * 1.0 = 1.0$$\n",
    "$$ etc $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor([1. 0. 0. 0.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([1. 1. 0. 0.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([1. 0. 1. 0.], shape=(4,), dtype=float32)\n",
      "WARNING:tensorflow:Calling GradientTape.gradient on a persistent tape inside its context is significantly less efficient than calling it outside the context (it causes the gradient ops to be recorded on the tape, leading to increased CPU and memory usage). Only call GradientTape.gradient inside the context if you actually want to trace the gradient in order to compute higher order derivatives.\n",
      "tf.Tensor(\n",
      "[[0. 0. 0. 0.]\n",
      " [1. 1. 1. 1.]\n",
      " [1. 1. 1. 1.]\n",
      " [0. 0. 0. 0.]], shape=(4, 4), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "@tf.function\n",
    "def gate_and(x, y):\n",
    "    return x * y\n",
    "\n",
    "x = tf.Variable([1, 0, 1, 0], dtype=tf.float32)\n",
    "y = tf.Variable([1, 1, 0, 0], dtype=tf.float32)\n",
    "\n",
    "with tf.GradientTape(persistent=True) as tape:\n",
    "    z = gate_and(x, y)\n",
    "    \n",
    "print(z)\n",
    "print(tape.gradient(z,x))\n",
    "print(tape.gradient(z,y))\n",
    "print(compute_hessian(gate_and, x, y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differentiable NOT gate\n",
    "\n",
    "Since, valid values are constrained between $0.0$ and $1.0$. We can negate by subtracting from $1.0$\n",
    "\n",
    "$$1.0 - 0.0 = 1.0$$\n",
    "$$1.0 - 1.0 = 0.0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor([0. 1. 0. 1.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([-1. -1. -1. -1.], shape=(4,), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "@tf.function\n",
    "def gate_not(x):\n",
    "    return 1 - x\n",
    "\n",
    "x = tf.Variable([1, 0, 1, 0], dtype=tf.float32)\n",
    "\n",
    "with tf.GradientTape(persistent=True) as tape:\n",
    "    z = gate_not(x)\n",
    "    \n",
    "print(z)\n",
    "print(tape.gradient(z,x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differentiable OR gate\n",
    "We compose the OR gate using [De Morgan's Law](https://en.wikipedia.org/wiki/De_Morgan's_laws)\n",
    "\n",
    "$$ X \\lor Y = \\overline{\\overline{X} \\land \\overline{Y}} $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor([1. 1. 1. 0.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([-0. -0.  1.  1.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([-0.  1. -0.  1.], shape=(4,), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "@tf.function\n",
    "def gate_or(x, y):\n",
    "    not_x = gate_not(x)\n",
    "    not_y = gate_not(y)\n",
    "    \n",
    "    x_not_y_not = gate_and(not_x, not_y)\n",
    "    x_or_y = gate_not(x_not_y_not)\n",
    "    \n",
    "    return x_or_y\n",
    "\n",
    "x = tf.Variable([1, 0, 1, 0], dtype=tf.float32)\n",
    "y = tf.Variable([1, 1, 0, 0], dtype=tf.float32)\n",
    "\n",
    "with tf.GradientTape(persistent=True) as tape:\n",
    "    z = gate_or(x, y)\n",
    "    \n",
    "print(z)\n",
    "print(tape.gradient(z,x))\n",
    "print(tape.gradient(z,y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differentiable XOR gate\n",
    "We compose XOR gate from definition\n",
    "$$ X \\oplus Y = (\\overline{X} \\land Y) \\lor (X \\land \\overline{Y}) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:Calling GradientTape.gradient on a persistent tape inside its context is significantly less efficient than calling it outside the context (it causes the gradient ops to be recorded on the tape, leading to increased CPU and memory usage). Only call GradientTape.gradient inside the context if you actually want to trace the gradient in order to compute higher order derivatives.\n",
      "tf.Tensor([0. 1. 1. 0.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([-1. -1.  1.  1.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([-1.  1. -1.  1.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([0. 0. 0. 0.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([-3. -1. -1. -3.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([-3. -1. -1. -3.], shape=(4,), dtype=float32)\n",
      "tf.Tensor([0. 0. 0. 0.], shape=(4,), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "@tf.function\n",
    "def gate_xor(x, y):\n",
    "    not_x = gate_not(x)\n",
    "    not_y = gate_not(y)\n",
    "    \n",
    "    x_and_not_y = gate_and(x, not_y)\n",
    "    y_and_not_x = gate_and(y, not_x)\n",
    "    \n",
    "    x_xor_y = gate_or(x_and_not_y, y_and_not_x)\n",
    "    \n",
    "    return x_xor_y\n",
    "\n",
    "x = tf.Variable([1, 0, 1, 0], dtype=tf.float32)\n",
    "y = tf.Variable([1, 1, 0, 0], dtype=tf.float32)\n",
    "\n",
    "with tf.GradientTape(persistent=True) as tape:\n",
    "    z = gate_xor(x, y)\n",
    "    grads = tape.gradient(z,[x,y])\n",
    "    \n",
    "print(z)\n",
    "print(grads[0])\n",
    "print(grads[1])\n",
    "hessian = tape.gradient(grads[0],[x,y])\n",
    "print(hessian[0])\n",
    "print(hessian[1])\n",
    "hessian = tape.gradient(grads[1],[x,y])\n",
    "print(hessian[0])\n",
    "print(hessian[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Circuit Graph\n",
    "\n",
    "### Candidate Encoding\n",
    "\n",
    "Candidate is the input to the circuit. The satisfiability problem is considered solved when we find a candidate that makes the circuit output $True$.\n",
    "\n",
    "Candidate $X_1 = 0$, $X_2 = 1$ and $X_3 = 0$ is encoded as tensor.\n",
    "\n",
    "\\begin{equation*}\n",
    "Candidate = \n",
    "\\begin{bmatrix}\n",
    "0.0 & 1.0 & 0.0 \\\\\n",
    "\\end{bmatrix}\\end{equation*}\n",
    "\n",
    "### Circuit Encoding\n",
    "\n",
    "Circuit is encoded as a 2D tensor. e.g.\n",
    "\n",
    "$$\n",
    "(X_1 \\lor X_2 \\lor \\overline{X_3}) \\land\n",
    "(X_1 \\lor \\overline{X_2} \\lor X_3) \\land \n",
    "(\\overline{X_1} \\lor X_2 \\lor X_3)\n",
    "$$\n",
    "\n",
    "is encoded as\n",
    "\n",
    "\\begin{equation*}\n",
    "Circuit = \n",
    "\\begin{bmatrix}\n",
    "0.0 & 0.0 & 1.0 \\\\\n",
    "0.0 & 1.0 & 0.0 \\\\\n",
    "1.0 & 0.0 & 0.0\n",
    "\\end{bmatrix}\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algorithm\n",
    "\n",
    "To make the best use of GPU we follow this method.\n",
    "\n",
    "First we broadcast/tile the candidate to match the dimensions of the problem.\n",
    "\n",
    "\\begin{equation*}\n",
    "Candidate = \n",
    "\\begin{bmatrix}\n",
    "0.0 & 1.0 & 0.0 \\\\\n",
    "\\end{bmatrix}\\end{equation*}\n",
    "\n",
    "becomes\n",
    "\n",
    "\\begin{equation*}\n",
    "Candidate_{broadcasted} = \n",
    "\\begin{bmatrix}\n",
    "0.0 & 1.0 & 0.0 & 0.0 & 1.0 & 0.0 & 0.0 & 1.0 & 0.0 \\\\\n",
    "\\end{bmatrix}\\end{equation*}\n",
    "\n",
    "Next we flatten the circuit\n",
    "\n",
    "\\begin{equation*}\n",
    "Circuit = \n",
    "\\begin{bmatrix}\n",
    "0.0 & 0.0 & 1.0 &\n",
    "0.0 & 1.0 & 0.0 &\n",
    "1.0 & 0.0 & 0.0\n",
    "\\end{bmatrix}\\end{equation*}\n",
    "\n",
    "The XOR of the broadcasted candidate and the flat circuit gives us the activations at each element. e.g.\n",
    "\n",
    "$$\n",
    "(X_1 \\lor X_2 \\lor \\overline{X_3}) \\land\n",
    "(X_1 \\lor \\overline{X_2} \\lor X_3) \\land \n",
    "(\\overline{X_1} \\lor X_2 \\lor X_3)\n",
    "$$\n",
    "$$=$$\n",
    "$$\n",
    "(0.0 \\lor 1.0 \\lor 1.0) \\land\n",
    "(0.0 \\lor 0.0 \\lor 0.0) \\land \n",
    "(1.0 \\lor 1.0 \\lor 0.0)\n",
    "$$\n",
    "\n",
    "We shall now reshape the activations into 2D for the sake of clarity\n",
    "\n",
    "\\begin{equation*}\n",
    "Activations = \n",
    "\\begin{bmatrix}\n",
    "0.0 & 1.0 & 1.0 \\\\\n",
    "0.0 & 0.0 & 0.0 \\\\\n",
    "1.0 & 1.0 & 0.0\n",
    "\\end{bmatrix}\\end{equation*}\n",
    "\n",
    "Now we parallel reduce the $\\lor$ (OR) operations. We get\n",
    "\n",
    "\\begin{equation*}\n",
    "Activations = \n",
    "\\begin{bmatrix}\n",
    "1.0 \\\\\n",
    "0.0 \\\\\n",
    "1.0\n",
    "\\end{bmatrix}\\end{equation*}\n",
    "\n",
    "On reducing the $\\land$ (AND) operations, we get the output of the circuit\n",
    "\n",
    "\\begin{equation*}\n",
    "Activations = \n",
    "\\begin{bmatrix}\n",
    "0.0\n",
    "\\end{bmatrix}\\end{equation*}\n",
    "\n",
    "Since, the whole thing is differentiable end to end. We can optimize using SGD.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(0.0, shape=(), dtype=float32)\n",
      "tf.Tensor([ 1. -1.  1.], shape=(3,), dtype=float32)\n",
      "tf.Tensor(\n",
      "[[ 0.  0.  0.]\n",
      " [ 1. -1.  1.]\n",
      " [ 0.  0.  0.]], shape=(3, 3), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "@tf.function\n",
    "def circuit_forward_pass(candidate, problem_encoding):\n",
    "    disjunction_count = tf.shape(problem_encoding)[0]\n",
    "    \n",
    "    broadcasted_candidate = tf.tile(candidate, [disjunction_count])\n",
    "    reshaped_problem = tf.reshape(problem_encoding, [-1])\n",
    "    \n",
    "    resolved_values = gate_xor(broadcasted_candidate, reshaped_problem)\n",
    "    reshaped_values = tf.reshape(resolved_values, [-1, disjunction_count])\n",
    "    reshaped_values = tf.transpose(reshaped_values, [1, 0])\n",
    "    \n",
    "    conjunctions = tf.scan(gate_or, reshaped_values)[-1]\n",
    "    conjunctions = tf.reshape(conjunctions, [-1,1])\n",
    "    disjunctions = tf.scan(gate_and, conjunctions)[-1][0]\n",
    "    \n",
    "    return disjunctions\n",
    "\n",
    "candidate = tf.Variable([0, 1, 0], dtype=tf.float32)\n",
    "problem_encoding = tf.Variable([\n",
    "   [0, 0, 1], \n",
    "   [0, 1, 0], \n",
    "   [1, 0, 0], \n",
    "], dtype=tf.float32)\n",
    "\n",
    "with tf.GradientTape(persistent=True) as tape:\n",
    "    z = circuit_forward_pass(candidate, problem_encoding)\n",
    "    \n",
    "print(z)\n",
    "print(tape.gradient(z,candidate))\n",
    "print(tape.gradient(z,problem_encoding))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constraint/Boundary Loss function\n",
    "\n",
    "Since, the tensors are of type float, we need to constrain them between 0.0 and 1.0 and make sure the values stay close to these two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.5625, 0.0, 0.03515625, 0.0625, 0.03515625, 0.0, 0.5625]\n",
      "[-3.0, 0.0, 0.1875, 0.0, -0.1875, 0.0, 3.0]\n"
     ]
    }
   ],
   "source": [
    "@tf.function\n",
    "def constraint_loss_fn(candidate):\n",
    "    x = candidate\n",
    "    a = (x ** 2)\n",
    "    b = (x - 1) ** 2\n",
    "    \n",
    "    return a * b\n",
    "\n",
    "candidate = tf.Variable([-0.5, 0, 0.25, 0.5, 0.75, 1, 1.5], dtype=tf.float32)\n",
    "with tf.GradientTape() as tape:\n",
    "    z = constraint_loss_fn(candidate)\n",
    "    print(list(z.numpy()))\n",
    "    print(list(tape.gradient(z,candidate).numpy()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x29c008b9e48>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXRc9Xn/8fczo9FmLZYsWZZkW7Jjg21BwAvChARoCA4QEkKTtJAATZvWcQL5Jac9bdKk5fTXvf21OQnQsGRpAiShaRZCiAkmkAQT8CKvIMnYxpJtWbKsxdZirTPz/P6YK0XII2tkzejOHT2vc+boztzr0XNnrM/c+d7v/X5FVTHGGON9PrcLMMYYEx8W6MYYkyIs0I0xJkVYoBtjTIqwQDfGmBSR5tYvLioq0srKSrd+vTHGeNKuXbvaVbU42jrXAr2yspKamhq3fr0xxniSiBydaJ01uRhjTIqwQDfGmBRhgW6MMSnCAt0YY1KEBboxxqSImANdRPwiskdEnomyTkTkfhE5LCL7RWRNfMs0xhgzmakcoX8WqJ9g3U3Acue2EXhomnUZY4yZopgCXUQWAu8DvjHBJrcCj2nENmCuiJTGqca3eONkD/+yuZ7ewWAint4YYxLqq788xMuH2hPy3LEeoX8F+CsgPMH6cuD4mPtNzmNvISIbRaRGRGra2tqmVOiI4519PPLSEQ60dF/QvzfGGLf0DQX5ygsHqTnamZDnnzTQReQW4JSq7jrfZlEeO2fmDFV9VFXXqeq64uKoV65Oqqo8D4DaZgt0Y4y31Lf0oApVZfkJef5YjtCvBj4gIo3Ak8C7ReSJcds0AYvG3F8INMelwnEW5GVSkB2gtrkrEU9vjDEJU+fk1qqyvIQ8/6SBrqp/raoLVbUSuB14UVXvHLfZ08DdTm+X9UCXqrbEv1wQEarK8u0I3RjjObXN3czNDlCWn5mQ57/gfugisklENjl3NwNHgMPA14FPx6G2CVWV5XGwtYeh4ERN+sYYk3xqm7upKstDJFor9fRNabRFVf018Gtn+eExjytwTzwLO59VZXkMh5TDp3oT9tXFGGPiaTgU5o2TPXz86sqE/Q5PXik6ckLB2tGNMV7xZlsvQ6EwVQk8CPVkoC8pmkNWwG/t6MYYz6g9EckrC/Rx/D5hZWkudRboxhiPqG3uJjPgY0lRTsJ+hycDHSLt6HUt3YTD53R3N8aYpFPb3MWKBXn4fYk5IQoeDvSqsnx6B4Mc6+xzuxRjjDkvVaWupTuhzS3g6UC3K0aNMd5wvLOfnoFgwq4QHeHZQL+oJBe/T6hrsZ4uxpjkNtIjz47QJ5AZ8LN8fo4doRtjkl5dSzd+n3DxgtyE/h7PBjpEToxaoBtjkl1tczfLinPIDPgT+ns8HehVZfm09QxyqmfA7VKMMWZCtc1dM3JVu8cD3U6MGmOSW3vvIK3dgwlvPwePB/rK0sgLZBcYGWOS1cgBpx2hTyI/K8CiwiwLdGNM0hrt4VKa2C6L4PFAh8iLZIN0GWOSVW1zNwsLssjPDiT8d3k/0MvyaOzoo2dg2O1SjDHmHHXNib9CdIT3A92ZY7S+pcflSowx5q16B4M0dpxl1Qw0t0Bsk0RnisgOEdknIrUi8n+jbHOdiHSJyF7ndl9iyj2XjY1ujElWB1q6nUmhZ+YIPZYZiwaBd6tqr4gEgJdF5FlV3TZuu62qekv8Szy/+bkZFOWkW9dFY0zSGcmlkZaERJs00J3p5XqduwHnljRj1ooIK0vzrKeLMSbp1DZ3UTgnnQV5iZkUeryY2tBFxC8ie4FTwPOquj3KZlc5zTLPikjVBM+zUURqRKSmra1tGmW/VVVZPodO2aTRxpjkkuhJoceLKdBVNaSqlwMLgWoRuWTcJruBClW9DHgAeGqC53lUVdep6rri4uLp1P0WVc6k0Qdb7cSoMSY5DAXDHGztmdGJ7KfUy0VVzwC/Bm4c93i3qvY6y5uBgIgUxavIyYyccLBmF2NMsjh0qofhkCZ8DPSxYunlUiwic53lLOA9wIFx2ywQ5zuFiFQ7z9sR/3Kjq5w3hznpfuvpYoxJGiMHmKtKZ+4IPZZeLqXAd0TETySof6Cqz4jIJgBVfRj4MPApEQkC/cDtzsnUGeHzRU6MWk8XY0yyqG3uJivgZ0nRnBn7nbH0ctkPrI7y+MNjlh8EHoxvaVNTVZbHD3c1EQ4rvgROwmqMMbGoa+5mZWluQieFHs/zV4qOqCrL5+xQiKM2abQxxmXh8Mik0DPXfg4pFOirRsdGt3Z0Y4y7jnX20TsYnLErREekTKAvL8khzSfWjm6Mcd3oFaJ2hH5hMtL8LC/JtUA3xriutrmLNJ+wvCRnRn9vygQ6RE6M1jV3MYMdbIwx5hx1Ld0sm5/4SaHHS7lAb+8d4lTPoNulGGNmsdrm7hm9QnREigV6pL3Krhg1xrjlVM8AbT2DM95+DikW6CtLcwHr6WKMcc/vTojaEfq05GYGqJiXbSdGjTGuGb3k3wJ9+qrKbAgAY4x7apu7WFyYTV5m4ieFHi8FAz2fY519dNuk0cYYF9Q2d8/ogFxjpVygr7KhdI0xLukZGOZoR58r7eeQgoFeNToEgAW6MWZm1bdEJtmZqTlEx0u5QJ+fm0lxboYdoRtjZtxIDzs3uixCCgY6jJwYta6LxpiZVdvcTVFOOvNzM1z5/Skb6IdP9TIYDLldijFmFolcIZo/Y5NCjxfLFHSZIrJDRPaJSK2I/N8o24iI3C8ih0Vkv4isSUy5sVlVmk8wrBw82etmGcaYWWQwGOJQa49rJ0QhtiP0QeDdqnoZcDlwo4isH7fNTcBy57YReCiuVU5RlY2NboyZYYdaewmG1bUuixBDoGvEyKFuwLmNH87wVuAxZ9ttwFwRKY1vqbFbXJhNTkaa9XQxxsyY350QTeJABxARv4jsBU4Bz6vq9nGblAPHx9xvch4b/zwbRaRGRGra2toutOZJ+XzCqtI86los0I0xM6OuuZs56X4q583cpNDjxRToqhpS1cuBhUC1iFwybpNoZwDOGZRcVR9V1XWquq64uHjq1U7BqrI86lu6CYVtbHRjTOLVNnezsjTP1Unqp9TLRVXPAL8Gbhy3qglYNOb+QqB5WpVNU1VZHn1DIRo7zrpZhjFmFgiHlfqWblebWyC2Xi7FIjLXWc4C3gMcGLfZ08DdTm+X9UCXqrbEvdopuKQ80rF/f9MZN8swxswCb7b1cnYoRFW5OxcUjYjlCL0U+JWI7Ad2EmlDf0ZENonIJmebzcAR4DDwdeDTCal2Ci4qySU3M40dDafdLsUYk+J2NHYCcEVloat1pE22garuB1ZHefzhMcsK3BPf0qbH7xOuqCxkp/NCG2NMouxs6KQ4N4PKedmu1pGSV4qOuKKykMOnemnvtTlGjTGJoapsb+ikurLQtStER6R0oFcviXz9qbGjdGNMgjSd7qela2A0b9yU0oF+aXk+mQEf2xss0I0xibHDyRcL9ARLT/OxZnHB6AtujDHxtrOxk7zMNC4uyXW7lNQOdIi0o9e3dNuUdMaYhNjR0En1kkJXLygakfKBfuWSQsIKu45a90VjTHyd6hngSPtZ17srjkj5QF+9uIA0n1izizEm7nY617kkQ/s5zIJAz0r38/aF+Rboxpi429HQQVbAP3pluttSPtABqpfMY3/TGQaGbQYjY0z87Gg8zdqKAgL+5IjS5KgiwaqXFDAcUvYcs3FdjDHx0dU3zIGT3UnT3AKzJNDXVhQigjW7GGPipuZoJ6ruj98y1qwI9PysACsX5LGjscPtUowxKWJHQycBv7B68Vy3Sxk1KwIdImehdx09zVAw7HYpxpgUsL2hk8sWziUz4He7lFGzKtAHhsO8bhNHG2OmqW8oyOsnupKq/RxmUaCPtHPttHZ0Y8w07Tl2hmBYvRfoIrJIRH4lIvUiUisin42yzXUi0iUie53bfYkp98IV52awtHiOnRg1xkzb9oZOfAJrKwrcLuUtJp3gAggCf6Gqu0UkF9glIs+rat247baq6i3xLzF+rlxSyDP7WwiFFX8SjLtgjPGmHQ0drCrLIzcz4HYpbzHpEbqqtqjqbme5B6gHyhNdWCJcUVlIz0CQN072uF2KMcajhoJh9hw7Q3XlPLdLOceU2tBFpJLIdHTbo6y+SkT2icizIlI1wb/fKCI1IlLT1tY25WKna6S9a0eDdV80xlyY106cYTAYTrr2c5hCoItIDvAj4HOq2j1u9W6gQlUvAx4Anor2HKr6qKquU9V1xcXFF1rzBVtYkE353Cx2NtrIi8aYCzMyYc4VlcnVfg4xBrqIBIiE+XdV9cfj16tqt6r2OsubgYCIFMW10jipXlLI9oZOIvNaG2PM1Oxo6GTZ/Bzm5WS4Xco5YunlIsA3gXpV/fIE2yxwtkNEqp3nTcp2jeolhbT3DtLQftbtUowxHhMKK7saTydlcwvE1svlauAu4DUR2es89kVgMYCqPgx8GPiUiASBfuB2TdJD4JH+6DsaOllanONyNcYYL6lv6aZnMMiVXg10VX0ZOG8fP1V9EHgwXkUl0tuK5zBvTjo7Gjq5vXqx2+UYYzxkx2j7eXIG+qy5UnSEiFC9pJAdjXaBkTFmanY0dLKoMIuyuVlulxLVrAt0iLSjN53u58SZfrdLMcZ4hKqys7EzaY/OYZYGuo3rYoyZqjfbztJxdihp289hlgb6ytI8cjPSRvuTGmPMZEbaz6uXJN8VoiNmZaD7fcK6ygK7YtQYE7MdDR0U52ZQOS/b7VImNCsDHSKfsm+2naW9d9DtUowxSU5V2d7QSfWSQpxLbpLSLA70SDtYjfV2McZMoul0Py1dA1Qn8QlRmMWBfml5PpkBn7WjG2MmtbNxpP3cAj0ppaf5WL2owCa8MMZMakdDJ3mZaVxckut2Kec1awMdIp+2dS3ddA8Mu12KMSaJ7XDaz31JPjHOrA70K5cUogq7bDhdY8wETvUMcKT9bFJfUDRiVgf66sUFpPnEhgEwxkyoxjngS/b2c5jlgZ6V7ufShfnWjm6MmdCOhk6yAn4uKc93u5RJzepAh8in7v6mM/QPhdwuxRiThLY3dLK2ooCAP/njMvkrTLArlxQyHFL2HLd2dGPMW3X1DXPgZLcnmlvAAp21FYX4BF45bMMAGGPe6tUjHah6o/0cYpuCbpGI/EpE6kWkVkQ+G2UbEZH7ReSwiOwXkTWJKTf+8rMCXFFZyPN1rW6XYoxJMs/XtZKfFWBtRfJNCB1NLEfoQeAvVHUlsB64R0RWjdvmJmC5c9sIPBTXKhNsQ9UC3mjtodHmGTXGOIKhMC8caOX6FfM90X4OMQS6qrao6m5nuQeoB8rHbXYr8JhGbAPmikhp3KtNkA2rSgDsKN0YM2pn42nO9A2zoarE7VJiNqWPHRGpBFYD28etKgeOj7nfxLmhj4hsFJEaEalpa2ubWqUJtKgwm5WleWypO+l2KcaYJLGl7iQZaT6uuajY7VJiFnOgi0gO8CPgc6raPX51lH+i5zyg+qiqrlPVdcXFyfUibVhVQs3R0zacrjEGVWVLbSvvWl5Ednqa2+XELKZAF5EAkTD/rqr+OMomTcCiMfcXAs3TL2/m3LCqBFV4od6aXYyZ7epaujlxpp8bVnmnuQVi6+UiwDeBelX98gSbPQ3c7fR2WQ90qWpLHOtMuKqyPMrnZrGl1gLdmNluS20rInD9Sm8FeizfJa4G7gJeE5G9zmNfBBYDqOrDwGbgZuAw0Af8cfxLTSwR4YZVJXxvxzHODgaZk+Gdr1nGmPjaUtfKuooCinIy3C5lSiZNLVV9meht5GO3UeCeeBXllg1VJXz7lUa2Hmrjxks800nHGBNHxzv7qG/p5ks3r3S7lCnzRufKGVJdWUh+VsCaXYyZxbY43Ze91n4OFuhvkeb3cf3K+bxw4BTDobDb5RhjXPB83UkuLsmlsmiO26VMmQX6OBtWLaCrf5idNqSuMbPO6bND7Gjo9NTFRGNZoI9zzUVFZKT5Rr92GWNmjxcOnCKskQM7L7JAHyc7PY13LS9mS+1JIud6jTGzxZbak5TmZ3JJeZ7bpVwQC/QoNlSV0Nw1QG3z+AtijTGpqn8oxEuH2tiwqoTI5TfeY4EexfUr5uOTyKe1MWZ22HqojYHhMBuqvNncAhboUc3LyWBdZaG1oxszi2ypayUvM80zk1lEY4E+gQ2rSjhwsodjHX1ul2KMSbBgKMwL9a1cv7LEM2OfR+PdyhNs5Cy3DalrTOqrOXqa033Do3MjeJUF+gQWz8tmxYJca3YxZhbYUttKusfGPo/GAv08Nqwqoaaxkw4bI92YlKWqbKk7yTuXFXl+UD4L9PPYULWAsEYuNjDGpKb6lh6aTvd7vrkFLNDPq6osj7L8TBusy5gUtqXupCfHPo/GAv08RIQNVQvYeqiNvqGg2+UYYxJgS20raxcXUJzrrbHPo7FAn8SGVSUMBsO8dLDd7VKMMXF2vLOPupZuzw7GNV4sU9B9S0ROicjrE6y/TkS6RGSvc7sv/mW654olzhjp1n3RmJTz/OjY5969OnSsWE7pfht4EHjsPNtsVdVb4lJRkgn4fVy/Yj4v1J8iGAqT5uGLDowxb/V8XSsXleSwxINjn0czaTqp6kvArB4cfENVCV39w+xonNUvgzEp5fTZIXY0dnp2qNxo4nW4eZWI7BORZ0WkaqKNRGSjiNSISE1bW1ucfnXiXXNRMRlpvtGvZ8YY73vxwClCYU2Z9nOIT6DvBipU9TLgAeCpiTZU1UdVdZ2qrisu9s4VWZEx0ovYUttqY6QbkyK21J1kQV4ml5bnu11K3Ew70FW1W1V7neXNQEBEiqZdWZLZsGoBJ870U9diY6Qb43X9QyF+c7CNDVXeHfs8mmkHuogsEOcVEZFq5zk7pvu8yeb6lZEx0p973Xq7GON1L42MfZ5C7ecQW7fF7wOvAheLSJOIfEJENonIJmeTDwOvi8g+4H7gdk3Bdol5ORlcvayI/93VRDAUdrscY8w0/M/O4xTnZnDlUu+OfR7NpN0WVfWOSdY/SKRbY8q7c30Fn3x8Fy8cOMV7PTyriTGz2fHOPn71xik+83vLPD32eTSptTcJdv2K+ZTlZ/LEtqNul2KMuUDf3X4Mnwh3XLnY7VLizgJ9CtL8Pj565WK2HmrnSFuv2+UYY6ZoYDjE/+w8xg0rSyjNz3K7nLizQJ+iP7hiEQG/8MS2Y26XYoyZos2vtXC6b5i7rqpwu5SEsECfovm5mdx4SSn/u+u4jcBojMc8vu0oS4vn8I63zXO7lISwQL8Ad19VQc9AkKf3NrtdijEmRq+f6GLPsTPceWVFSvU9H8sC/QKsqyhgxYJcHnv1qF05aoxHPP7qUbICfj60dqHbpSSMBfoFEBHuXF9BXUs3e46fcbscY8wkuvqG+em+E3xwdRn5WQG3y0kYC/QLdNvqcnIy0nj8VevCaEyy++HuJgaGw9y5PjVPho6wQL9AczLS+NCacn6+v4WO3kG3yzHGTCAcVp7YdpS1FQVUlaXOQFzRWKBPw53rKxgKhflBTZPbpRhjJvDbN9tpaD/LXSl+dA4W6NOyvCSXq5bO44ltRwmF7eSoMcno8VePUjgnnZsuTf3hOizQp+muqyo4caafX79xyu1SjDHjNJ/p55f1rfzhFYvISPO7XU7CWaBP0w2rSijJy+BxG9/FmKTzve3HUOBjKThuSzQW6NMU8Pu4o3oxvznYxtGOs26XY4xxDAXDPLnzGNevmM/Cgmy3y5kRFuhxcEf1YnwifHe7je9iTLL4Re1J2nuHUr6r4lixTHDxLRE5JSKvT7BeROR+ETksIvtFZE38y0xuJXmZvLeqhB/UHGdgOOR2OcYY4IlXj1IxL5trlntn/uLpiuUI/dvAjedZfxOw3LltBB6aflnec+f6Cs70DfPM/ha3SzFm1jtwspsdjZ3ceWUFPl9qjtsSzaSBrqovAZ3n2eRW4DGN2AbMFZHSeBXoFVctncey+Tk8/mqj26UYM+s9/upRMtJ8fDiFx22JJh5t6OXA8TH3m5zHziEiG0WkRkRq2tra4vCrk4eIcNf6CvY1dbHPxncxxjU9A8P8ZM8J3n9ZGQVz0t0uZ0bFI9CjfZ+JepWNqj6qqutUdV1xceq1a922ppzsdL9NUWeMi36y5wR9QyHuTtFJLM4nHoHeBCwac38hMCsHCs/LDHDb6nKe3tfM6bNDbpdjzKyjqjz26lEuW5jP2xfOdbucGRePQH8auNvp7bIe6FLVWXtm8M71FQwGw/xwl43vYsxM23akk8OnemdVV8WxYum2+H3gVeBiEWkSkU+IyCYR2eRsshk4AhwGvg58OmHVesDK0jyqlxTyjZePcHbQpqgzZqaoKg+8eIiC7ADvv6zM7XJckTbZBqp6xyTrFbgnbhWlgM/feDEfeuhVHvnNm/z5hovdLseYWeH5ulZeebODv7+1isxA6o/bEo1dKZoAaysKef9lZTzy0hFOnOl3uxxjUt5gMMQ/ba5n+fwcPlo9O8ZticYCPUG+cNMKAP712QMuV2JM6vvOK40c7ejjb25ZRZp/9sba7N3zBCufm8Unr1nKz/Y1s+vo+a7LMsZMR3vvIA+8cJh3r5jPtRelXnfoqbBAT6BPXvs2SvIy+Puf1RG2CTCMSYj/3HKQ/uEQX7x5pduluM4CPYHmZKTx+RtXsK+pi5/sOeF2OcaknLrmbv5n5zHuuqqCZfNz3C7HdRboCfbBy8u5bGE+//7cAevGaEwcqSr/8EwdeVkBPnf9RW6XkxQs0BPM5xPue/8qWrsHeeQ3b7pdjjEpY0tdK68e6eDPb7iI/OyA2+UkBQv0GbC2opAPON0Ym073uV2OMZ43GAzxz9ZN8RwW6DPk8zetQAT+7RdvuF2KMZ737d9Guin+7SzvpjievRIzpHxuFhuveRs/29dMTaN1YzTmQrX1DPLAi5FuitfM8m6K41mgz6BN1y6NdGN8xroxGnOhvvz8GwwMh/jS+6yb4ngW6DMoOz3SjXG/dWM05oLUNnfx5M7j3H1VJW8rtm6K41mgz7APXl7OZYvm8m+/sG6MxkzFSDfF/KwAn71+udvlJCUL9Bnm8wn33bKKUz2DPGzdGI2J2XO1rWw70mndFM/DAt0FaysKuPXyMh61bozGxMS6KcbGAt0ln78x0o3xnzfXExlS3hgzkW9sbeBYp3VTnExMr4yI3Cgib4jIYRH5QpT114lIl4jsdW73xb/U1FI2N4vPvHs5m187yX//ttHtcoxJWq8cbufLzx/kxqoF1k1xEpPOWCQifuC/gBuITAi9U0SeVtW6cZtuVdVbElBjyvrUtW9j7/Ez/OPP61g2P8f+sxozztGOs3z6e7tZWjSH//eRt7tdTtKL5Qi9GjisqkdUdQh4Erg1sWXNDj6f8JU/vJyLSnK593u7OdLW63ZJxiSNnoFhPvGdGgC+8UfryM20E6GTiSXQy4HjY+43OY+Nd5WI7BORZ0WkKtoTichGEakRkZq2trYLKDf1zMlI4+t3ryPN7+NPv1NDV/+w2yUZ47pQWPnsk3tpaD/L1z66hop5c9wuyRNiCXSJ8tj4s3i7gQpVvQx4AHgq2hOp6qOquk5V1xUXW/PCiEWF2Tz0sTUc6+zjM9/fQzAUdrskY1z1788d4MUDp/i796/iHcuK3C7HM2IJ9CZg0Zj7C4HmsRuoareq9jrLm4GAiNi7MAVXLp3HP3zwEl462Ma/2DykZhb78e4mHvnNET525WLuuqrS7XI8JZZA3wksF5ElIpIO3A48PXYDEVkgIuIsVzvP2xHvYlPdHdWL+fg7Kvnmyw38YOfxyf+BMSlmz7HTfOHHr7F+aSF/94GoLbfmPCbt5aKqQRG5F3gO8APfUtVaEdnkrH8Y+DDwKREJAv3A7Wqdqy/I37xvJW+29fKlp15jafEc1lUWul2SMTOipaufjY/vYkFeJg99bC0B628+ZeJW7q5bt05rampc+d3JrqtvmA9+7bd09w/z03uvZmFBttslGZNQ/UMhPvLIKzS29/HjT7+Di0py3S4paYnILlVdF22dfQQmofzsAF+/ex1DoTB/9tguG8TLpDRV5S9/uI/a5m6+evvlFubTYIGepJbNz+GBO1bzxslu/uIH+2z8dJOyHnzxMM/sb+Gv3ruC61eWuF2Op1mgJ7HrLp7PF29eyS9qT/KVXx50uxxj4u4Xr7fwn88f5LbV5Wy6dqnb5XjepCdFjbs+8c4lHGzt4f4XD9PZN8R9t1SRnmafw8bbVJVvv9LIP/28nssXzeVffv9SnI5yZhos0JOciPDPt11KwZx0HvnNEWqbu3noY2tZkJ/pdmnGXJC+oSB//ePX+OneZm5YVcJ//sFlZAb8bpeVEuxQzwPS/D7++qaVfO1jazh4sodbHtjKtiPWzd94z9GOs/z+117h6X3N/OV7L+aRO9eSZ2O0xI0FuofcfGkpT91zNXmZAT72je18Y+sRG0vdeMaLB1p5/wMvc7J7gG//cTX3/N4yfD5rZoknC3SPWV6Sy0/vvZr3rJzPP/68nv/z5F76hqxbo0le4bDylV8e5E++XcOiwmx+du87udaGik4IC3QPys0M8PCda/mrGy/m5/ubue2/XqGh/azbZRlzjq6+Yf70sRq+8stDfGjNQn70qXewqNAulEsUC3SPEhE+fd0yvvMn1bT2DPCBB1/ml3WtbpdlzKj6lm7e/+DLbD3Uxj988BL+4yNvt5OfCWaB7nHvWl7Mz+59JxXzsvnTx2r4j+feYGA45HZZZhYLh5Uf7mritq/9lsFgiCc3XsVd6yusW+IMsLFcUsTAcIi/fep1/ndXE0U56Xz8HZXcub6CudnpbpdmZonBYIin9pzgkZeOcKTtLNVLCnnwo6uZn2tdbOPpfGO5WKCnEFXl1Tc7eOSlI/zmYBvZ6X5uv2Ixn3jXEsrnZrldnklRXf3DfG/7Mb712wbaegapKstj4zVLed+lpaTZiIlxZ4E+C9U1d/P1rUf42b5mFPjAZWVsvGYpK0vz3C7NpIiWrn6+9XID399xnN7BIO9aXsQnr3kbVy+bZ80rCWSBPoudODPyR3eMvqEQ115UzCevXcpVS+2PzjsHO8cAAAiiSURBVFyYN0728OhLR/jp3hMocMvbS/mzdy3lkvJ8t0ubFSzQDV19wzyx/Sj//dtG2nsHuaQ8j+tXlLCmooDLF80lP8u7V+sNDIfo7h+my7l1DzjLfcN09QdHH+sZGGYoGGY4pAyFwgw7t9HHgpH7wbDi9wnpfh8BvxDw+yK3NB/pfiE9zTf6WE5GGvlZAfKyApGfmZH7+VkB8rMDo8tZAb9nP0AHhkO8fqKL3cdOs/VQO1sPtZMV8POHVyziE+9cYt0QZ9i0A11EbgS+SmTGom+o6r+OWy/O+puBPuDjqrr7fM9pge6OgeEQP9lzgie2HaWupRtVEIFlxTmsWVzA6sVzWVNRwLLiHNeu4lNVzg6F6OgdpL13kLaeIdpHlwedZeexnkHODp2/V8+cdD/5WQFyMwNkBEbCOBLU6W8Jax/paYLfJ4TCjAb+SOgPhZTh4O8eGwyG6R2MfGD0DJz/4q6MNB9FORkU5WZQnJMeWc7JoCgnnaLcjNH7xTkZ5GWluRb+qkrT6X72HD/D7qOn2XPsNHUt3QyHIjmxuDCbD69dyF3rKyiYYyfc3TCtQBcRP3AQuIHIhNE7gTtUtW7MNjcDnyES6FcCX1XVK8/3vBbo7usdDLLP+cPdfew0e46f4UzfMAC5mWlcvmguqxcXUJafSW5mgNzMNOcWORLNzQyQGfBNGD7DoTD9wyEGhkL0DYXoH47cuvuH6egdouPsoPNziI7eQednJKgHg+Goz1mQHfhdGOY6gZiTMXqEnD/uSDkvKzAjU5mFwkrPwDDdzjeC8bfTfUO09wzSNvqhNETn2UGiDXMf8AuFc9KZNyeDeTnpzJuTzrycyHKR89jc7ACZAT9ZAT/Z6WlkBfxkpkc+lCZ6PwaDIXoGgs5tePRn90CQ9t7ByP+FY2do6xkEICvg5+0L81lTUcCaxZFvcsW5GYl8GU0MzhfosYy2WA0cVtUjzpM9CdwK1I3Z5lbgMWce0W0iMldESlW1ZZq1mwTKyUjj6mVFXL2sCIgcnTW0n2X3sTORgD92hgdfPBQ1dEak+WQ05IHRAO8fDhGMYVKOdL8vElo5kQBbNj+HopwMJ9DSnSPaDIpzI48l6zyTfp8wNzt9St1EQ2GNBH3vIO1jvoWMfsD1DtF+dojGjrN09A7RN8k3EQCfRII4K91PZsBPmk/oHQzSPRBkaIIPyREV87J557Ii1iyOfJCvWJBrvVQ8JpZALwfGTkHfROQofLJtyoG3BLqIbAQ2AixevHiqtZoEExGWFuewtDiHD69dCETmejzdNzTmaC5yZNc97iivZyCIAFnO0WJWui9y1OiES5ZzNJmZ7ic3I40i54gzJ8O95gW3+X0y+m2DBZNv3zcUHP1Gc6ZviIHhMAPOt55+50N0YMxy/1CI4bCSk5HmfKNKG/NN63ffsnIz05ibHRj9UDbeFUugR/trG3/oFcs2qOqjwKMQaXKJ4Xcbl2Wl+8lKtz7sySA7PY3swjQ7CWkmFMv3qSZg0Zj7C4HmC9jGGGNMAsUS6DuB5SKyRETSgduBp8dt8zRwt0SsB7qs/dwYY2bWpE0uqhoUkXuB54h0W/yWqtaKyCZn/cPAZiI9XA4T6bb4x4kr2RhjTDQxzSmqqpuJhPbYxx4es6zAPfEtzRhjzFRYnyRjjEkRFujGGJMiLNCNMSZFWKAbY0yKcG20RRFpA45e4D8vAtrjWI6bbF+SU6rsS6rsB9i+jKhQ1eJoK1wL9OkQkZqJBqfxGtuX5JQq+5Iq+wG2L7GwJhdjjEkRFujGGJMivBroj7pdQBzZviSnVNmXVNkPsH2ZlCfb0I0xxpzLq0foxhhjxrFAN8aYFOGJQBeRj4hIrYiERWTCrj4icqOIvCEih0XkCzNZY6xEpFBEnheRQ87Pggm2axSR10Rkr4gkzeSrk73GzhDK9zvr94vIGjfqjEUM+3KdiHQ578FeEbnPjTonIyLfEpFTIvL6BOu99J5Mti+eeE8ARGSRiPxKROqd/PpslG3i+96oatLfgJXAxcCvgXUTbOMH3gSWAunAPmCV27VHqfPfgS84y18A/m2C7RqBIrfrneprTGQY5WeJzGK1Htjudt3T2JfrgGfcrjWGfbkGWAO8PsF6T7wnMe6LJ94Tp9ZSYI2znAscTPTfiyeO0FW1XlXfmGSz0cmsVXUIGJnMOtncCnzHWf4O8EEXa5mqWF7j0QnDVXUbMFdESme60Bh45f/LpFT1JaDzPJt45T2JZV88Q1VbVHW3s9wD1BOZa3msuL43ngj0GE00UXWyKVFnNifn5/wJtlNgi4jscibXTgaxvMZeeR9irfMqEdknIs+KSNXMlBZ3XnlPYuW590REKoHVwPZxq+L63sQ0wcVMEJFfEn3u8y+p6k9jeYooj7nSJ/N8+zKFp7laVZtFZD7wvIgccI5e3BS3CcOTQCx17iYybkaviNwMPAUsT3hl8eeV9yQWnntPRCQH+BHwOVXtHr86yj+54PcmaQJdVd8zzadImomqz7cvItIqIqWq2uJ8tTo1wXM0Oz9PichPiDQRuB3oqTRh+KR1jv3jU9XNIvI1ESlSVa8NEOWV92RSXntPRCRAJMy/q6o/jrJJXN+bVGpyiWUy62TwNPBHzvIfAed8+xCROSKSO7IMbACinvWfYak0Yfik+yIiC0REnOVqIn8vHTNe6fR55T2ZlJfeE6fObwL1qvrlCTaL73vj9pngGM8W30bkk2wQaAWecx4vAzaPO2N8kEjvhS+5XfcE+zIPeAE45PwsHL8vRHpe7HNutcm0L9FeY2ATsMlZFuC/nPWvMUGvpGS4xbAv9zqv/z5gG/AOt2ueYD++D7QAw87fySc8/J5Mti+eeE+cWt9JpPlkP7DXud2cyPfGLv03xpgUkUpNLsYYM6tZoBtjTIqwQDfGmBRhgW6MMSnCAt0YY1KEBboxxqQIC3RjjEkR/x98oNqCadpG0wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs = tf.range(-1,2.1,0.1)\n",
    "ys = constraint_loss_fn(xs)\n",
    "xs, ys = xs.numpy(), ys.numpy()\n",
    "plt.plot(xs, ys)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization loss\n",
    "Loss to make sure that the output is equal to 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(8.25, shape=(), dtype=float32)\n",
      "tf.Tensor(-1.0, shape=(), dtype=float32)\n",
      "tf.Tensor([-12.   0.  12.], shape=(3,), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "@tf.function\n",
    "def loss_fn(candidate, output):\n",
    "    optimization_loss = tf.math.squared_difference(output, 1.0)\n",
    "    constraint_loss = constraint_loss_fn(candidate)\n",
    "    constraint_loss = tf.reduce_sum(constraint_loss)\n",
    "    return optimization_loss + constraint_loss\n",
    "\n",
    "candidate = tf.Variable([-1, 1, 2], dtype=tf.float32)\n",
    "output = tf.Variable(0.5, dtype=tf.float32)\n",
    "with tf.GradientTape(persistent=True) as tape:\n",
    "    z = loss_fn(candidate, output)\n",
    "    \n",
    "print(z)\n",
    "print(tape.gradient(z,output))\n",
    "print(tape.gradient(z,candidate))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training/Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1 [0.000300008513 0.9997 0.000300008513]\n",
      "1000 0.372614384 [0.202071369 0.809448719 0.202071369]\n",
      "2000 0.302438021 [0.317299336 0.773605466 0.317299306]\n",
      "3000 0.217134744 [0.467763871 0.900396705 0.467763722]\n",
      "4000 0.101519562 [0.686918437 1.02080131 0.686917782]\n",
      "5000 0.012053174 [0.915573061 1.00370216 0.915572584]\n",
      "6000 2.50427984e-05 [0.996464252 1.00000632 0.996464193]\n",
      "7000 2.62632721e-10 [0.999989033 1.00000441 0.999989033]\n",
      "8000 2.00159576e-11 [0.999997437 1.00000262 0.999997437]\n",
      "9000 7.20489597e-12 [0.99999845 1.00000155 0.99999845]\n"
     ]
    }
   ],
   "source": [
    "candidate = tf.Variable([0, 1, 0], dtype=tf.float32)\n",
    "problem_encoding = tf.Variable([\n",
    "   [0, 0, 1], \n",
    "   [0, 1, 0], \n",
    "   [1, 0, 0], \n",
    "], dtype=tf.float32)\n",
    "\n",
    "opt = tf.keras.optimizers.Adam(learning_rate=3e-4)\n",
    "\n",
    "@tf.function\n",
    "def train_step():\n",
    "    with tf.GradientTape() as tape:\n",
    "        z = circuit_forward_pass(candidate, problem_encoding)\n",
    "        loss = loss_fn(candidate, z)\n",
    "    grads = tape.gradient(loss, candidate)\n",
    "    opt.apply_gradients(zip([grads], [candidate]))\n",
    "    return loss\n",
    "\n",
    "for i in range(10000):\n",
    "    loss = train_step()\n",
    "    if i % 1000 == 0:\n",
    "        tf.print(i, loss, candidate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Result\n",
    "$X_1 = 1.0$, $X_2 = 1.0$ and $X_3 = 1.0$ is a valid candidate which we can verify manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor([1. 1. 1.], shape=(3,), dtype=float32)\n",
      "tf.Tensor(1.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "print(tf.round(candidate))\n",
    "z = circuit_forward_pass(tf.round(candidate), problem_encoding)\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "f1d83f9e8f5c8837e86d6304fd35fdd4149e22cdc219dfc99f4444b7e631426b"
  },
  "kernelspec": {
   "display_name": "Python 3.7.6 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}