tutorials/clustering_euclidian.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustering with Logic Tensor Networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a basic example in which we use LTN for unsupervised clustering in a 2d space of real numbers. We define a theory that encodes the following facts\n",
    "* every point in the sample set should be assigned to a cluster\n",
    "* clusters should be to be disjoint\n",
    "* if two points are close (in terms of euclidian distance), they should belong to the same cluster\n",
    "\n",
    "We start generating normally distributed samples from $\\mathrm{nr\\_of\\_clusters}$ clusters. We use this to setup the theory and then optimize a loss function defined on the theory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging; logging.basicConfig(level=logging.INFO)\n",
    "import tensorflow as tf\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# adding parent directory to sys.path (ltnw is in parent dir)\n",
    "import sys, os\n",
    "sys.path.insert(0,os.path.normpath(os.path.join(os.path.abspath(''),os.path.pardir)))\n",
    "import logictensornetworks_wrapper as ltnw\n",
    "ltnw.ltn.set_universal_aggreg(\"min\")\n",
    "ltnw.ltn.set_existential_aggregator(\"max\")\n",
    "ltnw.ltn.set_tnorm(\"prod\")\n",
    "ltnw.ltn.LAYERS = 4\n",
    "import logictensornetworks_library as ltnl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate $\\mathrm{nr\\_of\\_points\\_x\\_cluster}$ samples for $\\mathrm{nr\\_of\\_clusters}$ clusters. Every sample is drawn from a normal distribution that is the underlying generator for each cluster. Samples from $\\mathbb{R}$^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "nr_of_clusters = 2\n",
    "nr_of_points_x_cluster = 10\n",
    "cov = np.array([[.01,0],[0,.01]])\n",
    "\n",
    "clusters=[]\n",
    "for i in range(nr_of_clusters):\n",
    "    mean = np.random.uniform([-1,-1],[1,1],2).astype(np.float32)    \n",
    "    clusters.append(np.random.multivariate_normal(mean=mean,cov=cov,size=nr_of_points_x_cluster).astype(np.float32) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every cluster is represented by a unary predicate $C_0,..C_{\\mathrm{nr\\_of\\_clusters}-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C_0\n",
      "C_1\n"
     ]
    }
   ],
   "source": [
    "[ltnw.predicate(\"C_\"+str(i),2) for i in range(nr_of_clusters)]\n",
    "print(\"\\n\".join(ltnw.PREDICATES.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the theory using a number of variables. First, we define two variables $x,y$ that encompass the complete sample set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data  = np.concatenate(clusters)\n",
    "ltnw.variable(\"?x\",data);\n",
    "ltnw.variable(\"?y\",data);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a function close takes an exponential envelope of the euclidian distance of two vectors. It is something like a euclidian similarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ltnw.predicate(\"close\",2,ltnl.equal_euclidian);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can define the theory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exists ?x: C_0(?x)\n",
      "exists ?x: C_1(?x)\n",
      "forall ?x,?y: (C_0(?x) & close(?x,?y)) -> C_0(?y)\n",
      "forall ?x,?y: (C_0(?x) & ~close(?x,?y)) -> (C_1(?y))\n",
      "forall ?x,?y: (C_1(?x) & close(?x,?y)) -> C_1(?y)\n",
      "forall ?x,?y: (C_1(?x) & ~close(?x,?y)) -> (C_0(?y))\n",
      "forall ?x: C_0(?x)|C_1(?x)\n",
      "forall ?x: ~(C_0(?x) & C_1(?x))\n"
     ]
    }
   ],
   "source": [
    "ltnw.axiom(\"forall ?x: %s\" % \"|\".join([\"C_%s(?x)\" % i for i in range(nr_of_clusters)]))\n",
    "for i in range(nr_of_clusters):\n",
    "    ltnw.axiom(\"exists ?x: C_%s(?x)\" % i)\n",
    "    ltnw.axiom(\"forall ?x,?y: (C_%s(?x) & close(?x,?y)) -> C_%s(?y)\" % (i,i))\n",
    "    ltnw.axiom(\"forall ?x,?y: (C_%s(?x) & ~close(?x,?y)) -> (%s)\" % (i,\"|\".join([\"C_%s(?y)\" % j for j in range(nr_of_clusters) if i!=j])))\n",
    "    \n",
    "    for j in range(i+1,nr_of_clusters):\n",
    "        ltnw.axiom(\"forall ?x: ~(C_%s(?x) & C_%s(?x))\" % (i,j))\n",
    "print(\"\\n\".join(sorted(ltnw.AXIOMS.keys())))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the knowledgebase and optimize. Notice that $\\mathrm{close}$ is not a very crisp predicate. Therefore reaching high sat levels is basically impossible. Hence, we never really reach sat levels near $1.0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:logictensornetworks_wrapper:Initializing knowledgebase\n",
      "INFO:logictensornetworks_wrapper:Initializing Tensorflow session\n",
      "INFO:logictensornetworks_wrapper:Initializing optimizer\n",
      "INFO:logictensornetworks_wrapper:Assembling feed dict\n",
      "INFO:logictensornetworks_wrapper:INITIALIZED with sat level = 0.6686284\n",
      "INFO:logictensornetworks_wrapper:TRAINING 0 sat level -----> 0.6686284\n",
      "INFO:logictensornetworks_wrapper:TRAINING 100 sat level -----> nan\n",
      "INFO:logictensornetworks_wrapper:TRAINING 200 sat level -----> nan\n",
      "INFO:logictensornetworks_wrapper:TRAINING 300 sat level -----> nan\n",
      "INFO:logictensornetworks_wrapper:TRAINING 400 sat level -----> nan\n",
      "INFO:logictensornetworks_wrapper:TRAINING 500 sat level -----> nan\n",
      "INFO:logictensornetworks_wrapper:TRAINING 600 sat level -----> nan\n",
      "INFO:logictensornetworks_wrapper:TRAINING 700 sat level -----> nan\n",
      "INFO:logictensornetworks_wrapper:TRAINING 800 sat level -----> nan\n",
      "INFO:logictensornetworks_wrapper:TRAINING 900 sat level -----> nan\n",
      "INFO:logictensornetworks_wrapper:TRAINING finished after 999 epochs with sat level nan\n"
     ]
    }
   ],
   "source": [
    "ltnw.initialize_knowledgebase(optimizer=tf.train.RMSPropOptimizer(learning_rate=0.01,decay=.9),\n",
    "                              initial_sat_level_threshold=.5)\n",
    "\n",
    "ltnw.train(max_epochs=1000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh0AAAF1CAYAAABML1hNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu4XWV57/3vj3CUsyRaIYSABksESuxS9HUrKEECbgNVigGpwdKm+IJ2S3ULLxYpbhRhVypbPKSKeChixEPTvmCKHLRVolkxEEi4IiEIJNCSQIIoAZK17v3HGDOMTOZhrHkcc6zf57rGteY4znvOtdZc93rG89yPIgIzMzOzbtuh3wGYmZnZ+OCkw8zMzHrCSYeZmZn1hJMOMzMz6wknHWZmZtYTTjrMzMysJ5x0WE9IukTSt3rwPCHpVd1+HjPbnqRvSzolx3FHSvp5L2Ky4nHSYT0naWqaHOzY5nXukPQXnYrLbLyTdIakYUm/k/SYpJsl/bd034cl/aekpyRdK2mXzHlHAn8E/HOz54iI5cAmSe/s2guxwnLSMU61+we/24oen1nZSDof+AfgU8DLgSnAF4CTJZ0AXAAcB0wFDgH+LnP6XwH/FPmrTf5Teo6NM046SkbSayUtk/S0pO9K+o6k/yXpWElrJX1M0n8CX0uP/0tJqyU9KWmhpP3T7S9qjci2LEg6S9J/SPrfkjZKelDSiZljD5b0kzSOW4CJmTB/mn7dlP5H9cb0ej+TdJWkJ4FLqm/JZGOSdBnwZuDz6TU+n7n+TEn3p3FdI0mdfZfNykXS3sClwLkR8f2I+H1EbImIf4mIjwJzga9GxIqI2Ah8Ejgrc4kTgZ9krvdFSTdm1j8j6dbM7+IdwHHZ1hIbH5x0lIiknYEfANcBLwW+DfxJ5pA/SLcfBMyT9Dbg08BpwCuAh4AbxvCURwOrSBKKK4CvZj5UrgeWpvs+SfKhVfGW9Os+EbFHRNyZud4a4GXAZY2eOCIuAv4dOC+9xnmZ3f8deB1Jc+9pwAljeE1m49EbgV1JPj9qeQ1wd2b9buDlkvaTtDtwMMlnQcXfAEem/0y8GTgbmFtpCYmIdcAW4NWdfRlWdG7CLpc3kHxPr05/ub8v6ZeZ/aPAJyLiOQBJ7wWujYhfpesXAhslTc35fA9FxD+m536dpCn25Wny8zpgZvpcP5X0Lzmu92hE/J/08dY2Giguj4hNJC0ptwNHAT9q9WJm48B+wIaI2Fpn/x7AU5n1yuM9SZIVgKcrOyPiGUlnkvzePQ18MCLWVl3zaWCfdgO3weKWjnLZH1hXdV/1kczj9RHxbNXxD1VWIuJ3wBPAATmf7z8z5z6TPtwjve7GiPh95tiHaO6R5oeMLS7gmTQmM6vvCWBig75UvwP2yqxXHj8NbEof75k9ISJ+SdJyKWBBjWvumTnXxgknHeXyGHBAVR+GAzOPqzt5PUpyqwWAtJl0P2AdUEkYXpI5/g/GEMe+6fUqpjSIo9723zd5fk+RbNYZdwLPAvWGvK4guV1Z8UfAf0XEE+k/Fw8Ah2ZPkHQusAvJ58z/rNq3P7Az29+SsXHASUe53AmMAOelnS1PBl7f4PjrgfdLOirt0PUp4BcR8ZuIWE+SfJwpaYKkPwdemSeIiHgIGAb+TtLO6ZC77PC49SS3eg5pcqm7gLdImpJ2dLuwav9/5biGmTUREU8BFwPXSDpF0ksk7STpRElXAN8AzpY0XdK+wMdJ+o5V3AQcU1mRdCjwv4AzgT8D/qekozLHHwvcVrnVa+OHk44SiYjngXeRdNraRPIL/69AzV/siLgV+FvgeyStE68E5mQO+UvgoyRNr68BxlLQ5wySjqFPAp8g+dCqPO8zJB1FfyZpk6Q31InvFuA7wHKSTqn/WnXI54BT01EqV48hNjOrEhGfBc4nSSjWk9zuPA/4YUT8iKSz+O0kt0ofIvm9rpgPvFeJHYFvAZ+JiLsj4n7g/wO+mRmt8l7gSz14WVYwyj+s2gaRpF8AX4qIr/U7FjMrL0nXAwsi4odNjjsCmB8Rb+xNZFYkTjpKRtIxJPdJN/DCfxOHRMRjfQ3MzMzGPd9eKZ9Xk4yhf4pkrPypTjjMzGys0nL3j0u6t85+Sbo6LTC5XNJrm17TLR1mZmZWTdJbSIZLfyMiDq+x/yTgg8BJJH34PhcRRze6pls6zMzM7EUi4qckgwHqOZkkIYmIWAzsI+kVja7ppMPMzMxacQDbF3VcS5PikoUrgz5x4sSYOnVqv8Mws9TSpUs3RMSkfscxFv4csbIay+/jCW/dPZ54cqT+tZY/t4KkKFzF/IiYP4Zwas1V0bDPRuGSjqlTpzI8PNzvMMwsJSlPCftC8eeIldVYfh83PLmVn/+ofsPDrvs/+GxEDLURzlq2r3o9maQCbV2+vWJmZlZCAWxlpO7SAQuB96WjWN4APNVstGThWjrMzMysfUEw0sYIVUnfJilZP1HSWpIqtDsBRMSXSMrfnwSsJplc8/3Nrumkw8x6QtIsktL1E4CvRMTlVfvPAc4lmT/od8C8iFgpaSpwHy9MDrY4Is7pVdxmgyqALYy2fn7E6U32B8nvbG7jN+lYvgBuvRSeWgt7T4bjLoYjT+t3VGalJGkCcA1wPMl94CWSFkbEysxh16f/PSFpNvBZYFa674GIyE4YZmZNBLAlWk86umF8Jh3LF8C/fAi2bE7Wn3okWQcnHmbd8XpgdUSsAZB0A8kY/21JR0T8NnP87jTpBW9mzRUr5RivHUlvvfSFhKNiy+Zku5l1Q67x/JLOlfQAyYymH8rsOljSMkk/kfTmWk8gaZ6kYUnD69ev72TsZgMpIni+wdIP4zPpeGrt2LabWbtyjeePiGsi4pXAx0imWAd4DJgSETNIpl6/XtJeNc6dHxFDETE0adJAlRUx64pAbGmw9MP4TDr2njy27WbWrrGO578BOAUgIp6LiCfSx0uBB4BDuxSnWWkkfTpUd+mH8Zl0HHcx7LTb9tt22i3ZbmbdsASYJulgSTsDc0jG+G8jaVpm9R3A/en2SWlHVCQdAkwD1vQkarMBFsAIqrv0w/jsSFrpLOrRK2Y9ERFbJZ0HLCIZMnttRKyQdCkwHBELgfMkzQS2ABuBuenpbwEulbSVZDjtORHRaBIqM6PS0lGstoXxmXRAkmA4yTDrmYi4iaSYUHbbxZnHf13nvO8B3+tudGblE4gtMaHfYWxn/CYdZmZmJRaI5510mJmZWbcFMFqwrptOOszMzEoowi0dZmZm1gPJ3CtOOszMzKzrxIhHr/SZJ3ozM7NxIBkyW6yWjrZSIEmzJK2StFrSBTX2T5F0ezpnwnJJJ7XzfG2rTPT21CNAvDDR2/IF2x9z1eFwyT7J1+w+MzOzAZEMmd2x7tIPLScdmamqTwSmA6dLml512MeBBemcCXOAL7T6fB3RbKK3PEmJmZnZABhNh8zWW/qhnZaObVNVR8TzJHMlnFx1TACViZn2pvFcC93XbKI3zz5rZmYlMsoOdZd+aKd9pdZU1UdXHXMJ8G+SPgjsDsysdSFJ84B5AFOmTGkjpCb2npy2YtTYDp591szMSiOieBVJ20l18kxVfTpwXURMBk4CvinpRc/Zsympm0305tlnzcysJAJ4Pnasu/RDO0lHnqmqzwYWAETEncCuwMQ2nrM9R54G77wa9j4QUPL1nVe/MHrFs8+amVlJVOZeqbf0QzupzrapqoF1JB1Fz6g65mHgOOA6SYeRJB3r23jO9jWa6M2zz5qZWUkEMFqWOh05p6r+G+AfJX2Y5PWfFRHVt2CKxbPPmplZCbQ7y6ykWcDnSP7GfyUiLq/aPwX4OrBPeswF6WzSdbV1UyfHVNUrgTe18xxmZmY2du10JM2UxTiepDvFEkkL07/rFZWyGF9MS2bcBExtdN1itbuYmZlZx4zEDnWXJrpSFmP8lUE3MzMbB3LcXpkoaTizPj8i5qePO1YWI8tJh5mZWQnlmHtlQ0QM1dk3lrIYfy/pjSRlMQ6PiNF6T+ikw8zMrIQCsbX1jqR5y2LMgqQshqRKWYzH613UfTrMzMxKKAJGQnWXJraVxZC0M0lZjIVVx1TKYpC3LMZgtnR4enozM7OGArF1tLWWjm6VxRi8pKMyE2xlYrbKTLBQzMTDCZKZmfVB0qej9Rsa3SiLMXi3VwZpJthKgvTUI0C8kCAtX9DvyMzMrPSSPh31ln4YvKRjkGaCHaQEyazLJM2StErSakkX1Nh/jqR7JN0l6T/SYkOVfRem562SdEJvIzcbTG326eiKwUs6Bmkm2EFKkMy6KFPd8ERgOnB6NqlIXR8RR0TEUcAVwGfTc6eTdGJ7DUlP+S+k1zOzBip9Ouot/TB4SccgzQQ7SAmSWXc1rW4YEb/NrO7OCzUBTgZuiIjnIuJBYHV6PTNrIICtsUPdpR8GryNpkWeCre40Ou3tcPf1299iKWqCZNZdeaobIulc4HxgZ+BtmXMXV517QI1z5wHzAKZMmdKRoM0GXWlmme2rIs4EW2tUzd3Xwx+dAff/W/ESJLPeylPdkIi4BrhG0hkkk0nNHcO584H5AENDQ8WezdqsByLUtxaNegYz6Siiep1G7/83+PC9/YnJrDjyVDfMugH4Yovnmhnp7ZVRJx39te0WyCOgCRAjsPeB27dAtFJbw51GzRrZVt0QWEfSMfSM7AGSpkXE/enqO4DK44XA9ZI+C+wPTAN+2ZOozQZY0pHUSUf/VN8CiZHka7bAGLRWfGzvyWk9jhrbzca5nNUNz5M0E9gCbCS5tUJ63AJgJbAVODei8strZo2M1rw72T/jK+modQukIls/o15tjUZJx3EXb5+sgDuNmmXkqG741w3OvQy4rHvRmZVPhG+v9FezWx2N9jc7t8ijaszMbNzx7ZV+q3cLJLsfWr9NUsRRNWZmNm6NFGz0SrGi6bZahcUqKrdCBqn4mJmZWR0RMBqqu/TD+Grp2O4WSIPRK9uO8W0SMzMbVGKkTLdXJM0CPkfSG/0rEXF5jWNOAy4hGTJ8d0ScUX1MT+W5BeLbJGZmNuACypN0ZCZwOp6keM8SSQsjYmXmmGnAhcCbImKjpJe1G7CZmZnlEMktliJpJwVqOoET8JfANRGxESAiHm/j+czMzCynIOlIWm/ph3aetdYETtWTMB0KHCrpZ5IWp7djXkTSPEnDkobXr1/fRkgtWL4ArjocLtkn+bp8QW+f38zMrCvEyGj9pR/a6dORZxKmHUlKFh9LMl/Cv0s6PCI2bXdSvyZqqjVJW57qo2ZmZgUXAaMF69PRTjR5JmFaC/xzRGyJiAeBVSRJSDHUm6StUpnUzMxsgBVtyGw7Sce2CZwk7UwygdPCqmN+CLwVQNJEktsta9p4zs7yJG1mZlZio6Oqu/RDy0lHRGwFKhM43QcsqEzgJGl2etgi4AlJK4HbgY9GxBPtBt0x9aqMepI2MzMbcEH9Vo48LR2SZklaJWm1pAvqHHOapJWSVki6vtk126rTkWMCpwDOT5fi8SRtZmZWVgHRYotGt8piFKuHSa8deRq88+qkIilKvr7zanciNTOzUohQ3aWJrpTFGF9l0Gtx9VEzMyuhgGZ9NyZKGs6sz09Hk0LtshhHV51/KICkn5FUJr8kIn7U6AmddJiZmZVR89srGyJiqM6+jpXFyBrft1fMzMzKLBosjXWlLIaTDjMzs1ISMVp/aaIrZTF8e8XMzKyM2hi9EhFbJVXKYkwArq2UxQCGI2Jhuu/taVmMEXKUxXDSYWZmVlZtVB7tRlkMJx1mZmZlVbCp7Z10mJmZlVEbt1e6xUmHmZlZWRUs6fDoFTPriWbzOEg6P53DYbmkWyUdlNk3IumudKnuQW9mtQRotP7SD27pMLOuyzOPA7AMGIqIZyR9ALgCeE+6b3NEHNXToM0GntrqSNoNbukws15oOo9DRNweEc+kq4tJihGZWTtGGyx94KTDzHqh1jwOBzQ4/mzg5sz6rpKGJS2WdEqtEyTNS48ZXr9+ffsRmw26IOnTUW/pA99eMbNeyDOPQ3KgdCYwBByT2TwlIh6VdAhwm6R7IuKB7S6WTFQ1H2BoaKhgAwXN+kMF+01wS4eZ9UKeeRyQNBO4CJgdEc9VtkfEo+nXNcAdwIxuBmtWGq3PvdIVTjrMrBeazuMgaQbwZZKE4/HM9n0l7ZI+ngi8Cch2QDWzOjSquks/+PaKmXVdznkcrgT2AL4rCeDhiJgNHAZ8WdIoyT9Kl1eNejGzWoK+dRitx0mHmfVEjnkcZtY57+fAEd2Nzqycitanw0mHmZlZWRWspaOtPh3NKgxmjjtVUkgaauf5zMzMLB9F8fp0tJx0ZCoMnghMB06XNL3GcXsCHwJ+0epzmZmZ2diVqQz6tgqDAJIqFQarO3h9kqSc8UfaeC6z0vvhsnVcuWgVj27azP777MZHT3g1p8xoVD/LzKyJgvXpaOf2StMKg+kQuAMj4l/beB6z0vvhsnVc+P17WLdpMwGs27SZC79/Dz9ctq7foZnZoCrghG/tJB0NKwxK2gG4Cvibphdy+WIb565ctIrNW0a227Z5ywhXLlrVp4jMrAzKlHQ0qzC4J3A4cIek3wBvABbW6kwaEfMjYigihiZNmtRGSGaD6dFNm8e03cxsELWTdDSsMBgRT0XExIiYGhFTSWaNnB0Rw21FbFZC+++z25i2m5k1VabbKxGxFahUGLwPWFCpMChpdqcCNBsPPnrCq9ltpwnbbdttpwl89IRX9ykiMyuFNqa270ZZjLaKgzWrMFi1/dh2nsusjLIjVvbebSd23WkHNj2zxaNXzKxtovUWjUxZjONJulMskbSwegqCsZbFcEVSsz6pjFipdCDdtHkLu+00gavec5STDTNrX7RVBr0rZTE8y6xZn3jEipl1Xeu3V7pSFsMtHWZ94hErZtZtTW6vTJSUHdwxPyLmV06tcXytshhnjSUeJx1mfbL/PruxrkaC4RErZtYR0TTp2BAR9Tp/jqUsBsAfkJTFaDhK1bdXzPrEI1bMrOuiwdJYV8piuKXDrE8qnUU934qZdUuro1ciYqukSlmMCcC1lbIYwHBELGx8hdqcdJj10SkzDnCSYWZdofZGr3SlLIaTDjMzs7Iq2CyzTjrMzMxKql/lzutx0mFmZlZGzUev9JyTDjMzs5Jy0mFmZma94T4dZuNTdnI3D481s64r4O0VFwcz64HK5G7rNm0mgHWbNnPh9+/hh8vW9Tu0nmk2Tbak8yWtlLRc0q2SDsrsmyvp/nSZ29vIzQZTZZbZeks/OOkw64HxPrlbZprsE4HpwOmSplcdtgwYiogjgRtJZq5E0kuBTwBHk8x8+QlJ+/YqdrNBptGou/SDkw6zHvDkbi9Mkx0RzwOVabK3iYjbI+KZdHUxyVwPACcAt0TEkxGxEbgFmNWjuM0GV6MS6H3q6+Gkw6wH6k3iNo4md2s6TXaVs4Gbx3KupHmShiUNr1+/vs1wzcrBt1fMxiFP7tZ4muztDpTOBIaAK8dybkTMj4ihiBiaNGlSy4GalYmTDrNx6JQZB/Dpdx3BAfvshoAD9tmNT7/riPE0eqXZNNkASJoJXEQyW+VzYznXzKrEC/Ov1Fr6wUNmzXqk35O79XnI7rZpsoF1JNNkn5E9QNIM4MvArIh4PLNrEfCpTOfRtwMXdj9ks8FWGb1SJE46zLqoW3/ox3rdypDdygiaypBdoCeJR85psq8E9gC+Kwng4YiYHRFPSvokSeICcGlEPNn1oM1KoF+jVOppK+mQNAv4HMmHyFci4vKq/ecDfwFsBdYDfx4RD7XznGaDolt/6Fu5bqMhu71q7Wg2TXZEzGxw7rXAtd2LzqyEAjTS/LBearlPRzvj7s3Gg27V5mjluh6yazZOlWjIbDvj7s1Kr1t/6Fu5rofsmo1DUa7iYO2Mu9+Ox9dbGXXrD30r1/WQXbPxqUxDZtsZd7/9SR5fbyXUrT/0rVzXQ3bNxh9F/VaOfrV0tNORdKzj7o/JjLs3K73KH/ROj15p9br9HrJrZr3Xr3oc9bSTdLQz7t5sXOjUH/paQ2R/dsHbOhChmZVZ0ep0tHx7JSK2ApVx9/cBCyrj7iXNTg/Ljru/S9LCtiM2G2cqQ2TXbdpM8MIQ2R8uW9fv0MysyAIYifpLE5JmSVolabWkC2rsP1/SSknLJd0q6aBm12yrTkc74+7NLJ8i1Ngws8HU6u2VTFmM40m6UyyRtDAiVmYOq5TFeEbSB0jKYryn0XU994pZwdUbCrtu02a3dphZQ210JO1KWQwnHWYF12gorG+zmFk9al6nY2KlXEW6zMuc3rGyGFlOOswKrtYQ2YpOVDg1s/LSSNRdgA2VchXpMj97ao3LtVQWI8sTvpkVXKXfxv/4zl0197uUuZnV1F65866UxXBLh9kAOGXGARzgUuZmNiZtFQfbVhZD0s4kZTG2G4GaKYsxO29ZDCcdZgPCpczNbEyi6e2V+qd2qSyGb6+YDYhuVTg1sxJro9x5N8piOOkwGyAuZW5mY6EoVh10Jx1mZmZlVKlIWiBOOszMzEpIBBot1uQrTjrMzMzKyrdXzMzMrOvS0StF4qTDrA9qTVXvDqJm1lkBvr1iNr5VpqqvzBxbmaoecOJhZp1TwI6kLg5m1mONpqovM0mzJK2StFrSBTX2v0XSryRtlXRq1b6RtPhQrgJEZpZQRN2lH9zSYdZj9eZKKfMcKpImANcAx5PM6bBE0sKIWJk57GHgLOAjNS6xOSKO6nqgZmUSwEixbq+4pcOsx+rNlVLyOVReD6yOiDUR8TxwA3By9oCI+E1ELAeK9SlpNrDSPh31lj5w0mHWY+N0DpUDgEcy62vTbXntKmlY0mJJp3Q2NLOSqrR01Fv6wLdXzNo01pEo43QOFdXYNpabylMi4lFJhwC3SbonIh7Y7gmkecA8gClTprQeqVlpBESxGg6ddJi1odWRKONwDpW1wIGZ9cnAo3lPjohH069rJN0BzAAeqDpmPjAfYGhoqFhd9s36oWx9OnL0Rt9F0nfS/b+QNLWd5zMrmvE6EqUFS4Bpkg6WtDMwB8g1CkXSvpJ2SR9PBN4ErGx8lpkB5enTkemNfiIwHThd0vSqw84GNkbEq4CrgM+0+nxmRTQeR6K0IiK2AucBi4D7gAURsULSpZJmA0h6naS1wJ8CX5a0Ij39MGBY0t3A7cDlVaNezKymSMqg11v6oJ3bK9t6owNIqvRGz34YnAxckj6+Efi8JEUUrBi8WYv232c31tVIMEo+EqUlEXETcFPVtoszj5eQ3HapPu/nwBFdD9CsbAIYGWl6WC+1c3slT2/0bcek/+k8BezXxnOaFco4HYliZgMhSjV6JU9v9Fw91t3r3AbVOB2JYmaDICAK1tLRTtKRpzd65Zi1knYE9gaerL6Qe53bIBuHI1HMbFAUrDdDO7dX8vRGXwjMTR+fCtzm/hxmZmY9EJH06ai39EHLSUee3ujAV4H9JK0GzgdeNKzWzMzMuiGIkZG6SzPdKIvRVnGwHL3RnyUZ/mZmZma91MbolZyTNG4riyFpDklZjPc0uq7nXjEzMyuhAGI06i5NNJ2kMV3/evr4RuA4SbUGkGzjpMPMzKyMoq3bK10pi1G4uVeWLl26QdJDXbj0RGBDF67bC469fwY5/k7FflAHrtFT6efI7ynG964oP0OO48WKEstY4sj9+/g0Gxf9eHTBxAaH7CppOLM+Px1NCh0si5FVuKQjIiZ147qShiNiqBvX7jbH3j+DHP8gx96uiJhUlNfvOIoZBxQnlm7FERGz2ji9Y2Uxsnx7xczMzKp1pSxG4Vo6zMzMrL8iYqukSlmMCcC1lbIYwHBELCQpi/HNtCzGkySJSUPjKemY3/yQwnLs/TPI8Q9y7J1QlNfvOLZXlDigOLEUJY7tdKMshlwg1MzMzHrBfTrMzMysJ0qbdEj6U0krJI1KqtsruFmZ136Q9FJJt0i6P/26b53jRiTdlS7VHXx6qhvlcnslR+xnSVqfea//oh9x1iLpWkmPS7q3zn5Jujp9bcslvbbXMXZCju/RVZnvz68lbUq3vzWz/S5Jz0o6Jd13naQHM/uO6kAcUyTdLmlZ+n6flNl3YXreKkkn5L1mp2ORdLykpZLuSb++LXPOHek1K+/Jy7oYx1RJmzPP9aXMOX+cxrc6/fltWHCqzTjeW/UzMlr5WejS+3GQpFvTGO6QNDmzb66Sz/37Jc3NbB/z+1FYEVHKBTgMeDVwBzBU55gJwAPAIcDOwN3A9ALEfgVwQfr4AuAzdY77Xb9jzfs+Av8v8KX08RzgO/2OewyxnwV8vt+x1on/LcBrgXvr7D8JuJlkPP0bgF/0O+ZufI+qjv8gSae36u0vJens9pJ0/Trg1A7/rMwHPpA+ng78JvP4bmAX4OD0OhNa/QxqM5YZwP7p48OBdZlz6n5ediGOqQ1+bn8JvDH9ub0ZOLFbcVQdcwSwpsvvx3eBuenjtwHfzPx8rkm/7ps+3reV96PIS2lbOiLivohY1eSwPGVe+yFbWvbrwCl9jCWPrpTL7ZGi/gzkEhE/pfG4+JOBb0RiMbCPpFf0JrqOGev36HTg2zW2nwrcHBHPdDGOAPZKH+/NC3UNTgZuiIjnIuJBYHV6vVZ//lqOJSKWRUQlrhUkBaJ2yfGcHY2jnvTnc6+IuDOSv7jfoPlnYKfiqPezk1eeOKYDt6aPb8/sPwG4JSKejIiNwC3ArBbfj8IqbdKRU54yr/3w8oh4DCD9Wq9Jb1dJw5IWV5qM+6Qr5XJ7JO/PwLvT5tAbJR1YY39RFfVnfCxyvwZJB5G0JNxWY/ccXvwH5bL0+3pVjj+8eeK4BDhT0lqSXv8fbHJuq9+fdmLJejewLCKey2z7Wnor4W9z/GPQbhwHp7c7fiLpzZlrrm1yzU7HUfEeXvwz0un3426S9x3gT4A9Je3X4NxW3o/CGuikQ9KPJd1bY8n7n+qYS7h2SgdiB5gSSRW7M4B/kPTKLoXbTFfK5fZInrj+BZgaEUcCP+aFFptBUNT3fSzG8hrmADdGxHYTS6T/LR5BUnOg4kLgD4HXkTRpf6wDcZwOXBcRk0lubX1T0g4Nzm31+9NOLMkFpNeQzAr6V5lz3hsRRwBvTpc39ozyAAAbIUlEQVQ/62Icj5F8hs0Azgeul7RXzmt2Mo7kAtLRwDMRke0f1Y334yPAMZKWAccA64CtDc4tw+/wNgNdpyMiZrZ5iTxlXruiUeyS/kvSKyLisfTD8vE616g0l66RdAfJvdoHuhFvE10pl9sjTWOPiCcyq/9I8kE9KPr2M95BY3kNc4Bza2w/DfhBRGypbKi0JgLPSfoayR+DduM4G5iVXv9OSbuSzKvR6NxWvj/txPJ42nnxB8D7ImLbZ0ZErEu/Pi3pepLbBd/oRhwR8TjwXLp9qaQHgEPTa07OnJ/nPWnr/Uj3v6glrBvvR/q5/S4ASXsA746Ip9IWmGOrzr2D1t6P4up3p5JuLzTuSLojSWedg3mh089rChDzlWzfkfSKGsfsC+ySPp4I3E+fOsHmeR9J/hBkO5Iu6Pf7PIbYX5F5/CfA4n7HXRXfVOp3yHsH23ck/WW/4+3G9yg97tXAb0jrD1XtWwy8tdb3NX1v/gG4vAM/KzcDZ6WPDyP54yDgNWzfkXQNSafDlj6D2oxln/T4d9e45sT08U4kfa/O6WIck4AJ6fZDSP7jf2m6viT9ea10nDypW3Gk6zuQ/HE/pAfvx0Rgh/TxZcCl6eOXAg+SfLbvmz5u6f0o8tL3ALr2wpI/DmtJMun/Ahal2/cHbsocdxLwa5IWgov6HXca034kHY3uT79WfvCGgK+kj/8f4J70h/oe4Ow2nu/TwP/IcdzLgftIk52qfS96H4FLgdnp411Jem2vJumJfUir8Xbh/W4W+6dJOtzdTdLx6w/7HXMm9m+TNFNvSX/ezwbOqXw4ph9S16Sv7R7G0BO/SEuz71G6fgk1EgeSpGxd5YM+s/229D25F/gWsEcHflamAz9Lf1buAt6eOfei9LxVZEYftPoZ1GoswMeB36fbKsvLgN2BpcDy9Of9c6RJQZfieHfm9+pXwDsz1xxKvy8PAJ+nRiLZ4e/NsVT9M9HF9+NUks/2XwNfIfN5Cvw5yWfkauD97bwfRV1ckbQkJJ1Bcl/0D4GnSX6pLgM2AX8P/DGwX0So6rxJ6bGviojNOZ7nC8B9EfF/OvsKzMys7Aa6I6klJJ1P0jz8KZLWiCnAF0iGYm0BFpD8B1zLWSQtP00TjtQ/sX3HMzMzs1zc0jHgJO1N0nT8/oj4boPjXgXcX6Ol4zaSQkrfStc/RnJr6r9FMsvgB4DzgD+OiGfTTqC/BQ6LiIe686rMzKyM3NIx+N5I0l/iBy2efwTJPeaKK4HngY9LmkbSenJmJLMJEkmNjdXAH7UcsZmZjUsDPWTWgKTT6YY0GWjFPiR9QACIiFFJ7yPp2PUekpEzy6rOeTo9z8zMLDe3dAy+J4CJ6W2PVmwE9sxuiIjfkIzSmEoy8qHaniQdVM3MzHJz0jH47gSepfVa/MtJCvJsk86++EaS4bpXVu3bEXgVybAzMzOz3Jx0DLiIeAq4GLhG0imSXiJpJ0knSrpCiV1JCtUgqXpyp5tISvGS7p8IfBX4C2Au8E5lpucmqcj3G3ciNTOzsfLolZKQ9F7gwySV9p4mKWpzGUnVvQerDn8oIqam500kqdMxLSI2S/o+8HhEnJPuP5EkCTkiIp6QdA2wKiKu7sHLMjOzEnHSYUj6FEmi8Q9NjnsZ8BNgRmU0i5mZWV5OOszMzKwn3KfDrIQkXSvpcUn31tkvSVdLWi1puaTX9jpGMxt/nHSYldN1pNN413EiMC1d5gFf7EFMZjbOOekwK6GI+CnwZINDTga+EYnFwD6SXtGb6MxsvHLSYTY+HQA8kllfm24zM+uawpVBnzhxYkydOrXfYZh13NKlSzdExKQ8x57w1t3jiSdH6l9r+XMrSIrCVcyPiPljCEc1trlXuZl1VeGSjqlTpzI8PNzvMMw6TlLugmobntzKz39Uv+Fh1/0ffDYihtoIZy1wYGZ9MklNFzOzrvHtFbMCCmArI3WXDlgIvC8dxfIG4KmIeKwTFzYzq6dwLR1mBkEw0kYNHUnfBo4lmQxwLfAJYCeAiPgSSfn7k4DVwDPA+9sM2cysqVxJh6RZwOeACcBXIuLyqv3nAOcCI8DvgHkRsVLSVOA+YFV66OJKeW0zqy+ALYy2fn7E6U32B8nvrJlZzzRNOiRNIJne/HiS+8BLJC2MiJWZw65P/3tC0mzgs7xQI+CBiDiqs2GblVsAW6L1pMPMrIjy9Ol4PbA6ItZExPPADSRj/LeJiN9mVnfHveDN2hIEWxosZmaDKE/SkWs8v6RzJT0AXAF8KLPrYEnLJP1E0ptrPYGkeZKGJQ2vX79+DOGblVTASIPFzGwQ5Uk6co3nj4hrIuKVwMeAj6ebHwOmRMQM4Hzgekl71Th3fkQMRcTQpEm5yhiYlVogtjRYzMwGUZ6kY6zj+W8ATgGIiOci4on08VLgAeDQ1kI1Gz+SPh2qu5iZDaI8SccSYJqkgyXtDMwhGeO/jaRpmdV3APen2yelHVGRdAjJ5FJrOhG4WZklSccOdRczs0HUdPRKRGyVdB6wiGTI7LURsULSpcBwRCwEzpM0E9gCbATmpqe/BbhU0laS4bTnRESjSajMjCTpGPFtFDMrmVx1OiLiJpJiQtltF2ce/3Wd874HfK+dAM3Go0BsiQn9DsPMrKNckdSsgALxvJMOMysZJx1mBRTAqKdGMrOScdJhVkARbukws/Jx0mFWQMncK046zKxcnHSYFVDSkdS/nmZWLv5UMyuoERcBM7OScdJhVkBu6TCzMvKnmlkBjXrIrJmVkJMOswKKwC0dZlY6/lQzKyQx6jLoZlYyTjrMCiiA593SYWYl4081swLy3CtmVkZOOswKKIBRT2FvZiXjTzWzAqq0dNRb8pA0S9IqSaslXVBj/xRJt0taJmm5pJM6/kLMzDLc0mFWQBHt3V6RNAG4BjgeWAsskbQwIlZmDvs4sCAivihpOnATMLX1qM3MGnPSYVZAAe326Xg9sDoi1gBIugE4GcgmHQHslT7eG3i0nSc0M2vGSYdZIYmRxn06JkoazqzPj4j5mfUDgEcy62uBo6uucQnwb5I+COwOzGw9XjOz5px0mBVQjpaODREx1GB/rSIfUbV+OnBdRPy9pDcC35R0eESMji1aM7N8nHSYFVAgtrZ3e2UtcGBmfTIvvn1yNjALICLulLQrMBF4vJ0nNjOrx6NXzAooAraM7lB3yWEJME3SwZJ2BuYAC6uOeRg4DkDSYcCuwPoOvgwzs+24pcOsgAK1VacjIrZKOg9YBEwAro2IFZIuBYYjYiHwN8A/Svowya2XsyKi+haMmVnHOOkwK6CkT0d7DZERcRPJMNjstoszj1cCb2rrSczMxsBJh1khtd2nw8yscHL9K5WjsuE5ku6RdJek/0gLDVX2XZiet0rSCZ0M3qysImAkVHcxMxtETVs6clY2vD4ivpQePxv4LDArTT7mAK8B9gd+LOnQiBjp8OswK5VAbB11S4eZlUuelo5tlQ0j4nmgUtlwm4j4bWZ1d16oB3AycENEPBcRDwKr0+uZWQMBbI0d6i5mZoMoT5+OPJUNkXQucD6wM/C2zLmLq849oMa584B5AFOmTMkTt1m5hVs6zKx88vzLlKeyIRFxTUS8EvgYyURSYzl3fkQMRcTQpEmTcoRkVm4BjKK6i5nZIMrT0pGnsmHWDcAXWzzXzEhvr+QrAmZmNjDyfKo1rWwoaVpm9R3A/enjhcAcSbtIOhiYBvyy/bDNyi3pSLpD3cXMbBA1benIWdnwPEkzgS3ARmBueu4KSQtIptPeCpzrkStmOQTuMGpmpZOrOFiOyoZ/3eDcy4DLWg3QbDwKYNT1OMysZFyR1KyAKrdXzMzKxEmHWUGN+PaKmZWMkw6zAorw7RUzKx8nHWaFJEZ8e8XMSsZJh1kBBTjpMLPScdJhVkTpLLNmZmXipMOsgAIIJx1mVjJOOswKSYyMOukws3Jx0mFWQBEw6j4dZlYyTjrMCsotHWZWNv5XyqygIlR3yUPSLEmrJK2WdEGdY06TtFLSCknXd/QFmJlVcUuHWQEFaqs4mKQJwDXA8cBaYImkhRGxMnPMNOBC4E0RsVHSy9oM28ysIbd0mBVRQIyq7pLD64HVEbEmIp4HbgBOrjrmL4FrImIjQEQ83tHXYGZWxUmHWUG1eXvlAOCRzPradFvWocChkn4mabGkWR0K3cysJt9eMSugAEYbt2hMlDScWZ8fEfMz67VOjqr1HYFpwLHAZODfJR0eEZvGHrGZWXNOOsyKKL290sCGiBhqsH8tcGBmfTLwaI1jFkfEFuBBSatIkpAlLURsZtaUb6+YFVL9/hw5+3QsAaZJOljSzsAcYGHVMT8E3gogaSLJ7ZY1HXwRZmbbcdJhVlTRYGl2asRW4DxgEXAfsCAiVki6VNLs9LBFwBOSVgK3Ax+NiCc6/TLMzCp8e8WsiJrfXml+iYibgJuqtl2ceRzA+eliZtZ1TjrMisoTvplZyTjpMCuiAFwG3cxKxkmHWUFFjr4bZmaDxEmHWVG5pcPMSibX6JVmE0dJOj+dNGq5pFslHZTZNyLprnSpHrJnZrUEaLT+YmY2iJq2dOSZOApYBgxFxDOSPgBcAbwn3bc5Io7qcNxmJSd3JDWz0snT0tF04qiIuD0inklXF5NUPzSzdow2WMzMBlCepCPPxFFZZwM3Z9Z3lTScTih1Sq0TJM1Ljxlev359jpDMSq4yeqXeYmY2gPJ0JM0zcVRyoHQmMAQck9k8JSIelXQIcJukeyLige0ulkxUNR9gaGjIffbNcN8NMyufPC0deSaOQtJM4CJgdkQ8V9keEY+mX9cAdwAz2ojXzMzMBlSepKPpxFGSZgBfJkk4Hs9s31fSLunjicCbgGwHVDOrQ6Oqu5iZDaKmt1ciYqukysRRE4BrKxNHAcMRsRC4EtgD+K4kgIcjYjZwGPBlSaMkCc7lVaNezKyWwB1Gzax0chUHyzFx1Mw65/0cOKKdAM3GK/fpMLOycUVSs6Jyl2ozKxknHWYFpMB9N8ysdJx0mBWUb6+YWdk46TArKt9eMbOScdJhVkThlg4zKx8nHWYF5aTDzMrGSYdZUfn2ipmVjJMOsyKKZASLmVmZ5CmDbmb90ObU9pJmSVolabWkCxocd6qkkDTUftBmZvU56TArIJH06ai3ND1fmgBcA5wITAdOlzS9xnF7Ah8CftHRF2BmVoOTDrMiivaSDuD1wOqIWBMRzwM3ACfXOO6TwBXAsx2L3cysDicdZkUVDRaYKGk4s8yrOvsA4JHM+tp02zbp7NAHRsS/ducFmJltzx1JzQqqSYvGhoho1AejVg31bV1TJe0AXAWc1UpsZmatcNJhVkTtFwdbCxyYWZ8MPJpZ3xM4HLhDEsAfAAslzY6I4bae2cysDt9eMSuqxrdXmlkCTJN0sKSdgTnAwm2XjngqIiZGxNSImAosBpxwmFlXOekwK6h2OpJGxFbgPGARcB+wICJWSLpU0uzuRm5mVptvr5gVkDpQHCwibgJuqtp2cZ1jj23v2czMmnPSYVZQnnvFzMrGSYdZUbkMupmVjJMOsyLy1PZmVkJOOswKykmHmZWNkw6zgnLSYWZlk2vIbLPZKiWdL2mlpOWSbpV0UGbfXEn3p8vcTgZvVlqNanS4r4eZDaimSUfO2SqXAUMRcSRwI8kEUkh6KfAJ4GiSCag+IWnfzoVvVk7tzjJrZlZEeVo6ms5WGRG3R8Qz6epikpLLACcAt0TEkxGxEbgFmNWZ0M3KTaNRdzEzG0R5ko6ms1VWORu4ucVzzQx8e8XMSilPR9KGs1Vud6B0JjAEHDOWc9NpuecBTJkyJUdIZuXn2yhmVjZ5WjqazVYJgKSZwEUkk0Y9N5ZzI2J+RAxFxNCkSZPyxm5Wau7TYWZlkyfpaDhbJYCkGcCXSRKOxzO7FgFvl7Rv2oH07ek2M2sknHSYWfk0vb0SEVslVWarnABcW5mtEhiOiIXAlcAewHclATwcEbMj4klJnyRJXAAujYgnu/JKzEpEtD/hm5lZ0eQqDtZstsqImNng3GuBa1sN0Gy88igVMysbVyQ1K6IAjfQ7CDOzznLSYVZQ7rthZmXjpMOsiAII314xs3Jx0mFWUG7pMLOycdJhVkAKlzs3s/Jx0mFWUB4ya2Zl46TDrKB8e8XMyiZPRVIz67UARqL+koOkWZJWSVot6YIa+8+XtFLSckm3Sjqo0y/DzCzLSYdZQbUztb2kCcA1wInAdOB0SdOrDlsGDEXEkcCNwBUdfglmZttx0mFWUIr6Sw6vB1ZHxJqIeB64ATg5e0BE3B4Rz6Sri0kmZDQz6xr36TArIEXTMugTJQ1n1udHxPzM+gHAI5n1tcDRDa53NnDzmAM1MxsDJx1mBaXGfTc2RMRQo9NrbKt5QUlnAkPAMfmjMzMbOycdZkUU1EkRclsLHJhZnww8Wn2QpJnARcAxEfFcW89oZtaE+3SYFVL9TqQ5i4YtAaZJOljSzsAcYGH2AEkzgC8DsyPi8Y6/BDOzKm7pMCuiaHp7pfHpEVslnQcsAiYA10bECkmXAsMRsRC4EtgD+K4kgIcjYnb7wZuZ1eakw6yo2iyDHhE3ATdVbbs483hmW09gZjZGTjrMCkqeZdbMSsZJh1kRVSqSmpmViJMOswISgUY9+YqZlYuTDrMickuHmZWQkw6zgnKfDjMrGycdZoUU4NsrZlYyTjrMisi3V8yshHJVJJU0S9IqSaslXVBj/1sk/UrSVkmnVu0bkXRXuiysPtfMalNE3cXMbBA1bemQNAG4BjieZD6HJZIWRsTKzGEPA2cBH6lxic0RcVQHYjUbPwIY8e0VMyuXPLdXXg+sjog1AJJuAE4GtiUdEfGbdJ8/Jc06wn06zKx88txeOQB4JLO+Nt2W166ShiUtlnTKmKIzG68qLR31FjOzAZSnpUM1to3lpvKUiHhU0iHAbZLuiYgHtnsCaR4wD2DKlCljuLRZWQWEkwszK5c8LR1rgQMz65OBR/M+QUQ8mn5dA9wBzKhxzPyIGIqIoUmTJuW9tFl5uaXDzEooT9KxBJgm6WBJOwNzgFyjUCTtK2mX9PFE4E1k+oKYWQOjo/UXM7MB1DTpiIitwHnAIuA+YEFErJB0qaTZAJJeJ2kt8KfAlyWtSE8/DBiWdDdwO3B51agXM6slAkZG6i9mZgMoV3GwiLgJuKlq28WZx0tIbrtUn/dz4Ig2YzQbn1yPw8xKxhVJzQop3HfDzErHSYdZEQWEb6OYWck46TArKt9eMbOScdJhVkSVjqRmZiWSa8I3M+u1IEZG6i555JiocRdJ30n3/0LS1A6/CDOz7TjpMCuioK0hs5mJGk8EpgOnS5peddjZwMaIeBVwFfCZzr4IM7PtOekwK6AAYjTqLjlsm6gxIp4HKhM1Zp0MfD19fCNwnKRa0x6YmXWEkw6zIoq2b6/kmahx2zFpEcCngP06EL2ZWU2F60i6dOnSDZJ+D2zodyzARBxHVlHigOLEMpY4Dsp70afZuOjHowsmNjhkV0nDmfX5ETE/s55nosZ2J3M0MxuTwiUdETFJ0nBEDPU7FsdRzDigOLF0K46ImNXmJfJM1Fg5Zq2kHYG9gSfbfF4zs7p8e8WsnPJM1LgQmJs+PhW4LcLFQcysewrX0mFm7YuIrZIqEzVOAK6tTNQIDEfEQuCrwDclrSZp4ZjTv4jNbDwoatIxv/khPeE4tleUOKA4sRQljhfJMVHjsyQzQ5uZ9YTcmmpmZma94D4dZmZm1hNdTzpylGK+StJd6fJrSZvS7W/NbL9L0rOSTkn3XSfpwcy+ozoQxxRJt0taJmm5pJMy+y5Mz1sl6YS81+x0LJKOl7RU0j3p17dlzrkjvWblPXlZF+OYKmlz5rm+lDnnj9P4Vku6Ok+xqTbieG/Vz8ho5WehS+/HQZJuTWO4Q9LkzL65ku5Pl7mZ7WN+P8zMSisiuraQdGB7ADgE2Bm4G5je4PgPknR4q97+UpKObi9J168DTu1kHCT35j+QPp4O/Cbz+G5gF+Dg9DoTxvraOhTLDGD/9PHhwLrMOXcAQz16T6YC99a57i+BN5LUgLgZOLFbcVQdcwSwpsvvx3eBuenjtwHfzPx8rkm/7ps+3reV98OLFy9eyrx0u6UjTynmrNOBb9fYfipwc0Q808U4Atgrfbw3L9Q0OBm4ISKei4gHgdXp9cb62tqOJSKWRUQlrhUkBaJ2yfGcHY2jHkmvAPaKiDsjIoBvAKf0KI56Pzt55YljOnBr+vj2zP4TgFsi4smI2AjcAsxq8f0wMyutbicdeUoxA0nTNUlLwm01ds/hxX9QLkubua/K8Yc3TxyXAGdKWkvS4/+DTc7N/do6GEvWu4FlEfFcZtvX0lsJf5ujGb/dOA5Ob3f8RNKbM9dc2+SanY6j4j28+Gek0+/H3STvO8CfAHtK2q/Bua28H2ZmpdXtpGMsZZbnADdGxHYTS6T/LR5BUm+g4kLgD4HXkTRpf6wDcZwOXBcRk4GTSOoX7NDg3FZLSLcTS3IB6TUkM4L+Veac90bEEcCb0+XPuhjHY8CUiJgBnA9cL2mvnNfsZBzJBaSjgWci4t7MOd14Pz4CHCNpGXAMsA7Y2uBclxk3M8vodtKRpxRzRa3WDIDTgB9ExJbKhoh4LBLPAV8jaRpvN46zgQXp9e8EdiWZV6PeuWN5bZ2KhbTz4g+A90XEA5UTImJd+vVp4Hq6+J6kt5qeSLcvJekLcWh6zcmZ8/O8J229H6kX/ex04/2IiEcj4l1psnVRuu2pBue28n6YmZVWt5OOPKWYkfRqkg54d9a4xovu1aetH6RN5qcA99Y4b6xxPAwcl173MJI/bOvT4+ZI2kXSwcA0ks6BuV5bJ2ORtA/w/wMXRsTPMu/HjpIqSclOwH/v5nsiaZKkCen2Q9L3ZE1EPAY8LekN6ffmfcA/dyuOdH0HkgJXN3T7/ZA0MdPCciFwbfp4EfB2SftK2hd4O7CoxffDzKy8ut1TlaQ5/Nck/w1flG67FJidOeYS4PIa504lacLeoWr7bcA9JH9IvgXs0W4cJJ0Ef0Zy3/4u4O2Zcy9Kz1tFZvRBrWt24j2pFwvwceD36bbK8jJgd2ApsJykg+nngAldjOPd6fPcDfwKeGfmmkPp9+UB4POkBei6+L05Flhcdb1uvR+nAvenx3wF2CVz7p+TdDJeDby/nffDixcvXsq6uCKpmZmZ9YQrkpqZmVlPOOkwMzOznnDSYWZmZj3hpMPMzMx6wkmHmZmZ9YSTDjMzM+sJJx1mZmbWE046zMzMrCf+L6Cg3i5n+miUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x432 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "prC = [ltnw.ask(\"C_%s(?x)\" % i) for i in range(nr_of_clusters)]\n",
    "n = 2\n",
    "m = (nr_of_clusters + 1) // n + 1\n",
    "\n",
    "fig = plt.figure(figsize=(3 * 3, m * 3))\n",
    "\n",
    "fig.add_subplot(m, n, 1)\n",
    "plt.title(\"groundtruth\")\n",
    "for c in clusters:\n",
    "    plt.scatter(c[:, 0], c[:, 1])\n",
    "data=np.concatenate(clusters)\n",
    "x0 = data[:, 0]\n",
    "x1 = data[:, 1]\n",
    "for i in range(nr_of_clusters):\n",
    "    fig.add_subplot(m, n, i + 2)\n",
    "    plt.title(\"C\" + str(i) + \"(x)\")\n",
    "    plt.scatter(x0, x1, c=prC[i].T[0])\n",
    "    plt.scatter(x0[:2], x1[:2], s=0, c=[0, 1])\n",
    "    plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}