{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "83778a8a",
   "metadata": {},
   "source": [
    "# Advanced Usage Example"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e52dafd",
   "metadata": {},
   "source": [
    "In this example we cover advanced usage of `BoolXAI.RuleClassifier`. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9275b4c7",
   "metadata": {},
   "source": [
    "## Input data\n",
    "\n",
    "We'll start with the same binarized data we used in the Basic Usage Example. In order to speed up the execution, we'll only use a subset of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "92481d9c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-08T22:51:41.201232Z",
     "iopub.status.busy": "2023-05-08T22:51:41.200861Z",
     "iopub.status.idle": "2023-05-08T22:51:42.049953Z",
     "shell.execute_reply": "2023-05-08T22:51:42.048537Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn import set_config\n",
    "\n",
    "set_config(transform_output=\"pandas\")\n",
    "\n",
    "from sklearn import datasets\n",
    "\n",
    "X, y = datasets.load_breast_cancer(return_X_y=True, as_frame=True)\n",
    "\n",
    "# Use a subset of the data to speed up execution.\n",
    "# For higher quality results, comment these lines out.\n",
    "X = X.iloc[:100, :100]\n",
    "y = y.iloc[:100]\n",
    "\n",
    "# Binarize the data\n",
    "from util import BoolXAIKBinsDiscretizer\n",
    "\n",
    "binarizer = BoolXAIKBinsDiscretizer(\n",
    "    n_bins=10, strategy=\"quantile\", encode=\"onehot-dense\"\n",
    ")\n",
    "X_binarized = binarizer.fit_transform(X)\n",
    "\n",
    "X_binarized.head();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16fe1ddc",
   "metadata": {},
   "source": [
    "## Number of starts\n",
    "\n",
    "Each training run consists of `num_starts` starts. The best rule and score for each start can be accessed in the `rules_` and `scores_` attributes (respectively). Next we'll train a classifier and inspect these attributes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "31f9b76c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-08T22:51:42.053567Z",
     "iopub.status.busy": "2023-05-08T22:51:42.053238Z",
     "iopub.status.idle": "2023-05-08T22:51:59.083081Z",
     "shell.execute_reply": "2023-05-08T22:51:59.082198Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best score and rule:\n",
      "  score=0.91 rule=Or(230, 50, 271, 222)\n",
      "\n",
      "Best score and rule for each start:\n",
      "  score=0.81 rule=Or(251, 120, 261, 160)\n",
      "  score=0.73 rule=AtMost1(~171, And(257, ~246), ~130)\n",
      "  score=0.74 rule=AtMost1(~50, ~30, 109)\n",
      "  score=0.86 rule=Or(20, 92, 72, 71, 231)\n",
      "  score=0.71 rule=Choose1(~225, And(9, ~37), 195)\n",
      "  score=0.78 rule=Choose1(~205, 28, 206, 256, 64)\n",
      "  score=0.76 rule=AtMost1(~261, ~202, 126)\n",
      "  score=0.71 rule=AtMost1(~101, ~270, 296)\n",
      "  score=0.74 rule=AtLeast2(220, 72, Choose1(138, ~129))\n",
      "  score=0.91 rule=Or(230, 50, 271, 222)\n"
     ]
    }
   ],
   "source": [
    "from boolxai import BoolXAI\n",
    "\n",
    "rule_classifier = BoolXAI.RuleClassifier(random_state=43)\n",
    "rule_classifier.fit(X_binarized, y)\n",
    "\n",
    "print(\"Best score and rule:\")\n",
    "print(f\"  score={rule_classifier.best_score_:.2f} rule={rule_classifier.best_rule_}\\n\")\n",
    "\n",
    "print(\"Best score and rule for each start:\")\n",
    "for score, rule in zip(rule_classifier.scores_, rule_classifier.rules_):\n",
    "    print(f\"  {score=:.2f} {rule=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1477c163",
   "metadata": {},
   "source": [
    "## Bagging\n",
    "\n",
    "Bagging or \"boostrap aggregation\" is commonly used to select only a partial sample of the dataset to avoid overfitting. It is often applied to single estimators that are then combined to form an ensemble estimator. In our case, our model is already highly regularized, so it's unlikely to overfit the data (unless it's very small/simple). For this reason, combined with the binary inputs, it's unnecessary to feed it huge amounts of data. So, in this case we use bagging not to mitigate overfitting, but rather to boost performance, since evaluation of rules scales linearly with the number of samples.\n",
    "\n",
    "In the rule classifier, one can control the level of bagging using the parameter `max_samples`. It controls the maximum number of samples that will be used in each start, with a default of 2,000. For large datasets, this can make a big difference. Note that the `best_rule_` attribute is populated with the best rule as evaluated over the whole (train) dataset, not only `max_samples`.\n",
    "\n",
    "Let's look at an example: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "18536ca5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-08T22:51:59.087101Z",
     "iopub.status.busy": "2023-05-08T22:51:59.086749Z",
     "iopub.status.idle": "2023-05-08T22:52:15.168716Z",
     "shell.execute_reply": "2023-05-08T22:52:15.167825Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Without bagging:\n",
      "score=0.91 rule=Or(230, 50, 271, 222)\n",
      "\n",
      "CPU times: user 21.3 ms, sys: 5.04 ms, total: 26.3 ms\n",
      "Wall time: 5.84 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "rule_classifier = BoolXAI.RuleClassifier(random_state=43)\n",
    "rule_classifier.fit(X_binarized, y)\n",
    "print(\"Without bagging:\")\n",
    "print(f\"score={rule_classifier.best_score_:.2f} rule={rule_classifier.best_rule_}\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e3822987",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-08T22:52:15.172888Z",
     "iopub.status.busy": "2023-05-08T22:52:15.172570Z",
     "iopub.status.idle": "2023-05-08T22:52:18.356418Z",
     "shell.execute_reply": "2023-05-08T22:52:18.355763Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With bagging:\n",
      "score=0.76 rule=Choose1(244, 226, ~168, 217, 219)\n",
      "\n",
      "CPU times: user 61 ms, sys: 6.37 ms, total: 67.4 ms\n",
      "Wall time: 2.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "rule_classifier = BoolXAI.RuleClassifier(\n",
    "    random_state=43, max_samples=len(X_binarized) // 5\n",
    ")\n",
    "rule_classifier.fit(X_binarized, y)\n",
    "print(\"With bagging:\")\n",
    "print(f\"score={rule_classifier.best_score_:.2f} rule={rule_classifier.best_rule_}\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7de06479",
   "metadata": {},
   "source": [
    "The wall clock time was much shorter, but the result was worse. That's expected for such a small dataset - using 20 samples for each start is not representative. However, for larger datasets we have observed that one can get a significant reduction in run time without a significant change in score. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abacb6ea",
   "metadata": {},
   "source": [
    "## Parallelization\n",
    "\n",
    "Since the runtime for training the classifier can be significant, it's useful to parallelize the computation. The number of parallel jobs (starts, in this case) is controlled by the argument `num_jobs`. We can probe the dependence of the runtime on `num_jobs` with a small experiment (results are dependent on the machine you are using):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ecc62ddc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-08T22:52:18.359819Z",
     "iopub.status.busy": "2023-05-08T22:52:18.359537Z",
     "iopub.status.idle": "2023-05-08T22:53:46.530637Z",
     "shell.execute_reply": "2023-05-08T22:53:46.529861Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "num_jobs=24, elapsed_time=7.77 sec, 0.834, AtLeast1(160, 21, 30, 261)\n",
      "num_jobs=12, elapsed_time=6.02 sec, 0.834, AtLeast1(160, 21, 30, 261)\n",
      "num_jobs=8, elapsed_time=4.22 sec, 0.834, AtLeast1(160, 21, 30, 261)\n",
      "num_jobs=4, elapsed_time=6.91 sec, 0.834, AtLeast1(160, 21, 30, 261)\n",
      "num_jobs=2, elapsed_time=11.32 sec, 0.834, AtLeast1(160, 21, 30, 261)\n",
      "num_jobs=1, elapsed_time=18.35 sec, 0.834, AtLeast1(160, 21, 30, 261)\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "seed = 42\n",
    "rule_classifier = BoolXAI.RuleClassifier(\n",
    "    num_starts=24, num_iterations=100, random_state=seed\n",
    ")\n",
    "\n",
    "all_num_jobs = [24, 12, 8, 4, 2, 1]\n",
    "times = []\n",
    "for num_jobs in all_num_jobs:\n",
    "    rule_classifier.num_jobs = num_jobs\n",
    "    start_time = time.time()\n",
    "    rule_classifier.fit(X_binarized, y)\n",
    "    elapsed_time = time.time() - start_time\n",
    "    times.append(elapsed_time)\n",
    "\n",
    "    print(\n",
    "        f\"{num_jobs=}, {elapsed_time=:.2f} sec, {rule_classifier.best_score_:.3}, {rule_classifier.best_rule_}\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0d8481b",
   "metadata": {},
   "source": [
    "It is clear that we can get a large speedup by using multiple cores. Note that the rule and score are the same, as expected, since we used the same seed. We see that there is a sweet spot located at the number of cores of the machine we are using (note: using the full data is advisable to get more representative results):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ca9d1aa5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-08T22:53:46.534046Z",
     "iopub.status.busy": "2023-05-08T22:53:46.533417Z",
     "iopub.status.idle": "2023-05-08T22:53:46.719207Z",
     "shell.execute_reply": "2023-05-08T22:53:46.718542Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAADZCAYAAADIdjAgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAmoUlEQVR4nO3de1xUZf4H8M+AMCjIHcGJS+IFFZG8paaSrSZQYV5Kc820rH6apqy3Ul+FohtKt02z/OW2aq2vbPuVlbuFrRckXaS8Ra7JCpKjAZIkVxWUeX5/PDsDIxdnYODM5fN+vc7LmTNnznyPox8Oz3nO86iEEAJERORQnJQugIiI2h/Dn4jIATH8iYgcEMOfiMgBMfyJiBwQw5+IyAEx/ImIHBDDn4jIAXVQuoC2ptPpUFBQgM6dO0OlUildDhFRqwkhUFFRAY1GAyenlp3D2334FxQUICQkROkyiIgs7sKFCwgODm7Re+0+/Dt37gxA/iV5enoqXA0RUeuVl5cjJCTEkG8tYffhr2/q8fT0vG34a7XA5ctNv+7vD4SGWrI6IqKWa01Ttt2Hv6m0WiAiArh+velt3NyAnBz+ACAi28fePv91+XLzwQ/I15v7zYCIyFYw/ImIHBDDn4jIATH8iYgcEMOfiMgBMfyJiBwQw5+IyAEx/P/L31/242+Om5vcjojI1vEmr/8KDZU3cOn78VdWAvfeKx/v3w94efEOXyKyHwz/ekJDjcM9LAw4fx5wcgIGDlSuLiIiS2OzTzOiouSfp04pWwcRkaUx/JuhD/8ff1S2DiIiS2P4N4PhT0T2iuHfjPrNPkIoWwsRkSUx/JsREQG4uACdOgG//qp0NURElsPePs1wcQEuXQJ8fJSuhIjIsnjmfxsMfiKyRwx/IiIHxPC/jZwcID4eGDdO6UqIiCyHbf630bEjkJYGdOgA1NQArq5KV0RE1Ho887+NkBA5rs/Nm/K3ACIie8Dwvw2VCujXTz7Ozla2FiIiS2H4m4B3+hKRvWH4m4DhT0T2huFvAoY/Edkbhr8J+vWTF33DwmSPHyIiW8eunibw8QGuXJEXf4mI7AHP/E3E4Ccie8LwN1NtrdIVEBG1HsPfRBkZQM+ewJgxSldCRNR6bPM3kbc3kJsrx/UXgs1ARGTbeOZvot695fg+ZWXAxYtKV0NE1DqKhn9GRgYSEhKg0WigUqnw+eefG70+a9YsqFQqoyUuLk6RWl1d5cxeAPv7E5HtUzT8q6qqEB0djU2bNjW5TVxcHAoLCw3LRx991I4VGuPNXkRkLxRt84+Pj0d8fHyz26jVagQFBbVTRc2LigJ27mT4E5Hts/o2//T0dHTp0gURERGYO3cuSkpKmt2+uroa5eXlRoul8MyfiOyFVYd/XFwcPvjgA+zbtw/r16/HwYMHER8fj9pmOtunpKTAy8vLsISEhFisnv79gYEDgWHDLLZLIiJFqIQQQukiAEClUmHXrl2YMGFCk9ucO3cO3bt3x969ezGmiQ731dXVqK6uNjwvLy9HSEgIysrK4OnpaemyiYjaXXl5Oby8vFqVaya1+W/YsMHsHT/55JPo3Lmz2e9rTnh4OPz9/ZGbm9tk+KvVaqjVaot+LhGRvTEp/BMTExEcHAxnZ2eTdnrhwgU89NBDFg//ixcvoqSkBF27drXofs114wZQWgoEBChaBhFRi5nc2+fo0aPo0qWLSduaGvqVlZXIzc01PM/Pz8fJkyfh6+sLX19frF69GpMnT0ZQUBDy8vKwbNky9OjRA7GxsaaWbXE7dgBPPgnExgK7dytWBhFRq5gU/klJSfDw8DB5pytWrICvr+9ttzt69Cjuu+8+w/NFixYBAGbOnIl3330X2dnZ2L59O0pLS6HRaDBu3DisWbNG0WadO+6QZ/7s8UNEtsxqLvi2FUtcGKmvpATw95ePy8oAXkMmovZmiVwzu6tnfn4+zp4922D92bNn8fPPP7eoCFvi5wdoNPLxqVPK1kJE1FJmh/+sWbPwr3/9q8H6rKwszJo1yxI1WT3e7EVEts7s8D9x4gRGjBjRYP2wYcNw8uRJS9Rk9Rj+RGTrzA5/lUqFioqKBuvLysqavfPWnjD8icjWmR3+MTExSElJMQr62tpapKSkYOTIkRYtzloNGQJMmQJMmqR0JURELWN2b5/Tp08jJiYG3t7eGDVqFADg22+/RXl5Ofbv349+/fq1SaEtZenePkRESlOkt0/fvn2RnZ2NKVOmoLi4GBUVFXjiiSdw5swZqwt+IiJqHPv5t5BOB+Tny8fdu1tst0REt6XImT8gm3kef/xx3HPPPfjll18AAB9++CEOHTrUoiJs0bp1QI8eQHKy0pUQEZnP7PD/9NNPERsbi44dO+L48eOG4ZPLysrwyiuvWLxAa9W3r/yTPX6IyBaZHf5r167F5s2bsWXLFri4uBjWjxgxAsePH7docdZM393z9Gng5k1layEiMpfZ4Z+Tk4OYmJgG6728vFBaWmqJmmxCt26AuztQXQ3UG5iUiMgmmB3+QUFBRsMw6x06dAjh4eEWKcoWODkBkZHyMZt+iMjWmB3+zzzzDBYuXIisrCyoVCoUFBRgx44dWLJkCebOndsWNVotfdNPdraydRARmcvkyVz0XnzxReh0OowZMwZXr15FTEwM1Go1lixZgueff74tarRaHOaBiGxVi/v519TUIDc3F5WVlejbt69Zk720p7a8wzc7G/jb34B77gEeeMCiuyYialK7TeDeGFdXV/Tt2xfl5eXYu3cvIiIi0KdPn5buzib17y8XIiJbY3ab/5QpU/D2228DAK5du4YhQ4ZgypQp6N+/Pz799FOLF0hERJZndvhnZGQYBnTbtWsXdDodSktLsWHDBqxdu9biBVq74mJgzx5e9CUi22J2+JeVlRkmZ09LS8PkyZPRqVMnPPjgg41O72jv1q8H4uKA999XuhIiItOZHf4hISHIzMxEVVUV0tLSMG7cOADAlStX4ObmZvECrR17/BCRLTL7gm9iYiKmT58ODw8PhIWFYfTo0QBkc1CUPgkdSP3wFwJQqZSth4jIFC3q6nns2DFotVrcf//9hi6e//jHP+Dt7d3o/L5KauvJXK5dAzw85BDPhYVAUJDFP4KIyIglco3j+VtARATwn/8A33wD3H9/m3wEEZFBu43nv2jRIlRVVZm80+XLl+O3335rUUG2RKsFjh8HgoPl86+/ls/1i1arbH1ERE0x6czf2dkZRUVFCAgIMGmnnp6eOHnypFUM9NZWZ/5arTzjv3696W3c3ICcHCA01GIfS0TUfnf4CiHQq1cvqEy8mmnObwm26vLl5oMfkK9fvszwJyLrY1L4b9261ewdBwYGmv0eIiJqHyaF/8yZM9u6DiIiakctmsCdiIhsG8OfiMgBMfyJiBwQw5+IyAG1OPxzc3OxZ88eXLt2DYDsDupI/P1lP/7muLnJ7YiIrI3ZA7uVlJRg6tSp2L9/P1QqFc6ePYvw8HDMnj0bPj4+eP3119uiTqsTGipv4Lp82Xh9UREwcSJQUwO89Rb7+BORdTL7zP8Pf/gDOnToAK1Wi06dOhnWT506FWlpaRYtztqFhgIDBxovDzwALFsmX3/rLTnSJxGRtTH7zP+bb77Bnj17EKwf0Oa/evbsifPnz1usMFu2bBlQUAAsX84hnonIOpkd/lVVVUZn/Hq//fYb1Gq1RYqydZ07c2YvIrJuZjf7jBo1Ch988IHhuUqlgk6nQ2pqKu677z6z9pWRkYGEhARoNBqoVCp8/vnnRq8LIfDyyy+ja9eu6NixI8aOHWuTU0UWFipdARGRMbPDPzU1Fe+99x7i4+NRU1ODZcuWoV+/fsjIyMD69evN2ldVVRWio6OxadOmJj9rw4YN2Lx5M7KysuDu7o7Y2Fhcv92IalaipgaYNQvo1g2wwZ9ZRGTHWjSZS1lZGd5++2388MMPqKysxMCBAzFv3jx07dq15YWoVNi1axcmTJgAQJ71azQaLF68GEuWLDF8bmBgILZt24bHHnvMpP22x2QuzXngATnO/8SJwGeftfvHE5EdarchnW/l5eWFlStXtugDTZWfn4+ioiKMHTvW6HOHDh2KzMzMJsO/uroa1dXVhufl5eVtWuftvPaanOFr1y7g4EHg3nsVLYeICEALw//69evIzs5GcXExdDqd0Wvjx4+3SGFFRUUAGg4NHRgYaHitMSkpKVi9erVFarCEvn2BZ54BNm8GFi0Cvv8ecOJ91USkMLPDPy0tDU888QQu33p3E2TTTW1trUUKa6nly5dj0aJFhufl5eUICQlRsCJg9Wpgxw45teOOHcCMGYqWQ0Rk/gXf559/Ho8++igKCwuh0+mMFksGf1BQEADg0qVLRusvXbpkeK0xarUanp6eRovSunQBVqyQj5cvB65eVbYeIiKzw//SpUtYtGhRm8/U1a1bNwQFBWHfvn2GdeXl5cjKysLw4cPb9LPbQmIiEBYmewCdPq10NUTk6Mxu9nnkkUeQnp6O7t27t/rDKysrkZuba3ien5+PkydPwtfXF6GhoUhMTMTatWvRs2dPdOvWDS+99BI0Go2hR5AtcXOTF33DwwEvL6WrISJHZ3ZXz6tXr+LRRx9FQEAAoqKi4OLiYvT6ggULTN5Xenp6ozeGzZw5E9u2bYMQAklJSXjvvfdQWlqKkSNH4p133kGvXr1M/gylu3oSEVmaJXLN7PB///33MWfOHLi5ucHPzw+qeoPXqFQqnDt3rkWFtBVrDH8hgP/7P6B3byAqSulqiMjWKNLPf+XKlVi9ejVefPFFOLHPYoskJQFr1gBjx8p7ADj4GxG1N7PTu6amBlOnTmXwt8KTTwKursDevfLuXyKi9mZ2gs+cORMff/xxW9TiMLp1AxYulI+XLAFu3lS2HiJyPGY3+9TW1iI1NRV79uxB//79G1zwfeONNyxWnD1bsQLYuhX46SdgyxZg7lylKyIiR2L2Bd/mhm1WqVTYv39/q4uyJGu84Ku3aRMwfz4QECBH/WQXUCIyhSK9fWyNNYf/jRtA//7AmTPACy8A69YpXRER2QJL5Bqv2irIxQV49VX5A+D++5WuhogciUlt/pMmTcK2bdvg6emJSZMmNbvtZxy03iwPPgjExwPOzkpXQkSOxKTw9/LyMtzM5cWGaYtSqYyDXwj2+yeitmdym39ycjKWLFnS6OTt1sya2/zru3YN+NOfgPR0IC2NPwCIqGntesHX2dkZhYWF6NKlS4s+SCm2Ev4FBUCvXkBVFfDxx8CUKUpXRETWql2Hd7DzTkGK02iAZcvk0A8vvACMHy9HAiUi+6LVAo3MhWXg7w+EhrZ9HWbd5KViW0SbWrwYeO894OefgY0bgaVLla6IiCxJqwUiIoDr15vexs0NyMlp+x8AZnX17NWrF3x9fZtdqOXc3YE//lE+XrsW+PVXZeshIsu6fLn54Afk6839ZmApZp35r169mr192tiMGcCGDXK+31Wr5F3ARGQfrGkcL5Mv+Do5OaGoqIgXfNtBejpw331Ax47y10R/f6UrIiJTVFQApaVASEjdumeflXfxa7XAhQuATnf7/Rw7Bgwc2PTr7XrBl+397Wf0aGD9emDiRAY/kTXatQvIzZWBfv68/FOrBa5cAQYMkL+56337rQx/a8PePlZq2TKlKyByLFVVdSFefzl/HvD0BL78sm7bFSuaDvTycuPnq1fLmzfDwoCyMiAuru2OwRwmh7/OlN9VqE38+CPQty+HgCBqKZ0OKC42DvTaWuOTrMGDmw50Pz/j5w8+CAwaJAM9NFQuYWGyuadzZ+Nt69+zU/83AqWZPZ4/ta+nnwbefx946SVgwoSGr7dXn2Aia3b9ugz1334Dhg2rW//UU0BGhmxrr6kxfo+fn3H4h4bKmy1vDXT94/pDr7z2WtsfU1tj+FsxrRbYvl0+XrNGLrdqrz7BREq5dbyrDz+UZ9D1m2aKi+Vrfn7G3SR/+QXIy5OPnZzkzZT1Q12nk+sB2ayjVrftsfj7y/+zt+vn3x7X+hj+Vuzy5dt3DdP3CWb4ky07f16GdP1mGf3jq1flmbveX/8KfPNNw324uwOBgfIM39VVrlu9Gli5Uoa9RiOHUW9KWwc/IP+f5uTY4B2+RETmEEI2xdwa6CUlchpTvWefbTzQ9aqqZLgDwOTJcg6MW5tlfHwaDohYvwnIWujrVRrDn4harKZGNq1otbK9fNq0uteefhrYuVMGd2Pefrsu0CMi5A+HWwNd/7hjx7r3Pfts2x2PI2H4E1GjhJBdE72969a9/z7wz3/WncEXFMjt9MaPrwt0laou+AMDGwZ6/fdt2NDmh0O3YPjbgVv7FROZ49Qp4IcfGm9vr6gAKivrAv3IETnkeH1qdV2wV1TUbbtypexNExxsfOZO1oHhbwfGjwdefhmYP5/DQFOd8vLGA12rlRMG6UP6rbeAP/+56f1cvCibZQDgkUeAPn2Mz+ADAup6zNR3550WPySyIIa/HaiokMM/b9wou4NOn84bwuxdbS1QWGgc6M8/X3eG/fzzsk29KRcuAL17y8cDBsghRRprbw8JAepP3hcbKxeyfQx/K2Zqn+C1a+UUkFotMHMm8Prrcmyg2FhOB2kp7T0Bh36ogR496ron/uUv8r4PrVaejd/aDXj8+LpA14+/6ONjHOj6UA8Kqnvfc8/JhRwLw9+KmdMn+Lnn5EWzlBQgOxuIj5dDQicltVu5dqstJ+A4ehQ4eNC4Web8edk9EgB++qku0H/5Rd6tqtehg2xP14d6h3r/mxcsABITGw41QKTH8LdypvYJ7thRTv/49NPAK68AmzcDjz1W9/qtd0mS6cyZgCM0tG6ogVsDXf/4q6/kfM0A8PXX8npNY7y8ZH94vQkTgJ496/5NdO3adPMep92g22H42xk/P9ns89JLxl305s2TvTJWruQw0W3t1VebDnRATtOpD/8hQ4Df/964WUa/3BrgUVFyIbIEhr+dqh/8+fnA//6vHMfkL38BXnwRWLjQ+EIeWU5oKODh0Xhbe2iovDtVLy7Oeob4Jcdi8kxetsoWZ/JqC998I5uFTp6Uz++4A0hOlheI2TOoecePy+F7b0c/+5JOJ5vY2MxGbcUSuWbWBO5ku8aNk+H04YfyDPSXX4DZs4HoaODf/1a6Ouvy44/y5qSPPmrZ+52cGPxk/Rj+DsTJCXj8cTlhxWuvyW6AhYXytwBHV1gor5XcdZdslnn11eb7yRPZOrb5OyA3N2DxYjnRxalTddcHhJDdQ2fMkP3LHcHOnXJ0yb176ybWdnEBHnpI/j0Q2Sue+TswHx9g1Ki6519+Ka8D9Okj7xDVT5BhT26djfSDD+T1EJ0OuOce4N13gaIi4LPPgIkT5Tb6m+2a014TcBBZCi/4ksHp03KYiK++ks89PGTb96JFdePA2KoffpDXO3buBA4fltc9ANnP/rvvZHNY9+5Nv7+97/Alao5Fck1YsaSkJAHAaImIiDBrH2VlZQKAKCsra6Mq7c/+/UIMHiyEbAgSIihIiM2bhbh5U+nKzHPxohDr1wsRFVV3LIBcR2TLLJFrVt/mHxkZib179xqed+hg9SXbvPvuA7KygE8+AVasAM6dkxc/n35a6cpMc+6cnPBj//66MeNdXYGEBNmOHx+vbH1E1sDqk7RDhw4Iqj8KFbULJydg6lTZ7r15s7wOoL8f4Pp1eb+AtUyRd/OmnFRE3+zSpQuQmSmDf+RIGfiPPiqvcRCRZPUXfM+ePQuNRoPw8HBMnz4dWq222e2rq6tRXl5utFDLubrKQcLuv79u3caNwPDhwKRJcjAzJQgBnDghr0eEhMizej0PDznJd14e8O238rcABj+RMasO/6FDh2Lbtm1IS0vDu+++i/z8fIwaNQoVFRVNviclJQVeXl6GJSQkpB0rdgzFxfI3g127gMhIYM4c2U++PVy8KIerjoqSd9O++absnfPLL/JPvYkTgfDw9qmJyBbZVG+f0tJShIWF4Y033sDs2bMb3aa6uhrV1dWG5+Xl5QgJCWFvHws7fVqOEbR7t3zeqZO8d2Dp0rYbRnjNGjlEtf5frFotx7B//HE5Po6ra9t8LpG1sURvH6tv86/P29sbvXr1Qm5ubpPbqNVqqNXqdqzKMfXtK+8L+PZbGfhZWTKctVpg27bW7//mTTlReJ8+ddMBDhgggz8mRrbjP/KI8QB2RGQ6mwr/yspK5OXlYQZvvbQao0bJi6uffiqHMX7xxbrXrl6V8wxcuGBaH3kh5CBqH34ox9UpLgaWL5fzEwByZrL8fM4NS2QJVh3+S5YsQUJCAsLCwlBQUICkpCQ4Oztj2rRpSpdG9ahU8ix88mTjAc0WLJCDyZ0+DdTUNP1+tVoOMb17t5y5Si8gwPjmMhcXBj+RpVh1+F+8eBHTpk1DSUkJAgICMHLkSBw5cgQBAQFKl0aNqB/8V67I+wRM6WxVXQ2kpsrHajXw8MOyWSc2tm7+WiKyLKsO/507dypdArWQj4/sBjp/vmwSup2BA+U8xI88wikIidqDVXf1JNsWFCTvEDbFli1yfgEGP1H7YPgTETkghj8RkQNi+BMROSCGPxGRA2L4U5viLFhE1smqu3qS7QsNlV0+OQsWkXWx+/DXj1vHoZ2V4+19+zF4+PUQmU6fZ60Zl9Puw18//DOHdiYie1NRUQGvFt4cY1NDOreETqdDQUEBhBAIDQ3FhQsXHGZoZ/1w1jxm++eIx+3Ix6zVaqFSqaDRaODk1LJLt3Z/5u/k5ITg4GDDr0menp4O8w9Fj8fsOBzxuB3xmL28vFp9zOztQ0TkgBj+REQOyGHCX61WIykpyaFm+eIxOw5HPG4ec+vY/QVfIiJqyGHO/ImIqA7Dn4jIATH8iYgcEMOfiMgBOUT4b9q0CXfeeSfc3NwwdOhQfPfdd0qX1KZWrVoFlUpltPTu3VvpsiwqIyMDCQkJ0Gg0UKlU+Pzzz41eF0Lg5ZdfRteuXdGxY0eMHTsWZ8+eVaZYC7ndMc+aNavB9x4XF6dMsRaSkpKCIUOGoHPnzujSpQsmTJiAnJwco22uX7+OefPmwc/PDx4eHpg8eTIuXbqkUMWtZ8oxjx49usF3PWfOHLM+x+7D/+OPP8aiRYuQlJSE48ePIzo6GrGxsSguLla6tDYVGRmJwsJCw3Lo0CGlS7KoqqoqREdHY9OmTY2+npqaig0bNmDz5s3IysqCu7s7YmNjcf369Xau1HJud8wAEBcXZ/S9f/TRR+1YoeUdPHgQ8+bNw5EjR/DPf/4TN27cwLhx41BVVWXY5g9/+AN2796NTz75BAcPHkRBQQEmTZqkYNWtY8oxA8Azzzxj9F2npqaa90HCzt19991i3rx5hue1tbVCo9GIlJQUBatqW0lJSSI6OlrpMtoNALFr1y7Dc51OJ4KCgsSrr75qWFdaWirUarX46KOPFKjQ8m49ZiGEmDlzpnj44YcVqae9FBcXCwDi4MGDQgj5vbq4uIhPPvnEsM1PP/0kAIjMzEylyrSoW49ZCCHuvfdesXDhwlbt167P/GtqanDs2DGMHTvWsM7JyQljx45FZmamgpW1vbNnz0Kj0SA8PBzTp0+HVqtVuqR2k5+fj6KiIqPv3cvLC0OHDrX77z09PR1dunRBREQE5s6di5KSEqVLsqiysjIAgK+vLwDg2LFjuHHjhtF33bt3b4SGhtrNd33rMevt2LED/v7+6NevH5YvX46rV6+atV+7Htjt8uXLqK2tRWBgoNH6wMBAnDlzRqGq2t7QoUOxbds2REREoLCwEKtXr8aoUaNw6tQpdO7cWeny2lxRUREANPq961+zR3FxcZg0aRK6deuGvLw8rFixAvHx8cjMzISzs7PS5bWaTqdDYmIiRowYgX79+gGQ37Wrqyu8b5kwwl6+68aOGQB+//vfIywsDBqNBtnZ2XjhhReQk5ODzz77zOR923X4O6r4+HjD4/79+2Po0KEICwvD3/72N8yePVvByqgtPfbYY4bHUVFR6N+/P7p374709HSMGTNGwcosY968eTh16pTdXb9qTlPH/OyzzxoeR0VFoWvXrhgzZgzy8vLQvXt3k/Zt180+/v7+cHZ2bnDl/9KlSwgKClKoqvbn7e2NXr16ITc3V+lS2oX+u3X07z08PBz+/v528b3Pnz8ff//733HgwAEEBwcb1gcFBaGmpgalpaVG29vDd93UMTdm6NChAGDWd23X4e/q6opBgwZh3759hnU6nQ779u3D8OHDFaysfVVWViIvLw9du3ZVupR20a1bNwQFBRl97+Xl5cjKynKo7/3ixYsoKSmx6e9dCIH58+dj165d2L9/P7p162b0+qBBg+Di4mL0Xefk5ECr1drsd327Y27MyZMnAcC877pVl4ttwM6dO4VarRbbtm0Tp0+fFs8++6zw9vYWRUVFSpfWZhYvXizS09NFfn6+OHz4sBg7dqzw9/cXxcXFSpdmMRUVFeLEiRPixIkTAoB44403xIkTJ8T58+eFEEKsW7dOeHt7iy+++EJkZ2eLhx9+WHTr1k1cu3ZN4cpbrrljrqioEEuWLBGZmZkiPz9f7N27VwwcOFD07NlTXL9+XenSW2zu3LnCy8tLpKeni8LCQsNy9epVwzZz5swRoaGhYv/+/eLo0aNi+PDhYvjw4QpW3Tq3O+bc3FyRnJwsjh49KvLz88UXX3whwsPDRUxMjFmfY/fhL4QQGzduFKGhocLV1VXcfffd4siRI0qX1KamTp0qunbtKlxdXcUdd9whpk6dKnJzc5Uuy6IOHDggADRYZs6cKYSQ3T1feuklERgYKNRqtRgzZozIyclRtuhWau6Yr169KsaNGycCAgKEi4uLCAsLE88884zNn+Q0drwAxNatWw3bXLt2TTz33HPCx8dHdOrUSUycOFEUFhYqV3Qr3e6YtVqtiImJEb6+vkKtVosePXqIpUuXirKyMrM+h0M6ExE5ILtu8yciosYx/ImIHBDDn4jIATH8iYgcEMOfiMgBMfyJiBwQw5+IyAEx/Mmh/fzzz1CpVIbb463BmTNnMGzYMLi5ueGuu+5qdJtZs2ZhwoQJJu8zPT0dKpWqwRg45LgY/qQo/dSD69atM1r/+eefQ6VSKVSVspKSkuDu7o6cnByjMWvqe+utt7Bt27b2LYzsCsOfFOfm5ob169fjypUrSpdiMTU1NS1+b15eHkaOHImwsDD4+fk1uo2Xl1eDMeyJzMHwJ8WNHTsWQUFBSElJaXKbVatWNWgC+dOf/oQ777zT8FzfFPLKK68gMDAQ3t7eSE5Oxs2bN7F06VL4+voiODgYW7dubbD/M2fO4J577oGbmxv69euHgwcPGr1+6tQpxMfHw8PDA4GBgZgxYwYuX75seH306NGYP38+EhMT4e/vj9jY2EaPQ6fTITk5GcHBwVCr1bjrrruQlpZmeF2lUuHYsWNITk6GSqXCqlWrGt3Prc0+1dXVWLBgAbp06QI3NzeMHDkS33//fYP3HT58GP3794ebmxuGDRuGU6dOGV47f/48EhIS4OPjA3d3d0RGRuKrr75q9PPJ9jH8SXHOzs545ZVXsHHjRly8eLFV+9q/fz8KCgqQkZGBN954A0lJSXjooYfg4+ODrKwszJkzB//zP//T4HOWLl2KxYsX48SJExg+fDgSEhIMUyCWlpbid7/7HQYMGICjR48iLS0Nly5dwpQpU4z2sX37dri6uuLw4cPYvHlzo/W99dZbeP311/Haa68hOzsbsbGxGD9+PM6ePQsAKCwsRGRkJBYvXozCwkIsWbLEpONetmwZPv30U2zfvh3Hjx9Hjx49EBsbi99++63Bcb7++uv4/vvvERAQgISEBNy4cQOAnDikuroaGRkZ+PHHH7F+/Xp4eHiY9Plkgyw+JB2RGepPOj5s2DDx1FNPCSGE2LVrl6j/z7OxSenffPNNERYWZrSvsLAwUVtba1gXEREhRo0aZXh+8+ZN4e7ubpjIPT8/XwAQ69atM2xz48YNERwcLNavXy+EEGLNmjVi3LhxRp994cIFAcAwUui9994rBgwYcNvj1Wg04o9//KPRuiFDhojnnnvO8Dw6OlokJSU1u5/6f2+VlZXCxcVF7Nixw/B6TU2N0Gg0IjU1VQhRNyLozp07DduUlJSIjh07io8//lgIIURUVJRYtWrVbY+B7APP/MlqrF+/Htu3b8dPP/3U4n1ERkbCyanun3VgYCCioqIMz52dneHn54fi4mKj99Wf+KNDhw4YPHiwoY4ffvgBBw4cgIeHh2Hp3bs3ANk+rzdo0KBmaysvL0dBQQFGjBhhtH7EiBGtOua8vDzcuHHDaL8uLi64++67G+y3/nH6+voiIiLCsM2CBQuwdu1ajBgxAklJScjOzm5xTWT9GP5kNWJiYhAbG4vly5c3eM3JyQniltHH9c0V9bm4uBg9V6lUja7T6XQm11VZWYmEhAScPHnSaDl79ixiYmIM27m7u5u8T2v09NNP49y5c5gxYwZ+/PFHDB48GBs3blS6LGojDH+yKuvWrcPu3buRmZlptD4gIABFRUVGPwAs2Tf/yJEjhsc3b97EsWPH0KdPHwDAwIED8e9//xt33nknevToYbSYE/ienp7QaDQ4fPiw0frDhw+jb9++La69e/fuhmsNejdu3MD333/fYL/1j/PKlSv4z3/+YzhOAAgJCcGcOXPw2WefYfHixdiyZUuL6yLr1kHpAojqi4qKwvTp07Fhwwaj9aNHj8avv/6K1NRUPPLII0hLS8PXX38NT09Pi3zupk2b0LNnT/Tp0wdvvvkmrly5gqeeegqAvBC6ZcsWTJs2DcuWLYOvry9yc3Oxc+dO/PnPf4azs7PJn7N06VIkJSWhe/fuuOuuu7B161acPHkSO3bsaHHt7u7umDt3rqFHU2hoKFJTU3H16lXMnj3baNvk5GT4+fkhMDAQK1euhL+/v6HXUGJiIuLj49GrVy9cuXIFBw4cMPrBQPaFZ/5kdZKTkxs0y/Tp0wfvvPMONm3ahOjoaHz33Xcm94Qxxbp167Bu3TpER0fj0KFD+PLLL+Hv7w8AhrP12tpajBs3DlFRUUhMTIS3t7fR9QVTLFiwAIsWLcLixYsRFRWFtLQ0fPnll+jZs2er6588eTJmzJiBgQMHIjc3F3v27IGPj0+D7RYuXIhBgwahqKgIu3fvhqurKwCgtrYW8+bNQ58+fRAXF4devXrhnXfeaVVdZL04jSORDZo2bRqcnZ3x17/+VelSyEbxzJ/Ihty8eROnT59GZmYmIiMjlS6HbBjDn8iGnDp1CoMHD0ZkZCTmzJmjdDlkw9jsQ0TkgHjmT0TkgBj+REQOiOFPROSAGP5ERA6I4U9E5IAY/kREDojhT0TkgBj+REQOiOFPROSA/h9NkoAx/lrbogAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4, 2))\n",
    "plt.plot(all_num_jobs, times, \"s--b\")\n",
    "plt.xlabel(\"Number of jobs\")\n",
    "plt.ylabel(\"Time [sec]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdd973cb",
   "metadata": {},
   "source": [
    "## Cross validation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c1dd241",
   "metadata": {},
   "source": [
    "Cross-validation is often used in order to quantify a model's generalization ability, as well as for hyperparameter optimization. Below we provide an example of running cross-validation for a rule classifier. Here we switch off the internal parallelization over starts, instead parallelizing over the cross-validation splits. Note that the binarization is data-dependent, so we make sure to run it separately on each split. This is done by instantiating an `sklearn`  `Pipeline` consisting of a binarizer and a rule classifier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3ea55b4d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-08T22:53:46.721972Z",
     "iopub.status.busy": "2023-05-08T22:53:46.721744Z",
     "iopub.status.idle": "2023-05-08T22:55:25.110908Z",
     "shell.execute_reply": "2023-05-08T22:55:25.110006Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train=0.92+-0.03, Time=36.68+-0.68sec\n",
      "Test=0.80+-0.08\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import cross_validate, StratifiedShuffleSplit\n",
    "from sklearn.pipeline import Pipeline, make_pipeline\n",
    "import numpy as np\n",
    "\n",
    "# Define the pipeline\n",
    "rule_classifier = BoolXAI.RuleClassifier(num_jobs=1, random_state=seed)\n",
    "binarizer = BoolXAIKBinsDiscretizer(\n",
    "    n_bins=10, strategy=\"quantile\", encode=\"onehot-dense\"\n",
    ")\n",
    "pipeline = make_pipeline(binarizer, rule_classifier)\n",
    "\n",
    "# Run cross-validation\n",
    "cv = StratifiedShuffleSplit(n_splits=8, train_size=0.7, random_state=seed)\n",
    "result = cross_validate(\n",
    "    pipeline,\n",
    "    X,\n",
    "    y,\n",
    "    cv=cv,\n",
    "    n_jobs=-1,  # Parallelization here, instead of inside classifier\n",
    "    return_train_score=True,\n",
    "    return_estimator=True,\n",
    "    error_score=\"raise\",\n",
    ")\n",
    "\n",
    "print(\n",
    "    f\"\"\"Train={np.mean(result[\"train_score\"]):.2f}+-{np.std(result[\"train_score\"]):.2f}, Time={np.mean(result[\"fit_time\"]):.2f}+-{np.std(result[\"fit_time\"]):.2f}sec\"\"\"\n",
    ")\n",
    "print(\n",
    "    f\"\"\"Test={np.mean(result[\"test_score\"]):.2f}+-{np.std(result[\"test_score\"]):.2f}\"\"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61e4db14",
   "metadata": {},
   "source": [
    "## Pareto frontier\n",
    "\n",
    "In most cases there isn't a clear reason to choose a specific `max_complexity` value. In general, we might expect that higher values will provide more expressivity and therefore higher metric scores, but this is not guaranteed. For these reasons, it's often a good idea to find the best rule for multiple `max_complexity` values. Viewing the resulting curve (the Pareto frontier) in the metric-complexity 2D space is instructive. We show an example of how to do this using cross-validation. Warning, this cell might take a few minutes to run!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6d35dd31",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-08T22:55:25.115231Z",
     "iopub.status.busy": "2023-05-08T22:55:25.114658Z",
     "iopub.status.idle": "2023-05-08T23:00:26.893350Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max_complexity=3, mean_complexity=3.0, train_score=0.78, test_score=0.71, fit_time=22.41sec\n",
      "max_complexity=6, mean_complexity=5.75, train_score=0.92, test_score=0.80, fit_time=30.02sec\n",
      "max_complexity=10, mean_complexity=8.125, train_score=0.91, test_score=0.80, fit_time=41.83sec\n",
      "max_complexity=15, mean_complexity=9.75, train_score=0.94, test_score=0.79, fit_time=60.93sec\n"
     ]
    }
   ],
   "source": [
    "# Run cross-validation\n",
    "train_scores = []\n",
    "test_scores = []\n",
    "fit_times = []\n",
    "complexities = []\n",
    "\n",
    "for max_complexity in [3, 6, 10, 15]:\n",
    "    # Use the pipeline we instantiated earlier, adjusting max_complexity\n",
    "    pipeline.set_params(ruleclassifier__max_complexity=max_complexity)\n",
    "\n",
    "    result = cross_validate(\n",
    "        pipeline,\n",
    "        X,\n",
    "        y,\n",
    "        cv=cv,\n",
    "        n_jobs=-1,  # Parallelization here, instead of inside classifier\n",
    "        return_train_score=True,\n",
    "        return_estimator=True,\n",
    "        error_score=\"raise\",\n",
    "    )\n",
    "\n",
    "    train_score = np.mean(result[\"train_score\"])\n",
    "    train_scores.append(train_score)\n",
    "\n",
    "    test_score = np.mean(result[\"test_score\"])\n",
    "    test_scores.append(test_score)\n",
    "\n",
    "    fit_time = np.mean(result[\"fit_time\"])\n",
    "    fit_times.append(fit_time)\n",
    "\n",
    "    pipelines = result[\"estimator\"]\n",
    "    mean_complexity = np.mean(\n",
    "        [pipeline[-1].best_rule_.complexity() for pipeline in pipelines]\n",
    "    )\n",
    "    complexities.append(mean_complexity)\n",
    "\n",
    "    print(\n",
    "        f\"{max_complexity=}, {mean_complexity=}, {train_score=:.2f}, {test_score=:.2f}, {fit_time=:.2f}sec\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9faf1d3",
   "metadata": {},
   "source": [
    "We can plot the results for easier comparison. Note that we make sure to use the mean of the complexities of the resulting rules on the x axis, and not the respective `max_complexity` (which is an upper bound on the former):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ae2e9197",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGwCAYAAABFFQqPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABcHUlEQVR4nO3deVhU1R8G8HdA2QUXZBMURNwKl1zINSsKN3JJxR01zco1bNEStfwpaaVo7kuaSy4papmihmmaJi64pbkrqICiAgICMnN/fxxBkQEZGObO8n6eZx5nLjOXL5M5L+ee8z0KSZIkEBEREZkgM7kLICIiIpILgxARERGZLAYhIiIiMlkMQkRERGSyGISIiIjIZDEIERERkcliECIiIiKTVU7uAvSRSqXC7du3UaFCBSgUCrnLISIiomKQJAkPHz6Em5sbzMyKN9bDIKTG7du34eHhIXcZREREVAJxcXFwd3cv1nMZhNSoUKECAPFG2tvby1wNERERFUdqaio8PDzyPseLg0FIjdzLYfb29gxCREREBkaTaS2cLE1EREQmi0GIiIiITBaDEBEREZkszhEqBaVSicePH8tdhsGysLAo9vJGIiKissAgVAKSJCEhIQHJyclyl2LQzMzM4OXlBQsLC7lLISIiE8UgVAK5IcjJyQk2NjZsulgCuU0r4+PjUb16db6HREQkCwYhDSmVyrwQVKVKFbnLMWhVq1bF7du3kZOTg/Lly8tdDhERmSBO0NBQ7pwgGxsbmSsxfLmXxJRKpcyVEBGRqWIQKiFeyik9vodERCQ3XhojIiKiMhUbCyQlFf51R0egenXd1fMsBiEiIiIqM7GxQJ06QGZm4c+xsgIuXJAnDDEI6Zg+p2JNeXp6YuzYsRg7dqzcpRARkZ5KSio6BAHi60lJDEJGT65U/KK5OJMnT8aUKVM0Pu/Ro0dha2tbwqqIiIjkxyCkQ3Kl4vj4+Lz7GzZswKRJk3DhwoW8Y3Z2dnn3JUmCUqlEuXIv/qtRtWpV7RVJREQkA64a0wJJAtLTX3x79Kh453v0qHjnk6Tinc/FxSXv5uDgAIVCkff4v//+Q4UKFbBz5040adIElpaWOHjwIK5cuYIuXbrA2dkZdnZ2aNasGf7444985/X09ER4eHjeY4VCgWXLlqFbt26wsbGBj48Pfv3112K+i0REZIz0fRMGBiEtyMgA7OxefGvdunjna926eOfLyNDezzB+/Hh88803OH/+PBo0aIC0tDR07NgRUVFRiImJQfv27REYGIjY2Ngiz/PVV1+hV69eOH36NDp27Ih+/frh/v372iuUiIgMwvnzwAcfAB07yl1J0RiECADw9ddf46233oK3tzcqV66Mhg0bYvjw4Xj55Zfh4+ODqVOnwtvb+4UjPIMGDUKfPn1Qq1YtTJ8+HWlpaYiOjtbRT0FERHKSJGDXLqBDB6B+fWDxYiArS+6qisY5QlpgYwOkpb34eSdPFm9U6OBBoFGj4n1fbWnatGm+x2lpaZgyZQp+//13xMfHIycnB48ePXrhiFCDBg3y7tva2sLe3h537tzRXqFERKR3MjKANWuAOXOAc+fEMYUC6NIF6NQJGDZM3vqKwiCkBQoFUJzFU9bWxTuftXXxzqdNz6/++uSTT7Bnzx589913qFWrFqytrdGjRw9kZ2cXeZ7n9wxTKBRQqVRar5eIiOR36xawYIEY+bl3TxyzswPeew8YPRqoWVOsmLayevGKaUdH3dT8PAYhUuvvv//GoEGD0K1bNwBihOj69evyFkVERHrh2DEgPBzYsAHIyRHHPD1F+BkyBHBwePrc6tVFWxh97aHHIKRDjo76nYqf5ePjg4iICAQGBkKhUCA0NJQjO0REJiwnB9i2DZg9G/j776fH27YFxo4F3nkHMDdX/9rq1fW3WTCDkA7peyp+1qxZszBkyBC0bNkSjo6O+Pzzz5Gamip3WUREpGMpKcCyZcAPPwA3bohj5csDvXsDY8YATZrIW19pKSSpuN1oTEdqaiocHByQkpICe3v7fF/LzMzEtWvX4OXlBSsrK5kqNA58L4mI9Nfly8DcucCKFU8XBDk6iiXxH34IuLnJW586RX1+F4YjQkRERARALH/ft0/M//ntt6eNe196SVz+6tev+At/DAWDEBERkYnLygLWrRMB6NSpp8c7dgQ+/hh4802xQtoYMQgRERGZqMREYNEisQQ+t+WbjQ0waJBYAVanjqzl6QSDEBERkYk5dUqM/vz8M5DbHs7dHRg1Chg6FKhcWdbydIpBiIiIyASoVMD27SIA/fnn0+N+fuLyV/fuYjWYqWEQIiIiMmJpaWLl19y5YiUYIPr99OghJkC/+qqs5cmOQYiIiMgIXb8OzJsnegClpIhjFSsC778PjBihHz3r9AGDEBERkZGQJODQIdH9ecsWcTkMAGrXFs0Pg4N1v5elvmMQIiIiMnDZ2cCmTSIAHTv29Li/v5j/0749YGYmX336jG+LXJRK0bVq3Trxp1JZZt9KoVAUeZsyZUqpzr1161at1UpERMV37x4QFgZ4eYlmh8eOAZaWYvf3M2eAPXtELyCGoMJxREgOERFijPLmzafH3N2BOXPEtH0ti4+Pz7u/YcMGTJo0CRcuXMg7Zmdnp/XvSUREZefcOfGRsXo18OiROObiIub+DB8OVK0qb32GhBlR1yIixFT9Z0MQANy6JY5HRGj9W7q4uOTdHBwcoFAo8h1bv3496tWrBysrK9StWxcLFizIe212djZGjhwJV1dXWFlZoUaNGggLCwMAeHp6AgC6desGhUKR95hMR2wscOJE4bfYWLkrJDIekgRERorLXC+9BCxZIkJQ48bAqlViQ9SJExmCNMURIW2QJCAj48XPUypFq051+9xKkuhfPmaMuKhrbv7i89nYlLrn+dq1azFp0iTMmzcPjRs3RkxMDIYNGwZbW1sEBwdj7ty5+PXXX7Fx40ZUr14dcXFxiIuLAwAcPXoUTk5OWLFiBdq3bw/z4tRMRiM2VnSdzcws/DlWVsCFC1ydQlQaGRli5GfOHOD8eXFMoQC6dhXL39u0Md7tL3SBQUgbMjIAbVxekiQxUuTgULznp6WVevr/5MmT8f3336P7k0tyXl5eOHfuHBYvXozg4GDExsbCx8cHrVu3hkKhQI0aNfJeW/XJrx0VK1aEi4tLqeogw5OUVHQIAsTXk5IYhIhK4tYtYP58YPFi4P59caxCBTH/Z9QooGZNeeszFgxCJiw9PR1XrlzBe++9h2HDhuUdz8nJgcOTMDZo0CC89dZbqFOnDtq3b4/OnTvj7bfflqtkIiKjd/So6P68cSOQkyOOeXmJCwaDBwP29rKWZ3QYhLTBxkaMzrzIX3+J6fsvsmMH0LZt8b5vKaQ9qXnp0qXw8/PL97Xcy1yvvPIKrl27hp07d+KPP/5Ar1694O/vj02bNpXqe5PpOHYMcHIC3Ny4coWoMDk5ou9PeLjoA5TrtdfE5a/AwOLNmCDNMQhpg0JRvEtUb78tVofduqV+npBCIb7+9ts6+Rvv7OwMNzc3XL16Ff369Sv0efb29ggKCkJQUBB69OiB9u3b4/79+6hcuTLKly8PZRku/SfDN3y4+NPaGvD2BmrVEjcfn6f33d0Zksg0JSeLzs8//PB0cUH58kCfPmIE6JVXZC3PJDAI6ZK5uZjt1qOHCD3PhqHcmW7h4TqN/V999RVGjx4NBwcHtG/fHllZWTh27BgePHiAkJAQzJo1C66urmjcuDHMzMzwyy+/wMXFBRUrVgQgVo5FRUWhVatWsLS0RKVKlXRWOxkGd3cgPl6sbjl7VtyeZ2n5NCQ9G5B8fMTr+ZswGZtLl8TeXytWAOnp4pijI/Dhh+Lm6ipvfaaEQUjXuncX7T/V9REKDy+TPkJFGTp0KGxsbPDtt9/i008/ha2tLXx9fTF27FgAQIUKFTBz5kxcunQJ5ubmaNasGXbs2AGzJ7++f//99wgJCcHSpUtRrVo1XL9+Xaf1k/7btg3w9RVLey9dEps+Xr789P61a0BWluiLcu5cwddbWIhJoc8HpFq1xCRshiQyFJIkdn0PDxe7wOf+Lvzyy+LyV9++YuSUdEshSequ0Zi21NRUODg4ICUlBfbPzUrLzMzEtWvX4OXlBSsrq5J/E6USOHBA/Krs6irWP5rYv+haey9JFkuXis0bX+T48aKH93NyREjKDUjPhqSrV4HHjwt/bfnyYhKpupBUowZQjr/qkR7IzBSbCISHA6dPPz3eqZPY/uKNN7j8XVuK+vwuDP+ZkIu5OdCundxVEJXImTPAuHEvfp6VlRjuL0q5cuKymLc3EBCQ/2tKpZg38XxAunwZuHJF7K908aK4qTuvl5f6OUmeniJEEZWlxERg4UJxu3NHHLOxAQYNEhcFateWtTx6gkGIiDQSGys62z58CDRtKqa9FTag5+hYuh5C5uYizHh5AW+9lf9rSqW4uvx8QLp0SYSkrCxx/9Il9ef19FQfkry8xOU4opI6eVKM/qxbJ8I6AHh4iN4/Q4cCnEqpXxiEiKjY7t0Toza3b4sW/7t3y/ePurm5uPxVowbw5pv5v6ZSicWZ6kLS5cviUsWVK+K2a1f+15qZiXOqm7jt5SUmdhM9T6kU837Cw8U+2rlefVVc/urWjaOQ+opBiIiKJSND9DL57z8xt3/nTv39zdbMTPwG7uEBvP56/q+pVGJqnrqAdPmy+DmvXRO3PXvyv1ahECNc6uYk1axZ+MgYGa+HD8XKr7lzRbAGREjv2VNc/nr1VXnroxdjECohzjEvPb6HhiMnR/Q1OXwYqFhRbPzo4SF3VSVjZgZUqyZuz0/TkyQgIaHwkJSWJiZ237gB/PFH/tcqFOI9yQ1Iz4Ykb2+uBjI216+L3j/LlgGpqeJYpUpiAcGIEYb7/4cpYhDSUPknY5sZGRmw5r9spZL95OI5N2vVb5IEfPQR8OuvYsTjt9/EZTFjpFCIRZyurgWbu0uSmPyqLiBduiRGBmJjxW3v3oLndndXPyfJ27vUWwaSjkgS8Pff4vLXli1idBEQmw+PGQMMHMj/loaIQUhD5ubmqFixIu48WQJgY2MDBdc9akylUuHu3buwsbFBOa5x1mtffSWWypuZicmfrVvLXZE8FArAxUXcnn8PJAm4e7fwkJSSIiZ237yZf/5ILjc39XOSvL21s58zlU52NvDLLyIAHTv29Phbb4n5PwEB7IxuyNhHSI0X9SGQJAkJCQlITk7WfXFGxMzMDF5eXrDgEh29tXgx8MEH4v7ChU/vU/FJkphk/nxAyn384EHRr3dxKTwkcfPNspWUJP4fmD9fzCsDxKjogAHA6NGiESLpl5L0EWIQUqO4b6RSqcTjorq9UZEsLCzyOlST/tm6FXj3XTH8P2mSGBki7bt/v/DVbffuFf1aJyf1E7dr1QIcHHRTvzH691/RFmL1arHCEBCXS0eMEHOAqlaVtz4qHIOQlpTkjSQyJgcPimH/zEzR92TJEna+lcODB+o7bl++LC7FFaVqVfUBqVYt/V3tJyeVSrRSCA8XbSFyvfKKuPzVqxf7SxkCBiEtYRAiU/bvv2IOTHKyWC4fEcGtKvRRcrJYrq0uJCUmFv3aKlXUhyQfH6ByZZ2UrzfS08XIz5w5ojUEIOb7dO0q9v9q3Zq/BBgSBiEtYRAiU3XzJtCixdM///hDbAlAhiU19WlIev6SW0JC0a+tVEl9QKpVSwQoYwkFN28C8+aJ0c7ceVoVKogR0FGjRPNMMjwMQlrCIESm6MEDsffvv/8CdeuKy2NVqshdFWlbWlrhIen27aJf6+BQ8DJb7uOqVQ0jJEVHA7Nni1VgSqU4VrOmWP4+aBAnoBs6BiEtYRAiU/PokVgCfOCAWMp96JDYZoJMS3o6cPWq+onbN28W/Vp7e/UBqVYtwNlZ3pCUkyP6/syeLZqC5nrtNTH/p3Nn0Q2aDB+DkJYwCJEpUSrFdgBbtojf+A8cAHx95a6K9M2jR/nnJD0bkuLiRJuAwtjZFT4nycWlZCEpNlYsby+Mo6O41LVsmegAHRcnjpcvD/TtK0aAGjfW/PuSfivJ5zenQBKZMEkCRo4UIcjCAti2jSGI1LO2Fn1z1PXOycwUI0nqJm7fuCEux508KW7Ps7EpPCS5uqpvVBgbK7o55y5tV8fcXPydfvRIPK5aFfjwQ3FzcSnJO0DGikGIyIRNmwYsWiR+I1+7VlwqINKUlRVQv764PS8rS2xgq66h5PXrYpPb06fF7XnW1qJx5PMB6dGjokMQIEY6Hz0SwX7sWDEKxE1xSR1eGlODl8bIFCxbBgwbJu7/8IMYGSLSpexsEYbUTdy+fv3pZOaSWrgQGD7cMCZxk3bw0hgRFctvv4kPCAD44guGIJKHhQVQu7a4Pe/xY3FZTd22JFevFi8kNW/OEEQvxiBEZGIOHwaCgkQn3UGDgP/9T+6KiAoqX/7p5bDnRUcDfn66r4mMEzd6IjIh//0nlgo/egR07MitM8gwsdM5aRODEJGJuH1b9Aq6f19cMti4UfzWTURkyhiEiExAcjLQvr1Ydly7NvD774CtrdxVERHJT/YgNH/+fHh6esLKygp+fn6Ijo4u9LmPHz/G119/DW9vb1hZWaFhw4aIjIws1TmJjF1mpthA8swZ0T8lMlI0myMyVI6OL14Kb2XFv+dUPLJead2wYQNCQkKwaNEi+Pn5ITw8HAEBAbhw4QKcnJwKPH/ixIlYs2YNli5dirp162LXrl3o1q0bDh06hMZPWoRqek4iY6ZUAgMGAPv3iy67O3dyM0kyfNWrAxcuvLizdPXququJDJesfYT8/PzQrFkzzJs3DwCgUqng4eGBUaNGYfz48QWe7+bmhi+//BIjRozIO/buu+/C2toaa9asKdE5ASArKwtZWVl5j1NTU+Hh4cE+QmTQJAkYPVrssG1hIUaCXn9d7qqIiMpOSfoIyXZpLDs7G8ePH4e/v//TYszM4O/vj8PP7or3jKysLFg9Nx5qbW2NgwcPlvicABAWFgYHB4e8m4eHR2l+NCK98M03IgQpFMDq1QxBRETqyBaEkpKSoFQq4ezsnO+4s7MzEhIS1L4mICAAs2bNwqVLl6BSqbBnzx5EREQgPj6+xOcEgAkTJiAlJSXvFpe7Ox+RgVq5UjRKBIDwcKBXLzmrISLSX7JPltbEnDlz4OPjg7p168LCwgIjR47E4MGDYaZuVz4NWFpawt7ePt+NyFDt2AEMHSruf/65uDxGRETqyRaEHB0dYW5ujsTExHzHExMT4VLI1sBVq1bF1q1bkZ6ejhs3buC///6DnZ0datasWeJzEhmTI0eAnj3FJOmBA4GwMLkrIiLSb7IFIQsLCzRp0gRRUVF5x1QqFaKiotCiRYsiX2tlZYVq1aohJycHmzdvRpcuXUp9TiJDd/Ei0KmT2M07IEBsqsqu0URERZN1+XxISAiCg4PRtGlTNG/eHOHh4UhPT8fgwYMBAAMHDkS1atUQ9uTX2iNHjuDWrVto1KgRbt26hSlTpkClUuGzzz4r9jmJjFF8vAg/9+4BTZsCmzaxazQRUXHIGoSCgoJw9+5dTJo0CQkJCWjUqBEiIyPzJjvHxsbmm/+TmZmJiRMn4urVq7Czs0PHjh2xevVqVKxYsdjnJDI2KSlAhw7A9etig8rffwfs7OSuiojIMMjaR0hflaQPAZEcsrLE5ql79wJOTmJn+SdT5oiITI5B9REiotJRqYDgYBGC7OxE12iGICIizTAIERkgSQJCQoANG8RcoIgI4JVX5K6KiMjwMAgRGaDvvgPmzBH3V64E3npL1nKIiAwWgxCRgVm9GshdKPn990DfvvLWQ0RkyBiEiAzIrl3AkCHi/rhx4vIYERGVHIMQkYE4dgx4910gJ0eMAs2cKXdFRESGj0GIyABcviyWyaenA/7+wIoVQCm32CMiIjAIEem9xETRNfruXbEyLCICsLCQuyoiIuPAIESkxx4+FCNBV6+KHkE7dgAVKshdFRGR8WAQItJT2dliTtCJE0DVqkBkJMCdYoiItItBiEgPqVRiddiePYCtrdg/zMdH7qqIiIwPgxCRHvr8c2DtWqBcObGTfLNmcldERGScGISI9MysWaJzNAD8+CPQvr289RARGTMGISI9sm6daJQIADNmAAMGyFsPEZGxYxAi0hN//CF2kweAMWOATz+Vtx4iIlPAIESkB2JigG7dgMePgV69xOUxhULuqoiIjB+DEJHMrl4FOnQA0tKA118HVq1i12giIl3hP7dEMrpzR3SNTkwEGjYEtmwBLC3lroqIyHQwCBHJJC0N6NxZ7CNWowawcyfg4CB3VUREpoVBiEgGjx8DPXsCR48CVaoAu3YBrq5yV0VEZHoYhIh0TJKAoUPFlhk2NqJrdJ06cldFRGSaGISIdOyLL8SEaHNzYONGwM9P7oqIiEwXgxCRDs2dC3zzjbi/dCnQqZO89RARmToGISId2bgRGDtW3J82DRg8WNZyiIgIDEJEOvHnn2K7DEkCRowAJkyQuyIiIgIYhIjK3KlTQNeuQHY28O67wJw57BpNRKQvGISIytD166JrdGoq0LYtsGaNmCRNRET6gUGIqIwkJQHt2wPx8cDLLwPbtgFWVnJXRUREz2IQIioD6emia/SFC4CHh+gZVLGi3FUREdHzGISItCwnBwgKAo4cASpVEl2jq1WTuyoiIlKHQYhIiyQJGD5cdIu2sgK2bwfq1ZO7KiIiKgyDEJEWTZoE/PgjYGYGbNgAtGwpd0VERFQUBiEiLVmwAPjf/8T9RYuAd96Rtx4iInoxBiEiLYiIAEaOFPe/+goYNkzeeoiIqHgYhIhK6a+/gL59n84PCg2VuyIiIiouBiGiUjhzRlwCy8oS3aPnz2fXaCIiQ8IgRFRCsbGia3RKCtCqFfDzz+waTURkaBiEiErg/n3RNfrWLaB+feDXXwFra7mrIiIiTTEIEWno0SMgMBA4f140SoyMBCpXlrsqIiIqCQYhIg3k5AC9ewOHDoktM3btEltoEBGRYWIQIiomSQJGjBCXwSwtxZ8vvSR3VUREVBoaB6HXXnsNq1atwqNHj8qiHiK99fXXwJIlomv0unVAmzZyV0RERKWlcRBq3LgxPvnkE7i4uGDYsGH4559/yqIuIr2yZAkwZYq4P38+0K2brOUQEZGWaByEwsPDcfv2baxYsQJ37txB27ZtUb9+fXz33XdITEwsixqJZLVtG/Dhh+J+aCjwwQfy1kNERNpTojlC5cqVQ/fu3bFt2zbcvHkTffv2RWhoKDw8PNC1a1fs3btX23USyeLvv8XkaJUKeO89sX0GEREZj1JNlo6OjsbkyZPx/fffw8nJCRMmTICjoyM6d+6MTz75RFs1Esni3DmxTD4zE+jcWWykyq7RRETGRSFJkqTJC+7cuYPVq1djxYoVuHTpEgIDAzF06FAEBARA8eRT4uDBg2jfvj3S0tLKpOiylpqaCgcHB6SkpMDe3l7uckgGN28CLVsCcXHAq68CUVGAjY3cVRERUVFK8vldTtNv4u7uDm9vbwwZMgSDBg1C1apVCzynQYMGaNasmaanJtILDx6IrTPi4oA6dYDt2xmCiIiMlcZBKCoqCm1esG7Y3t4ef/75Z4mLIpJLZibQpQtw9izg6ioaJlapIndVRERUVjSeI+Tu7o5Lly4VOH7p0iVcv35dGzURyUKpBPr1Aw4cAOztxdYZNWrIXRUREZUljYPQoEGDcOjQoQLHjxw5gkGDBmmjJiKdkyRg9GggIgKwsBBL5hs0kLsqIiIqaxoHoZiYGLRq1arA8VdffRUnT57URk1EOjd9OrBggVgVtmYN0K6d3BUREZEuaByEFAoFHj58WOB4SkoKlEqlVooi0qUffwQmThT358wBevaUtx4iItIdjYNQ27ZtERYWli/0KJVKhIWFoXXr1lotjqisbd8OvP++uD9hAjBqlLz1EBGRbmm8amzGjBlo27Yt6tSpk7d67MCBA0hNTWVHaTIo//wD9OolJkkPGgRMmyZ3RUREpGsajwjVr18fp0+fRq9evXDnzh08fPgQAwcOxH///YeXX365LGok0rr//gM6dQIePQI6dhSbqrJrNBGR6dG4s7QpYGdp43b7tugafeMG0Lw5sHcvYGsrd1VERFRaOuksnSsjIwOxsbHIzs7Od7wB1xyTHktJEV2jb9wAatcGfv+dIYiIyJRpHITu3r2LwYMHY+fOnWq/zpVjpK9yu0afPg24uIiGiY6OcldFRERy0niO0NixY5GcnIwjR47A2toakZGR+Omnn+Dj44Nff/21LGokKjWlEhgwANi/H6hQAdi5E/DykrsqIiKSm8YjQnv37sW2bdvQtGlTmJmZoUaNGnjrrbdgb2+PsLAwdOrUqSzqJCoxSQLGjgU2bQLKlwe2bgUaNZK5KCIi0gsajwilp6fDyckJAFCpUiXcvXsXAODr64sTJ05otzoiLZgxA5g3T9xfvRp44w156yEiIv2hcRCqU6cOLly4AABo2LAhFi9ejFu3bmHRokVwdXXVuID58+fD09MTVlZW8PPzQ3R0dJHPDw8PR506dWBtbQ0PDw98/PHHyMzMzPv6lClToFAo8t3q1q2rcV1kHFauFI0SASA8HAgKkrMaIiLSNxpfGhszZgzi4+MBAJMnT0b79u2xdu1aWFhYYOXKlRqda8OGDQgJCcGiRYvg5+eH8PBwBAQE4MKFC3mjTs/6+eefMX78ePz4449o2bIlLl68iEGDBkGhUGDWrFl5z3vppZfwxx9/PP0hy5V4cRwZsJ07gaFDxf3PPgPGjJG3HiIi0j+l7iOUkZGB//77D9WrV4ejhktw/Pz80KxZM8x7ct1CpVLBw8MDo0aNwvjx4ws8f+TIkTh//jyioqLyjo0bNw5HjhzBwYMHAYgRoa1bt2q0AWxWVhaysrLyHqempsLDw4N9hAxYdDTw+utARoaYJL1yJWCm8fgnEREZkpL0EdLoo+Hx48fw9vbG+fPn847Z2NjglVde0TgEZWdn4/jx4/D3939ajJkZ/P39cfjwYbWvadmyJY4fP553+ezq1avYsWMHOnbsmO95ly5dgpubG2rWrIl+/fohNja2yFrCwsLg4OCQd/Pw8NDoZyH9cvGi6BqdkQEEBADLlzMEERGRehp9PJQvXz7ffJzSSEpKglKphLOzc77jzs7OSEhIUPuavn374uuvv0br1q1Rvnx5eHt7o127dvjiiy/ynuPn54eVK1ciMjISCxcuxLVr19CmTRs8fPiw0FomTJiAlJSUvFtcXJxWfkbSvYQEEX6SkoCmTZ+uFCMiIlJH49+TR4wYgRkzZiAnJ6cs6inSvn37MH36dCxYsAAnTpxAREQEfv/9d0ydOjXvOR06dEDPnj3RoEEDBAQEYMeOHUhOTsbGjRsLPa+lpSXs7e3z3cjwpKaKrtHXrwPe3qJrtJ2d3FUREZE+03gW8dGjRxEVFYXdu3fD19cXts/tTxAREVGs8zg6OsLc3ByJiYn5jicmJsLFxUXta0JDQzFgwAAMfTID1tfXF+np6Xj//ffx5ZdfwkzN9Y+KFSuidu3auHz5crHqIsOUlQV07w6cPAk4OQG7dok/iYiIiqLxiFDFihXx7rvvIiAgAG5ubvnm1jg4OBT7PBYWFmjSpEm+ic8qlQpRUVFo0aKF2tdkZGQUCDvm5uYAgMLmfKelpeHKlSslWtpPhkGlAgYNAqKixAjQjh1iRIiIiOhFNB4RWrFihda+eUhICIKDg9G0aVM0b94c4eHhSE9Px+DBgwEAAwcORLVq1RAWFgYACAwMxKxZs9C4cWP4+fnh8uXLCA0NRWBgYF4g+uSTTxAYGIgaNWrg9u3bmDx5MszNzdGnTx+t1U36Q5KAceOA9euBcuWAiAigSRO5qyIiIkMha4OdoKAg3L17F5MmTUJCQgIaNWqEyMjIvAnUsbGx+UaAJk6cCIVCgYkTJ+LWrVuoWrUqAgMDMW3atLzn3Lx5E3369MG9e/dQtWpVtG7dGv/88w+qVq2q85+Pyt7334tGiYBYIv/WW3JWQ0REhkbjPkJeXl5QKBSFfv3q1aulLkpuJelDQLq3Zo3oEQQA330nRoaIiMh0leTzW+MRobFjx+Z7/PjxY8TExCAyMhKffvqppqcjKpHdu4EnV1AREsIQREREJVOiLTbUmT9/Po4dO1bqgohe5NgxsUIsJwfo0wf49lu5KyIiIkOltX67HTp0wObNm7V1OiK1Ll8GOnYE0tMBf39unUFERKWjtY+QTZs2oXLlyto6HVEBiYlA+/bA3btA48bA5s2AhYXcVRERkSHT+NJY48aN802WliQJCQkJuHv3LhYsWKDV4ohyPXwo9g+7cgXw8hK9gjiPnYiISkvjINS1a9d8j83MzFC1alW0a9cOdevW1VZdRHmys4EePYDjxwFHR9E1upDm40RERBrROAhNnjy5LOogUkulAoYMEavEbGzESJCPj9xVERGRsdB4jtCOHTuwa9euAsd37dqFnTt3aqUoolzjxwNr14qu0Zs3A82ayV0REREZE42D0Pjx46FUKgsclyQJ48eP10pRRAAwe/bTpfHLl4uJ0kRERNqkcRC6dOkS6tevX+B43bp1ucM7ac369aJRIgB88w0wcKC89RARkXHSOAg5ODio3Ubj8uXLsLW11UpRZNqiop4Gn9Gjgc8+k7ceIiIyXhoHoS5dumDs2LG4cuVK3rHLly9j3LhxeOedd7RaHJmemBigWzfg8WOgVy9xeayIre2IiIhKReMgNHPmTNja2qJu3brw8vKCl5cX6tWrhypVquC7774rixrJRFy7BnToIHoGtWsHrFrFrtFERFS2NF4+7+DggEOHDmHPnj04deoUrK2t0aBBA7Rt27Ys6iMTcfcuEBAgukc3aABs3QpYWspdFRERGTuFJEmS3EXom9TUVDg4OCAlJQX2bF9c5tLSgDfeAI4eBWrUAA4dAtzc5K6KiIgMTUk+vzW+8DB69GjMnTu3wPF58+Zh7Nixmp6OTFzuXKCjR4EqVUTXaIYgIiLSFY2D0ObNm9GqVasCx1u2bIlNmzZppSgyDZIEDBsG7NwJWFsD27cDderIXRUREZkSjYPQvXv34ODgUOC4vb09kpKStFIUmYYvvwR++gkwNwc2bgRefVXuioiIyNRoHIRq1aqFyMjIAsd37tyJmjVraqUoMn4//ACEhYn7S5YAnTvLWw8REZkmjVeNhYSEYOTIkbh79y7eeOMNAEBUVBS+//57hIeHa7s+MkK//AKMGSPu/+9/YlNVIiIiOWgchIYMGYKsrCxMmzYNU6dOBQB4enpi4cKFGMh9EOgF9u0D+vcX84M++gj44gu5KyIiIlNWquXzd+/ehbW1Nezs7AAA9+/fR+XKlbVWnFy4fL5snD4NtGkDpKYC3buLeUHm5nJXRURExkIny+efVbVqVdjZ2WH37t3o1asXqlWrVprTkRG7cUPsHp+aCrRtC6xdyxBERETyK3EQunHjBiZPngxPT0/07NkTZmZmWLVqlTZrIyNx757oGh0fD7z8MrBtG2BlJXdVREREGs4Rys7ORkREBJYtW4a///4b/v7+uHnzJmJiYuDr61tWNZIBy8gQK8IuXAA8PETPoIoV5a6KiIhIKPaI0KhRo+Dm5oY5c+agW7duuHnzJn777TcoFAqY8xoHqZGTAwQFAf/8A1SqBERGAu7ucldFRET0VLFHhBYuXIjPP/8c48ePR4UKFcqyJjICkgR88IHoFm1lBfz2G1C/vtxVERER5VfsEaHVq1cjOjoarq6uCAoKwvbt26FUKsuyNjJgkycDy5cDZmbA+vWAml1ZiIiIZFfsINSnTx/s2bMHZ86cQd26dTFixAi4uLhApVLh3LlzZVkjGZiFC4EnLaawcCHQpYu89RARERWmxH2EJEnC7t27sXz5cvz6669wdHRE9+7d1e5Mb2jYR6jkIiKAHj3EpbEpU8TIEBERkS6U5PO7VA0Vc92/fx+rVq3CihUrcOrUqdKeTnYMQiVz4ADw1ltAVhbw/vvAokWAQiF3VUREZCpkC0LGhkFIc2fPiq7RycniUtimTUA5jTdwISIiKjmdd5YmAoC4ONE1OjkZaNkSWLeOIYiIiAwDgxCVyv37omv0rVtAvXpimby1tdxVERERFQ+DEJXYo0fAO+8A588D1aqJholGsOcuERGZEAYhKpGcHKBPH+DvvwEHBxGCqleXuyoiIiLNFGsmx+nTp4t9wgYNGpS4GDIMkgSMGCE2T7W0BH79VWymSkREZGiKFYQaNWoEhUIBSZKgeMF6aHabNn5TpwJLloil8T//DLRtK3dFREREJVOsS2PXrl3D1atXce3aNWzevBleXl5YsGABYmJiEBMTgwULFsDb2xubN28u63pJZkuXPm2SOH8+0L27vPUQERGVRrFGhGrUqJF3v2fPnpg7dy46duyYd6xBgwbw8PBAaGgounbtqvUiST/8+qvYSBUAJk4EPvxQ3nqIiIhKS+PJ0mfOnIGXl1eB415eXtxzzIgdOgQEBQEqFTBkCPD113JXREREVHoaB6F69eohLCwM2dnZeceys7MRFhaGevXqabU40g/nzwOdOwOZmeLPxYu5dQYRERkHjfv/Llq0CIGBgXB3d89bIXb69GkoFAr89ttvWi+Q5HXrlmiY+OAB8OqrwIYN7BpNRETGo0R7jaWnp2Pt2rX477//AIhRor59+8LW1lbrBcqBe40Jycli/7CzZ4E6dYCDBwFHR7mrIiIiUq8kn98l+t3e1tYW77//fkleSgYiM1Nsnnr2LODqCuzaxRBERETGp0SdpVevXo3WrVvDzc0NN27cAADMnj0b27Zt02pxJA+lEujfH/jrL8DeXnSNfmbhIBERkdHQOAgtXLgQISEh6NChAx48eJDXQLFSpUoIDw/Xdn2kY5IEjBkDbN4MWFiI7tFsFk5ERMZK4yD0ww8/YOnSpfjyyy9R7plZs02bNsWZM2e0WhzpXliYaJSoUABr1gDt2sldERERUdnROAhdu3YNjRs3LnDc0tIS6enpWimK5LFiBfDll+L+nDlAz57y1kNERFTWNA5CXl5eOHnyZIHjkZGR7CNkwH7/HRg2TNyfMAEYNUreeoiIiHRB41VjISEhGDFiBDIzMyFJEqKjo7Fu3TqEhYVh2bJlZVEjlbEjR8Toj1IJBAcD06bJXREREZFuaByEhg4dCmtra0ycOBEZGRno27cv3NzcMGfOHPTu3bssaqQydOEC0KkT8OgR0KGD2FSVXaOJiMhUlKihYq6MjAykpaXByclJmzXJzlQaKt6+DbRsCdy4ATRrBvz5J2AkPTGJiMgE6aSh4rVr15CTkwMfHx/Y2NjAxsYGAHDp0iWUL18enp6emp6SZJCSIkaAbtwAfHzEHCGGICIiMjUaT5YeNGgQDh06VOD4kSNHMGjQIG3URGUsKwvo2hU4fRpwcRFdo6tWlbsqIiIi3dM4CMXExKBVq1YFjr/66qtqV5ORflGpgAEDgH37gAoVgB07AC8vuasig6dUir9U69aJP580WiUi0ncaXxpTKBR4+PBhgeMpKSl5XaZJP0kSMHYs8MsvQPnywJYtgJqWUESaiYgQ7chv3nx6zN1dNKPq3l2+uoiIikHjEaG2bdsiLCwsX+hRKpUICwtD69attVocadfMmcAPP4j7q1YBb74pbz1kBCIigB498ocgALh1SxyPiJCnLjIdHI2kUtJ41di5c+fQtm1bVKxYEW3atAEAHDhwAKmpqdi7dy9efvnlMilUl4xx1dhPPwG5U7hmzxYjQ0SlolQCnp4FQ1AuhUKMDF27Bpib67Q0MhEcjaTnlOTzu0TL52/fvo158+bh1KlTsLa2RoMGDTBy5EhUrlxZ46L1kbEFoZ07gcBA8bn16adiZIio1PbtA15//cXPq14dsLcHzMxEIDI3L9794j7P0O+/6HlmZmzupU7uaOTzH2G579WmTQxDJkhnQcjYGVMQio4Wn1UZGUD//mJkyEzjC6JEaqxbB/TtK3cVpiE3EMkd2vQlzEqSGOJOSlL/fnE00mTppI8QACQnJyM6Ohp37tyBSqXK97WBAwdqdK758+fj22+/RUJCAho2bIgffvgBzZs3L/T54eHhWLhwIWJjY+Ho6IgePXogLCwMVlZWJT6nsbp0SXSNzsgA3n4bWL6cIYi0yNW1eM+bPRvw9RVLFpVKcTOU+7r4HsX5XVSlErecnNL9NzMVkgTExQFdugB+fmJU0sPj6e2ZzwsijYPQb7/9hn79+iEtLQ329vZQPDNkq1AoNApCGzZsQEhICBYtWgQ/Pz+Eh4cjICAAFy5cUNut+ueff8b48ePx448/omXLlrh48SIGDRoEhUKBWbNmleicxiohAQgIEL8wNWkiRoktLOSuioxGTo645lqU3N/KR43ib+VFkaSnwcgYgp0u7t+9K7rBvsjvv4vb86pWfRqOnv/Tw0OEfP6dNRkaXxqrXbs2OnbsiOnTp+d1lS4pPz8/NGvWDPPmzQMAqFQqeHh4YNSoURg/fnyB548cORLnz59HVFRU3rFx48bhyJEjOHjwYInOqY6hXxpLTQXatQNiYgBvb+DvvwFnZ7mrIqMRHw/06QPs3//0mEKRf2SD8zSoLBV3flpwsAg0sbFihCguTgyRv0i5coCbm/qQlHu/cmXO3dJDOrk0duvWLYwePbrUISg7OxvHjx/HhAkT8o6ZmZnB398fhw8fVvuali1bYs2aNYiOjkbz5s1x9epV7NixAwMGDCjxOQEgKysLWVlZeY9TU1NL9bPJKTtbfO7ExABOTqJrNEMQac3evSIE3bkD2NkBy5aJplTqVu6EhzMEUdlo00b8Hbt1S/2lxdzRyOXL84/sSBJw/74IRLnh6Pk/b90SI56xseJWGBubwkNS7p+l/Jwk3dA4CAUEBODYsWOoWbNmqb5xUlISlEolnJ/7lHZ2dsZ///2n9jV9+/ZFUlISWrduDUmSkJOTgw8++ABffPFFic8JAGFhYfjqq69K9fPoA5VKzB+MihL7hu3YIUaEiEpNpQKmTwcmTxb3fX1FZ846dcTXu3QBDhwQo0WuruKDipcWqKyYm4sl8j16FD4aGR5e8O+gQgFUqSJujRqpP7dSKf4eqwtJuffv3hUjSxcuiFthqlQpOiy5uYnRJ5KVxv8FOnXqhE8//RTnzp2Dr68vypcvn+/r77zzjtaKe96+ffswffp0LFiwAH5+frh8+TLGjBmDqVOnIjQ0tMTnnTBhAkJCQvIep6amwsPDQxsl69Snn4qFPOXKiZWlTZrIXREZhaQksS9LZKR4PHgwMG9e/t92zc3F9VgiXeneXVx61fZopLm5OIe7O9CihfrnPHokvqe6kJT7Z1oacO+euBW2/ZSZmQhDRYUlR0degitjGgehYcOGAQC+/vrrAl9TKBTF3mbD0dER5ubmSExMzHc8MTERLi4ual8TGhqKAQMGYOjQoQAAX19fpKen4/3338eXX35ZonMCgKWlJSwtLYtVt776/nvgyXxxrFghVokRldrhw0CvXuIffSsrYMECEYSI9EH37vKMRlpbAz4+4qaOJAEpKYWPKMXFif+nHj8Wf968Kf5fU8fKKv+KN3XzlipUKLuf1QRoHISeXy5fUhYWFmjSpAmioqLQtWvXvHNHRUVh5MiRal+TkZEBs+fWf5s/+QsvSVKJzmkM1q4FPvlE3P/2W9EviKhUJEn8Vv3ZZ2K+hI+P+O27QQO5KyPKTx9HIxUKoGJFcSvs/xmVCkhMLHxEKS5OLP/NzBS9UC5dKvz7VayoPiTl/lmtGpcNF0HWi5MhISEIDg5G06ZN0bx5c4SHhyM9PR2Dn/zGOXDgQFSrVg1hYWEAgMDAQMyaNQuNGzfOuzQWGhqKwMDAvED0onMam927n26d8fHHwLhxspZDxiAlBRgy5Ok+YT17iknRBriCkkhvmZmJESxXV9HrSJ2sLDF5u6jJ3SkpQHKyuJ05o/48CgXg4lJ0WHJy0k2jOaVS7+YTligIpaenY//+/YiNjUV2dna+r40ePbrY5wkKCsLdu3cxadIkJCQkoFGjRoiMjMyb7BwbG5tvBGjixIlQKBSYOHEibt26hapVqyIwMBDTpk0r9jmNyfHjwLvvil/Ye/cGvvuOl5KplGJixATUq1fFarBZs4ARI/gXi0gOlpZAzZriVpjU1MJHlHJvWVkieMTHi+0G1LGwEPOiigpLDg6l+3n0dG84jfsIxcTEoGPHjsjIyEB6ejoqV66MpKQk2NjYwMnJCVevXi2rWnXGEPoIXbkCtGwpVjG/+aboGWbg05xITpIELF0KjB4t/tGsUQPYuBEwwY7sREZFksQqt8JGlOLigNu3i9fh3N6+8LlK1auLUFPYB5GO9obTyV5j7dq1Q+3atbFo0SI4ODjg1KlTKF++PPr3748xY8aguxH0DdH3IHTnjghBV66IFaD79/OqBZVCejrwwQfAmjXicadOwKpVomEcERm/x49FGCpqcvf9+8U7l5NTwZDk7i5+ybpzR/1rtLg3nE6CUMWKFXHkyBHUqVMHFStWxOHDh1GvXj0cOXIEwcHBRfbrMRT6HITS0sS8wOPHAS8v4NAhcemXqETOnxe/pZ07J/4BmjZN9GHgpnRE9Kz09KLbBcTFibYCpfHnn6We+K6TztLly5fPm7fj5OSE2NhY1KtXDw4ODoiLi9P0dKSB7GwxJ+j4cdFaIjKSIYhK4eefgfffF//AuboC69cDbdvKXRUR6SNbW6BuXXFTR5JEzyR1YenECeDixRd/j/h47dZcTBoHocaNG+Po0aPw8fHBa6+9hkmTJiEpKQmrV6/Gyy+/XBY1EsRKy/feE6vEbGzEnKDateWuigxSZiYwdiyweLF4/MYbIhQZ4YICItIRhUL8hu7oCDRunP9rxd0bztW1TEp7EY3Hv6dPnw7XJ8VOmzYNlSpVwocffoi7d+9iyZIlWi+QhAkTxBQOc3Mxp4xzWKlErl4FWrUSIUihAEJDRbpmCCKispK7N1xhq08VCjGXqE0b3daV++01nSNkCvRtjlB4uOgRBAArV4oNlYk0tnWraDqVkiL2QFqzBmjfXu6qiMgU5K4aA9TvDSfjqjHOiNRz69c/DUFhYQxBVAKPH4tOm926iRDUooXoF8QQRES6krs3XLVq+Y+7u2stBJVUsUaEGjduDEUxG6qdOHGi1EXJTV9GhPbuFZ9Vjx8Do0aJnlPsa0cauXkTCAoSywsBkapnzBDNEomIdK2MO0uX2aqx3H27SHdOngS6dhUhqGdPYPZshiDS0O7dQL9+Yvd4e3txXbVbN7mrIiJTpod7w3GOkBpyjwhduyYaJiYkiL8vO3eKDYiJikWpBL76Cvjf/8S1+MaNgV9+Aby95a6MiKhM6aSPEJWtu3eBgAARgho0EPNbGYKo2BITxShQVJR4PHy4mG3Pv0RERGppHISUSiVmz56NjRs3qt109X5x23BTAenpQOfOwKVLojP5zp2l3+OOTMiBA2I+UHy8aDa1eDHQv7/cVRER6TWNV4199dVXmDVrFoKCgpCSkoKQkBB0794dZmZmmDJlShmUaBoePwZ69RIbA1euDOzaBbi5yV0VGQSVCpg5UzQsi48H6tUDjh5lCCIiKgaNg9DatWuxdOlSjBs3DuXKlUOfPn2wbNkyTJo0Cf/8809Z1Gj0JEnsdLBjB2BtDWzfXngXc6J87t8HunQBPv9czA3q10+k6fr15a6MiMggaByEEhIS4OvrCwCws7NDSkoKAKBz5874/ffftVudiZg4USzoMTcHNm4UbV6IXujoUeCVV0RytrQUl8JWrwbs7OSujIjIYGgchNzd3RH/ZGM0b29v7N69GwBw9OhRWFpaarc6EzBvHjB9uri/eLGYI0RUJEkSf3FatQJu3ABq1hR9gt5/nz0WiIg0pHEQ6tatG6KerEgZNWoUQkND4ePjg4EDB2LIkCFaL9CYbdoEjB4t7k+dKjZVJSrSw4dAnz6iw+bjx6Iv0PHjYmSIiIg0Vuo+QocPH8bhw4fh4+ODwMBAbdUlK233EYqNFT3tnnXsGDBiBJCTI+a0rlrFX+bpBc6cEXv1XLwIlCsnOkR//DH/4hARPVGSz282VFRDm0EoNhaoUwfIzCz8OVZWwIULYsk8kVorVwIffQQ8eiT25tmwQXTdJCKiPDppqHjv3j1UqVIFABAXF4elS5fi0aNHeOedd9CmTRtNT2f0kpKKDkGA+HpSEoMQqZGRIS6D/fijeBwQIHaNd3SUty4iIiNR7DlCZ86cgaenJ5ycnFC3bl2cPHkSzZo1w+zZs7FkyRK8/vrr2Lp1axmWSmRiLl4EXn1VhCAzMzGRbMcOhiAiIi0qdhD67LPP4Ovri7/++gvt2rVD586d0alTJ6SkpODBgwcYPnw4vvnmm7Kslch0/PIL0LSpmBfk5ATs2SP6LJhpvL6BiIiKUOxLY0ePHsXevXvRoEEDNGzYEEuWLMFHH30Esyf/MI8aNQqvvvpqmRVKZBKys4FPPgF++EE8btsWWLeObcaJiMpIsX+9vH//PlxcXACIRoq2traoVKlS3tcrVaqEhw8far9CIlNx4wbQps3TEPT552LzVIYgIqIyo9FkacVzy3Sff0xEJfT778CAAcCDB0ClSqKfArtrEhGVOY2C0KBBg/K6R2dmZuKDDz6Ara0tACArK0v71REZu5wcIDQUyJ1f16yZ2GfF01PWsoiITEWxg1BwcHC+x/3V7Gw9cODA0ldkZBwdRZ+gF/UR4kIgExQfD/TuDfz1l3g8ciTw3Xdi3zAiItIJNlRUQxedpZ/l6MgeQiZn716xVcadO2KT1GXLgKAguasiIjJoOmmoSJqrXp1Bh55QqcQuu5Mni/u+vmLTudq15a6MiMgkMQgR6UpSkthYbtcu8XjwYLGLvI2NvHUREZkwBiEiXTh8GOjVC7h5U0wKW7BABCEiIpIV29QSlSVJAmbPFo0Rb94Ul8COHGEIIiLSExwRIiorycnAkCHAli3ica9ewNKlgBYm4BMRkXYwCBGVhZgYoEcP4OpVoHx5YNYsYMQIgE1IiYj0CoMQkTZJkhj1GT0ayMoCatQQDRKbN5e7MiIiUoNzhIi0JT0dGDgQGD5chKDOnYETJxiCiIj0GIMQkTacOycCz5o1gLm52DJj2zagcmW5KyMioiLw0hhRaa1dC7z/PpCRAbi6AuvXi1ViRESk9zgiRFRSmZnABx+IJokZGcAbb4hJ0gxBREQGg0GIqCSuXAFatgQWLxYrwUJDgd27AWdnuSsjIiIN8NIYkaa2bBENEVNSgCpVxKWxgAC5qyIiohLgiBBRcT1+DIwbB3TvLkJQy5biUhhDEBGRweKIEFFxxMUBQUFizzAACAkRK8PKl5e3LiIiKhUGIaIX2bUL6NcPuHdPbI+xciXQrZvcVRERkRbw0hhRYZRKYNIkoEMHEYIaNxYNEhmCiIiMBkeEiNRJTBSjQFFR4vHw4UB4OGBlJWtZRESkXQxCRM/76y+gd28gPh6wsRFL5Pv3l7sqIiIqA7w0RpRLpQJmzBCNEePjgXr1gKNHGYKIiIwYR4SIAOD+fSA4GNi+XTzu1w9YtAiws5O3LiIiKlMMQkRHjwI9ewI3bgCWlsDcucCwYaJjNBERGTVeGiPTJUnAvHlAq1YiBNWsCRw6JDZQZQgiIjIJHBEi0/TwITB0KLBxo3jcrRuwYgXg4CBvXUREpFMcESLTc+YM0LSpCEHlygGzZgGbNzMEERGZII4IkWlZuRL46CPg0SPA3R3YsEHsGUZERCaJI0JkGjIygCFDxK7xjx6JjVJjYhiCiIhMHIMQGb+LF4FXXxVzgMzMgKlTgR07AEdHuSsjIiKZ8dIYGbeNG4H33gPS0gAnJ2DdOtEwkYiICBwRImOVlQWMGgUEBYkQ1LatuBTGEERERM9gECLjc+MG0KaN6BEEAOPHi81T3dzkrYuIiPQOL42Rcdm+HRg4EHjwAKhUCVi1CujcWe6qiIhIT3FEiIxDTo4Y+QkMFCGoWTPgxAmGICIiKhJHhMjw3b4N9OkD/PWXeDxqFPDdd4CFhbx1ERGR3tOLEaH58+fD09MTVlZW8PPzQ3R0dKHPbdeuHRQKRYFbp06d8p4zaNCgAl9v3769Ln4U0rW9e4HGjUUIsrMTDRLnzmUIIiKiYpF9RGjDhg0ICQnBokWL4Ofnh/DwcAQEBODChQtwcnIq8PyIiAhkZ2fnPb537x4aNmyInj175nte+/btsWLFirzHlpaWZfdDkO6pVMC0acCUKeK+ry+waRNQu7bclRERkQGRfURo1qxZGDZsGAYPHoz69etj0aJFsLGxwY8//qj2+ZUrV4aLi0vebc+ePbCxsSkQhCwtLfM9r1KlSrr4cUgXkpKAjh2BSZNECBo8GPjnH4YgIiLSmKxBKDs7G8ePH4e/v3/eMTMzM/j7++Pw4cPFOsfy5cvRu3dv2Nra5ju+b98+ODk5oU6dOvjwww9x7969Qs+RlZWF1NTUfDfSU4cOiUthu3YB1tbAjz+Km42N3JUREZEBkjUIJSUlQalUwtnZOd9xZ2dnJCQkvPD10dHROHv2LIYOHZrvePv27bFq1SpERUVhxowZ2L9/Pzp06AClUqn2PGFhYXBwcMi7eXh4lPyHorIhSWKX+NdeA27eFKM/R46I0SAiIqISkn2OUGksX74cvr6+aN68eb7jvXv3zrvv6+uLBg0awNvbG/v27cObb75Z4DwTJkxASEhI3uPU1FSGIX2SnCw2TN2yRTzu1QtYuhSwt5e1LCIiMnyyjgg5OjrC3NwciYmJ+Y4nJibCxcWlyNemp6dj/fr1eO+99174fWrWrAlHR0dcvnxZ7dctLS1hb2+f70Z6IiYGaNJEhKDy5UW36PXrGYKIiEgrZA1CFhYWaNKkCaKiovKOqVQqREVFoUWLFkW+9pdffkFWVhb69+//wu9z8+ZN3Lt3D66urqWumXREkoDFi4EWLYCrV4EaNYC//wZGjAAUCrmrIyIiIyH7qrGQkBAsXboUP/30E86fP48PP/wQ6enpGPxk7sfAgQMxYcKEAq9bvnw5unbtiipVquQ7npaWhk8//RT//PMPrl+/jqioKHTp0gW1atVCQECATn4mKqW0NGDAAOCDD8TmqZ07iy7RzZrJXRkRERkZ2ecIBQUF4e7du5g0aRISEhLQqFEjREZG5k2gjo2NhZlZ/rx24cIFHDx4ELt37y5wPnNzc5w+fRo//fQTkpOT4ebmhrfffhtTp05lLyFDcO4c0KMHcP48YG4OTJ8OfPIJYCZ7ZiciIiOkkCRJkrsIfZOamgoHBwekpKRwvpAurV0LvP8+kJEBuLqKuUBt28pdFRERGYiSfH7z12ySX2YmMHw40L+/CEFvvikmSTMEERFRGWMQInlduQK0bAksWSImQYeGimaJz/WWIiIiKguyzxEiE7Zli2iImJICODoCa9YAnNBOREQ6xBEh0r3Hj4Fx44Du3UUIatlSXApjCCIiIh1jECLdiosT22TMmiUejxsH7NsHuLvLWhYREZkmXhoj3dm1C+jXD7h3D3BwAFauBLp2lbsqIiIyYRwRorKnVAKTJgEdOogQ1LgxcPw4QxAREcmOI0JUthITgb59gb17xePhw4HwcMDKStayiIiIAAYhKkt//QX07g3ExwM2NmKJfL9+cldFRESUh5fGSPtUKmDGDOCNN0QIql8fOHqUIYiIiPQOR4RIu+7fB4KDge3bxeP+/YFFiwBbW3nrIiIiUoNBiLQnOhro1Qu4cQOwtATmzgWGDRMdo4mIiPQQL41R6UkSMG8e0Lq1CEHe3sDhw2IDVYYgIiLSYxwRotJJTRWjPhs3isfdugErVog+QURERHqOI0JUcqdPA02bihBUrpzoFr15M0MQEREZDI4IUcmsWAF89BGQmSm2x9i4EWjRQu6qiIiINMIRIdJMRgYwZIi4ZWaKjVJjYhiCiIjIIDEIUfFdvAi8+qoYDTIzA6ZOBXbsABwd5a6MiIioRHhpjIpn40bgvfeAtDTAyQlYt040TCQiIjJgHBGiomVlAaNGAUFBIgS1bQucPMkQRERERoFBiAp3/TrQpo3oEQQA48cDUVGAq6usZREREWkLL42Retu3AwMHAg8eAJUqAatXA506yV0VERGRVnFEiPLLyREjP4GBIgQ1by5WhTEEERGREeKIED11+zbQuzdw4IB4PGoU8N13gIWFvHURERGVEQYhEqKigL59gTt3gAoVgOXLgZ495a6KiIioTPHSmKlTqUQ/oLfeEiGoQQPg2DGGICIiMgkcETJld+8CAwYAu3aJx0OGiBVi1tby1kVERKQjDEKm6tAhoFcv4NYtEXwWLAAGDZK7KiIiIp3ipTFTI0lil/jXXhMhqHZt4MgRhiAiIjJJHBEyJcnJ4vLXli3icVAQsHSpmBxNRERkghiETMWJE2IC9NWrQPnywOzZwEcfAQqF3JURERHJhkHI2EkSsGQJMGaM2DesRg3gl1+AZs3kroyIiEh2nCNkzNLSxKqwDz4QIahzZzEyxBBEREQEgEHIeJ07J7bHWLsWMDcHZswAtm0DKleWuzIiIiK9wUtjxmjNGmD4cCAjQ+wUv3490Lat3FURERHpHY4IGZPMTBGABgwQIejNN4GTJxmCiIiICsEgZCyuXAFathQToxUKYNIk0THayUnuyoiIiPQWL40Zgy1bgMGDgZQUwNFRzAt6+225qyIiItJ7HBEyZI8fAyEhQPfuIgS1bAnExDAEERERFRODkKGKixPbZMyeLR6PGwfs2we4u8taFhERkSHhpTFDFBkJ9O8P3LsHODgAK1cCXbvKXRUREZHB4YiQIVEqgdBQoGNHEYIaNxYNEhmCiIiISoQjQoYiMRHo2xfYu1c8/uADcVnMykreuoiIiAwYg5Ah+OsvoHdvID4esLUVS+T79pW7KiIiIoPHS2P6TKUCvvkGeP11EYLq1weOHmUIIiIi0hKOCOmr+/eB4GBg+3bxuH9/YNEiMSJEREREWsEgpI+io4FevYAbNwBLS+CHH4ChQ0XHaCIiItIaXhrTJ5IkQk/r1iIEeXsDhw8Dw4YxBBEREZUBjgjpi9RUMerzyy/icbduwIoVok8QERERlQmOCOmD06eBpk1FCCpXTiyL37yZIYiIiKiMcURIl5RK4MABsQLM1RVo0wZYtQr46CMgM1Nsj7FxI9CihdyVEhERmQQGIV2JiADGjAFu3nx6zMYGyMgQ99u3B1avFrvHExERkU4wCOlCRATQo4eYDP2s3BDUpw+wZg1gxiuVREREusRP3rKmVIqRoOdD0LMOHiz660RERFQmGITK2oED+S+HqRMXJ55HREREOsUgVNbi47X7PCIiItIaBqGy5uqq3ecRERGR1jAIlbU2bcSy+MI6QysUgIeHeB4RERHpFINQWTM3B+bMEfefD0O5j8PDxfOIiIhIpxiEdKF7d2DTJqBatfzH3d3F8e7d5amLiIjIxLGPkK507w506VKwszRHgoiIiGTDIKRL5uZAu3ZyV0FERERP6MWlsfnz58PT0xNWVlbw8/NDdHR0oc9t164dFApFgVunTp3yniNJEiZNmgRXV1dYW1vD398fly5d0sWPQkRERAZE9iC0YcMGhISEYPLkyThx4gQaNmyIgIAA3LlzR+3zIyIiEB8fn3c7e/YszM3N0bNnz7znzJw5E3PnzsWiRYtw5MgR2NraIiAgAJmZmbr6sYiIiMgAKCRJ3r0d/Pz80KxZM8ybNw8AoFKp4OHhgVGjRmH8+PEvfH14eDgmTZqE+Ph42NraQpIkuLm5Ydy4cfjkk08AACkpKXB2dsbKlSvRu3fvF54zNTUVDg4OSElJgb29fel+QCIiItKJknx+yzoilJ2djePHj8Pf3z/vmJmZGfz9/XH48OFinWP58uXo3bs3bG1tAQDXrl1DQkJCvnM6ODjAz8+v0HNmZWUhNTU1342IiIiMn6xBKCkpCUqlEs7OzvmOOzs7IyEh4YWvj46OxtmzZzF06NC8Y7mv0+ScYWFhcHBwyLt5eHho+qMQERGRAZJ9jlBpLF++HL6+vmjevHmpzjNhwgSkpKTk3eLi4rRUIREREekzWYOQo6MjzM3NkZiYmO94YmIiXFxcinxteno61q9fj/feey/f8dzXaXJOS0tL2Nvb57sRERGR8ZM1CFlYWKBJkyaIiorKO6ZSqRAVFYUWLVoU+dpffvkFWVlZ6N+/f77jXl5ecHFxyXfO1NRUHDly5IXnJCIiItMie0PFkJAQBAcHo2nTpmjevDnCw8ORnp6OwYMHAwAGDhyIatWqISwsLN/rli9fjq5du6JKlSr5jisUCowdOxb/+9//4OPjAy8vL4SGhsLNzQ1du3bV1Y9FREREBkD2IBQUFIS7d+9i0qRJSEhIQKNGjRAZGZk32Tk2NhZmZvkHri5cuICDBw9i9+7das/52WefIT09He+//z6Sk5PRunVrREZGwsrKqlg15XYU4OoxIiIiw5H7ua1JZyDZ+wjpo5s3b3LlGBERkYGKi4uDu7t7sZ7LIKSGSqXC7du3UaFCBSgUCq2eOzU1FR4eHoiLizPpSdl8H57ie/EU34un+F48xfdC4PvwVGHvhSRJePjwIdzc3ApcTSqM7JfG9JGZmVmxk2RJcXWawPfhKb4XT/G9eIrvxVN8LwS+D0+pey8cHBw0OodB9xEiIiIiKg0GISIiIjJZDEI6ZmlpicmTJ8PS0lLuUmTF9+EpvhdP8b14iu/FU3wvBL4PT2nzveBkaSIiIjJZHBEiIiIik8UgRERERCaLQYiIiIhMFoMQERERmSwGIR1YuHAhGjRokNf4qUWLFti5c6fcZemFb775Jm+jXFMzZcoUKBSKfLe6devKXZZsbt26hf79+6NKlSqwtraGr68vjh07JndZOufp6Vng74VCocCIESPkLk2nlEolQkND4eXlBWtra3h7e2Pq1Kka7SFlTB4+fIixY8eiRo0asLa2RsuWLXH06FG5yypzf/31FwIDA+Hm5gaFQoGtW7fm+7okSZg0aRJcXV1hbW0Nf39/XLp0SaPvwSCkA+7u7vjmm29w/PhxHDt2DG+88Qa6dOmCf//9V+7SZHX06FEsXrwYDRo0kLsU2bz00kuIj4/Pux08eFDukmTx4MEDtGrVCuXLl8fOnTtx7tw5fP/996hUqZLcpenc0aNH8/2d2LNnDwCgZ8+eMlemWzNmzMDChQsxb948nD9/HjNmzMDMmTPxww8/yF2aLIYOHYo9e/Zg9erVOHPmDN5++234+/vj1q1bcpdWptLT09GwYUPMnz9f7ddnzpyJuXPnYtGiRThy5AhsbW0REBCAzMzM4n8TiWRRqVIladmyZXKXIZuHDx9KPj4+0p49e6TXXntNGjNmjNwl6dzkyZOlhg0byl2GXvj888+l1q1by12GXhozZozk7e0tqVQquUvRqU6dOklDhgzJd6x79+5Sv379ZKpIPhkZGZK5ubm0ffv2fMdfeeUV6csvv5SpKt0DIG3ZsiXvsUqlklxcXKRvv/0271hycrJkaWkprVu3rtjn5YiQjimVSqxfvx7p6elo0aKF3OXIZsSIEejUqRP8/f3lLkVWly5dgpubG2rWrIl+/fohNjZW7pJk8euvv6Jp06bo2bMnnJyc0LhxYyxdulTusmSXnZ2NNWvWYMiQIVrfAFrftWzZElFRUbh48SIA4NSpUzh48CA6dOggc2W6l5OTA6VSCSsrq3zHra2tTXYUGQCuXbuGhISEfJ8jDg4O8PPzw+HDh4t9Hm66qiNnzpxBixYtkJmZCTs7O2zZsgX169eXuyxZrF+/HidOnDCJ69tF8fPzw8qVK1GnTh3Ex8fjq6++Qps2bXD27FlUqFBB7vJ06urVq1i4cCFCQkLwxRdf4OjRoxg9ejQsLCwQHBwsd3my2bp1K5KTkzFo0CC5S9G58ePHIzU1FXXr1oW5uTmUSiWmTZuGfv36yV2azlWoUAEtWrTA1KlTUa9ePTg7O2PdunU4fPgwatWqJXd5sklISAAAODs75zvu7Oyc97XiYBDSkTp16uDkyZNISUnBpk2bEBwcjP3795tcGIqLi8OYMWOwZ8+eAr/dmJpnf7Nt0KAB/Pz8UKNGDWzcuBHvvfeejJXpnkqlQtOmTTF9+nQAQOPGjXH27FksWrTIpIPQ8uXL0aFDB7i5ucldis5t3LgRa9euxc8//4yXXnoJJ0+exNixY+Hm5maSfydWr16NIUOGoFq1ajA3N8crr7yCPn364Pjx43KXZvB4aUxHLCwsUKtWLTRp0gRhYWFo2LAh5syZI3dZOnf8+HHcuXMHr7zyCsqVK4dy5cph//79mDt3LsqVKwelUil3ibKpWLEiateujcuXL8tdis65uroW+KWgXr16JnupEABu3LiBP/74A0OHDpW7FFl8+umnGD9+PHr37g1fX18MGDAAH3/8McLCwuQuTRbe3t7Yv38/0tLSEBcXh+joaDx+/Bg1a9aUuzTZuLi4AAASExPzHU9MTMz7WnEwCMlEpVIhKytL7jJ07s0338SZM2dw8uTJvFvTpk3Rr18/nDx5Eubm5nKXKJu0tDRcuXIFrq6ucpeic61atcKFCxfyHbt48SJq1KghU0XyW7FiBZycnNCpUye5S5FFRkYGzMzyf0SZm5tDpVLJVJF+sLW1haurKx48eIBdu3ahS5cucpckGy8vL7i4uCAqKirvWGpqKo4cOaLRHFxeGtOBCRMmoEOHDqhevToePnyIn3/+Gfv27cOuXbvkLk3nKlSogJdffjnfMVtbW1SpUqXAcWP3ySefIDAwEDVq1MDt27cxefJkmJubo0+fPnKXpnMff/wxWrZsienTp6NXr16Ijo7GkiVLsGTJErlLk4VKpcKKFSsQHByMcuVM85/pwMBATJs2DdWrV8dLL72EmJgYzJo1C0OGDJG7NFns2rULkiShTp06uHz5Mj799FPUrVsXgwcPlru0MpWWlpZvlPzatWs4efIkKleujOrVq2Ps2LH43//+Bx8fH3h5eSE0NBRubm7o2rVr8b+JFle2USGGDBki1ahRQ7KwsJCqVq0qvfnmm9Lu3bvlLktvmOry+aCgIMnV1VWysLCQqlWrJgUFBUmXL1+WuyzZ/Pbbb9LLL78sWVpaSnXr1pWWLFkid0my2bVrlwRAunDhgtylyCY1NVUaM2aMVL16dcnKykqqWbOm9OWXX0pZWVlylyaLDRs2SDVr1pQsLCwkFxcXacSIEVJycrLcZZW5P//8UwJQ4BYcHCxJklhCHxoaKjk7O0uWlpbSm2++qfH/NwpJMtE2nURERGTyOEeIiIiITBaDEBEREZksBiEiIiIyWQxCREREZLIYhIiIiMhkMQgRERGRyWIQIiIiIpPFIEREREQmi0GIiEzS9evXoVAocPLkSa2ds127dhg7dqzWzkdEZY9BiIjKVEJCAkaNGoWaNWvC0tISHh4eCAwMzLdRorGIiIjA1KlT8x57enoiPDxcvoKI6IVMczc/ItKJ69evo1WrVqhYsSK+/fZb+Pr64vHjx9i1axdGjBiB//77T+4Stapy5cpyl0BEGuKIEBGVmY8++ggKhQLR0dF49913Ubt2bbz00ksICQnBP//8AwCIjY1Fly5dYGdnB3t7e/Tq1QuJiYl555gyZQoaNWqEH3/8EdWrV4ednR0++ugjKJVKzJw5Ey4uLnBycsK0adPyfW+FQoGFCxeiQ4cOsLa2Rs2aNbFp06Yi6z179iw6dOgAOzs7ODs7Y8CAAUhKSgIA7Nu3DxYWFjhw4EDe82fOnAknJ6e8ep+9NNauXTvcuHEDH3/8MRQKBRQKBdLT02Fvb1+gjq1bt8LW1hYPHz4s2RtNRCXGIEREZeL+/fuIjIzEiBEjYGtrW+DrFStWhEqlQpcuXXD//n3s378fe/bswdWrVxEUFJTvuVeuXMHOnTsRGRmJdevWYfny5ejUqRNu3ryJ/fv3Y8aMGZg4cSKOHDmS73WhoaF49913cerUKfTr1w+9e/fG+fPn1dabnJyMN954A40bN8axY8cQGRmJxMRE9OrVC8DTkDNgwACkpKQgJiYGoaGhWLZsGZydnQucLyIiAu7u7vj6668RHx+P+Ph42Nraonfv3lixYkW+565YsQI9evRAhQoVNHqPiUgLNNqrnoiomI4cOSIBkCIiIgp9zu7duyVzc3MpNjY279i///4rAZCio6MlSZKkyZMnSzY2NlJqamrecwICAiRPT09JqVTmHatTp44UFhaW9xiA9MEHH+T7fn5+ftKHH34oSZIkXbt2TQIgxcTESJIkSVOnTpXefvvtfM+Pi4uTAEgXLlyQJEmSsrKypEaNGkm9evWS6tevLw0bNizf81977TVpzJgxeY9r1KghzZ49u8D7Ym5uLt2+fVuSJElKTEyUypUrJ+3bt6/Q94mIyg5HhIioTEiS9MLnnD9/Hh4eHvDw8Mg7Vr9+fVSsWDHfyI2np2e+0RJnZ2fUr18fZmZm+Y7duXMn3/lbtGhR4HFhI0KnTp3Cn3/+CTs7u7xb3bp1AYgRKQCwsLDA2rVrsXnzZmRmZmL27Nkv/Bmf17x5c7z00kv46aefAABr1qxBjRo10LZtW43PRUSlx8nSRFQmfHx8oFAotDIhunz58vkeKxQKtcdUKlWJv0daWhoCAwMxY8aMAl9zdXXNu3/o0CEA4tLf/fv31V72e5GhQ4di/vz5GD9+PFasWIHBgwdDoVCUuHYiKjmOCBFRmahcuTICAgIwf/58pKenF/h6cnIy6tWrh7i4OMTFxeUdP3fuHJKTk1G/fv1S15A7IfvZx/Xq1VP73FdeeQX//vsvPD09UatWrXy33LBz5coVfPzxx1i6dCn8/PwQHBxcZPiysLCAUqkscLx///64ceMG5s6di3PnziE4OLgUPyURlQaDEBGVmfnz50OpVKJ58+bYvHkzLl26hPPnz2Pu3Llo0aIF/P394evri379+uHEiROIjo7GwIED8dprr6Fp06al/v6//PILfvzxR1y8eBGTJ09GdHQ0Ro4cqfa5I0aMwP3799GnTx8cPXoUV65cwa5duzB48GAolUoolUr0798fAQEBGDx4MFasWIHTp0/j+++/L/T7e3p64q+//sKtW7fyVp8BQKVKldC9e3d8+umnePvtt+Hu7l7qn5WISoZBiIjKTM2aNXHixAm8/vrrGDduHF5++WW89dZbiIqKwsKFC6FQKLBt2zZUqlQJbdu2hb+/P2rWrIkNGzZo5ft/9dVXWL9+PRo0aIBVq1Zh3bp1hY40ubm54e+//4ZSqcTbb78NX19fjB07FhUrVoSZmRmmTZuGGzduYPHixQDE5bIlS5Zg4sSJOHXqlNpzfv3117h+/Tq8vb1RtWrVfF977733kJ2djSFDhmjlZyWiklFIxZnRSERkYBQKBbZs2YKuXbvKXYpaq1evxscff4zbt2/DwsJC7nKITBYnSxMR6VBGRgbi4+PxzTffYPjw4QxBRDLjpTEiIh2aOXMm6tatCxcXF0yYMEHucohMHi+NERERkcniiBARERGZLAYhIiIiMlkMQkRERGSyGISIiIjIZDEIERERkcliECIiIiKTxSBEREREJotBiIiIiEzW/wGgyTx+8Z19kQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(complexities, train_scores, \"s-b\", label=\"Train\")\n",
    "plt.plot(complexities, test_scores, \"o-r\", label=\"Test\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Complexity\")\n",
    "plt.ylabel(\"Balanced Accuracy\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2aa3f823",
   "metadata": {},
   "source": [
    "## Rule class diagram\n",
    "\n",
    "For the next two sections, it will be useful to refer to the rule class diagram:\n",
    "\n",
    "![title](rule_class_diagram.png)\n",
    "\n",
    "Note that only the leaves are concrete classes - specific operators, literals, or trivial rules. The rest of the classes are abstract classes. The `Operator` class is not actually implemented, it is only pictured for clearer presentation. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16f110b9",
   "metadata": {},
   "source": [
    "## Changing the allowed operators\n",
    "\n",
    "The rule classifier has an argument `operators` which by default includes all the operators provided in `boolxai`. But what if we want to allow only some of the operators? As an example, we might want to only include the unparameterized operators `And` and `Or`. We can do that like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3c26cbd6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from boolxai import Operator\n",
    "\n",
    "rule_classifier = BoolXAI.RuleClassifier(\n",
    "    random_state=43, operators=[Operator.And, Operator.Or]\n",
    ")\n",
    "rule_classifier.fit(X_binarized, y);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7054b17d",
   "metadata": {},
   "source": [
    "We can inspect the best rule seen, as well as the best rule for each start, to make sure that they indeed include only those operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ed64809d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best score and rule:\n",
      "  score=0.93 rule=Or(252, 222, 220, 201, 33)\n",
      "\n",
      "Best score and rule for each start:\n",
      "  score=0.71 rule=And(Or(~205, 42), ~109, ~267)\n",
      "  score=0.87 rule=Or(71, 170, 250, 200, 62)\n",
      "  score=0.93 rule=Or(252, 222, 220, 201, 33)\n",
      "  score=0.80 rule=Or(232, 190, 0)\n",
      "  score=0.76 rule=And(~38, Or(230, 23, 41))\n",
      "  score=0.89 rule=Or(220, 222, 272, 221)\n",
      "  score=0.80 rule=And(Or(192, 30, 272), ~9)\n",
      "  score=0.81 rule=Or(251, 130, 260, 261)\n",
      "  score=0.71 rule=Or(And(~126, 1, ~296), 270)\n",
      "  score=0.92 rule=Or(230, 150, 232, 252, 202)\n"
     ]
    }
   ],
   "source": [
    "print(\"Best score and rule:\")\n",
    "print(f\"  score={rule_classifier.best_score_:.2f} rule={rule_classifier.best_rule_}\\n\")\n",
    "\n",
    "print(\"Best score and rule for each start:\")\n",
    "for score, rule in zip(rule_classifier.scores_, rule_classifier.rules_):\n",
    "    print(f\"  {score=:.2f} {rule=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "983e6765",
   "metadata": {},
   "source": [
    "## Adding custom operators\n",
    "\n",
    "We can even add custom operators. To do so, we follow the template of the included operators in `boolxai.rules.operators`. In particular, we must subclass `UnparametrizedOperator` or `ParametrizedOperator`, and implement the `_evaluate()` method. Let's add the `AllEqual()` operator, which returns `True` only if all included literals are true:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "857dbc62",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import List, Optional, Union\n",
    "\n",
    "from boolxai import Literal\n",
    "from boolxai import UnparametrizedOperator\n",
    "\n",
    "\n",
    "class AllEqual(UnparametrizedOperator):\n",
    "    \"\"\"Returns True if all subrules are equal.\"\"\"\n",
    "\n",
    "    def _evaluate(self, state):\n",
    "        sum_subrules = np.sum(\n",
    "            np.column_stack([subrule.evaluate(state) for subrule in self.subrules]),\n",
    "            axis=1,\n",
    "        )\n",
    "        return np.logical_or(sum_subrules == 0, sum_subrules == len(self.subrules))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb0f2d23",
   "metadata": {},
   "source": [
    "Let's try our new operator to see if it works as expected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e33eca91",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([False, False,  True, False,  True])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.array([[0, 0, 1], [0, 0, 0], [0, 1, 1], [1, 1, 1], [1, 0, 0]])\n",
    "\n",
    "rule = AllEqual([Literal(0, negated=True), Literal(1), Literal(2)])\n",
    "rule.evaluate(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "426c3ea3",
   "metadata": {},
   "source": [
    "Indeed, the result we were expecting! Now we can include this operator class in the training of the rule classifier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "1542e443",
   "metadata": {},
   "outputs": [],
   "source": [
    "rule_classifier = BoolXAI.RuleClassifier(\n",
    "    random_state=43, operators=[Operator.And, AllEqual]\n",
    ")\n",
    "rule_classifier.fit(X_binarized, y);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa69b11c",
   "metadata": {},
   "source": [
    "Once again, we can inspect the rules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "bfe57a8c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best score and rule:\n",
      "  score=0.85 rule=AllEqual(And(~230, ~21), 232)\n",
      "\n",
      "Best score and rule for each start:\n",
      "  score=0.82 rule=AllEqual(~28, ~205, ~265, ~249, ~206)\n",
      "  score=0.74 rule=AllEqual(72, ~70)\n",
      "  score=0.76 rule=AllEqual(~261, 260)\n",
      "  score=0.85 rule=AllEqual(And(~230, ~21), 232)\n",
      "  score=0.81 rule=AllEqual(~239, ~225, ~178, ~88, ~7)\n",
      "  score=0.78 rule=AllEqual(~62, AllEqual(270, ~61))\n",
      "  score=0.81 rule=AllEqual(And(~230, ~262, ~263), 212)\n",
      "  score=0.83 rule=AllEqual(109, 66, 235, 78, 219)\n",
      "  score=0.80 rule=AllEqual(~274, ~179, ~299, ~138, ~207)\n",
      "  score=0.76 rule=AllEqual(AllEqual(79, ~136), 76, 108)\n"
     ]
    }
   ],
   "source": [
    "print(\"Best score and rule:\")\n",
    "print(f\"  score={rule_classifier.best_score_:.2f} rule={rule_classifier.best_rule_}\\n\")\n",
    "\n",
    "print(\"Best score and rule for each start:\")\n",
    "for score, rule in zip(rule_classifier.scores_, rule_classifier.rules_):\n",
    "    print(f\"  {score=:.2f} {rule=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbca2b87",
   "metadata": {},
   "source": [
    "We see that the only operators appearing in these rules are `And` and our new operator `AllEqual`, as expected. In this example we added an unparameterized operator. Adding parameterized operators follows the same pattern, provided that the operators have one integer parameter $k$ for which $1 \\leq k \\leq \\tilde{m}$ where $\\tilde{m}$ is the number of subrules of the operator. The value of $k$ is set on the fly by the rule classifier so that it is always in that range. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1794e1d",
   "metadata": {},
   "source": [
    "## Optimizing part of a rule\n",
    "\n",
    "We can also find the best subrule, given a pre-determined base rule. For example, we might already know that two features `f1` and `f2` are well correlated with the target (labels), and we'd like them to definitely be included. Then, we can formulate the problem as optimizing over `Or(f1, f2, *)`, where `*` is the part to be optimized over, which we refer to as the \"wildcard node\". We'll show how to do this step by step below.\n",
    "\n",
    "First, let us construct a rule containing two chosen features. We use two literals that appear in the best solution found at the top of this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a6b98194",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Or(230, 50)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.metrics import balanced_accuracy_score\n",
    "\n",
    "rule = Operator.Or([Literal(230), Literal(50)])\n",
    "rule"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7cefbca",
   "metadata": {},
   "source": [
    "We can check the score that this rule obtains on the dataset directly by using the rule's `evaluate()` method, followed by an application of the metric (in this case balanced accuracy):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "aee3c838",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "score=0.74\n"
     ]
    }
   ],
   "source": [
    "y_predict = rule.evaluate(X_binarized.values)\n",
    "score = balanced_accuracy_score(y, y_predict)\n",
    "print(f\"{score=:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a70f8bc",
   "metadata": {},
   "source": [
    "Now we'd like to add the additional wildcard node to this rule and optimize over that part of the rule. This node will be replaced by the result of the optimization we will run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d9b2af4d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rule before Or(230, 50, *)\n"
     ]
    }
   ],
   "source": [
    "from boolxai import Wildcard\n",
    "\n",
    "base_rule = Operator.Or([Literal(230), Literal(50), Wildcard()])\n",
    "\n",
    "print(\"Rule before\", base_rule)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a318c0ef",
   "metadata": {},
   "source": [
    "We instantiate a new `RuleClassifier` as done previously, but this time we pass in the `base_rule` containing the wildcard node:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9f866430",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rule after Or(230, 50, 222, 291)\n"
     ]
    }
   ],
   "source": [
    "rule_classifier = BoolXAI.RuleClassifier(random_state=43, base_rule=base_rule)\n",
    "rule_classifier.fit(X_binarized, y)\n",
    "print(\"Rule after\", rule_classifier.best_rule_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b2ca59e",
   "metadata": {},
   "source": [
    "Now we can evaluate the new rule, with the optimized subtree:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "f941fa9b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "score=0.90\n"
     ]
    }
   ],
   "source": [
    "y_predict = rule_classifier.best_rule_.evaluate(X_binarized.values)\n",
    "score = balanced_accuracy_score(y, y_predict)\n",
    "print(f\"{score=:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bc64be1",
   "metadata": {},
   "source": [
    "Or we can just obtain the score from the classifier, and of course they match:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ddc25e5d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best score and rule:\n",
      "  score=0.90 rule=Or(230, 50, 222, 291)\n"
     ]
    }
   ],
   "source": [
    "print(\"Best score and rule:\")\n",
    "print(f\"  score={rule_classifier.best_score_:.2f} rule={rule_classifier.best_rule_}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a7119e7",
   "metadata": {},
   "source": [
    "As before, we can also print the rule and score found in each start. In this case, we can verify that all the rules contain our base rule: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "b01d2698",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best score and rule for each start:\n",
      "  score=0.87 rule=Or(230, 50, 1, 232)\n",
      "  score=0.83 rule=Or(230, 50, Choose1(197, 272))\n",
      "  score=0.85 rule=Or(230, 50, 222)\n",
      "  score=0.85 rule=Or(230, 50, 222)\n",
      "  score=0.85 rule=Or(230, 50, 222)\n",
      "  score=0.85 rule=Or(230, 50, 222)\n",
      "  score=0.86 rule=Or(230, 50, AtMost1(~142, ~272))\n",
      "  score=0.88 rule=Or(230, 50, 171, 62)\n",
      "  score=0.85 rule=Or(230, 50, 222)\n",
      "  score=0.90 rule=Or(230, 50, 222, 291)\n"
     ]
    }
   ],
   "source": [
    "print(\"Best score and rule for each start:\")\n",
    "for score, rule in zip(rule_classifier.scores_, rule_classifier.rules_):\n",
    "    print(f\"  {score=:.2f} {rule=}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "boolxai_test",
   "language": "python",
   "name": "boolxai_test"
  },
  "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}