{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2c9533a6",
   "metadata": {},
   "source": [
    "# **Solving a Sudoku problem**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cebb2eb7",
   "metadata": {},
   "source": [
    "***sudoku.ipynb***\n",
    "\n",
    "Sudoku: place numbers from 1 to 9 into a 9x9 grid such that no number repeats in any row, in any column, and in any 3x3 sub-grid. \n",
    "\n",
    "© Copyright 2025 Fair Isaac Corporation\n",
    "\n",
    "Licensed under the Apache License, Version 2.0 (the \"License\"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.\n",
    " \n",
    "Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an \"AS IS\" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.\n",
    "\n",
    "This example uses FICO® Xpress software. By running it, you agree to the Community License terms of the [Xpress Shrinkwrap License Agreement](https://www.fico.com/en/shrinkwrap-license-agreement-fico-xpress-optimization-suite-on-premises) with respect to the FICO® Xpress software. See the [licensing options](https://www.fico.com/en/fico-xpress-trial-and-licensing-options) overview for additional details and information about obtaining a paid license."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c858d7e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install the xpress package\n",
    "%pip install -q xpress"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "413d5489",
   "metadata": {},
   "source": [
    "## Problem description and formulation\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12f8b922",
   "metadata": {},
   "source": [
    "In classic sudoku, the objective is to fill a $n \\times n$ grid with digits so that each column, each row, and each of the nine $q \\times q$ subgrids that compose the grid (also called \"boxes\", \"blocks\", or \"regions\") contains all of the digits from 1 to $n$. The puzzle setter provides a partially completed grid, which for a well-posed puzzle has a single solution.\n",
    "\n",
    "Several formulations exist for this problem, where choosing the right variables is a fundamental step. Although all cells must contain integer numbers, using integer decision variables would make it hard to guarantee that they are different within a given block (row/column) with a mathematical programming formulation. Alternatively, a Constraint Programming formulation can be applied, as shown in the blog post [Solving Sudoku: A Constraint Programming Play](https://community.fico.com/s/blog-post/a5Q4W000001V7gdUAC/fico2486).\n",
    "\n",
    "In this example, we use binary variables $assign_{i,j,k}$ that indicate whether a value $k \\in \\{1,..,9\\}$ is assigned to a given cell $i,j \\in \\mathcal{N}$ of the grid (=1) or not (=0). Also, no objective function is needed: this is a __feasibility__ problem not an __optimization__ problem, subject to the following constraints:\n",
    "\n",
    "* Each cell can only have one value: \n",
    "$$\\sum_{k \\in \\mathcal{N}} assign_{i,j,k} = 1, \\qquad \\forall i,j \\in \\mathcal{N}$$\n",
    "\n",
    "* Assign values already in grid ($g_{i,j}$ has a positive value): \n",
    "$$assign_{i,j,k} = 1, \\qquad \\forall i,j \\in \\mathcal{N}, k = g_{i,j},  g_{i,j} > 0$$\n",
    "\n",
    "* Every number must appear once on every row:\n",
    "$$\\sum_{j \\in \\mathcal{N}} assign_{i,j,k} = 1, \\qquad \\forall i,k \\in \\mathcal{N} $$\n",
    "\n",
    "* Every number must appear once on every column:\n",
    "$$\\sum_{i \\in \\mathcal{N}} assign_{i,j,k} = 1, \\qquad \\forall j,k \\in \\mathcal{N} $$\n",
    "\n",
    "* Every number must appear once in every $q \\times q$ block:\n",
    "$$\\sum_{i,j \\in \\mathcal{Q}: n = i+q.h, m = j+q.l} assign_{n,m,k} = 1, \\qquad \\forall h,l \\in \\mathcal{Q}, \\forall k \\in \\mathcal{N} $$\n",
    "\n",
    "where $\\mathcal{N}=\\{1,..,9\\}$ and $\\mathcal{Q}=\\{1,..,3\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01a45272",
   "metadata": {},
   "source": [
    "## Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22d335b2",
   "metadata": {},
   "source": [
    "The input is a starting grid where the unknown numbers are replaced by zero. We start with a $9 \\times 9$ grid, and provide a $16 \\times 16$ option at the bottom of this notebook. \n",
    "\n",
    "|   |   |   |   |   |   |   |   |   |\n",
    "|---|---|---|---|---|---|---|---|---|\n",
    "| 8 |   |   |   |   |   |   |   |   |\n",
    "|   |   | 3 | 6 |   |   |   |   |   |\n",
    "|   | 7 |   |   | 9 |   | 2 |   |   |\n",
    "|   | 5 |   |   |   | 7 |   |   |   |\n",
    "|   |   |   |   | 4 | 5 | 7 |   |   |\n",
    "|   |   |   | 1 |   |   |   | 3 |   |\n",
    "|   |   | 1 |   |   |   |   | 6 | 8 |\n",
    "|   |   | 8 | 5 |   |   |   | 1 |   |\n",
    "|   | 9 |   |   |   |   | 4 |   |   |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5ab9befe",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xpress as xp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "import math\n",
    "\n",
    "# Initial data for the sudoku solver\n",
    "grid3x3 = \\\n",
    " [[8, 0, 0, 0, 0, 0, 0, 0, 0],\n",
    "  [0, 0, 3, 6, 0, 0, 0, 0, 0],\n",
    "  [0, 7, 0, 0, 9, 0, 2, 0, 0],\n",
    "  [0, 5, 0, 0, 0, 7, 0, 0, 0],\n",
    "  [0, 0, 0, 0, 4, 5, 7, 0, 0],\n",
    "  [0, 0, 0, 1, 0, 0, 0, 3, 0],\n",
    "  [0, 0, 1, 0, 0, 0, 0, 6, 8],\n",
    "  [0, 0, 8, 5, 0, 0, 0, 1, 0],\n",
    "  [0, 9, 0, 0, 0, 0, 4, 0, 0]]\n",
    "\n",
    "grid = grid3x3\n",
    "q = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc2bceef",
   "metadata": {},
   "source": [
    "## Model implementation and results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36488945",
   "metadata": {},
   "source": [
    "When passing sets, lists, or range objects to [prob.addVariables()](https://www.fico.com/fico-xpress-optimization/docs/latest/solver/optimizer/python/HTML/problem.addVariables.html), the result is a Python dictionary of variables, whose keys are tuples of indices. Variables $x$ are created this way.\n",
    "\n",
    "The constraints are then created and added to the problem directly by passing the corresponding expression as a list comprehension. Note that no objective is set as this is a feasibility problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aeeb2976",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main dimensions of the problem: q is the size of the qxq block (3x3 in the classic Sudoku game).\n",
    "\n",
    "n = q**2       # the size must be the square of the size of the subgrids\n",
    "N = range(n)   # set of numbers from 0 to n-1\n",
    "Q = range(q)   # set of numbers from 0 to q-1\n",
    "\n",
    "# Create a model\n",
    "prob = xp.problem()\n",
    "\n",
    "assign = prob.addVariables(N, N, N, vartype=xp.binary)\n",
    "\n",
    "# Constraint 1: each cell can only have one value\n",
    "prob.addConstraint(xp.Sum(assign[i,j,k] for k in N) == 1 for i in N for j in N)\n",
    "\n",
    "# Constraint 2: fix the cells in the starting grid\n",
    "prob.addConstraint(assign[i,j,grid[i][j] - 1] == 1 for i in N for j in N if grid[i][j] > 0)\n",
    "\n",
    "# Constraint 3a: Every number must appear once on every row\n",
    "prob.addConstraint(xp.Sum(assign[i,j,k] for j in N) == 1 for i in N for k in N)\n",
    "\n",
    "# Constraint 3b: ... and on every column\n",
    "prob.addConstraint(xp.Sum(assign[i,j,k] for i in N) == 1 for j in N for k in N)\n",
    "\n",
    "# Constraint 3c: Every number must appear once in every qxq block\n",
    "prob.addConstraint(xp.Sum(assign[i+q*h,j+q*l,k] for i in Q for j in Q) == 1 for h in Q for l in Q for k in N)\n",
    "\n",
    "prob.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36ecae50",
   "metadata": {},
   "source": [
    "Now we use *matplotlib* to visualize the solution in a $n\\times n$ grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fa6ea6d2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAGKCAYAAAASfgYQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQ4VJREFUeJzt3XlsTWvbP/CrhrYHxzE8pVVap6aiiOkVhJIiERoioaaklP8qipAYIoj5HyEhqEgJjhJDqYTWLIKoMUTSowg1Rd682h5Tq6frl++V326qPOup7j7d1z79fpL1svtKn+vc91r3dU973QGO4zhCRET0b9T7d/8PIiIiYKIgIiJXTBREROSKiYKIiFwxURARkSsmCiIicsVEQURErpgoiIjIVQOpouLiYr08ysrK5P/+7/+kZcuWEhAQUNVfQ0REBuC71n/99Ze0adNG6tWrVzOJYv369bJq1aqaiI+IiIzIz8+Xtm3buv6bgKq+wqPyiKKwsFAiIiLk5s2bEhYWJhaVlpbK+fPnJS4uTho0qHJOrFX+ECN6Hd26dZNHjx7Jr7/+Khb5Qzn6Q4ys67oT45s3b+R//ud/pKCgQH777TfXf1vl/4KgoCC9KkOS+E/ZyFe+fv0q//rXvyQ8PFwaNmwoFvlDjEVFRfonYmzatKlY5A/l6A8xsq7rToweVVk64GI2ERG5YqIgIiJXTBREROSKiYKIiFwxURARkSsmCiIicsVEQURErpgoiIjIFRMFERG5YqIgIiJXTBREROSKiYKIiGwnir///luWL18uv//+u/zyyy/SoUMHWb16tb4r3ZIrV65IfHy8vrsdL9HKyMgQS/Aa+P79++sbP1u1aiXjx4+X3NxcsaR9+/ZadpWv5ORksWLlypXfxRcdHS2WbN++XXr27Kkv7cM1cOBAOX36tFjz6tUrmT59up5Zg2e7R48ecuvWLbFqw4YNWt/z5s0Ta2/0nTdvnkRGRmo5Dho0SHJycmo1Bp+//3bjxo164+/du1e6d++uN9LMmTP1tbdz584VKz5+/Ci9evWSpKQkmTBhglhz+fJlbXCRLPCK46VLl8qoUaP0ddGNGzcWC3Bzo2Pg8fDhQxk5cqRMnDhRLMF9eO7cufLP1l4Tjbc1o1Hr1KmTdqjw7IwbN07u3r2rsVvw/v17GTx4sAwfPlyTWEhIiDx+/FiaN28uFuHe3LlzpyZga2bPnq3Pyr59+7Sjun//fhkxYoQ+23g7bW3w+RNw7do1vcnHjBlT3us8ePCgnnNhyejRo/Wy6syZM9983rNnj44sbt++LUOHDhUL0FhUhMYOI8jY2FixBIkhNDRUrMLItqK1a9dqZ+vGjRtmEgU6gO3atZO0tLTyn2HWwKIPHz7ItGnTZNeuXbJmzRqx5PPnz3L06FE5ceJE+XOMUW9mZqbWeW3F6/OpJwyjcMDHn3/+qZ/v378vV69eNd0o+wMcLAUtWrQQi0pKSrRnhBGataN00fNFzy0qKkobkBcvXohVGKGlp6friBdTUFacPHlS+vXrp6NFdFh69+6tDbFFGImjo4peujWlpaVax8HBwd/8HFNQaCdri89HFIsXL9bDUjAPXL9+fS0U9JDwgFL14DxzzGli6B8TEyMWYY0HJ2vNmDFDLBkwYICOxrp06aIngOH43yFDhujQ39KJbw8ePNDE8OXLF2nSpIkcP35cT6az4unTp9rjXbBggU6DYmoHU8mBgYGSmJgoViDJ3rlzp9bn/KsK9xzqGeu2Xbt2ldatW+uMy/Xr16Vjx45SZxLF4cOH5cCBA/LHH3/osPnevXvayKFHZ+mG8ifoIaFhq80ex8/avXu3jhpRz5ZUHMlivhqJA4uIuE9nzZolViCR4VnByPHIkSP6rGCdykqyQGcFI4p169bpZ4wocE/u2LHDzHONs6JTUlLk7Nmz3/XYLdm3b5+OvLEegc50nz59ZMqUKTqtXGcSxaJFi3RUMXnyZP2MnRHPnz/XXTxWbih/MmfOHDl16pTu0rJ6RC3qF4vFx44dE+uaNWsmnTt3lry8PLEEPXNPj7Jv377aI96yZYsuyFqAI5IrJy30iDHfbgUa2nfv3mnD64EZDTw7W7duleLiYm2YfQ3reOgEYHoRsy8o24SEBJ0arS0+X6P49OmT1Kv3bRioHPRIqOqw+wVJAlMQFy5cMLtwCFjgxLy1ZwODZVjofPLkiT6cluF5QcNmBaY9K2/PxjokRmdWxMXF6RQeRmaeC6MgTHvj7xaSREXYvYj7EDvKsrKydBNQnRlRYAcH1iQiIiJ06glb/DZt2qRDLWsNRsVe5bNnz/RmwmIxYrcw3YTpO+yOwLzm27dv9efYZoyFL0sNGhIFRovWtp3CwoUL9Z5Eg/b69WtZsWKFNhgY6luxZMkSnSLDfYc99qj3S5cuaeNhxfz583WjCqaeJk2apLsYU1NT9bICz0nlNTw0xvjeh6W1vaysLO0IYroRbRBmYbCmi68R1BqnmgoLC/GNOCc/P9/xRlFRkZOSkuJEREQ4wcHBTlRUlLNs2TKnuLjY8VZJSYmTkZGhf3rr4sWL+t9b+UpMTDQR449iw5WWluZ4y1PX+NNbWVlZ+rtyc3OdmlRT5ZiQkOCEhYU5gYGBTnh4uH7Oy8szFWNSUpITGRmpMYaEhDhxcXFOdnZ2jcRYk3WdmZnpxMTEOEFBQU50dLSTmppaIzHW5HNdWWxsrLZHlmI8dOiQtouo79DQUCc5OdkpKCjw+vei7a5qXTewkNU3b96sl2XDhg0z923xiizHVhG+BGg5VuyCsQ4bAfzB2LFj9fInGJlZM2nSJL18yedrFEREZBsTBRERuWKiICIiV0wURETkiomCiIhcMVEQEZErJgoiInLFREFERK6YKIiIyBUTBRERuWKiICIiV0wURETkiomCiIhcMVEQEZErr18zXlpaKl+/fhWLPHFZjc9fYkQd4/Aj1vU/P0bWdd2JsbS0tMr/NgCHUvzML9+2bZteOFsWRxvidK1GjRpVJ04iIvLhMdRTp06VwsJCadq0ac0mCg8c8o1jNnEkaHh4uFiEbH727FkZOXKkNGzYUCzyhxhx3Gbbtm3l5cuXetCURf5Qjv4QI+u67sT46tUr+f3336uUKLyeesK5x1YLwgPxMUbv6vjz58+s6zoQI+u67sTY4CfOrOdiNhERuWKiICIiV0wURETkiomCiIhcMVEQEZErJgoiInLFREFERK6YKIiIyBUTBRERuWKiICIiV0wURETkiomCiIjsJwq8sXLevHkSGRmp78IfNGiQ5OTkiFUbNmyQgIAAjdmK7du3S8+ePfUtkLgGDhwop0+fFotvrJw+fbq0bNlS67pHjx5y69YtsQKvz1++fLm+VRPxdejQQVavXi3VfMnyf8X69eulf//++nbXVq1ayfjx4yU3N1esuXLlisTHx0ubNm30ecnIyBBr2rdvr7FVvpKTk8WKlStXfhdfdHR0rcbg9dtja8Ls2bPl4cOHsm/fPr2p9u/fLyNGjJBHjx6Ze4U5EtjOnTu1UbYEr4ZGAuvUqZM2anv37pVx48bJ3bt3pXv37mLB+/fvZfDgwTJ8+HBNYiEhIfL48WNp3ry5WLFx40ZNuig/lBuS2MyZM/WV+nPnzhULLl++rA0ZkgUOn1m6dKmMGjVKn5fGjRuLFR8/fpRevXpJUlKSTJgwQSzC84zOgQfaIbwafOLEiWJJ9+7d5dy5c9V68+s/IlHglcZHjx6VEydOyNChQ8szaGZmpj6wa9asESs+fPgg06ZNk127dpmKC9Bzq2jt2rVafjdu3DCTKNAIt2vXTtLS0sp/hp67JdeuXdMEO2bMmPIe58GDB+XmzZtixZkzZ775vGfPHh1Z3L59u/wZsmD06NF6WYbOSkXobGEUGRsbK5Y0aNBAQkND6+7UE3pEyOjBwcHf/BzD/qtXr4ol6MWhAcFoxzKUZ3p6uvboMAVlxcmTJ6Vfv37aW0PD1rt3b026lmDa8/z583p6I9y/f1/vQ8sNHg6egRYtWvg6FL9WUlKisxkYAWF6x5LHjx/rbEtUVJR2Vl+8eFG3RhSYZ0Vjhnngrl27SuvWrbUHd/36denYsaNYgYb3zp07ptdOHjx4oGX55csXadKkiRw/fly6desmVjx9+lRHOQsWLNDpEpQlpnMCAwMlMTFRLFi8eLGe3og54Pr162vSxegMD6dFZWVlulaGKb2YmBhfh+PXsIZSUFAgM2bMEEsGDBigo8YuXbrImzdvZNWqVTJkyBCdJqutUwh9nigAaxPI4liPwMPZp08fmTJlig6lLcjPz5eUlBQ92rDyyMcS3Ej37t3THuaRI0e08cV8tpVkgUYNI4p169bpZ4wocLPv2LHDTKI4fPiwHDhwQM+Cx5QdyhMNMXpzVmKsPMpFGVobffuj3bt368gRdW3J6AqjWayNInFg4w/u1VmzZtWdRIE5QTRomCpBby4sLEwSEhJ0mGUBEta7d+80gXmgp4ldHVu3bpXi4mJNcL6GnrlnFNa3b1/tsW/ZskUX3y1AvVZOWhhFYo3KikWLFumoYvLkyfoZu7KeP3+uO42sJYo5c+bIqVOn9D7EZgaqPtQxFouPHTsm1jVr1kw6d+4seXl5dWeNoiLs2EBjgt0xWVlZuqhoQVxcnE7roHfpudAzxnQE/m4hSfy7HjySmBWYHqm8jRNrAegdWfHp0yepV+/bxwL1i7K0ArvakCQwtXjhwgVzGwL8ETZYYN3Ms4nBsg8fPsiTJ0+0rawtJkYUSAq4+TF1giyJXh3miLEt0QLMA1ae/0VSw3cBrMwLL1myRIeoERER+r0UTJ1cunRJy9aK+fPn62Ixpp4mTZqkO4lSU1P1srR7DGsSKEdMPWF78aZNm3Rq1NJ0E+oXOwVxb759+1Z/ji282ARiqUGr2Ot99uyZdqyw6I7ytQKdACQKjBhre9tpVSxcuFDvS3SoXr9+LStWrNDOC6bna41TTYWFhfgGkpOfn+9469ChQ05UVJQTGBjohIaGOsnJyU5BQYHXv7ekpMTJyMjQP2tabGysk5KSYibGpKQkJzIyUsswJCTEiYuLc7Kzs52a4Klr/OmtzMxMJyYmxgkKCnKio6Od1NTUGomxpsqxqKhI6zUiIsIJDg7W+3LZsmVOcXGxmRhRFz+60tLSTNX1xYsXfxhnYmKiqec6KytL48rNzXVqSkkNxpiQkOCEhYXpsx0eHq6f8/LyvP69aLurWtcm0id6l7j8CXrr1hbi/MHYsWP1sgo99M2bN+tllaVvibsZNmyYX8SKLytajjM9Pd3XIdhaoyAiInuYKIiIyBUTBRERuWKiICIiV0wURETkiomCiIhcMVEQEZErJgoiInLFREFERK6YKIiIyBUTBRERuWKiICIiV0wURETkiomCiIhcef2a8dLSUvn69atY5InLanz+EiPqGAfisK7/+TGyrutOjKWlpVX+twE4lOJnfvm2bdv0wpnROMYSJ201atSoOnESEZEPj/2dOnWqFBYWStOmTWs2UXgUFRXp0Ys43jA8PFwsQjY/e/asjBw5Uho2bCgW+UOMOFq1bdu28vLlSz3YxyJ/KEd/iJF1XXdifPXqlZ63XpVE4fXUE86YtVoQHoiPMXpXx58/f2Zd14EYWdd1J8YGP3E+OBeziYjIFRMFERG5YqIgIiJXTBREROSKiYKIiFwxURARkSsmCiIicsVEQURErpgoiIjIFRMFERG5YqIgIiJXTBRERGQ/UeAthtOnT5eWLVvqu/B79Oght27dEivat28vAQEB313JyclixcqVK7+LLzo6WizbsGGDxjlv3jyx9gZVxBQZGan346BBgyQnJ0csuXLlisTHx0ubNm20DDMyMsSS9evXS//+/fUNtK1atZLx48dLbm6uWIKjEpYvX65vUEU9d+jQQVavXi3VfKH2f8327dulZ8+e+oZXXAMHDpTTp09LbfL67bHeev/+vQwePFiGDx+u//EhISHy+PFjad68uViBRgI3lcfDhw/19cETJ04US7p37y7nzp2r1tshfVGmO3fu1AfAmtmzZ2sd79u3Txvi/fv3y4gRI+TRo0dmXqn/8eNH6dWrlyQlJcmECRPEmsuXL2tHCskCB+QsXbpURo0apWXYuHFjsWDjxo3aCO/du1efHXROZ86cqccnzJ07V6xo27atdqo6deqkSQzxjhs3Tu7evatx14YGFiqrXbt2kpaWVv4zZHhLkLwqQqWh9xEbGyuWIDGEhoaKdR8+fJBp06bJrl27ZM2aNWIJXrF99OhROXHihAwdOrR8tJaZmamNipV4R48erZdVZ86c+ebznj17dGRx+/bt8nL1tWvXrmmDO2bMmPKZg4MHD8rNmzfFkvj4+G8+r127Vu/FGzdu1Fqi8PnU08mTJ6Vfv37aO8eN1Lt3b21ArCopKdEeJnpyGPJbgpEYesBRUVHaEL948UIsQk8TDyd66dag94vRY3Bw8Dc/x9TE1atXfRaXv8PhONCiRQuxAlOK58+f15M64f79+1rHlhPw33//Lenp6TqixBRUbfH5iOLp06eaHRcsWKDDU0xJYNgXGBgoiYmJYg3mggsKCmTGjBliyYABA7TX1qVLF3nz5o2sWrVKhgwZolMolk4qw01+584dc3P+HigrPICYq+7atau0bt1ae5nXr1+Xjh07+jo8v1RWVqZrPphijomJESsWL16sJ3ViLa9+/fraCKO3jk6WNQ8ePND78suXL9KkSRM5fvy4dOvWre4kCtxEGFGsW7dOP2NEgcZtx44dJhPF7t27tceBnrslFXtBmPdH4sBi7OHDh2XWrFliQX5+vqSkpOgRkZV77JZgbQIjRqxHoAHp06ePTJkyRadNqHojSDzT1kZkeDYOHDggf/zxh07h3Lt3TxManm1rbU+XLl00PozMjhw5ovFhHai2koXPE0VYWNh3/7HoyWGe2Jrnz5/rYvGxY8fEumbNmknnzp0lLy9PrEBD++7dO214PdCLww6erVu3SnFxsTbMvob1JzyEGN6jx4l7NCEhQaf06OfMmTNHTp06pXWMRVlLFi1apKOKyZMn62fstsQzjh1b1hJFYGBg+Yi2b9++OiLfsmWLbgipE2sUGI5W3jaHOUP0hq3BgjvWUTyLX9YXjJ88eaKNnBVxcXE6hEbPyHNhNImhPv5uIUlUhN05KD/szMvKytKFT6oa7M5BksAUyYULF8xtUIFPnz5JvXrfNoG4BzHLYV1ZWZl2rGqLz0cU8+fP10UlTD1NmjRJdxykpqbqZa1ikCjQ07C47XThwoW6OwIJ9vXr17JixQq96TFlYmn+v/IcNRpjfH/G0tw1kgIaOgz3MSJDzxPz2Ng6aakjUHG0+OzZM022WCyOiIgQC9NNmNLB7jHU+9u3b/Xn2HqKjQEW4HnBmgTKC1NP2G66adMmnXa0ZMmSJTq1jDjxHR+U66VLl/Q+rTVONRUWFuJbKU5+fr7jrczMTCcmJsYJCgpyoqOjndTUVKcmlJSUOBkZGfqnt7KysvS/Nzc316lJNRVjQkKCExYW5gQGBjrh4eH6OS8vr0Zi9NQ1/qxpsbGxTkpKiqm6PnTokBMVFaVlGRoa6iQnJzsFBQWmYrx48aLWSeUrMTHRRF3/KDZcaWlpjpVyLCoq0nsvIiLCCQ4O1jpftmyZU1xcbCZGSEpKciIjI/V+DAkJceLi4pzs7GzHW2i7q1rXJrrGY8eO1csyfFnI2jc2K+8m8kfoGVmDkS0uy4YNG2b6frQcmwdGOps3b9bLst27d/s6BN+vURARkW1MFERE5IqJgoiIXDFREBGRKyYKIiJyxURBRESumCiIiMgVEwUREblioiAiIldMFERE5IqJgoiIXDFREBGRKyYKIiJyxURBRESuqvyacZymVPFEJRwRCThIw/N3a75+/aqnWCG+hg0bikX+EKOnfq3Ws7+Uoz/EyLquOzH+9ddfVf63ATiUoir/cOXKlbJq1Spv4iIiImMKCwuladOmNZMofjSiaNeunTx69EjCw8PFalbPzs7WQ4esZnV/iNFT1/n5+f/xhvIVfyhHf4iRdV13Ynz16pV069atSomiylNPQUFBev3olCjLN1SjRo00PquV5Q8xeiBG1vU/O0YP1vU/P8ain5he5GI2ERG5YqIgIiJXTBREROSKiYKIiFwxURARkSsmCiIicsVEQURErpgoiIjIFRMFERG5YqIgIiJXTBREROSKiYKIiOwniitXrkh8fLy0adNGAgICJCMjQyzBK9YRV8UrOjpaLNuwYYPGOW/ePLFi/fr10r9/f32RZKtWrWT8+PGSm5srlmzfvl169uxZ/lK8gQMHyunTp8WS9u3bf3c/4kpOThYr/v77b1m+fLn8/vvv8ssvv0iHDh1k9erVUsWXVdfqmQx4RiIjIzXOQYMGSU5Ojlh7y+v06dOlZcuWGmOPHj3k1q1btRpDld8e+9/08eNH6dWrlyQlJcmECRPEou7du8u5c+fKPzdoYKLofgg3+s6dO7XBs+Ty5cvamCFZlJaWytKlS/U1zHhVfePGjcWCtm3bapLt1KmTNmp79+6VcePGyd27d/UesFK/aIg9Hj58KCNHjpSJEyeKFRs3btSki/JDuaFhmzlzpvz2228yd+5csWL27Nlafvv27dOO6v79+2XEiBFmjk94//69DB48WIYPH64dlpCQEHn8+LE0b968VuMw0dqNHj1aL8uQGEJDQ8W6Dx8+yLRp02TXrl2yZs0aseTMmTPffN6zZ4+OLG7fvi1Dhw4VCzCyrWjt2rXa4N24ccNMokBjURESG3rssbGxYsW1a9c0wY4ZM6Z8FHTw4EG5efOmWPH582c5evSonDhxovz+w+xBZmam1rmF52fjxo16PkhaWlr5zzBKq5NTT/4AWRw9jqioKG2IX7x4IRahx46HE70i63BgCrRo0UIsQq89PT1dR7yYgrKopKREe8EYjWP6yQpM4Zw/f17+/PNP/Xz//n25evWqqQ4hRrWo4+Dg4G9+jukdxGrByZMnpV+/fjpaRKeqd+/e2gmsbSZGFNYNGDBAe79dunSRN2/e6JGwQ4YM0SEr5tutQKN2584dc3OsP1JWVqZzwxhWx8TEiCUPHjzQxPDlyxdp0qSJHD9+XE8CswjreQUFBTJjxgyxZPHixXowDtby6tevrw0yRmfoZFmBZxf1jLWTrl27SuvWrXXUc/36denYsaNY8PTpUx3dLFiwQKdq8Wxj6i4wMFASExNrLQ4miiqo2AvCvD8SBxa/Dh8+LLNmzRILcHRlSkqKnD179rsektWRDxKtlZ5bRegQ3Lt3T0c8R44c0QcS6ysWk8Xu3bv1/sRo1xI8GwcOHJA//vhDp+xQnugYIM7abOD+E6xNYDSG9QgktD59+siUKVN0OtRKh6pfv36ybt06/YwRBZ6bHTt2MFFY16xZM+ncubPk5eWJFbix3717pze6B3px2FG2detWPe8cD4IFc+bMkVOnTmlsWDy2Br01T4+yb9++2ovbsmWLbhCw5Pnz57rB4tixY2LNokWLdFQxefJk/YydOogXO98sJQqs7aATgOlFjIDCwsIkISFBp5gtCAsL+66DgtEP1lZqE9coqrlg/OTJE61EK+Li4nTKBD03z4WeCIb6+LuFJIFdREgSmMq5cOGCTxblqturQ6K1BgucmLf2LBhb8unTJ6lX79vmBfcgytIi7LrD84xdRllZWboQb8HgwYO/20KOdR/MaNSmBlYa3oq982fPnmnjhkXOiIgI8bWFCxfqbhhUzuvXr2XFihV602OIamm+tfJcP25+7L22sgaA6SZMRWCXCeJ9+/at/hxbJrGAaMGSJUt0Kgf3HfbYI95Lly5p42EJGlwkCvTOLW7VxvOCNQmUI6aesL1406ZNOs1jCeoVHRhMN6INwkgI6yrYymvB/PnzdWMApp4mTZqku8ZSU1P1qlVONRUWFuKbM05+fr7jrYsXL+rvqnwlJiZ69XtLSkqcjIwM/dMbCQkJTlhYmBMYGOiEh4fr57y8PK9+Z03H+COxsbFOSkqK17/HU9f40xs/qmNcaWlpZsoxKSnJiYyM1LoOCQlx4uLinOzsbK/jq8kYISsrS8suNzfXqUk1VddFRUV670VERDjBwcFOVFSUs2zZMqe4uNhUOR46dEhjQ32HhoY6ycnJTkFBgakYMzMznZiYGCcoKMiJjo52UlNTnZqAtruqdW2iKzJs2DBz39isvJvIH6EnbInlOq64OOwP8EVFy+WJEePmzZv1sgy9dFyWjR07Vi9f4hoFERG5YqIgIiJXTBREROSKiYKIiFwxURARkSsmCiIicsVEQURErpgoiIjIFRMFERG5YqIgIiJXTBREROSKiYKIiFwxURARkSsmCiIicuX1a8ZLS0vl69evYpEnLqvx+UuMqGMcLMS6/ufHyLquOzGWlpZW+d8G4FCKn/nl27Zt0wvnMeNIPpwA1qhRo+rESUREPjyudurUqVJYWChNmzat2UThgYPIcYQlji0NDw8Xi5DNz549KyNHjpSGDRuKRf4QI44Ebdu2rbx8+VIPpLHIH8rRH2JkXdedGF+9eqXn1lclUXg99YTzeq0WhAfiY4ze1fHnz59Z13UgRtZ13YmxwU+ctc7FbCIicsVEQURErpgoiIjIFRMFERG5YqIgIiJXTBREROSKiYKIiFwxURARkSsmCiIicsVEQURErpgoiIjIFRMFERHZThTbt2+Xnj176tsLcQ0cOFBOnz4tFt+0OH36dGnZsqW+r79Hjx5y69Ytsfbmz3nz5klkZKTGOGjQIMnJyREr8Gr65cuX6xsrEV+HDh1k9erVUs0XGNeKDRs2SEBAgJarFStXrtSYKl7R0dFiTfv27b+LE1dycrJYsH79eunfv7++JbdVq1Yyfvx4yc3NFWuuXLki8fHx0qZNGy2/jIyMWo/B67fHeguvNMbD2KlTJ20w9u7dK+PGjZO7d+9K9+7dxYL379/L4MGDZfjw4ZrEQkJC5PHjx9K8eXOxZPbs2fLw4UPZt2+f3lT79++XESNGyKNHj0y8Cn7jxo3aMUAdo26RaGfOnKmvq587d65YgyS7c+dO7chYg/I7d+5ctd4EWpvlh86BB+5NvHZ74sSJYsHly5c1aSFZ4BCfpUuXyqhRo/R5ady4sVjx8eNH6dWrlyQlJcmECRN8EoPP7y5kyorWrl2rjcmNGzfMJAo0cO3atZO0tLTyn6FXbAleDX306FE5ceKEDB06tLznmZmZqeW5Zs0aX4co165d007AmDFjynucBw8elJs3b4o1Hz58kGnTpsmuXbtMlF1lSAyhoaFiGTpUFaFDiFFkbGysWHDmzJlvPu/Zs0dHFrdv3y5/hiwYPXq0XnV66qki9D7S09M1g2IKyoqTJ09Kv379tCeEG6l3797agFiCHhHKLzg4+JufY4rn6tWrYgGmws6fP68nI8L9+/c1Nl8/BD+CniYSGkZkFmFEi1FjVFSUJrQXL16IZSUlJTrCRa8Y0ycW4QAfaNGiha9DMcfnIwp48OCBJoYvX75IkyZN5Pjx49KtWzex4unTp9orX7BggQ5PMaTGVElgYKAkJiaKBZhnRRlizr9r167SunVr7a1fv35dOnbsKBYsXrxYT0bEfHr9+vU1sWEEiYbOEnRW7ty5Y2p9p6IBAwZo77dLly7y5s0bWbVqlQwZMkSndqyeSod59YKCApkxY4ZYVFZWputQmGKOiYnxdTjmmEgUuOHv3bunGf3IkSPa+GL+0EqywE2EEcW6dev0M0YUeCh37NhhJlEA1ibQY8N6BBriPn36yJQpU3QobcHhw4flwIEDes46phVR53g40TO2Uo75+fmSkpKix1hWHp1ZUXEEhvUTJA5sYED5zpo1SyzavXu3xo26tggjSDzTVkbf1phIFOiZe3q9ffv21Z7cli1bdCHRgrCwsO+SFnrtWBOwBPO/SLCYukPPHXEnJCTo9IQFixYt0lHF5MmT9TN2jj1//lx3n1hJFEiq79690yTrgZEPdp5s3bpViouLNQlb0qxZM+ncubPk5eWJRahjLLwfO3ZMLJozZ46cOnVK6xiba8j4GkXFHjweSCswHK28bQ7z7OjFWYQdG0gS2K2VlZWlC8gWfPr0SerV+/aWQ6OL+rYiLi5Op0Ix2vFcGE1iegx/t5YkPAvvT5480Tq3CJtAsLbn2cRgBXZZIklgqvvChQvmNqhY4vMRxZIlS3RIGhERod8DwLTEpUuXtIGzYv78+boQi6mnSZMm6S6d1NRUvSxBmeHmx1QeepfowWM9AFtQLcAON6xJoK4x9YQt0Js2bdLpMiswx195jhqJF9+fsTJ3vXDhQi1LdFRev34tK1as0ASGaUZr0AlAosCI0doWXkw3ob3BTkHU+9u3b/Xn2K6NTSCWOgJ5FUaLz549004LFt3xLNUKp5oKCwvxLSknPz/f8UZSUpITGRnpBAYGOiEhIU5cXJyTnZ3t1ISSkhInIyND//RWZmamExMT4wQFBTnR0dFOamqquRgPHTrkREVFaVmGhoY6ycnJTkFBgde/11PX+NMbRUVFTkpKihMREeEEBwdrrMuWLXOKi4tNlWNlsbGxGreVGBMSEpywsDCt5/DwcP2cl5fn1ISaqmuPrKws/X25ublOTampckRcP7rS0tLMxAgXL178YZyJiYmON9B2V7WuG1hY5PIHY8eO1csyjHZwWYVe2+bNm/XyJxjhWtuV5S/wBTar37y3Gldlw4YN83msJtcoiIjIDiYKIiJyxURBRESumCiIiMgVEwUREblioiAiIldMFERE5IqJgoiIXDFREBGRKyYKIiJyxURBRESumCiIiMgVEwUREblioiAiIldev2a8tLRUvn79KhZ54rIan7/EiDrGQS6s639+jKzruhNjaWlplf9tAA6l+Jlfvm3bNr1wjjCOA8UJUY0aNapOnERE5MOjiadOnSqFhYXStGnTmk0UHkVFRXpkII7lCw8PF4uQzc+ePSsjR46Uhg0bikX+ECOOqMWh8y9fvtTDhyzyh3L0hxhZ13UnxlevXuk54VVJFF5PPeEcXKsF4YH4GKN3dfz582fWdR2IkXVdd2Js8BNnmHMxm4iIXDFREBGRKyYKIiJyxURBRESumCiIiMgVEwUREblioiAiIldMFERE5IqJgoiIXDFREBGRKyYKIiJyxURBRES2E8X69eulf//++qbKVq1ayfjx4yU3N1cswSvVly9frm9axLv6O3ToIKtXr5Zqvnj3v6J9+/YSEBDw3ZWcnCyWXLlyReLj46VNmzYaX0ZGhlh8q+b06dOlZcuWWt89evSQW7duibW3vM6bN08iIyM1xkGDBklOTo5YsXLlyu/uxejoaLFm+/bt0rNnT317Kq6BAwfK6dOnxaoNGzZoWaLua5PXb4/11uXLl7UxQ7LAQRpLly6VUaNGyaNHj6Rx48ZiwcaNG/WG2rt3r3Tv3l0bjZkzZ+pr1ufOnSsWoJFAQvN4+PChvuJ44sSJYsnHjx+lV69ekpSUJBMmTBBr3r9/L4MHD5bhw4drgxESEiKPHz+W5s2biyWzZ8/WOt63b58m3f3798uIESP0ubHy2n88K+fOnavW20prC16pjsa3U6dO2vHDMz5u3Di5e/euxm9JTk6O7Ny5UxNbbfN5zZ05c+abz3v27NGRxe3bt2Xo0KFiwbVr1/TmGTNmTHnv/eDBg3Lz5k2xAg1aRbj5MfKJjY0VS0aPHq2XVegUtGvXTtLS0sp/hpGkJXgN+NGjR+XEiRPlzwh68JmZmdqhWbNmjViAxBAaGiqWYXRb0dq1a7UMb9y4YSpRfPjwQaZNmya7du3ySf36fOqpMhyiAS1atBArMKw/f/68nugH9+/fl6tXr5pt8EpKSrSHiV47hqlUdSdPnpR+/frpSAwdlt69e+vDaQlG3hg9BgcHf/NzTEHhvrQCIzGMdqKiorSRe/HihViGMk1PT9dRL6agLElOTtaOKkaNvuDzEUVFZWVlOveGoX9MTIxYsXjxYj3RD3Os9evX1xsKPQ/c/BZh3r+goEBmzJjh61D8ztOnT7VHuWDBAp0GxXAf04uBgYGSmJgoFmA9Dw0Z1sm6du0qrVu31hHu9evXpWPHjmLBgAEDdHagS5cu8ubNG1m1apUMGTJEp8usnZz34MEDLc8vX75IkyZN5Pjx49KtWzexIj09Xe7cuePTNagG1rImbiRLvSI4fPiwHDhwQM8Hx3D03r17mtDQW7LSeFS0e/duHe0gPvr5zgpGFOvWrdPPGFHgntyxY4epusbaBEaMWI9A56VPnz4yZcoUnbK1oOJoG3PqSBxYeMezNGvWLLEEyQzPNGYzjhw5ovWMtVMLySI/P19SUlL0WNXKI8g6mSjmzJkjp06d0l0xWGCyZNGiRTqqmDx5sn7GLpjnz5/rji1LjQcgLiwgHjt2zNeh+KWwsLDvGgj02rEmYAnWn9CYYZoEo13EnZCQoNM8FjVr1kw6d+4seXl5Yg1Gi56RWN++fbXnvmXLFl049rXbt2/Lu3fvtCPggRkNtJNbt26V4uJi7Sj849cosNMASQLDvQsXLphbOIRPnz5JvXrfFhUqB71Pa7AIi7l1z8I7/RxMe1beno21KfSGLcLOQCQJ7NbKysrSTRcWYTH2yZMnGqt1eK7RAFsQFxenU2MY8XgujHgx7Y2/10aSMDGiwHQTpnSwgwNzl2/fvtWfY+spFues7IzAmkRERIROPWHr3KZNm3Tob+0GR6LAKMfiVkRPg1GxV/ns2TO94bF5AeXra/Pnz9fNC5h6mjRpku5sS01N1csSJAV0sjBtgvLEqBdraNi2bcHChQv1uUGCff36taxYsUIbNUyPWbJkyRKdJsO9h++moC26dOmSlq8Fv/7663frtegc4Ds+tbqO61RTYWEhvm3m5OfnO97A7/jRlZaW5nirpKTEycjI0D+9UVRU5KSkpDgRERFOcHCwExUV5SxbtswpLi42EyNkZWVp2eXm5jo1yVPX+NNbFy9e/GF9JyYmminHzMxMJyYmxgkKCnKio6Od1NRUr39nTcd46NAhvQ8DAwOd0NBQJzk52SkoKDBT1wkJCU5YWJjGFx4erp/z8vIca+WYlJTkREZGapwhISFOXFyck52dbSrGymJjY7U98hba7qrWtc+7nZa+3eyW1Tdv3qyXZfiiovXyHDZsmPkYx44dq5dlGO3gsgo7dfwBNn74m0uXLtX6/6bP1yiIiMg2JgoiInLFREFERK6YKIiIyBUTBRERuWKiICIiV0wURETkiomCiIhcMVEQEZErJgoiInLFREFERK6YKIiIyBUTBRERuWKiICIiV1V+zThOfKp46hOOXwQc9uH5uzVfv37V0+kQX8OGDcUif4jRU79W69lfytEfYmRd150Y//rrryr/2wAcSlGVf7hy5UpZtWqVN3EREZExhYWF0rRp05pJFD8aUbRr104ePXok4eHhYjWrZ2dn64E+VrO6P8Toqev8/Pz/eEP5ij+Uoz/EyLquOzG+evVKunXrVqVEUeWpp6CgIL1+dPqb5RuqUaNGGp/VyvKHGD0QI+v6nx2jB+v6nx9j0U9ML3Ixm4iIXDFREBGRKyYKIiJyxURBRESumCiIiMgVEwUREblioiAiIldMFERE5IqJgoiIXDFREBGRKyYKIiJyxURBRES2EwVeXx4QEPDNFR0dLZasX79e+vfvry9AbNWqlYwfP15yc3PFmitXrkh8fLy0adNGyzEjI0MsvgN/3rx5EhkZKb/88osMGjRIcnJyxIr27dt/dz/iSk5OFqs2bNigMaJcrdi+fbv07Nmz/OWCAwcOlNOnT4vFN6hOnz5dWrZsqfdjjx495NatW2LF33//LcuXL5fff/9d4+vQoYOsXr1aqvjS7xpT5bfH/jd1795dzp07V/65QQMTYZW7fPmyNhRIFqWlpbJ06VJ9fTBesd64cWOx4uPHj9KrVy9JSkqSCRMmiEWzZ8+Whw8fyr59+zSh7d+/X0aMGGHmdfVIWng4PRDryJEjZeLEiWIR4t25c6c2ypa0bdtWE1inTp20Udu7d6+MGzdO7t69q8+7Be/fv5fBgwfL8OHDNYmFhITI48ePpXnz5mLFxo0bNemi/FBuSGIzZ86U3377TebOnVtrcZhokZEYQkNDxaozZ85883nPnj06srh9+7YMHTpUrBg9erReVn3+/FmOHj0qJ06cKC83jCgzMzP1YVizZo2vQ9TGoiI0dujFxcbGijUfPnyQadOmya5du0yUXUUY2Va0du1areMbN26YSRRohHH2RlpaWvnP0HO35Nq1a5pgx4wZUz7iPXjwoNy8ebNuTT0Bsjh6l1FRUXrjv3jxQizDQR/QokULX4fiVzAaQ289ODj4m59jSH316lWxpqSkREc8GKFhascajHLRgGBEZhnqPD09XUe8mIKy4uTJk9KvXz8dLaLj17t3b026lgwaNEjOnz8vf/75p36+f/++Piu13SH0+YhiwIAB2kPv0qWLvHnzRo9bHTJkiA75sSZgTVlZmc4FY8gaExPj63D8CuoTDQXmWLt27SqtW7fW3tH169elY8eOYg3WeAoKCmTGjBliDRreO3fumFrfqezBgwda31++fJEmTZrI8ePH9UQ1K54+faqjnAULFuh0MsoS0zmBgYGSmJgoFixevFgPGMK6bf369TXpYnSGDnWdShQVMyPmWZE4sNB5+PBhmTVrlljsxSGJWewB+wOsTaCHjvUI3Ph9+vSRKVOm6DSeNbt379b7E6NdS3BMaUpKipw9e/a70Zkl6Pzdu3dPR+BHjhzRxhfrfVaSBTp9GFGsW7dOP2NEgWd7x44dZhLF4cOH5cCBA/LHH3/olB3KEx1V3JO1GaPPE0VlzZo1k86dO0teXp5YM2fOHDl16pTuLsJiHf08zPejscA0BHpKYWFhkpCQoNOOljx//lw3WBw7dkysQVJ99+6dJlkP9DRxX27dulXPtkcS9jX0zD0jxb59+2qPfcuWLbr4bgHuvcpJCyNdrKNZsWjRIh1VTJ48WT9jVxbuTezErM1EYWKNovIC3ZMnT7QSrcCuDSQJDJ0vXLhgbsHLH2G3GOoYO0+ysrJ0wc4SLHBi3tqziGhJXFycTuugd+m50DPGdAT+biFJ/LsePJKYFZg+rrzNHWsBmNGw4tOnT1Kv3rfNNOoXZVmbfD6iWLhwoe6QQOW8fv1aVqxYoQWB6QhL000Y+mG3DubZ3759qz/HFjUsxFpKshVHYs+ePdOGA4vuERERYgGSAhIvpiUQK3pMmH/Flj8r8BAiUaDHZm2rNuAerLw+hsSL7wJYWTdbsmSJTtvhvsN3Z/D8XLp0Sevfivnz5+tiMaaeJk2apDuJUlNT9bIiPj5e1yRQjph6wvbiTZs26fRtrXKqqbCwEN/4cPLz8x1vJCQkOGFhYU5gYKATHh6un/Py8pyaUFJS4mRkZOif3sB/54+utLQ0MzHCxYsXfxhnYmKiV7/XU9f401uHDh1yoqKitL5DQ0Od5ORkp6CgwFQ5ZmVl6X9vbm6uU5NqMsbKYmNjnZSUFK9/T03VdVJSkhMZGan1HBIS4sTFxTnZ2dmOtXLMzMx0YmJinKCgICc6OtpJTU01FWNRUZHWa0REhBMcHKzPzrJly5zi4mKvY0TbXdW6bmBh94Z1tf0tyOoaNmyY+VjRc8NlGb5Mab0cK0Nv3dpGAH8wduxYvaz69ddfZfPmzXr5krk1CiIisoWJgoiIXDFREBGRKyYKIiJyxURBRESumCiIiMgVEwUREblioiAiIldMFERE5IqJgoiIXDFREBGRKyYKIiJyxURBRESumCiIiMiV168ZLy0tla9fv4pFnrisxucvMaKOcUAT6/qfHyPruu7EWFpaWuV/G4BDKX7ml2/btk0vnNGLYwNxclWjRo2qEycREfnwmNWpU6dKYWGhNG3atGYThUdRUZEeBYrjNsPDw8UiZPOzZ8/KyJEjpWHDhmKRP8SIoyzbtm0rL1++1INULPKHcvSHGFnXdSfGV69eye+//16lROH11BPOFLZaEB6IjzF6V8efP39mXdeBGFnXdSfGBj9xHjwXs4mIyBUTBRERuWKiICIiV0wURETkiomCiIhcMVEQEZErJgoiInLFREFERK6YKIiIyBUTBRERuWKiICIiV0wURETkX4liw4YNEhAQIPPmzRMrtm/fLj179tQ3LOIaOHCgnD59WizBa9+XL1+ub4PEeQIdOnSQ1atXSzVfDvxfsXLlSq3bild0dLRYsn79eunfv7++ObVVq1Yyfvx4yc3NFYtv/pw+fbq0bNlS67tHjx5y69YtsfYmWjzHkZGRGuOgQYMkJydHLLly5YrEx8dLmzZt9H7MyMgQa9q3b//dc4MrOTm51mLw+u2xNQk30c6dO7VRtgSvXUYC69Spkza8e/fulXHjxsndu3ele/fuYsHGjRs1oSE2xIRGY+bMmfoq+Llz54oViO3cuXPVeoNlbbh8+bI+gEgWONhl6dKlMmrUKHn06JE0btxYLHj//r0MHjxYhg8frh2WkJAQefz4sTRv3lwsmT17tjx8+FD27dunDfH+/ftlxIgRWpZWjib4+PGj9OrVS5KSkmTChAliUU5OjnYEPVCmeH35xIkTay0GM0/phw8fZNq0abJr1y5Zs2aNWIIeR0Vr167VRvnGjRtmEsW1a9c0eY0ZM6a8F3Lw4EG5efOmWILEEBoaKladOXPmm8979uzRkcXt27dl6NChYqVT0K5dO0lLSyv/GUaSluBV5UePHpUTJ06UlxtGlJmZmfrsWHnGR48erZdlISEh33xGpxUzBrGxsXVv6gm9ODRy6HFYhsyenp6uPRFMQVmBYf358+f11EG4f/++XL161dxDgJ4vepdRUVHaMXjx4oVYhkNdoEWLFmLFyZMnpV+/ftqjRBLr3bu3drAswWgMz0pwcPA3P8cUFO5Lqp6SkhIdmWEEhOmnOjWiQMN7584dc/OXFT148EATw5cvX6RJkyZy/Phx6datm1ixePFiPXUQc/7169fXhxQjHzTGVgwYMEB76F26dJE3b97IqlWrZMiQITqUtniaWllZmc6xY5onJiZGrHj69Kn2yhcsWKBTY3huML0YGBgoiYmJYgHqE88L1sm6du0qrVu31hHu9evXpWPHjr4Oz29lZGRIQUGBzJgxo1b/d32eKPLz8yUlJUWPDazc+7AEjdu9e/e0h3nkyBF9IDGfbSVZHD58WA4cOKBnmGM6DLGikUPv3UrjUXF0g3UoJA4sdCL2WbNmiTUY5SKJWesBI4FhRLFu3Tr9jBEF4tyxY4eZugasTaDni/UIdF769OkjU6ZM0Wk8qp7du3frc4Tnuk4lCtw0796905vIA71h7EbYunWrFBcX603ma+iteXpCffv21V7cli1bdPHdgkWLFumoYvLkyfoZu2CeP3+uu3gsNR4VNWvWTDp37ix5eXlizZw5c+TUqVN6H2IzgyVhYWHfdVDQa8eagCWYR0dnCtO0GO0i7oSEBJ12pJ+H5xkbQY4dOya1zedrFHFxcTqtgx6w50JvCVMm+LuFJPHvenVIYlZ8+vRJ6tX7tjpRdojTKmxgePLkiTYgVmBXG5IEphYvXLhgbpEYMBVWecsu1qYwOrMIu8VQx9itlZWVpZsu6Odh8wLWpDwbVurUiAJzmZXnf3FjYX+4lXnhJUuW6HAvIiJC94ZjeufSpUt601vamYU1CcSIqSds3d20aZMO/a1YuHChxokG7fXr17JixQpNZpiOsDTdhPrFbh3cm2/fvtWfY5sxFmItmD9/vm5ewNTTpEmTdGdbamqqXpbg+UDixbQtRo0Y9WINDdu2LXVWKo5onz17ph1UbF7As2RFWVmZJgrMDvhkS7lTTYWFhfgml5Ofn+/UtNjYWCclJcXr31NSUuJkZGTon95ISkpyIiMjncDAQCckJMSJi4tzsrOzvY6vJmMsKirSMouIiHCCg4OdqKgoZ9myZU5xcbHXMXrqGn96IyEhwQkLC9NyDA8P1895eXmOpXLEf+ePrrS0NDMxQmZmphMTE+MEBQU50dHRTmpqqlMTaqqu4dChQ3ofor5DQ0Od5ORkp6CgwFQ5Xrx48Yf1nZiYaCZGyMrK0rhyc3OdmoK2u6p17fMRxY+gt25tAck69H43b96sl1XY3WadpW+yuxk7dqxelmG0g8uyYcOG+UWdjxo1yqdx+nyNgoiIbGOiICIiV0wURETkiomCiIhcMVEQEZErJgoiInLFREFERK6YKIiIyBUTBRERuWKiICIiV0wURETkiomCiIhcMVEQEZErJgoiInLl9WvGS0tL5evXr2KRJy6r8flLjKhjHNrDuv7nx8i6rjsxlpaWVvnfBuBQip/55du2bdML51rj+EWcBtaoUaPqxElERD48Pnnq1KlSWFgoTZs2rdlE4YHD0nE8JI4ODA8PF4uQzc+ePSsjR46Uhg0bikX+ECOOf23btq28fPlSD0iyyB/K0R9iZF3XnRhfvXqlZ8JXJVF4PfWE81utFoQH4mOM3tXx58+fWdd1IEbWdd2JscFPnL3NxWwiInLFREFERK6YKIiIyBUTBRERuWKiICIiV0wURETkiomCiIhcMVEQEZErJgoiInLFREFERK6YKIiIyBUTBRER2U4U7du3l4CAgO+u5ORksfZWzXnz5klkZKS+r3/QoEGSk5Mj1t4GOX36dGnZsqXG2KNHD7l165ZYsX37dunZs6e+qRLXwIED5fTp02LNlStXJD4+Xtq0aaP3YkZGhliCV/wvX75c3/yJeu7QoYOsXr1aqvki6FqxYcMGLUs8Q5asXLnyu7YnOjpaLFm/fr30799f3+bbqlUrGT9+vOTm5tZqDF6/PdZbaGxx43s8fPhQX807ceJEsWT27Nka2759+7QB2b9/v4wYMUIePXpk4jXr79+/l8GDB8vw4cO18Q0JCZHHjx9L8+bNxQq8vhoNRqdOnbRR27t3r4wbN07u3r0r3bt3Fys+fvwovXr1kqSkJJkwYYJYs3HjRk26KD+UGzoDM2fO1Nf+z507V6zBM75z507tJFiEMjx37ly13qpaGy5fvqwdZyQLHDa0dOlSGTVqlLY9jRs3rpUYfF4iaNAqQkOCHlJsbKxYgdcuHz16VE6cOCFDhw4t74lkZmbqA7tmzRoTjUe7du0kLS2t/GfocVqCXnpFa9eu1fK7ceOGqUQxevRovay6du2aJtgxY8aUj8oPHjwoN2/eFGs+fPgg06ZNk127dpl4Tn4EiSE0NFSsOnPmzDef9+zZoyOL27dvl7dH//ipp4pKSkq0p46eHIaAViCLY9QTHBz8zc8x7L969apYcPLkSenXr5+OxHAT9e7dWx9Oq1Ce6enp2nvHFBRVHaY9z58/rydMwv379/U+tJjc0BNGQsPo2yqMvDFLEBUVpUntxYsXYllhYaH+2aJFi1r73/T5iKIizAUXFBTIjBkzxBLMDaIxwzxw165dpXXr1tqDu379unTs2FEsePr0qfbOFyxYoENTDPcxDREYGCiJiYlixYMHD7Qsv3z5Ik2aNJHjx49Lt27dfB2WX1m8eLGeMIm59Pr162vSxegMjZwl6AjcuXPH3FpeRQMGDNAeepcuXeTNmzeyatUqGTJkiE4zWzzhr6ysTNd5MM0cExNTNxPF7t27tVeE7G4N1iYw0sF6BB7OPn36yJQpU3T4Z+UGwohi3bp1+hkjCtzsO3bsMJUo8EDeu3dPe0VHjhzR2DAHy2RRdYcPH5YDBw7oefWYskN5ovHAc2OlrvPz8yUlJUWPA608Erek4igMayhIHNiwgjKeNWuWWByhPXz4sNZnMswkiufPn+uC0rFjx8QirJugQcNUCXpzYWFhkpCQoMNVCxBP5cYWox+srViCEY5nFNa3b1/tbW7ZskUXO6lqFi1apKOKyZMn62fsbsPzg90xVhIFOlDv3r3TDpUHRj7YUbZ161YpLi7WDpc1zZo1k86dO0teXp5YM2fOHDl16pSWITaG1CYzaxRYhMXcumeBzirsMkCjjF1GWVlZuqhoAYailbfMYQ4bvSPLMBJCo0FV9+nTJ6lX79tHF40uytKKuLg4nWbEaMdzYcSL6TH83WKS8Cy+P3nyRJ9xKxzH0SSBadoLFy74ZJOKiREFbnAkCvSGrG1N80BSQIVh6gS9DfTqMEeMbYkWzJ8/Xxc5MfU0adIk3QGTmpqqlxVLlizRoX5ERIR+LwVTJ5cuXdKytdZYVOxRPnv2TBs3LB4idgu7x7AmgVgw9YTtxZs2bdKpUSswv195Dh2dLHzHpzbn1v+ThQsXanmiQ/X69WtZsWKFJjFMK1uabvrjjz901yXK9e3bt/pzbIfGhppa4VRTYWEhvt3j5OfnO97KysrS35Wbm+vUpJKSEicjI0P/9NahQ4ecqKgoJzAw0AkNDXWSk5OdgoICUzFmZmY6MTExTlBQkBMdHe2kpqY6NcFT1/jTG0lJSU5kZKSWYUhIiBMXF+dkZ2fXSIw1WY4XL17U/97KV2JiookYi4qKnJSUFCciIsIJDg7W+3LZsmVOcXGxY6WufyQ2NlbjtlTXCQkJTlhYmN6T4eHh+jkvL89UjPKDexFXWlqaV78XbXdV69pE9x1fHrH8rVJALx2XZWPHjtXLKmxW8AfDhg0zfT+iV7l582a9/AlGj9ZgZ5Z1joF70cwaBRER2cREQURErpgoiIjIFRMFERG5YqIgIiJXTBREROSKiYKIiFwxURARkSsmCiIicsVEQURErpgoiIjIFRMFERG5YqIgIiJXTBREROSqyq8ZxylkFU8iw5nHgAPJrSotLZX//d//lVevXpk9EMkfYsQhQ4AYcQysRf5Qjv4QI+u67sT45v+33VV6jXlVD7lYsWLFvz1AgxcvXrx4iV9eT548+Y/tfwD+T3VGFAUFBXp84IsXL/RIPqv69+8vOTk5Ypn1GNGzbNeuneTn50vTpk3FKuvl6A8xsq7rToyFhYV6nO779++lWbNmrv+2ymOioKAgvSpDkrB8Q+H8W8vx+UuMgBgtx+kP5egPMQLrum7ECPXq/eel6n/8YjYOJrfOH2L0B/5Qjv4Qoz/wh3JM9oMYq6rKU08/GqJiNIHhiz9kTao+1nXdwbquO4p+oq6rPaLANNSKFSt+OB1F/yys67qDdV13BP1EXVd7REFERHXDP36NgoiIvMNEQURErpgoiIjIFRMFERG5YqIgIiJXTBREROSKiYKIiFwxURARkbj5f2uRbqx4za7vAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Visualize solution\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "min_val, max_val, diff = 0, n, 1\n",
    "\n",
    "N_points = (max_val - min_val) / diff\n",
    "ind_array = np.arange(min_val, max_val, diff)\n",
    "\n",
    "# This is used to visualize the Sudoku solution with the 16x16 grid below.\n",
    "encode = {1: '1', 2: '2', 3: '3', 4: '4', 5: '5',\n",
    "          6: '6', 7: '7', 8: '8', 9: '9', 10: 'A',\n",
    "          11: 'B', 12: 'C', 13: 'D', 14: 'E', 15: 'F', 16: 'G'}\n",
    "\n",
    "xsol = prob.getSolution(assign)\n",
    "\n",
    "for i1 in Q:\n",
    "    for i2 in Q:\n",
    "        for j1 in Q:\n",
    "            for j2 in Q:\n",
    "                c = encode[math.floor(1 + sum(xsol[i1*q + i2,j1*q + j2,k]*k for k in N) + 0.5)]\n",
    "                ax.text (j1*q + j2, n - 1 - (i1*q + i2), c, va='center', ha='center')\n",
    "\n",
    "# Set up the plot dimensions\n",
    "ax.set_aspect('equal', 'box')\n",
    "ax.set_xlim(min_val-diff/2, max_val-diff/2)\n",
    "ax.set_ylim(min_val-diff/2, max_val-diff/2)\n",
    "\n",
    "# Hide the axis labels\n",
    "ax.set_xticklabels([])\n",
    "ax.set_yticklabels([])\n",
    "\n",
    "# Draw the sudoku grid\n",
    "ax.yaxis.set_minor_locator(ticker.FixedLocator(np.arange(-0.5, q*q, 1)))\n",
    "ax.xaxis.set_minor_locator(ticker.FixedLocator(np.arange(-0.5, q*q, 1)))\n",
    "ax.xaxis.set_major_locator(ticker.FixedLocator(np.arange(-0.5, q*q, q)))\n",
    "ax.yaxis.set_major_locator(ticker.FixedLocator(np.arange(-0.5, q*q, q)))\n",
    "ax.grid(which='minor')\n",
    "ax.grid(which='major', color='black')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "703fbd28",
   "metadata": {},
   "source": [
    "## A variant with a $16\\times 16$ grid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d34f0be",
   "metadata": {},
   "source": [
    "Below is a sudoku on a $16\\times 16$ grid. Just run this cell and restart from the first cell under the \"*Model implementation and results*\" section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b7d47eeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid4x4 = \\\n",
    " [[ 0, 0,12, 0, 0, 2, 0, 0, 0, 7, 3, 0,13,15, 0, 0],\n",
    "  [15, 0, 0, 0, 0, 3, 0, 0, 9, 0, 0, 0,12, 0, 0,10],\n",
    "  [ 0, 0, 0, 0, 9, 0, 6, 0, 0, 0,12, 0, 0, 0, 2, 5],\n",
    "  [ 6,11, 1, 0, 0,10, 5, 0, 0, 2, 0,15, 0, 0, 0, 0],\n",
    "  [ 4, 6, 3, 0, 0, 0,13,14, 0, 0, 0, 0, 0, 7, 0, 0],\n",
    "  [ 0,15,11, 0, 7, 0, 9, 0, 0, 0, 0, 0, 0, 0, 1, 0],\n",
    "  [ 0, 1, 0,10,15, 0, 0, 0,11, 3,14, 0, 6, 0, 0, 0],\n",
    "  [13, 0, 8, 7, 0, 5, 0, 0, 0, 1, 9,12, 0, 0, 0, 0],\n",
    "  [ 0, 0, 0, 6, 3, 7,15, 4, 0, 0, 0, 0, 0,14, 0, 0],\n",
    "  [ 0, 8, 0, 0, 0, 0, 0, 0, 0,11, 7, 0, 4, 0, 0, 0],\n",
    "  [ 0, 0, 0, 0, 0, 0, 0, 0,13, 0, 0, 6, 9, 0, 3, 0],\n",
    "  [ 0, 0, 0, 0, 2, 8,14, 0, 3, 0, 0,10, 0, 0,13, 7],\n",
    "  [ 0, 0, 0, 8, 0, 0, 0, 7,10, 0, 0, 0, 0, 0, 5, 1],\n",
    "  [ 0, 4,10, 1, 6, 0, 0, 0, 0,12, 0,14, 7, 3, 9,15],\n",
    "  [ 3, 0,15, 0, 0, 0, 0, 8, 0, 0, 1, 0,14,12, 0, 0],\n",
    "  [ 2, 0, 0, 9,12, 0, 0, 1, 0, 0, 0, 0, 0, 6, 8, 0]]\n",
    "\n",
    "grid = grid4x4\n",
    "q = 4"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}