OverallVersion_FullyLearned.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# === Header ===\n",
    "#\n",
    "# @date: 05 / 07 / 2019 (Xi'an)\n",
    "# @author: ZHE LI\n",
    "# @title: Learned KD-Tree - Fully Learned\n",
    "#\n",
    "import math\n",
    "import time\n",
    "import random\n",
    "#import pysnooper\n",
    "import statistics\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "from numpy import genfromtxt\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import Rectangle\n",
    "from sklearn import preprocessing\n",
    "from sklearn.preprocessing import normalize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# === Data Loading ===\n",
    "original_domains = [\\\n",
    "    [1,1.20000000e+07],\\\n",
    "    [1,4.00000000e+05],\\\n",
    "    [1,2.00000000e+04],\\\n",
    "    [1,7.00000000e+00],\\\n",
    "    [1,5.00000000e+01],\\\n",
    "    [900.99,1.04899500e+05],\\\n",
    "    [0,1.00000000e-01],\\\n",
    "    [0,8.00000000e-02]\\\n",
    "]\n",
    "original_domains = np.asarray(original_domains)\n",
    "original_dataset = genfromtxt('/Users/lizhe/Desktop/LearnedKDTree/DataAndWorkload/SyntheticData/TPCH_12M_8Field.csv', delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# === Global Configuration ===\n",
    "\n",
    "# dimension for query and build learned KD-Tree\n",
    "used_dimensions = [1,2]\n",
    "\n",
    "# testing set percentage\n",
    "testing_set_percentage = 0.2\n",
    "\n",
    "# percentage of each dimension\n",
    "percentage_ = [1, 1]\n",
    "\n",
    "# maximum expected ranges\n",
    "ranges_ = [10000, 500]\n",
    "\n",
    "# maximum expected sigma\n",
    "sigmas_ = [1000, 50]\n",
    "\n",
    "# clusters amount\n",
    "cluster_amount_ = 100\n",
    "\n",
    "# queries in each cluster\n",
    "queries_each_cluster_ = 10\n",
    "\n",
    "# page block size\n",
    "block_size = 32000\n",
    "\n",
    "domains_ = np.copy(original_domains)\n",
    "domains_ = domains_[used_dimensions]\n",
    "dataset = np.copy(original_dataset)\n",
    "dataset = dataset[:,used_dimensions]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# # measure the differnet calculation method of median index\n",
    "# np.random.shuffle(dataset)\n",
    "\n",
    "# # === calculated by sorting ===\n",
    "# start_time = time.time()\n",
    "# total_size = len(dataset)\n",
    "# dataset = dataset[dataset[:, 1].argsort()]\n",
    "# median = dataset[int(total_size / 2), 1]\n",
    "# end_time = time.time()\n",
    "# print(\"find median by sorting(s): \", end_time-start_time)\n",
    "# print(\"median: \", median)\n",
    "\n",
    "# np.random.shuffle(dataset)\n",
    "\n",
    "# start_time = time.time()\n",
    "# mid = np.median(dataset[:,1])\n",
    "# end_time = time.time()\n",
    "# print(\"find median by fast median(s): \", end_time-start_time)\n",
    "# print(\"median: \", mid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# print('variance for the first query dimension: ',np.var(dataset[:,0]))\n",
    "# print('variance for the second query dimension: ',np.var(dataset[:,1]))\n",
    "# print('variance for the third query dimension: ',np.var(dataset[:,2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Try to improve query generation:\n",
    "# 1. the number of queries in each cluster should be random (under a maximum number, maybe under zipfian distribution)\n",
    "# 2. the number of clusters should be random (under a maximum number)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# === Query Generation ===\n",
    "#\n",
    "# by default, use Gaussian distribution to generate the synthetic query workload\n",
    "#\n",
    "# @cluster_amount: the cluster amount of the query; integer\n",
    "# @queries_in_each_cluster: the number of queries in each cluster; interger\n",
    "# @dimensions: the domain of each dimension; numpy object\n",
    "# @maximum_expected_range: the maximum expected range of each dimension; array object\n",
    "# @maximum_sigma: the maximum sigma of Gaussian distribution of each dimension; array object\n",
    "#\n",
    "# return @temp_queries: array object\n",
    "#\n",
    "def generate_query_clustered(cluster_amount, queries_in_each_cluster, dimensions, maximum_expected_range, maximum_sigma, percentage):\n",
    "    random.seed()\n",
    "    num_dims = len(dimensions)\n",
    "    \n",
    "    queries = []\n",
    "    for j in range(num_dims):\n",
    "        queries.append([])\n",
    "    \n",
    "    # for each cluster\n",
    "    for i in range(cluster_amount):\n",
    "        \n",
    "        # determin the cluster center on each dimension\n",
    "        cluster_center = []\n",
    "        \n",
    "        # for each dimension\n",
    "        for j in range(num_dims):\n",
    "            \n",
    "            # to assure the distribution do not exceed the border\n",
    "            dim_range = dimensions[j][1]-dimensions[j][0]\n",
    "            center1D = random.randrange(dim_range)\n",
    "            cluster_center.append(center1D)\n",
    "            \n",
    "            query_range = random.randrange(maximum_expected_range[j])\n",
    "            range_lower = center1D - query_range/2\n",
    "            range_upper = center1D + query_range/2\n",
    "            \n",
    "            if range_lower <= dimensions[j][0]:\n",
    "                range_lower = dimensions[j][0]\n",
    "            if range_upper >= dimensions[j][1]:\n",
    "                range_upper = dimensions[j][1]\n",
    "            \n",
    "            sigma = random.randrange(maximum_sigma[j])\n",
    "            \n",
    "            # for each query\n",
    "            for k in range(queries_in_each_cluster):\n",
    "                \n",
    "                # determine if to let this query dimension invalid by percentage (i.e., cover the entire domain)\n",
    "                possibility = random.uniform(0, 1)\n",
    "                if possibility > percentage[j]:\n",
    "                    lower1D = dimensions[j][0]\n",
    "                    upper1D = dimensions[j][1]\n",
    "                    queries[j].append([lower1D, upper1D])\n",
    "                else:\n",
    "                    #determin the lower bound\n",
    "                    lower1D = random.gauss(range_lower, sigma)\n",
    "                    upper1D = random.gauss(range_upper, sigma)\n",
    "                    if lower1D <= dimensions[j][0]:\n",
    "                        lower1D = dimensions[j][0]\n",
    "                    if lower1D >= dimensions[j][1]:\n",
    "                        lower1D = dimensions[j][1]\n",
    "                    if upper1D <= dimensions[j][0]:\n",
    "                        upper1D = dimensions[j][0]\n",
    "                    if upper1D >= dimensions[j][1]:\n",
    "                        upper1D = dimensions[j][1]\n",
    "                    if lower1D > upper1D:\n",
    "                        lower1D, upper1D = upper1D, lower1D\n",
    "                    queries[j].append([lower1D, upper1D])\n",
    "                \n",
    "    # transform the query structure\n",
    "    temp_queries = []\n",
    "    values = []\n",
    "    for i in range(len(queries[0])):\n",
    "        for j in range(num_dims):\n",
    "            values.append(queries[j][i])\n",
    "        temp_queries.append(values)\n",
    "        values=[]\n",
    "            \n",
    "    # return queries\n",
    "    return temp_queries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Example of usage (Generate Query)\n",
    "\n",
    "#@pysnooper.snoop()\n",
    "def query_seperation(query_collection, training_set_percentage = 0.8):\n",
    "    np.random.shuffle(query_collection)\n",
    "    training_set_size = int(training_set_percentage * len(query_collection))\n",
    "    query_collection = np.asarray(query_collection) # as the shuffle make np to array, we need to make it back\n",
    "    training_set = query_collection[0:training_set_size,:]\n",
    "    testing_set = query_collection[training_set_size:-1,:]\n",
    "    return training_set, testing_set\n",
    "\n",
    "query_collection = generate_query_clustered(cluster_amount_, queries_each_cluster_, domains_, ranges_, sigmas_, percentage_)\n",
    "\n",
    "training_set, testing_set = query_seperation(query_collection, 1-testing_set_percentage)\n",
    "# training_set = query_collection\n",
    "# testing_set = query_collection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAD8CAYAAABU4IIeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3X2sXHd95/H3Z+8loeUpDtxErp3c\nG5AJTarKJLMhFQvNspvESREOq3bX0aqxaCTzkEiw7WpxirQxsCtRWgqKygYMeJNsISE8bSxENrhe\nKCttIL4G4ziE4Jtwb3OxZRvMQySqdK/z3T/mN/G543meMzPnzHxe0tGc+c05c35zZuZ8z+/h/I4i\nAjMzs7z8s1FnwMzMxosDi5mZ5cqBxczMcuXAYmZmuXJgMTOzXDmwmJlZrtoGFkkXSPqGpMclPSbp\n3Sn9XEl7JB1Oj2tSuiTdIWlB0kFJl2Xea2ta/rCkrZn0yyU9mta5Q5IG8WHNzGzwOimxrAB/FhG/\nDVwJ3CLpEmA7sDciNgB703OA64ANadoG3AnVQATcDrwOuAK4vRaM0jLbMutt6v+jmZnZKLQNLBFx\nNCK+m+afAR4H1gGbgbvTYncDN6T5zcA9UfVt4BxJa4FrgT0RcTIifg7sATal114aEQ9H9WrNezLv\nZWZmJTPdzcKS5oDXAt8Bzo+Io1ANPpLOS4utA57OrLac0lqlLzdIb7T9bVRLNrzoRS+6/DWveU03\n2Tczm3j79+//aUTMDHIbHQcWSS8GvgS8JyJ+1aIZpNEL0UP6mYkRO4GdAJVKJebn59tl2wZhbg6W\nllov46GCzApJUps/b/866hUm6QVUg8pnI+LLKflYqsYiPR5P6cvABZnV1wNH2qSvb5BuRdUuqJjZ\nROukV5iAzwCPR8RfZ17aDdR6dm0FHsik35R6h10J/DJVmT0EXCNpTWq0vwZ4KL32jKQr07ZuyryX\nFV3E6cnMjM6qwl4P/DHwqKQDKe3PgQ8B90u6GfgH4I/Sa18DrgcWgF8DbwOIiJOSPgjsS8t9ICJO\npvl3AncBvwE8mCYzMyshlXXYfLexjFC2fS37+2mWbmaFIWl/RFQGuQ1feW/Nzc1Vg0X9lNUs3cwm\nVlfdjQup/oDmM+X81Brp6/epSyZm1oJLLGZmlqvyB5apqVHnwMzMMsofWFZW3OV1UGZnq4/N2lgc\n1M1Go1HbZ6u20CErfxuLDc7i4qhzYGb15uZGnYO2yl9iMTOzQnFg6UStaDntAp5Z1wpSPTM2SlCT\n4CNlN06dGnUOzMwK357sEksnZmerU8G/TDOzInCJpRMlKHqamRWFSyxm1rvp6fZtkLOz8LKXucQ/\nQVxiMbPe7Nixut2xWRukS/yD16hjxAgDuUssZtabHTtGnQNrZYTfj0ssZu1MT3fWI3B2dvLOzmdn\nTw9WWhupwYavYNWMDixm7bibeXOTFkitI64KM2unkzPxSSytmDXhwGKdmeT69MXF1QOdNpocVMye\n1zawSNol6bikQ5m0z0s6kKZFSQdS+pykf8y89onMOpdLelTSgqQ7pGo3BknnStoj6XB6XDOID1pY\nBR6h9Hlzc/D+9xcnP2ZWaJ2UWO4CNmUTIuLfRcTGiNgIfAn4cublJ2uvRcQ7Mul3AtuADWmqved2\nYG9EbAD2pueTp/4MuEhqjbNmZh1oG1gi4lvAyUavpVLHvwXubfUektYCL42IhyMigHuAG9LLm4G7\n0/zdmfTJUsSSiplZD/ptY3kDcCwiDmfSLpL0PUl/L+kNKW0dsJxZZjmlAZwfEUcB0uN5feapXNrV\n3ReBu5GaWRf67W58I6tLK0eBCyPiZ5IuB/6npEuBRqfgXR81JW2jWp3GhRde2EN2rSdumDazLvRc\nYpE0Dfwb4PO1tIh4NiJ+lub3A08Cr6ZaQlmfWX09cCTNH0tVZbUqs+PNthkROyOiEhGVmZmZXrNu\nZkWVHXssO5Xgrol2Wj9VYf8a+GFEPF/FJWlG0lSafyXVRvqnUhXXM5KuTO0yNwEPpNV2A1vT/NZM\nutlkm5sr/L3Nc1d/MWrtdhXj3oGkWUCtn0pys8FOuhvfCzwMXCxpWdLN6aUtnNlo/0bgoKTvA18E\n3hERtYb/dwKfBhaolmQeTOkfAq6WdBi4Oj03s+zBtKjtb4O2tOTSSlZJRoFQlPQHWqlUYn5+ftTZ\nMBucTksjJf0PN9TqM4/T52xkbq59ySyHER4k7Y+ISl9v0kY5ylVmkyg7wKONvzHqJOMhXcxGKduO\nUl9/3slQMuN+Fm+l5MBiNkrZEklJ6s9HwtdSlYoDi5kVR6MA4pGjS2f82lia3ZTJP06z4vN/dCyM\nV2Bp1aPEjaBmVq/+mOE2q1y4KsxslKamGs/b4I3bxaUFMl4lFrOyWVkZdQ7McjeegaW+OOszEzOz\noRnPwOJAYmbtuD1lYNzGYmZmuRqvEovPQMzMRs4lFjMzy5UDi5mZ5cqBxczMcjVebSxmNlyTeM8Y\na8uBxcyGLxuQHHTGjqvCzGw0HFDGlkssZta7foKDL2QeW21LLJJ2STou6VAmbYekn0g6kKbrM6/d\nJmlB0hOSrs2kb0ppC5K2Z9IvkvQdSYclfV7SWXl+QDMrIN8Bc6x1UhV2F7CpQfpHI2Jjmr4GIOkS\nYAtwaVrnv0makjQFfBy4DrgEuDEtC/AX6b02AD8Hbu7nA5mZ2Wi1DSwR8S3gZIfvtxm4LyKejYgf\nAwvAFWlaiIinIuKfgPuAzZIEvAn4Ylr/buCGLj+DmZkVSD+N97dKOpiqytaktHXA05llllNas/SX\nA7+IiJW69IYkbZM0L2n+xIkTfWS9B9LqyczMGuo1sNwJvArYCBwFPpLSGx1xo4f0hiJiZ0RUIqIy\nMzPTXY7NzGwoeuoVFhHHavOSPgV8NT1dBi7ILLoeOJLmG6X/FDhH0nQqtWSXNzOzfoyodqWnEouk\ntZmnbwVqPcZ2A1sknS3pImAD8AiwD9iQeoCdRbWBf3dEBPAN4A/T+luBB3rJUweZbj2ZmY2zIfbA\na1tikXQvcBXwCknLwO3AVZI2Uq22WgTeDhARj0m6H/gBsALcEhGn0vvcCjwETAG7IuKxtIn3AvdJ\n+i/A94DP5Pbp8jQ1BadOnZ633szNwdJS49dmZ2FxcZi5MZscQzyBVpS0H3mlUon5+flRZ8O61e7H\nXdLfo3Vgevr0yVkjPrHIX4MTOcH+iKgMcrO+8t5GpxZEXBVpNhiNAvUQ/m8eK8xGx+1bk2Vlpflr\nU1MurYwRBxaoFhfbNe5Pu3DXk/r96PapyZYdyiU7tQo6Vjo+WkLzxuSsVnXD1rlTp9yOYjbmXGKB\naqNhL8u4y7KZ2RlcYgHX7Zp1InvyNDXl6itryiWWfnno7+50Ujq04pmbW/3cVcPWggOLDVbE6gb7\npSWPfGA25lwV1i8fFNvz2W35LS6W77feaX6HXePQKF9jVuvhwGLDVf8HKtvBapKN2cHPBseBpR/+\no5lZr8Z45AkHFhuuMfwTWYmN8oLdMf4vuPHeBq9VT7DZWZf8utXJSBFjfNDqWLOr/H3F/8C5xGKD\n5+uE8tXJSBF2Wm1U5aKMnjwBJ1IOLGZl1ewA5dLKarVeiQ7IQ+PAYlZWDiBWUG5jMSsbjxBtBecS\ni1nZuMG5O7Xbins4oaFpW2KRtEvScUmHMml/KemHkg5K+oqkc1L6nKR/lHQgTZ/IrHO5pEclLUi6\nQ6qW4yWdK2mPpMPpcc0gPqiZTaiVFfj93y9Gw/2E6KQq7C5gU13aHuB3IuJ3gR8Bt2VeezIiNqbp\nHZn0O4FtwIY01d5zO7A3IjYAe9NzM7P8fPObo85Bb6anS9mdvG1giYhvASfr0r4eEbXy+LeB9a3e\nQ9Ja4KUR8XBEBHAPcEN6eTNwd5q/O5NuZjbZsuPslWgk9Twa7/8EeDDz/CJJ35P095LekNLWAcuZ\nZZZTGsD5EXEUID2e12xDkrZJmpc0f+LEiRyybmZmeesrsEh6H7ACfDYlHQUujIjXAn8KfE7SS4FG\n5bauQ29E7IyISkRUZmZmes22mVk5ZHsAlqAKrKbnwCJpK/Bm4N+n6i0i4tmI+Fma3w88Cbyaagkl\nW122HjiS5o+lqrJaldnxXvNkZta1dkPk1N/kbBCabXtlpXFvtoL3cOspsEjaBLwXeEtE/DqTPiNp\nKs2/kmoj/VOpiusZSVem3mA3AQ+k1XYDW9P81ky6mdngtbsif5hX7DdqR1lcPHOMs4L3cOuku/G9\nwMPAxZKWJd0M/A3wEmBPXbfiNwIHJX0f+CLwjoioNfy/E/g0sEC1JFNrl/kQcLWkw8DV6bmZWXv1\npY1eShftzv4LXjooIkVJehnUq1QqMT8/P+psmE2e2sG7CGfNzdocynZcq32O+nu0DOBzSNofEZXc\n3zjDV96Pu9rIrlC+P5sVT/b3ND3tUQDyVpLG+XYcWMad7zdvecr+nvzbys+YnfR5EEozK6+I020g\ntZvGjdlBuoxcYjGzcitCW4+t4hKLmXUue8Geh++3JlxiGXe1IcPN8uDGeuuASyzjrnbl7qD64re6\nankYVyyPevvjLrt/vT+tQw4sk2BxcXD10K2uSl5aGtwBqXbAa7f9oh0MWw0dUrQD+Nzc6v27tNR8\nGHezDAeWSZLHVcqtRDSud68FmDy3ubR0Zu+fRj2ChjkcR16KkudG+XC1qnXAbSyTotFZ5SAOYK0O\nPBE+u61p1CW2TPsmm/8y5duGwoHF8tPuACN53KUaH4xtjLkqzPrTaaCoXbw2yGsOxqW+vwzB1+0r\n1oJLLNafol+cVrSDdNmvCp+dheXl1VWevp7F6rjEYqPRb0+o2dn2Z8tTU8UPfEXWKCgvLla7sGfv\nDeJrW6yOSyyTouhnyt12JHDAGDzvY+uRSyw2Oh4w0GwsObCYmVmuHFiGrf4ixUnmfWA2PHNz8MIX\nDmVTHQUWSbskHZd0KJN2rqQ9kg6nxzUpXZLukLQg6aCkyzLrbE3LH5a0NZN+uaRH0zp3SGN8tCnK\nVdVmrTQbuqUow81Y95aW4Nlnh7KpTkssdwGb6tK2A3sjYgOwNz0HuA7YkKZtwJ1QDUTA7cDrgCuA\n22vBKC2zLbNe/bZs3GR7FbmtpViytx+u5xMj60BHgSUivgWcrEveDNyd5u8Gbsik3xNV3wbOkbQW\nuBbYExEnI+LnwB5gU3rtpRHxcEQEcE/mvcaP+/xbr7Ilh+kBduisBRUHfutRP20s50fEUYD0eF5K\nXwc8nVluOaW1Sl9ukH4GSdskzUuaP3HiRB9ZH6H6Pv9lGSW2k1F5i5z/cePBIK3ABtF43+joEj2k\nn5kYsTMiKhFRmZmZ6SOLBVDms8FG1Vhl+wxmk2aI/9N+AsuxVI1Fejye0peBCzLLrQeOtElf3yDd\nbPBa3SjMJTCznvQTWHYDtZ5dW4EHMuk3pd5hVwK/TFVlDwHXSFqTGu2vAR5Krz0j6crUG+ymzHuN\nr2YHMPe86U2jANHJfnRj9Jlq7YCjDLLZXmlWOp12N74XeBi4WNKypJuBDwFXSzoMXJ2eA3wNeApY\nAD4FvAsgIk4CHwT2pekDKQ3gncCn0zpPAg/2/9EKrF01UvbGWA40nWkUIBw0erOy0ryTyTA6n8zN\nrW5DGvXvv9FJyyA7T4wBRUnrxiuVSszPz486G4PR6Cxt1N9Tp2eOo8pns/y1y0/97XezZmc9XtYo\nFO333+tvq6Ak7Y+IyiC34bBbBkUY+r2kf6K2HDisU74DasfGY0iXcWt0ra9u8MGvvUbBtwgB2cqv\n9jsq+3FliFxiKSLf36J7Dr7jY3Z2dfXkqE8QFhfPDCqjzlPBjUdgGddqmnHSqi0D/B3aaUU8SfDv\nsyvjURVmZmaF4cBiw1HEs1AzG4jxqAqzcnB1gtlEcInFzMxy5cBiZsWTvdp91FfeW9ccWMyseLI9\nCJeWPIRKyTiwmFlx1drlfP+ZUnFgMTOzXDmwmFlxeRiVUnLFZb3sD9ndY80Gr5Pg4SFUSsWBxcyK\npX4UYZ/glY6rwhrxD9nMrGcOLGZmlisHFjMzy1XPgUXSxZIOZKZfSXqPpB2SfpJJvz6zzm2SFiQ9\nIenaTPqmlLYgaXu/H6pv7oliNjrZ/1/9Te+sFHpuvI+IJ4CNAJKmgJ8AXwHeBnw0Iv4qu7ykS4At\nwKXAbwF/J+nV6eWPA1cDy8A+Sbsj4ge95s3MSsRtmmMnr15h/wp4MiKW1PxsfzNwX0Q8C/xY0gJw\nRXptISKeApB0X1p2NIHFP3Izs77k1cayBbg38/xWSQcl7ZK0JqWtA57OLLOc0pqln0HSNknzkuZP\nnDiRU9bNzCxPfQcWSWcBbwG+kJLuBF5FtZrsKPCR2qINVo8W6WcmRuyMiEpEVGZmZvrKt5nVyY4o\nnJ08urB1KY8Sy3XAdyPiGEBEHIuIUxHxHPApTld3LQMXZNZbDxxpkW5mw5QdUbg+3cHFupBHYLmR\nTDWYpLWZ194KHErzu4Etks6WdBGwAXgE2AdskHRRKv1sScuaWTvNShm1qZfh5iNOTzXNgk4zLu1M\ntL4a7yX9JtXeXG/PJH9Y0kaq1VmLtdci4jFJ91NtlF8BbomIU+l9bgUeAqaAXRHxWD/5MpsYS0vV\ncbSaHfhHPdx8twHJxoKipL2gKpVKzM/PjzobNunm5hofPGdnYXFx8Ntv1gszG2w6/Y+3un6r28/j\nsb4KS9L+iKgMchu+8t6sH63aJUYp7+33GiQ9KvFE8ujGZnnInpUPc+SGVtVgtdc7lWfJwqWUiebA\nYlZmw6huM+uSA4tZHjy+nNnz3MZi1g8Pkmh2BpdYzPqxsjLqHJgVjkssZmaWK5dYzCx/rdqc3GNs\n7LnEYmb5ajeMSy/DzIyL7HA7Y8yBxczyVX9dTf24Y6MeZsYGzoHFzAZrAs7QB67ZIKM7dow6Zw05\nsJhZvtpd7e8u2vmYmnJgMbMJ0W40AHfR7l7EmQF7/frR5KUDE9yKZmYD455f+SvR8D0usZjZeGp1\nE7RR3YCs0U3UxpADi5mNp1ajPo/6tgZjzoHFrJ1mZ76TctvdubnJ+ayWC7exmLXS7A6RUE2fnh7v\nxmh3E7Ye9F1ikbQo6VFJByTNp7RzJe2RdDg9rknpknSHpAVJByVdlnmfrWn5w5K29psvs6GYpIv9\nXGqxDuVVFfYvI2Jj5j7K24G9EbEB2JueA1wHbEjTNuBOqAYi4HbgdcAVwO21YGQ2UouLra/L8K13\ni8vf28gMqipsM3BVmr8b+Cbw3pR+T0QE8G1J50ham5bdExEnASTtATYB9w4of2adK1E3z4Eq6n5o\nVV1ZM+a9sIomjxJLAF+XtF/StpR2fkQcBUiP56X0dcDTmXWXU1qz9FUkbZM0L2n+xIkTOWR9QCa1\nkbcTEzII39goQ/fYbFDJ5rfo+R5jeQSW10fEZVSruW6R9MYWyzY6mkSL9NUJETsjohIRlZmZmd5y\nOwru2midmPTeZzY2+g4sEXEkPR4HvkK1jeRYquIiPR5Piy8DF2RWXw8caZE+Hlyfe1ptnKhG40XV\nDqyTeiBt1fvMmsv+vxoFZhu6vgKLpBdJekltHrgGOATsBmo9u7YCD6T53cBNqXfYlcAvU1XZQ8A1\nktakRvtrUlo51RfFi1o3PQorK9V9Ut9FN1tP7gOpdaNdBwsPejl0/Tbenw98RdWzgmngcxHxvyTt\nA+6XdDPwD8AfpeW/BlwPLAC/Bt4GEBEnJX0Q2JeW+0CtIb9Qpqfbdy91nW5vHEyqB8BJ6r6cJ5+8\nFYqipAfCSqUS8/Pzw91oJ8XqXvZno4A1OztZf5ZsiWXSPntWs5OXkv5Px9qOHfD+9zd/XYLnnhta\ndjolaX/m0pCB8JAuo9bsQDJpZ/CLi646hNNVhe7ZVHw7drSugrvwwqFlpWg8pItZGfhajWKa5JOg\nFlxi6cYg7ozX7AzVBwnL6qQEO+3zRCsG/xK74bMTKzI3/FtBuMRiVgadlIZ9vZQVhEssZmUwzkPz\n29hxicXMzHLlwDJppqc9BpWZDZQDy6TJNvBO2rUyZjYUDixmZpYrBxYzM8uVA4uZmeXKgWXS1K51\nmJ311f1mNhC+jqWdRoNETk2V97oCjx5gZgPmEks7jYbJ8NAZZmZNObCY2WSo3fra13ENnAOLmU2G\n+uu2fB3XwDiwtOP7ZZuNF3daGbieA4ukCyR9Q9Ljkh6T9O6UvkPSTyQdSNP1mXVuk7Qg6QlJ12bS\nN6W0BUnb+/tITfRaDF5ZOXPUWI8ia1Zendxi3PrST6+wFeDPIuK7kl4C7Je0J7320Yj4q+zCki4B\ntgCXAr8F/J2kV6eXPw5cDSwD+yTtjogf9JG31Rr9kJaWqj2+OundNU49qTq5EyFM9n3nJ1mtF2SZ\nez42MzW1uuONTxAHpufAEhFHgaNp/hlJjwPrWqyyGbgvIp4FfixpAbgivbYQEU8BSLovLZtfYLHT\nOq1Xdv3zZKodeMex5+O4BcoCy6WNRdIc8FrgOynpVkkHJe2StCalrQOezqy2nNKapeen2a1/J/2H\n5tshm9kA9B1YJL0Y+BLwnoj4FXAn8CpgI9USzUdqizZYPVqkN9rWNknzkuZPnDjRb9Yt2+ZUm2pc\nTWD16tsp8+i62+w93RW41PoKLJJeQDWofDYivgwQEcci4lREPAd8itPVXcvABZnV1wNHWqSfISJ2\nRkQlIiozMzP9ZH1ydRIw3L4yuVqVXltVj/ZSdSo1X89VsaXWT68wAZ8BHo+Iv86kr80s9lbgUJrf\nDWyRdLaki4ANwCPAPmCDpIsknUW1gX93r/myNhYXm1cN1iYHFbOqVqW0+lK+Pa+fXmGvB/4YeFTS\ngZT258CNkjZSrc5aBN4OEBGPSbqfaqP8CnBLRJwCkHQr8BAwBeyKiMf6yJeZlUl96ahIB2uXnHqi\nKGmDbaVSifn5+VFnw8x6VQsgzQJLEY5NnQS5IuSzC5L2R0RlkNvwlfe9crHYLB9F/g+5E0tPPGy+\nmRVPUYZScntjTxxY+pUtBhftbMs60+p7G8cr0IuiZFVI1jlXhdlka3e9xDhegW42YC6x9MullHJz\nrx+z3LnEYpPNjbNmuXOJpVeuHx4Pbpw1y50Di5nZOBtBdb2rwszMJsEQa1kcWMzMLFcOLGaDdtVV\no86B2VCrxNzGYjZItT+z5A4fNhojGOTTJRYzM8uVA4uZmeXKVWFmZlntqopcpdmWSyxmZpYrBxYz\ns6xWw/x4CKCOuCps3NWK9R7+fTRmZ+EXv6hOVg4e5qdvDiyTwsO/j4YPUjaBClMVJmmTpCckLUja\nPur8mJlZbwoRWCRNAR8HrgMuAW6UdMloc2VmZr0oRGABrgAWIuKpiPgn4D5g84jzNF6Kcg9xMxt7\nRWljWQc8nXm+DLyufiFJ24Bt6emzkg4NIW/9egXw01FnglOn2vXPL0Y+WytDHsH5zJvzma+LB72B\nogSWRke8M65CioidwE4ASfMRURl0xvrlfOanDHkE5zNvzme+JM0PehtFqQpbBi7IPF8PHBlRXszM\nrA9FCSz7gA2SLpJ0FrAF2D3iPJmZWQ8KURUWESuSbgUeAqaAXRHxWJvVdg4+Z7lwPvNThjyC85k3\n5zNfA8+nwgOqmZlZjopSFWZmZmPCgcXMzPIVEaWagE3AE8ACsH2I210EHgUOAPMp7VxgD3A4Pa5J\n6QLuSHk8CFyWeZ+tafnDwNZM+uXp/RfSuuowX7uA48ChTNrA89VsG13mcwfwk7RPDwDXZ167LW3z\nCeDadt8/cBHwnZSfzwNnpfSz0/OF9PpcizxeAHwDeBx4DHh3Efdni3wWbX++EHgE+H7K5/t7fe+8\n8t9lPu8CfpzZnxsL8D+aAr4HfLWI+/L5dTo5eBVlSjv1SeCVwFnph3DJkLa9CLyiLu3DtS8A2A78\nRZq/Hngw/QCvBL6T+RE9lR7XpPnaQeoR4PfSOg8C13WYrzcCl7H6gD3wfDXbRpf53AH8xwbLXpK+\n27PTj/rJ9N03/f6B+4Etaf4TwDvT/LuAT6T5LcDnW+RxLekgAbwE+FHKS6H2Z4t8Fm1/Cnhxmn8B\n1YPTld2+d5757zKfdwF/2GD5Uf6P/hT4HKcDS6H25fP57OTgVZQpfTEPZZ7fBtw2pG0vcmZgeQJY\nm+bXAk+k+U8CN9YvB9wIfDKT/smUthb4YSZ91XId5G2O1Qfsgeer2Ta6zOcOGh8IV32vVHsL/l6z\n7z/9WX8KTNf/TmrrpvnptFynpcEHgKuLuj8b5LOw+xP4TeC7VEfU6Oq988x/l/m8i8aBZSTfO9Xr\n+/YCbwK+2sv3NKx9WbY2lkZDv6wb0rYD+Lqk/WloGYDzI+IoQHo8r00+W6UvN0jv1TDy1Wwb3bpV\n0kFJuySt6TGfLwd+ERErdemr3iu9/su0fEuS5oDXUj17Lez+rMsnFGx/SpqSdIBqNegeqmfF3b53\nnvnvKJ8RUduf/zXtz49KOrs+nx3mJ6/v/WPAfwKeS897+Z4Gvi+hfI33HQ39MiCvj4jLqI7AfIuk\nN7ZYtlk+u03PW9HydSfwKmAjcBT4SErPM59dfwZJLwa+BLwnIn7VatEc89m1Bvks3P6MiFMRsZHq\n2fYVwG/38N4D38/1+ZT0O1TP2F8D/HOq1VvvzTmfHZP0ZuB4ROzPJrd435HtSyhfYBnZ0C8RcSQ9\nHge+QvVPckzSWoD0eLxNPlulr2+Q3qth5KvZNjoWEcfSH/o54FNU92kv+fwpcI6k6br0Ve+VXn8Z\ncLJZniS9gOrB+rMR8eU2n3Vk+7NRPou4P2si4hfAN6m2SXT73nnmv9N8boqIo1H1LPDf6X1/5vG9\nvx54i6RFqqO/v4lqCaaY+7KTeseiTFTrCp+i2uhUa2C6dAjbfRHwksz8/6Xag+IvWd3w9uE0/wes\nbtx7JKWfS7WXyZo0/Rg4N722Ly1ba9y7vov8zbG67WLg+Wq2jS7zuTYz/x+A+9L8paxuYHyKauNi\n0+8f+AKrGxjfleZvYXUj5v0t8ifgHuBjdemF2p8t8lm0/TkDnJPmfwP4P8Cbu33vPPPfZT7XZvb3\nx4APFeR/dBWnG+8LtS+fz2OnB6+iTFR7ZPyIal3t+4a0zVemHV3rjvi+lP5yqo1ph9Nj7Uckqjcu\ne5JqF8NK5r3+hGp3vgXgbZlrQ/jzAAAAuklEQVT0CnAorfM3dN7AfC/Vao//R/Ws4+Zh5KvZNrrM\n5/9I+ThIdWy47IHxfWmbT5DpIdfs+0/f0SMp/18Azk7pL0zPF9Lrr2yRx39BtZh/kEyX3aLtzxb5\nLNr+/F2qXWMPps/8n3t977zy32U+/3fan4eAv+V0z7GR/Y/SOldxOrAUal/WJg/pYmZmuSpbG4uZ\nmRWcA4uZmeXKgcXMzHLlwGJmZrlyYDEzs1w5sJiZWa4cWMzMLFf/H/xK3sqfkXqvAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1bdb21b160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# === Query Visualization ===\n",
    "#\n",
    "# plot the queries, only in 2D\n",
    "#\n",
    "# @query_collection: could be more than 2 dimensions, but only plot the first 2 dimensions\n",
    "# @domains: the domain of the plot dimensions, shoud be corresponding with the order of plotdim\n",
    "#\n",
    "def plot_queries_2d(query_collection, domains):\n",
    "    \n",
    "    fig, ax = plt.subplots(1)\n",
    "    \n",
    "    plt.xlim(domains[0][0],domains[0][1])\n",
    "    plt.ylim(domains[1][0],domains[1][1])\n",
    "    \n",
    "    # transform the queries into boxes\n",
    "    query_size = len(query_collection)\n",
    "    # the number of boxed that will be plotted\n",
    "    boxes = []\n",
    "    for i in range(query_size):\n",
    "        lower1 = query_collection[i][0][0]\n",
    "        upper1 = query_collection[i][0][1]\n",
    "        lower2 = query_collection[i][1][0]\n",
    "        upper2 = query_collection[i][1][1]\n",
    "        rect = Rectangle((lower1,lower2),upper1-lower1,upper2-lower2,fill=False,edgecolor='r',linewidth=1)\n",
    "        ax.add_patch(rect)\n",
    "    \n",
    "    plt.show()\n",
    "\n",
    "plot_queries_2d(query_collection, domains_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# === Assistant Functions ===\n",
    "#\n",
    "# this works for one dimension only !!! An implementation of query bounding.\n",
    "\n",
    "def getoverlap(al, au, bl, bu):\n",
    "    return max(0, min(au,bu)-max(al,bl))\n",
    "\n",
    "# currently not used.\n",
    "def bounding_union(query_collection):\n",
    "    \n",
    "    # should keep it ordered first by the lower interval !!!!!!\n",
    "    query_collection = query_collection[query_collection[:,0].argsort()]\n",
    "    \n",
    "    remaining_query = query_collection\n",
    "    bounded_intervals = []\n",
    "    \n",
    "    while len(remaining_query) != 0:\n",
    "        \n",
    "        initial_interval = [remaining_query[0][0], remaining_query[0][1]]\n",
    "        temp_interval = []\n",
    "        \n",
    "        for i in range(len(remaining_query)-1):\n",
    "            \n",
    "            overlap = getoverlap(initial_interval[0],initial_interval[1],remaining_query[i+1][0], remaining_query[i+1][1])\n",
    "            \n",
    "            # there is no overlap\n",
    "            if overlap == 0:\n",
    "                temp_interval.append([remaining_query[i+1][0], remaining_query[i+1][1]])\n",
    "            else: # update interval border\n",
    "                initial_interval[0] = min(initial_interval[0], remaining_query[i+1][0])\n",
    "                initial_interval[1] = max(initial_interval[1], remaining_query[i+1][1])\n",
    "                \n",
    "        bounded_intervals.append(initial_interval)\n",
    "        remaining_query = temp_interval\n",
    "    \n",
    "    return bounded_intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 893,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "training time for Learned KD-Tree(s):  33.436580657958984\n",
      "learned KD-Tree leaf nodes:  512\n"
     ]
    }
   ],
   "source": [
    "# === Learned KD-Tree Split ===\n",
    "#\n",
    "# asssumption: the query boundings will not overlap. divide the KD-Tree recursively\n",
    "#\n",
    "# @dataset: contains the data only in this subnode, ordered in original load order;\n",
    "# @query: contains all the queries (bounded or original); ordered in original generated order;\n",
    "# @domains: the current domain of the node of every dimension [first lower, second upper],[]...; array object\n",
    "# @threshold: maximum page size\n",
    "# @level: the current tree depth\n",
    "#\n",
    "# return @kdnodes: contains the domain of each node and the correpsonding records amount, notice the domain is\n",
    "# ordered by the original load order as dataset\n",
    "\n",
    "def LearnedResuriveDivide(dataset, query, domains, threshold, level, current_dim = 0):\n",
    "    # print(\"level: \",level)\n",
    "    # print(\"dataset size: \",len(dataset))\n",
    "\n",
    "    # check if the threshold is already satisfied\n",
    "    total_size = len(dataset)\n",
    "    if total_size <= threshold:\n",
    "        # the kdnodes should be an global object outside the function\n",
    "        kdnodes = []\n",
    "        kdnodes.append([domains, total_size])\n",
    "        return kdnodes\n",
    "\n",
    "    split_distance_each_dim = []\n",
    "    split_position_each_dim = []\n",
    "    split_value_each_dim = []\n",
    "    # for each dimension, we calculated the distance from median to its first non-cross split\n",
    "    for D in range(len(dataset[0])):\n",
    "\n",
    "        # sort the dataset in this dimension\n",
    "        dataset = dataset[dataset[:, D].argsort()]\n",
    "\n",
    "        # median\n",
    "        median = dataset[int(total_size / 2), D]\n",
    "        median_low = domains[D][0]\n",
    "        median_up = domains[D][1]\n",
    "\n",
    "        # split position\n",
    "        split_distance = 0\n",
    "        split_position = int(total_size / 2)\n",
    "        split_low = 0\n",
    "        split_up = total_size\n",
    "\n",
    "        # that is the place where we need query bounding !\n",
    "        if len(query) == 0:\n",
    "            split_distance_each_dim.append(split_distance)\n",
    "            split_position_each_dim.append(split_position)\n",
    "            split_value_each_dim.append(median)\n",
    "            continue\n",
    "\n",
    "        query = np.asarray(query)\n",
    "        query_in_this_dim = query[:, D]\n",
    "\n",
    "        # generate bounded intervals in this dim\n",
    "        query_in_this_dim = query_in_this_dim[query_in_this_dim[:, 0].argsort()]\n",
    "        remaining_query = query_in_this_dim\n",
    "        bounded_intervals = []\n",
    "\n",
    "        while len(remaining_query) != 0:\n",
    "\n",
    "            initial_interval = [remaining_query[0][0], remaining_query[0][1]]\n",
    "            temp_interval = []\n",
    "\n",
    "            for i in range(len(remaining_query) - 1):\n",
    "\n",
    "                overlap = getoverlap(initial_interval[0], initial_interval[1], remaining_query[i + 1][0],\n",
    "                                     remaining_query[i + 1][1])\n",
    "                # there is no overlap\n",
    "                if overlap == 0:\n",
    "                    temp_interval.append([remaining_query[i + 1][0], remaining_query[i + 1][1]])\n",
    "                else:  # update interval border\n",
    "                    initial_interval[0] = min(initial_interval[0], remaining_query[i + 1][0])\n",
    "                    initial_interval[1] = max(initial_interval[1], remaining_query[i + 1][1])\n",
    "\n",
    "            bounded_intervals.append(initial_interval)\n",
    "            remaining_query = temp_interval\n",
    "        query_bound = bounded_intervals\n",
    "        \n",
    "        # check if the split position intersect some query boundings in this dim\n",
    "        for i in range(len(query_bound)):\n",
    "\n",
    "            # if intersect some query bounds\n",
    "            if median > query_bound[i][0] and median < query_bound[i][1]:\n",
    "\n",
    "                # check if the two end already exceeds domain\n",
    "                if query_bound[i][0] < domains[D][0] and query_bound[i][1] > domains[D][1]:\n",
    "                    split_distance = int(total_size / 2)\n",
    "                    break;\n",
    "\n",
    "                # if not exceeds, determine which side is closer to the median\n",
    "                else:\n",
    "                    if query_bound[i][0] >= domains[D][0]:\n",
    "                        # get the number of records from median to the end\n",
    "                        for j in range(split_position - 1, -1, -1):\n",
    "                            if dataset[j][D] <= query_bound[i][0]:\n",
    "                                split_low = j\n",
    "                                median_low = dataset[split_low, D]\n",
    "                                break\n",
    "\n",
    "                    if query_bound[i][1] <= domains[D][1]:\n",
    "                        # get the number of records from median to the end\n",
    "                        for j in range(split_position, total_size, 1):\n",
    "                            if dataset[j][D] >= query_bound[i][1]:\n",
    "                                split_up = j\n",
    "                                median_up = dataset[split_up, D]\n",
    "                                break\n",
    "\n",
    "                # if not exceeds then choose the one that is closest from the median (in terms of #records!)\n",
    "                if (total_size / 2) - split_low < (split_up - total_size / 2) and split_low != 0:\n",
    "                    split_position = split_low\n",
    "                    median = median_low\n",
    "                    split_distance = int(total_size / 2) - split_position\n",
    "                elif (total_size / 2) - split_low >= (split_up - total_size / 2) and split_up != total_size:\n",
    "                    split_position = split_up\n",
    "                    median = median_up\n",
    "                    split_distance = split_position - int(total_size / 2)\n",
    "\n",
    "        split_distance_each_dim.append(split_distance)\n",
    "        split_position_each_dim.append(split_position)\n",
    "        split_value_each_dim.append(median)\n",
    "\n",
    "    # aftern calculating the distance from median to its first non-cross split\n",
    "\n",
    "    split_distance_each_dim = np.asarray(split_distance_each_dim)\n",
    "    \n",
    "    split_dimension = 0\n",
    "    split_position = int(total_size / 2)\n",
    "    split_value = 0\n",
    "    \n",
    "    # degradation mechansim (if no valid split position, then using round robin)\n",
    "    # if the median do not cross any historical query, split round robin to enhance robustness\n",
    "    if min(split_distance_each_dim) == int(total_size / 2) or max(split_distance_each_dim) == 0:\n",
    "        # using variance is not good\n",
    "#         dim_variance = []\n",
    "#         for D in range(len(dataset[0])):\n",
    "#            dim_variance.append(np.var(dataset[:,D]))\n",
    "#         split_dimension = np.argmax(dim_variance)\n",
    "#         #split_dimension = random.randrange(len(dataset[0])) # using random is not good\n",
    "#         split_position = int(total_size / 2)\n",
    "#         dataset = dataset[dataset[:, split_dimension].argsort()]\n",
    "#         split_value = dataset[split_position, split_dimension]\n",
    "        split_dimension = current_dim + 1\n",
    "        if split_dimension >= len(domains):\n",
    "            split_dimension %= len(domains)\n",
    "        split_position = int(total_size / 2)\n",
    "        dataset = dataset[dataset[:, split_dimension].argsort()]\n",
    "        split_value = dataset[split_position, split_dimension]\n",
    "    else:\n",
    "        split_dimension = np.argmin(split_distance_each_dim)  # get the split dimension\n",
    "        #print(\"level: \",level,\" split_dimension: \",  split_dimension, \" split_distances: \", split_distance_each_dim)\n",
    "        # split_dimension =  split_distance_each_dim.argsort()[0] # get the split dimension\n",
    "        dataset = dataset[dataset[:, split_dimension].argsort()]\n",
    "        split_position = split_position_each_dim[split_dimension]  # get the split position\n",
    "        split_value = split_value_each_dim[split_dimension]\n",
    "\n",
    "    # split the dataset according to the split position\n",
    "    sub_dataset1 = dataset[0:split_position, :]\n",
    "    sub_dataset2 = dataset[split_position:-1, :]\n",
    "\n",
    "    # change the domains\n",
    "    sub_domains1 = np.copy(domains)\n",
    "    sub_domains1[split_dimension][1] = split_value\n",
    "    sub_domains2 = np.copy(domains)\n",
    "    sub_domains2[split_dimension][0] = split_value\n",
    "\n",
    "    # filter the queries for each sub node\n",
    "    sub_query1 = []\n",
    "    sub_query2 = []\n",
    "    for i in range(len(query)):\n",
    "        if query[i][split_dimension][0] < split_value:\n",
    "            sub_query1.append(query[i])\n",
    "        if query[i][split_dimension][1] > split_value:\n",
    "            sub_query2.append(query[i])\n",
    "\n",
    "    # used to see the current depth\n",
    "    level += 1\n",
    "\n",
    "    # recursion\n",
    "    kdnodes = []\n",
    "    kdnodes.extend(LearnedResuriveDivide(sub_dataset1, sub_query1, sub_domains1, threshold, level, split_dimension))\n",
    "    kdnodes.extend(LearnedResuriveDivide(sub_dataset2, sub_query2, sub_domains2, threshold, level, split_dimension))\n",
    "\n",
    "    return kdnodes\n",
    "\n",
    "start_time = time.time()\n",
    "kdnodes = LearnedResuriveDivide(dataset, training_set, domains_, block_size, 0)\n",
    "end_time = time.time()\n",
    "print(\"training time for Learned KD-Tree(s): \", end_time-start_time)\n",
    "print('learned KD-Tree leaf nodes: ',len(kdnodes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 922,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "training time for fast Learned KD-Tree(s):  9.001773834228516\n",
      "fast learned KD-Tree leaf nodes:  513\n"
     ]
    }
   ],
   "source": [
    "# === Learned KD-Tree Split (Fast Version) ===\n",
    "#\n",
    "# asssumption: the query boundings will not overlap. divide the KD-Tree recursively\n",
    "#\n",
    "# @dataset: contains the data only in this subnode, ordered in original load order;\n",
    "# @query: contains all the queries (bounded or original); ordered in original generated order;\n",
    "# @domains: the current domain of the node of every dimension [first lower, second upper],[]...; array object\n",
    "# @threshold: maximum page size\n",
    "# @level: the current tree depth\n",
    "#\n",
    "# return @kdnodes: contains the domain of each node and the correpsonding records amount, notice the domain is\n",
    "# ordered by the original load order as dataset\n",
    "def FastLearnedResuriveDivide(dataset, query, domains, threshold, level, current_dim = 0):\n",
    "\n",
    "    # check if the threshold is already satisfied\n",
    "    total_size = len(dataset)\n",
    "    if total_size <= threshold:\n",
    "        kdnodes = []\n",
    "        kdnodes.append([domains, total_size])\n",
    "        return kdnodes\n",
    "    \n",
    "    split_distance_each_dim = []\n",
    "    split_position_each_dim = []\n",
    "    split_value_each_dim = []\n",
    "    \n",
    "    # for each dimension, we calculated the distance from median to its first non-cross split\n",
    "    for D in range(len(dataset[0])):\n",
    "\n",
    "        # median, with fast median algorithm\n",
    "        median = np.median(dataset[:,D])\n",
    "        median_low = domains[D][0]\n",
    "        median_up = domains[D][1]    \n",
    "\n",
    "        # split position\n",
    "        split_distance = 0\n",
    "        split_position = int(total_size / 2)\n",
    "\n",
    "        # that is the place where we need query bounding !\n",
    "        if len(query) == 0:\n",
    "            split_distance_each_dim.append(split_distance)\n",
    "            split_position_each_dim.append(split_position)\n",
    "            split_value_each_dim.append(median)\n",
    "            continue\n",
    "\n",
    "        query = np.asarray(query)\n",
    "        query_in_this_dim = query[:, D]\n",
    "        \n",
    "        # bound the projected queries in this dimension\n",
    "        query_bound = bounding_union(query_in_this_dim)\n",
    "        \n",
    "        # check if the split position intersect some query boundings in this dim\n",
    "        for i in range(len(query_bound)):\n",
    "\n",
    "            # if intersect some query bounds (only possible to intersect one bounded query)\n",
    "            if median > query_bound[i][0] and median < query_bound[i][1]:\n",
    "\n",
    "                # check if the two end already exceeds domain, if yes, split from the middle\n",
    "                if query_bound[i][0] <= domains[D][0] and query_bound[i][1] >= domains[D][1]:\n",
    "                    split_distance = int(total_size / 2)\n",
    "                    break;\n",
    "                # if not exceeds, determine which side is closer to the median\n",
    "                else:\n",
    "                    # for the left side\n",
    "                    if query_bound[i][0] > domains[D][0]:\n",
    "                        median_low = query_bound[i][0]         \n",
    "                    # for the right side\n",
    "                    if query_bound[i][1] < domains[D][1]:\n",
    "                        median_up = query_bound[i][1]\n",
    "\n",
    "                # if not exceeds then choose the one that is closest from the median (in terms of #records!)\n",
    "                number_of_records_from_low_to_median = len(dataset[(dataset[:,D]>=median_low) & (dataset[:,D] < median)])\n",
    "                number_of_records_from_up_to_median = len(dataset[(dataset[:,D]<=median_up) & (dataset[:,D] > median)])\n",
    "                \n",
    "                if number_of_records_from_low_to_median <= number_of_records_from_up_to_median:\n",
    "                    median = median_low\n",
    "                    split_distance = number_of_records_from_low_to_median\n",
    "                else:\n",
    "                    median = median_up\n",
    "                    split_distance = number_of_records_from_up_to_median\n",
    "                    \n",
    "        # for each dimension, record its result\n",
    "        split_distance_each_dim.append(split_distance)\n",
    "        split_value_each_dim.append(median)\n",
    "\n",
    "    # aftern calculating the distance from median to its first non-cross split\n",
    "    split_distance_each_dim = np.asarray(split_distance_each_dim)\n",
    "    split_dimension = 0\n",
    "    split_value = 0\n",
    "    \n",
    "    # degradation mechansim (if no valid split position, then using round robin)\n",
    "    # if the median do not cross any historical query, split round robin to enhance robustness\n",
    "    if min(split_distance_each_dim) >= int((total_size / 2)-10) or max(split_distance_each_dim) <= 10:\n",
    "        split_dimension = current_dim + 1\n",
    "        if split_dimension >= len(domains):\n",
    "            split_dimension %= len(domains)\n",
    "        split_value = np.median(dataset[:,split_dimension])\n",
    "    else:\n",
    "        split_dimension = np.argmin(split_distance_each_dim)  # get the split dimension\n",
    "        split_value = split_value_each_dim[split_dimension]\n",
    "\n",
    "    # split the dataset according to the split position\n",
    "    sub_dataset1 = dataset[dataset[:,split_dimension] <= split_value]\n",
    "    sub_dataset2 = dataset[dataset[:,split_dimension] > split_value]\n",
    "\n",
    "    # change the domains\n",
    "    sub_domains1 = np.copy(domains)\n",
    "    sub_domains1[split_dimension][1] = split_value\n",
    "    sub_domains2 = np.copy(domains)\n",
    "    sub_domains2[split_dimension][0] = split_value\n",
    "\n",
    "    # filter the queries for each sub node\n",
    "    sub_query1 = query[query[:,split_dimension,0] < split_value]\n",
    "    sub_query2 = query[query[:,split_dimension,1] > split_value]\n",
    "\n",
    "    # used to see the current depth\n",
    "    level += 1\n",
    "\n",
    "    # recursion\n",
    "    kdnodes = []\n",
    "    kdnodes.extend(FastLearnedResuriveDivide(sub_dataset1, sub_query1, sub_domains1, threshold, level, split_dimension))\n",
    "    kdnodes.extend(FastLearnedResuriveDivide(sub_dataset2, sub_query2, sub_domains2, threshold, level, split_dimension))\n",
    "\n",
    "    return kdnodes\n",
    "\n",
    "start_time = time.time()\n",
    "fast_kdnodes = FastLearnedResuriveDivide(dataset, training_set, domains_, block_size, 0)\n",
    "end_time = time.time()\n",
    "print(\"training time for fast Learned KD-Tree(s): \", end_time-start_time)\n",
    "print('fast learned KD-Tree leaf nodes: ',len(fast_kdnodes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 813,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# === Performance Evaluation ===\n",
    "#\n",
    "# evaluate the blocks of data to be fetched, when physical data of kdnodes are seperate !!!\n",
    "# @queries: a collection of queries contains the lower and upper value in all dimensions; numpy object\n",
    "# @kdnodes: the kdnodes generated above; array object\n",
    "\n",
    "def Query(queries, kdnodes):\n",
    "    \n",
    "    counts = []\n",
    "    count_single_query = 0;\n",
    "    \n",
    "    # number of dimensions\n",
    "    dims = int(len(queries[0]))\n",
    "    \n",
    "    # for each query\n",
    "    for i in range(len(queries)):\n",
    "        \n",
    "        count_single_query = 0\n",
    "        \n",
    "        # check for intersection for each kdnode\n",
    "        for j in range(len(kdnodes)):\n",
    "            \n",
    "            # for each dimension\n",
    "            intersection_tag = True\n",
    "            for k in range(dims):\n",
    "                \n",
    "                # an intersection holds if it intersecs in all dimensions\n",
    "                if queries[i][k][0] >= kdnodes[j][0][k][1] or queries[i][k][1] <= kdnodes[j][0][k][0]:\n",
    "                    intersection_tag = False\n",
    "                    break\n",
    "                \n",
    "            # if the query intersect with this kdnode\n",
    "            if intersection_tag:\n",
    "                count_single_query += 1\n",
    "            \n",
    "        counts.append(count_single_query)\n",
    "    \n",
    "    #print(\"blocks IO: \", counts)\n",
    "    print(\"blocks IO(average): \", statistics.mean(counts))\n",
    "\n",
    "# === Performance Evaluation (dense) ===\n",
    "#\n",
    "# evaluate the blocks of data to be fetched, when physical data of kdnodes are dense !!!\n",
    "# @queries: a collection of queries contains the lower and upper value in all dimensions; numpy object\n",
    "# @kdnodes: the kdnodes generated above; array object\n",
    "#\n",
    "def QueryDense(queries, dense_kdnodes):\n",
    "    \n",
    "    counts = []\n",
    "    count_single_query = 0;\n",
    "    \n",
    "    # number of dimensions\n",
    "    dims = int(len(queries[0]))\n",
    "    \n",
    "    # for each query\n",
    "    for i in range(len(queries)):\n",
    "        \n",
    "        pages = []\n",
    "        \n",
    "        # check for intersection for each kdnode\n",
    "        for j in range(len(dense_kdnodes)):\n",
    "            \n",
    "            # for each dimension\n",
    "            intersection_tag = True\n",
    "            for k in range(dims):\n",
    "                \n",
    "                # an intersection holds if it intersecs in all dimensions\n",
    "                if queries[i][k][0] >= dense_kdnodes[j][0][k][1] or queries[i][k][1] <= dense_kdnodes[j][0][k][0]:\n",
    "                    intersection_tag = False\n",
    "                    break\n",
    "                \n",
    "            # if the query intersect with this kdnode\n",
    "            if intersection_tag:\n",
    "                pages.extend(dense_kdnodes[j][1]) # remember to remove repeated\n",
    "        \n",
    "        counts.append(len(set(pages)))\n",
    "    \n",
    "    #print(\"blocks IO: \", counts)\n",
    "    print(\"blocks IO(average): \", statistics.mean(counts))    \n",
    "\n",
    "# === Physical Storage ===\n",
    "#\n",
    "# Reform the kdnodes from seperate physcial stroage to dense, continous physical storage\n",
    "#\n",
    "# @kdnodes: the previous generated kdnodes, already in order; array\n",
    "# @threshold: the maximum page size \n",
    "#\n",
    "# return @dense_kdnodes: array\n",
    "\n",
    "def DenseKDNodes(kdnodes, threshold):\n",
    "    \n",
    "    dense_kdnodes = []\n",
    "    previous_records = 0\n",
    "    current_records = 0\n",
    "    page_count = 0\n",
    "    \n",
    "    for i in range(len(kdnodes)):\n",
    "        previous_records = current_records\n",
    "        current_records += kdnodes[i][1]\n",
    "        if current_records > threshold:\n",
    "            # determine how many pages exceeds\n",
    "            remaining = kdnodes[i][1] - (threshold - previous_records)\n",
    "            num_pages = math.ceil(remaining / threshold) # num of new pages required\n",
    "            pages = [i+page_count for i in range(num_pages+1)]\n",
    "            page_count += num_pages\n",
    "            current_records = remaining % threshold\n",
    "            dense_kdn = [kdnodes[i][0],pages]\n",
    "            dense_kdnodes.append(dense_kdn)\n",
    "        else:\n",
    "            dense_kdn = [kdnodes[i][0],[page_count]]\n",
    "            dense_kdnodes.append(dense_kdn)\n",
    "                         \n",
    "    return dense_kdnodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 814,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blocks IO(average):  6.889447236180905\n",
      "blocks IO(average):  9.025125628140703\n"
     ]
    }
   ],
   "source": [
    "# Example of usage (Performance Evaluation Learned)\n",
    "Query(testing_set, kdnodes)\n",
    "\n",
    "# Example of usage (Physical Storage)\n",
    "dense_kdnodes = DenseKDNodes(kdnodes, block_size)\n",
    "\n",
    "# Example of usage (Performance Evaluation (dense))\n",
    "QueryDense(testing_set, dense_kdnodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 815,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "traditional KD-Tree leaf nodes:  512\n",
      "traditional variance KD-Tree leaf nodes:  512\n"
     ]
    }
   ],
   "source": [
    "# === Compare With Existing === \n",
    "#\n",
    "# TraditionalKDTree (sparse), split by a round robin order\n",
    "# The parameters are just like the learned kdtree\n",
    "#\n",
    "# return: @kdnodes the traditional kd-tree nodes\n",
    "#\n",
    "def TraditionalKDTree(dataset, currentDim, domains, threshold, level):\n",
    "    \n",
    "    #print(\"level: \",level)\n",
    "    #print(\"dataset size: \",len(dataset))\n",
    "    \n",
    "    # check if the threshold is already satisfied\n",
    "    total_size = len(dataset)\n",
    "    if total_size <= threshold:\n",
    "        # the kdnodes should be an global object outside the function\n",
    "        kdnodes = []\n",
    "        kdnodes.append([domains,total_size])\n",
    "        return kdnodes\n",
    "    \n",
    "    # the current dimension\n",
    "    # divideDim = Dimorder[currentDim]\n",
    "    divideDim = currentDim # the dataset is already ordered as the Dimorder\n",
    "    \n",
    "    # sort according to the current dimension\n",
    "    dataset = dataset[dataset[:,divideDim].argsort()]\n",
    "    \n",
    "    # find the medium\n",
    "    medium = dataset[int(total_size/2),divideDim]\n",
    "\n",
    "    # start check split position from the medium\n",
    "    split_position = int(total_size/2)\n",
    "            \n",
    "    # split the dataset according to the split position\n",
    "    sub_dataset1 = dataset[0:split_position,:]\n",
    "    sub_dataset2 = dataset[split_position:-1,:]\n",
    "    \n",
    "    # change the domains\n",
    "    sub_domains1 = np.copy(domains)\n",
    "    sub_domains1[divideDim][1] = medium\n",
    "    sub_domains2 = np.copy(domains)\n",
    "    sub_domains2[divideDim][0] = medium\n",
    "    \n",
    "    # change the divideDim\n",
    "    currentDim += 1\n",
    "    if currentDim >= len(domains):\n",
    "        currentDim %= len(domains)\n",
    "    \n",
    "    # used to see the current depth\n",
    "    level += 1\n",
    "    \n",
    "    # recursion\n",
    "    kdnodes = []\n",
    "    kdnodes.extend(TraditionalKDTree(sub_dataset1, currentDim, sub_domains1, threshold, level))\n",
    "    kdnodes.extend(TraditionalKDTree(sub_dataset2, currentDim, sub_domains2, threshold, level))\n",
    "    \n",
    "    return kdnodes\n",
    "\n",
    "\n",
    "def TraditionalVarianceKDTree(dataset, domains, threshold, level):\n",
    "    \n",
    "    #print(\"level: \",level)\n",
    "    #print(\"dataset size: \",len(dataset))\n",
    "    \n",
    "    # check if the threshold is already satisfied\n",
    "    total_size = len(dataset)\n",
    "    if total_size <= threshold:\n",
    "        # the kdnodes should be an global object outside the function\n",
    "        kdnodes = []\n",
    "        kdnodes.append([domains,total_size])\n",
    "        return kdnodes\n",
    "    \n",
    "    # calculate the min-max scaling\n",
    "    min_max_scaler = preprocessing.MinMaxScaler()\n",
    "    scaled_dataset = min_max_scaler.fit_transform(dataset)\n",
    "    \n",
    "    # calculate the variance\n",
    "    variances = []\n",
    "    for i in range(len(scaled_dataset[0])):\n",
    "        variances.append(np.var(scaled_dataset[:,i]))\n",
    "        \n",
    "    # use the max normalized variance dimension to split\n",
    "    divideDim = np.argmax(variances)\n",
    "    \n",
    "    # sort according to the current dimension\n",
    "    dataset = dataset[dataset[:,divideDim].argsort()]\n",
    "    \n",
    "    # find the medium\n",
    "    medium = dataset[int(total_size/2),divideDim]\n",
    "\n",
    "    # start check split position from the medium\n",
    "    split_position = int(total_size/2)\n",
    "            \n",
    "    # split the dataset according to the split position\n",
    "    sub_dataset1 = dataset[0:split_position,:]\n",
    "    sub_dataset2 = dataset[split_position:-1,:]\n",
    "    \n",
    "    # change the domains\n",
    "    sub_domains1 = np.copy(domains)\n",
    "    sub_domains1[divideDim][1] = medium\n",
    "    sub_domains2 = np.copy(domains)\n",
    "    sub_domains2[divideDim][0] = medium\n",
    "    \n",
    "    # used to see the current depth\n",
    "    level += 1\n",
    "    \n",
    "    # recursion\n",
    "    kdnodes = []\n",
    "    kdnodes.extend(TraditionalVarianceKDTree(sub_dataset1, sub_domains1, threshold, level))\n",
    "    kdnodes.extend(TraditionalVarianceKDTree(sub_dataset2, sub_domains2, threshold, level))\n",
    "    \n",
    "    return kdnodes\n",
    "\n",
    "traditional_kdnodes = TraditionalKDTree(dataset, 0, domains_, block_size, 0)\n",
    "traditional_variance_kdnodes = TraditionalVarianceKDTree(dataset, domains_, block_size, 0)\n",
    "print('traditional KD-Tree leaf nodes: ', len(traditional_kdnodes))\n",
    "print('traditional variance KD-Tree leaf nodes: ', len(traditional_variance_kdnodes))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 816,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# variance after standard score will be 1\n",
    "\n",
    "# a = [1,2,3,4,5]\n",
    "# a = np.asarray(a)\n",
    "# print(np.var(a))\n",
    "# b = stats.zscore(a)\n",
    "# print(b)\n",
    "# print(np.var(b))\n",
    "\n",
    "# c = [20000,4,500000,8,10000000,19,27999,10000000000000]\n",
    "# print(np.var(c))\n",
    "# d = stats.zscore(c)\n",
    "# print(d)\n",
    "# print(np.var(d))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 817,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blocks IO(average):  8.859296482412061\n",
      "blocks IO(average):  14.522613065326633\n",
      "=========\n",
      "blocks IO(average):  13.125628140703517\n",
      "blocks IO(average):  16.22613065326633\n"
     ]
    }
   ],
   "source": [
    "Query(testing_set, traditional_kdnodes)\n",
    "dense_traditional_kdnodes = DenseKDNodes(traditional_kdnodes, block_size)\n",
    "QueryDense(testing_set, dense_traditional_kdnodes)\n",
    "print(\"=========\")\n",
    "Query(testing_set, traditional_variance_kdnodes)\n",
    "dense_traditional_variance_kdnodes = DenseKDNodes(traditional_variance_kdnodes, block_size)\n",
    "QueryDense(testing_set, dense_traditional_variance_kdnodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 818,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# === KDTree Visualization ===\n",
    "#\n",
    "# restricted to 2D\n",
    "#\n",
    "# @kdnodes: the generated kdnodes, domains in them are ordered according to the learned split dimension order,\n",
    "#   should only generated on the domains for visualization?\n",
    "# @query: the query in the query dimension order, should be the same as the order of 'kdnode domains'\n",
    "# @domains: the domains for visualization, should be in the origin order\n",
    "# @realtive_order: whether to adjust the order of domains in kdnodes to corresponding the 'domains' parameteer\n",
    "#\n",
    "def visualize_kdnodes_and_query(kdnodes, query, domains):\n",
    "    \n",
    "    fig, ax = plt.subplots(1)\n",
    "    \n",
    "    plt.xlim(domains[0][0],domains[0][1])\n",
    "    plt.ylim(domains[1][0],domains[1][1])\n",
    "    # transform the queries into boxes\n",
    "    \n",
    "    # first plot all the kdnodes\n",
    "    for i in range(len(kdnodes)):\n",
    "        \n",
    "        kd_domains = kdnodes[i][0]\n",
    "        kd_domains = np.asarray(kd_domains)\n",
    "        \n",
    "        lower1 = kd_domains[0][0]\n",
    "        upper1 = kd_domains[0][1]\n",
    "        lower2 = kd_domains[1][0]\n",
    "        upper2 = kd_domains[1][1]\n",
    "        \n",
    "        rect = Rectangle((lower1,lower2),upper1-lower1,upper2-lower2,fill=False,edgecolor='b',linewidth=1)\n",
    "        ax.add_patch(rect)\n",
    "    \n",
    "    # plot the queries\n",
    "    for i in range(len(query)):\n",
    "        \n",
    "        lower1 = query[i][0][0]\n",
    "        upper1 = query[i][0][1]\n",
    "        lower2 = query[i][1][0]\n",
    "        upper2 = query[i][1][1]\n",
    "        \n",
    "        rect = Rectangle((lower1,lower2),upper1-lower1,upper2-lower2,fill=False,edgecolor='r',linewidth=3)\n",
    "        ax.add_patch(rect)\n",
    "        \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 923,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAD8CAYAAABU4IIeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztnX/QJVWZ3z/fBWVdf4RBhBoZ2MHU\nuLtgpRAmyJaJO3EjDMQSTWkCqZIppTLqQkUTUxF0K/MiSRW6urpkLXRcJ0DigmypyxQlwZGV0qqI\nMKOEHyI7A6K8zhQjOyqk3GID++SPPnem3/veH933nnP7nHufT1W/b9/Tp59++pzT5zk/nj4tM8Nx\nHMdxYvFrXSvgOI7jzBduWBzHcZyouGFxHMdxouKGxXEcx4mKGxbHcRwnKm5YHMdxnKiMNSySTpb0\nTUkPS3pI0vtD+HGSdknaG/6vCeGSdK2kfZLul3RmTdaWEH+vpC218LMkPRDOuVaSUtys4ziOk54m\nPZbngA+a2e8A5wCXSToNuAK408w2AHeG3wDnAxvCthW4DipDBGwDXgecDWzrGaMQZ2vtvM3T35rj\nOI7TBWMNi5kdMLPvhf1ngIeBk4ALgRtCtBuAt4b9C4EbreJu4FhJa4HzgF1mdsjMfg7sAjaHYy8z\ns+9Y9bbmjTVZjuM4TmEc3SaypPXAa4HvAiea2QGojI+kE0K0k4Anaqcth7BR4csDwgddfytVzwbp\nxWeZ/XYb9R3HcRz2PGVmr0h5hcaGRdJLgC8DHzCzp0dMgww6YBOErw402w5sr/TZaGa7x6ndCgli\nrnATW14OMpvEjRVnXPymYZPq0UR+2/MnJUW+p5DfdflMJTNFmW5LvDzSj6eXMppGXmGSXkBlVL5o\nZl8JwU+GYSzC/4MhfBk4uXb6OmD/mPB1A8Idx3GcAmniFSbgC8DDZvbHtUM7gZ5n1xbg1lr4JcE7\n7Bzgl2HI7A7gXElrwqT9ucAd4dgzks4J17qkJis/pGab4zjOgtJkKOz1wDuBByTdF8I+DFwD3CLp\nUuAnwDvCsa8BFwD7gF8B7wIws0OSrgbuDfE+amaHwv77gOuBFwG3h81xHCc+UjXWPnAQ3ld7j4FK\nXTa/szmWpr0RsyzGhmPL9DkWn2OZhZykMkc9w0MuNGdzLHvMbOP0kobTyivMcRxnFbWKekVPoNBG\nqzM9blja4g+L04KY022pp+4mkW8Dx5MmlzcLWSmf4BzzqAuKNiwpEjm2zBJ0dNKRoh2S1XDTmLKY\n5fBawucn99cVenJTU7RhyebhmpG8lDIdx3Fi4asbO47jOFFxw+I4juNEpeihMCcePhzmRMUs6VzG\nVCR6FcA5ghsWB2j3HkvMeE7heO3sDMANi5OEJi+TLQKp7rMUb8NYMnPQrUn8ErxKZ4EbFqczSqp0\nJ6WUBn3OQ0OL7K3p7sYdkEMrZtbyUsnsgrZLr3Qt11kwmhSOXK1xxxRtWBatFZNLC8Yr4+aUlFYl\n6erkTdGGxekGnz9pTq4N2pyHvvrx8lQe/h6L4ziOExXvsThOQnJubeesWz/dzH026NLV5Ph86hHc\nsDhOQnIdbiptKCz379qU8i2WnqzUuGHJjRG5fvhbF6XUCI7jLCRFG5Z5dA1uYjJK7R7PAkPjlxKZ\noWHOOa9y1s0pm7GT95J2SDoo6cFa2Jck3Re2xyXdF8LXS/rb2rHP1s45S9IDkvZJulaqirWk4yTt\nkrQ3/F/TVHmzeFtseatk/uZ6DA3efnP9inPG3vcgGVPq6aQhdpnKubyn1NUpiyZeYdcDm+sBZvav\nzewMMzsD+DLwldrhR3vHzOy9tfDrgK3AhrD1ZF4B3GlmG4A7w2/HWYW0chsW5jhOt4w1LGb2LeDQ\noGOh1/GvgJtGyZC0FniZmX3HzAy4EXhrOHwhcEPYv6EW7jgraNLydhyne6adY/mnwJNmtrcWdqqk\n7wNPA39oZt8GTgKWa3GWQxjAiWZ2AMDMDkg6YUqdsuJIK/rx4ZF+zOF5gSZ1o0Ks+vfGvbU+O/rT\nelTa55wvOevWT0xdS12jrqT8mtawXMzK3soB4BQz+xtJZwF/Kel0Bk+ntm5fStpKNZwGnNJa2S5o\n3IoeVGrMVoZbz6D0zpngOkMu3abQxlrlNZcHpa1hGJTWg9xBS3HpLUHPEnSEMvTM2t1Y0tHAvwTO\n6oWZ2bPAs2F/j6RHgVdT9VDW1U5fB+wP+09KWht6K2uBg8OuaWbbge3V9Tdmnn0JSFgimjwMTR+a\nNg9XG5mx6JfV1FDE1sNx5pVpeiz/HPihmR0e4pL0CuCQmT0v6VVUk/SPmdkhSc9IOgf4LnAJ8N/C\naTuBLcA14f+tU+jkLDjCxvYcrO9YUn0KMUQl6FmCjkAjRdV+wKYoxhoWSTcBm4DjJS0D28zsC8BF\nrJ60fwPwUUnPAc8D7zWz3sT/+6g8zF4E3B42qAzKLZIuBX4CvKOp8vP0HsugYiY1Hy8s5qFbMHIf\nFoFyhm9y1xFaPO8d3sss6gpZCbk1AGmjme1OIHf6TJ9IxpjB/JQPVuwhrlRDYU2GrNrOdTSVO0l4\nSZVh7nqWoCP0GoN5f8dF0h4z25jyGkW/eT9XlPDUOK0ppSdZgp4l6NiUebqXQbhhcZyElNBeKKE3\nUIKOEAxGA0W7vJVZGDX/HovjOI4TFTcsjuM4TlSKHgrL+Q3aHFZKdrqnlDwrQc8SdIRy9ExJ0YYl\nxphrirHbkj7605PnpKGUeYHc9SxBRyhDz6zfvHcWE++JOY4zDjcsXTGqNs24yZNi+ZVZL+kyS0rR\nuwQ9S9ARytEzJW5YHCchGbcRgHyHbkoZoi5FZr/81LhXmOM4jhOVonsssSxvF/MGoxokJS0575SP\nlyUnNkUblhy7oI1lTDHFksIrrIQFPUukhHTIdSjMKZeiDUvRhKc5lzHu2BPobRahnGdyyNtRzHv6\nO93ghsVxHGcK3Divxg2L4ySkhEonVx1LeWcqxZB8rnnSFDcsMVH4LlzTQmFWfbuhRXynLHLOslyG\nYQeRYh4xxb2mMgCluxsXbVhya9HM8hktvUXjOKMY2+DK1SI6QOGGpZMvPY4UGFHWGHytsDLIPW1z\n1c/NRtkUbVgcJwbDKlcJtrF0+PdVYX9U/H5ybljnPBQ2y0aaE5+x37yXtAN4M3DQzF4TwpaAfwv8\nLET7sJl9LRy7ErgUeB74d2Z2RwjfDPwJcBTwZ2Z2TQg/FbgZOA74HvBOM/u7sYpH+OZ9DssxTKJD\nV+PPnX3zfuzbpkcExPrm/cBrm8W/tw7JWr8Wed5UXKo5lq7rkPby8/jm/fXAnwI39oV/ysw+UQ+Q\ndBpwEXA68ErgG5JeHQ5/BngTsAzcK2mnmf0A+FiQdbOkz1IZpesmvB/H6YaBFaFlO9TUo6l+1qAL\noagDWOPeEm4vMefvN81C5iwZa1jM7FuS1jeUdyFws5k9C/xI0j7g7HBsn5k9BiDpZuBCSQ8DbwT+\nTYhzA7CEGxZnTojV8kzZE27Waxwvp01PsS0leImVJDM108yxXC7pEmA38EEz+zlwEnB3Lc5yCAN4\noi/8dcDLgV+Y2XMD4q9C0lZga/XrlLn50uMkOnS1BEvseE3ith0pGySv9Bag45TEpIblOuBqqmf+\nauCTwLsZ3K4xBq+ibCPiD8TMtgPboTfH0kzZruYGSmrBDJU5wXdjoqdNy+H2tvk4kt/7vWaCHMc5\nzESGxcye7O1L+jxwW/i5DJxci7oO2B/2B4U/BRwr6ejQa6nHd2ZJzi5NQybNZzL5fNddiS/gOPPH\nRN9jkbS29vNtwINhfydwkaRjgrfXBuAe4F5gg6RTJb2QaoJ/p1V9jm8Cbw/nbwFunUQnx3EcJw/G\n9lgk3QRsAo6XtAxsAzZJOoNq2Opx4D0AZvaQpFuAHwDPAZeZ2fNBzuXAHVTuxjvM7KFwiQ8BN0v6\nL8D3gS9Eu7txrGqlZ9JCd5xh9JfZWfUqc+m9OkUw9j2WXGnzHsvQIZO+h1TY4s6xDJrIMMtjjmVI\n/GEyUr5rMjD+mImaWG64ZgOuNUUeTBK3DbEdJkrwCiuDPN5jcRaBQhsYJRAjaVNXWrkvpljKx+hK\naUimxg1LYnJwZ55O5ohSXUwLzRmH91jiyHQq3LAkppQWjC9qWQjes3QKwA2L40xLhHkOx5kn3LA4\nk5Pzuy+O43TGYhuWVW5G3ajhdMssluBxnEViYQzLwMZ1nyWx2l+nbNpU/qk6WbE9mVIawBJW/i3B\nK6wUmalZGMMy+D2WmauxMHT9MEy1VljjteTHXySWq3HKb/yU4hCS60rJqeSW8n7RIBbGsDQldiuv\nlBbMJDJHvXPa5AXJro2P4zhpcMPSR+5v3qeQP/Fb2hO8rN8fd+zFVp7RTLDjOJ1StGHxFm/HLIoH\nWO8+cyhwDXTwuUKna4o2LNG+ueE4JeAF2SmEog3L1PRZpmoxRmcuaFwJN/wCnOM4jVlsw+I4iSnF\nJuXiEDIreSW4V6eQNyvcsDgLz6qV6Cc8bxAlrGwM7m6co1x3N3aGMm0mDvxOCke+8VFqi6YRtacq\n5TDl6oe32ZVaO7U5zoLghiUxU7c4hlROPbfeUlrEWWJW1kuybXy5HadD3LA488m8ukI3uC93QnG6\n5tfGRZC0Q9JBSQ/Wwv5I0g8l3S/pq5KODeHrJf2tpPvC9tnaOWdJekDSPknXSlXTS9JxknZJ2hv+\nr0lxo1FZWmIbS7C01LUmjuM42TH2m/eS3gD8X+BGM3tNCDsX+Csze07SxwDM7EOS1gO39eL1ybkH\neD9wN/A14Fozu13Sx4FDZnaNpCuANWb2obGKx/jm/YTxVgxHjDghylDVsKEPs+7fvE8cZ1z8Ud+8\nb0PKRSh9rbD85PnkfQbfvDezbwWDUQ/7eu3n3cDbR8mQtBZ4mZl9J/y+EXgrcDtwIbApRL0BuAsY\na1g6Zds2lq6CpW1dK7JYDDIYI+zuqngxjFBbYslPvbS/uxvnKbfUKbUYcyzvBr5U+32qpO8DTwN/\naGbfBk4ClmtxlkMYwIlmdgDAzA5IOmHYhSRtBbZWv05Jsix4s3hLAFx1FXBVnOsOY1iDpSe31II3\nCcMMwyAj0oaUaTjrFmzs3mdTclgXL8W998cbdV4p6Zm9u7GkjwDPAV8MQQeAU8zsbySdBfylpNMZ\n7HvTOtnMbDuwvbr2xnGjeBPRVeEYHnd1YIqHrItVnXMxiqkq4lzuz3FmzcSGRdIW4M3A71uo4s3s\nWeDZsL9H0qPAq6l6KOtqp68D9of9JyWtDb2VtcDBSXVyJqfJ0NG0c1XTynQcpwzGeoUNQtJmqnmQ\nt5jZr2rhr5B0VNh/FbABeCwMdT0j6ZzgDXYJcGs4bSewJexvqYU7juM4BTK2xyLpJqrJ9eMlLQPb\ngCuBY4BdwWv4bjN7L/AG4KOSngOeB95rZoeCqPcB1wMvopq0vz2EXwPcIulS4CfAO6LcmeM4jtMJ\nY92Nc6WNu3E7ubnNsUwWd1p5OQ6FTSOzzfltdYgVP4Y8n7yPJ29QvPmYvM/A3dhxUmJoiGtHmQ2e\nfnKeGxqa9gMjz0d+pKYEl+hZULRhKcUX3VlcjnzGeUShalhpRy2XLYX5M9GMMnosaeTWKdqwlOBu\n7Dhj6S8oGfYOmgxDJX0faFTvKsP0WnSKNiyO4+RD0vrdG2lF4YbFcZzxjOmOWO2v47hhcZxZ0O9a\n5DhzjBsWpxVT14lDBKi/tet170DiLZVjQ79OmiPC8ppKKWBerEvcsDitSLX8SpPJYSfuezRSg8o6\no4SPvf5c7is6Z5T0rSnasJTgbpyqoJe0yGPupKyMBsUvqC7PirhGNU0jKffv2vRkpqZow5L6Lehp\nM3Xat8onje+r8LZj5m/ej0pvs7GfSnCc3JloEcri6Tndj9oWiG0sHbnvHD63PO95YbZqE9bJOP3Y\nRwFbtU0sa4rNKYuieyzOgiDhrqxpmMyWDT6pvshAicM3rfDJ+pEsZo/FcRzHScbCrG68ohXVpPlj\ntlBzLFmtbjwgf4a5m06re5u0ip1PbeN2XU5iXXcWMmNO3nexunFamb66seM4OTOgEWCw2kGh0AZs\nW3J3YZ4VbljGUGKmOk5TpinfbV6wXJTnqIweS1x5g3DDMoaJ3k4OJWFRHianXKaqtFqU72krR3+W\nysIn7x3HcZyoNDIsknZIOijpwVrYcZJ2Sdob/q8J4ZJ0raR9ku6XdGbtnC0h/l5JW2rhZ0l6IJxz\nrZS4fTLgPYJVm9Mdnh8zY1bvlvh7LItF0x7L9cDmvrArgDvNbANwZ/gNcD6wIWxbgeugMkTANuB1\nwNnAtp4xCnG21s7rv1ZeeMXnzAlN2lhD214tIk9zHX+8yqORYTGzbwGH+oIvBG4I+zcAb62F32gV\ndwPHSloLnAfsMrNDZvZzYBewORx7mZl9xyrf5xtrshzHcZzCmGby/kQzOwBgZgcknRDCTwKeqMVb\nDmGjwpcHhK9C0laqng1wyhSqT4n3zR1nfqk93/2u04d/ezdqJCm8wgbVujZB+OpAs+3AdqhekJxk\nnHcW1K/T9Jqx7yWGvKZhdQ6vJjUgXi+8f72padNoYPiqQGt3fksdYsRvErfrctIls7r3JiYjRllq\nSkl51GMaw/KkpLWht7IWOBjCl4GTa/HWAftD+Ka+8LtC+LoB8ceS8q3iw+dP7m0cXafYcmLKa5Ne\nk16rlZ59euT8tnSK+KXInCR+8vIfoQxPkp4x3upvet3UTONuvBPoeXZtAW6thV8SvMPOAX4Zhszu\nAM6VtCZM2p8L3BGOPSPpnOANdklNVvfkPNO4tHTEbSaHVYkdx3Fo2GORdBNVb+N4SctU3l3XALdI\nuhT4CfCOEP1rwAXAPuBXwLsAzOyQpKuBe0O8j5pZzyHgfVSeZy8Cbg9bPJqYaB8zjU9/c8xxnIWg\nkWExs4uHHPr9AXENuGyInB3AjgHhu4HXNNHFKQg3Jk6JeCNzauZ2deMVY40T9FhyGHtPLSemvCNz\nLOl6h+3mWJpNsuSQz6XMh+RwX5PInbW8pjK7m2Px1Y07JWcPj9i6xZAnTedR0/QaQOM13A57oI2I\nnl3HapRChTYEcye7MlA4blhG0HVLNpWMtDLH9wwmvdTKXmhDbRK2sEfJcdISO41TPFOLjBsWx3GK\nI/ZQmBOXxTAsvVI4qAT50EL59Oeh1xSO0ymLYVh6uBFxSsDLqVM4i2VYnPkk1x5KQ70qZzo3JkCj\nNLPa37mkoUdjzrhhcZwMaOtkUAol6RoNrV4vr399vHmnaMOSeoG9XN2Nc9VrFtc44m6c5poLWREG\nUiyCOa8Gsy2D0iFV2cwhHYs2LKl6iDm79OasWyp5g+Q3vcCoWKn0zMVIpzYCXS52mON7XD1SFf1Y\nZXUW5bNow+I4i0qKt+RLIuuGUBdpmdl6iNOsbuw4zhBKq6izoUHlt2jzFSXiPRbHSUBVPzarAJsu\nhbMwjO2OzUaNzijQC6wfNyyOM4dY9anOhpHLr8iywiz5POO46wOddpt9KMxxHMeJStE9lpK8TmLK\nzFm3VPJSyfe5kDIprbzOlAxupmjDsqqrGemN1RJcemPKLNHdONe07Mly0pJ7eZ114y/15yraUrRh\ncRynAb5I51ww0vi18DaeRfZPPMci6bck3Vfbnpb0AUlLkn5aC7+gds6VkvZJekTSebXwzSFsn6Qr\npr0px3FqSCs3Z/4wG7/NkIl7LGb2CHAGgKSjgJ8CXwXeBXzKzD5Rjy/pNOAi4HTglcA3JL06HP4M\n8CZgGbhX0k4z+8GkujmO4zjdEWso7PeBR83sxxreIroQuNnMngV+JGkfcHY4ts/MHgOQdHOI296w\nuNuk4zhO58RyN74IuKn2+3JJ90vaIWlNCDsJeKIWZzmEDQtfhaStknZL2g0/i6S648wfosHQSAdD\nJM5iMHWPRdILgbcAV4ag64CrqRwVrgY+CbybwdNLxmDjNrC0m9l2YHt13Y22qO7GKWSW5r6Zc1pO\nIitF/HnWoYTyOkn6jDqnpOmxGENh5wPfM7MnAXr/ASR9Hrgt/FwGTq6dtw7YH/aHhY9kRWOrTar3\ntdJKcI8tQcdUcmfxFnNsd+NSVt7OOd9TyEohN3f390GyUxNjKOxiasNgktbWjr0NeDDs7wQuknSM\npFOBDcA9wL3ABkmnht7PRSGu4ziOUyBT9Vgk/QaVN9d7asEfl3QG1XDW471jZvaQpFuoJuWfAy4z\ns+eDnMuBO4CjgB1m9tA0ejmO4zjdISt08k7aaGa76wHNT/ahsJnJL+HeY1/Dh8LKSM/FHQrTHjPb\nmEZ6hS9C6TiO40TFl3QJlODFVIKOqeTm8rnfLmSllp1zvqeQlUJu7l6Ks6Zow7Iy4a36BsW4c7CB\njs+lDTPlOiRSwr2mkhlbdmlDgKlk5p6eud7vKPmpKdqwrE788bkxKEbJLQPHcZzc8DkWZzb0L4Q4\nbosp38mLTZtW5s+mTV1r5ETGDYvjOI4TFTcsjuM4TlSKnmNxHKdA7rqraw2cxBRtWEpwZ0wpP1e3\n00Fy2jq5tFmcb5z8YeeW4BJcmpt1Kpm5p2eu99sVRRuW3N0kY8mcC7fTlg/KIF1G6jhCfpv7GnaN\nNunTJG6sONOQa3lPJTOH9GyrQ6r0TE3RhsUpiNSWsSb/8MPYe4KaPEmFLm3kODnihsVxuqLf4Llx\nc+YE9wpzHMdxouI9FscZh1Q5B9Q7GN67cJyhuGFxFoueQSjd7cbJm4GukZE8EJrQccOnaMNSgptk\nLJnudrqYuBt8XJkllkeJRgvs9p/TJUUbllzdGl1mxm6SXT9xLenaPdZltpOZgsrDcYJzhuDuxo6T\nA2Zp3oHweRpnTpnaK0zS45IekHSfpN0h7DhJuyTtDf/XhHBJulbSPkn3SzqzJmdLiL9X0pZp9XIc\nx3G6IZa78T8zszNq31G+ArjTzDYAd4bfAOcDG8K2FbgOKkMEbANeB5wNbOsZI8eZGLPmm+PkTJuy\nnEF5TvUey4XADWH/BuCttfAbreJu4FhJa4HzgF1mdsjMfg7sAjYn0s1xHCctGVb2sySGYTHg65L2\nSNoawk40swMA4f8JIfwk4InaucshbFj4CiRtlbS7GnL7WQTVHcdxnNjEmLx/vZntl3QCsEvSD0fE\nHeSPMMznYfWHh822A9sBpI2Ws1ujy+xWx2muHWM15CZxY8WZhhLK0bzJbHvdwhwZgQiGxcz2h/8H\nJX2Vao7kSUlrzexAGOo6GKIvAyfXTl8H7A/hm/rC7xp/7Wm1X0nK1U9TyY4pt6msXFZozSE9u1rJ\nOFf38ByeodgrUXd97dbX37YEV101OMK2bWjIoZhMNRQm6cWSXtrbB84FHgR2Aj3Pri3ArWF/J3BJ\n8A47B/hlGCq7AzhX0powaX9uCOueUd9S92+rO47jrGLaHsuJwFdVVapHA39uZv9L0r3ALZIuBX4C\nvCPE/xpwAbAP+BXwLgAzOyTpauDeEO+jZnZoSt0cx3GcDpAV6q0gbTSz3ZFlDuhutumJjEjLXIZu\nYsjqetgmlcxJ5PpQWBo508he+KGwsWVNe2qvhiTB37x3nH5aNCYOr3o8ZW1qaPSyHYU2AJ3FxA3L\nOPyBdjIlVy8+n250ijYsSQpwQ6Fa7Q0dS3Q2cmNSirtoSSziUJhTBkUbliQPVsPC2/baucwJjJOV\nipLmWBzHmY6iDYtTKG1r71kPR7a4XsoWuuOUihuWUfTXGN6cdRIhzA2UMzekWoTScRzHWVC8xzIK\n76EMxZPGcZxhuGFxJqLty3wpyEXuvBjZEtyNUyxCGltmimunun4qijYsKfz4m9aXk1w7l4qwC+o6\ntp1KGHZ/pcxJNF1KrmuDlrvnXlcrOHT5Jn+b+G2unZqiDUuaiqXZ980nqRxLeHBTFboVekZwCptK\nz0jL9LShhArWcWJRtGFx4jOTFs+Yi6RuebWpk73CdZz2uGFximXiVvuMOyxunJxFww2L4zjlMspq\nlzIJN4e4YXGKZfKeQJtZ00mv4aSgP89H5aT3FLtj/g3LhBO1pRTKUvRMQe7OEHW5i0pyz7UpOyyL\nnDcpKdqwNCkUk07Upqi0UhTi2N5GThpyW+LeVtXI5uXTiUbRhqVRoZ2gwzITl9sI+INWDtm5Gw8o\nO14+nVhMvFaYpJMlfVPSw5IekvT+EL4k6aeS7gvbBbVzrpS0T9Ijks6rhW8OYfskXTHdLTmOM5IB\nNf7qHkwhmA3fnM6YpsfyHPBBM/uepJcCeyTtCsc+ZWafqEeWdBpwEXA68ErgG5JeHQ5/BngTsAzc\nK2mnmf1gCt2O4AXMcZxZIB35VPVY5rtemtiwmNkB4EDYf0bSw8BJI065ELjZzJ4FfiRpH3B2OLbP\nzB4DkHRziBvHsDhZMNIrtO9J7H84J/lap+M43RFl2XxJ64HXAt8NQZdLul/SDklrQthJwBO105ZD\n2LBwpwN6Tgajtqbx6vGnGbEYdk4bHVLcVxu5Thy6yPOm13WOMPXkvaSXAF8GPmBmT0u6DriaquF5\nNfBJ4N0M7iAag43bwOpG0lZga/XrlKSZmUJ2CTJjL4zXi5eCmLq2nRAfFH+YjBQVzzTyht3mopXP\nJnneqly0vNcUvfhcDNxUhkXSC6iMyhfN7CsAZvZk7fjngdvCz2Xg5Nrp64D9YX9Y+ArMbDuwvZK9\n0WKvItqmokxRAbchlUwnDbHyKk6+2+rMNos64OjlczRj02bMvU6TtrNIx2m8wgR8AXjYzP64Fr62\nFu1twINhfydwkaRjJJ0KbADuAe4FNkg6VdILqSb4d06ql+M4DaiNK/ocViR6aTlq3HdBPNam6bG8\nHngn8ICk+0LYh4GLJZ1B1Xt7HHgPgJk9JOkWqkn554DLzOx5AEmXA3cARwE7zOyhKfRy5o0UzV/H\ncZKhRuNJGVINhe1uEM+HwlLJjJZeTfrmfQJi65p6jiWvobDFkDnzOZYJ4o8UNIopLiJpj5ltnFhA\nA4p+896ZD4RVz8k8DaI7zhQcfiYKxQ1LwXg97Dj5Eeu5LPn5dsMyI1IUEve6KYeYaVuCW3AqmSUQ\n6+NwqXoss8gXNywzwo3AYuNzLHFkOmXghsXJn0IGmw0NfP/Aan+dBWFpafD+guBeYRPEda+wdvFj\ne2RNcs4svMJGv0odJ7NyyPc1Oj8vAAAHuUlEQVRSZHbqFVYvCxMkRNqhMPcKczLHhyccZwDbtnWt\nQae4YXGmwj//6jgDWMDhrzrzbVg05PsIHQz/LbrXzThdJ7mXtueUlF7ObGhSJroqZyWX16INy9iX\nUyc4r2lmts30rsebm8pMRVfzXCviDbtBS/O99349cpQ17zJjvsCeYk4mxjUnkZmaog3LpCuEDjsv\n1SRySS2PknSNTQ6GtSldLROTxMGj4ZI+sfVc5LKemqINixOftpXGPLHI9+44MZlvw1KoK7Xj5Mqw\nd3VGneEsHvNtWBynTuKGxrCFAyWvXp3Fwg2L4zjd4iMLIylx2NUNy4zIwUPG6Y65zf9hXbThP+ee\n2PdbolfYxJ8mzgEp7paSJl8rbbMNlbltieoDqQO2bUuL/LXUTplZ/k8oq4R7L4XY+VwiRfdYSrTk\njjNX+EPjDKDoHovjOI6TH9msbixpM/AnwFHAn5nZNaPjN1vduJ0O7XtBXbyNm4vMFKtBJ33zPkK8\nUfGHycghr6LJatlDifmJ3RQvSMYu72W8eZ9+deMseiySjgI+A5wPnAZcLOm0brVyHKcfueO004As\nDAtwNrDPzB4zs78DbgYu7Fgnx3EGsQizz85U5DJ5fxLwRO33MvC6/kiStgJbw89nJT0YW5EEq+we\nDzxVgLvpRHqmWLRzRNzjgadSXTvGvYewFXpOInvS67fksJ6R8idK/AEk0zPyvSR51hOUnd+KLrGP\nXAzLoKRb1dwxs+3AdgBJu1OPE8bA9YxHCTqC6xkb1zMukuJOTg8gl6GwZeDk2u91wP6OdHEcx3Gm\nIBfDci+wQdKpkl4IXATs7Fgnx3EcZwKyGAozs+ckXQ7cQeVuvMPMHhpz2vb0mkXB9YxHCTqC6xkb\n1zMuyfXM5j0Wx3EcZz7IZSjMcRzHmRPcsDiO4zhxMbOiNmAz8AiwD7hihtd9HHgAuA/YHcKOA3YB\ne8P/NSFcwLVBx/uBM2tytoT4e4EttfCzgvx94Vw11GsHcBB4sBaWXK9h12ip5xLw05Cm9wEX1I5d\nGa75CHDeuPwHTgW+G/T5EvDCEH5M+L0vHF8/QseTgW8CDwMPAe/PMT1H6Jlbev46cA/wf4KeV00q\nO5b+LfW8HvhRLT3PyOA5Ogr4PnBbjml5+JwmlVcuW0jUR4FXAS8MBeG0GV37ceD4vrCP9zIAuAL4\nWNi/ALg9FMBzgO/WCtFj4f+asN+rpO4BfjeccztwfkO93gCcycoKO7lew67RUs8l4D8OiHtayNtj\nQqF+NOT90PwHbgEuCvufBd4X9v8A+GzYvwj40ggd1xIqCeClwF8HXbJKzxF65paeAl4S9l9AVTmd\n01Z2TP1b6nk98PYB8bt8jv4D8OccMSxZpeVhPZtUXrlsIWPuqP2+ErhyRtd+nNWG5RFgbdhfCzwS\n9j8HXNwfD7gY+Fwt/HMhbC3ww1r4ingNdFvPygo7uV7DrtFSzyUGV4Qr8pXKW/B3h+V/eFifAo7u\nLye9c8P+0SFe097grcCbck3PAXpmm57AbwDfo1pRo5XsmPq31PN6BhuWTvKd6v2+O4E3ArdNkk+z\nSsvS5lgGLf1y0oyubcDXJe0JS8sAnGhmBwDC/xPG6DkqfHlA+KTMQq9h12jL5ZLul7RD0poJ9Xw5\n8Asze26AnofPCcd/GeKPRNJ64LVUrdds07NPT8gsPSUdJek+qmHQXVSt4rayY+rfSE8z66Xnfw3p\n+SlJx/Tr2VCfWPn+aeA/AX8ffk+ST8nTEsqbvG+09EsiXm9mZ1KtwHyZpDeMiDtMz7bhsclNr+uA\nfwicARwAPhnCY+rZ+h4kvQT4MvABM3t6VNSIerZmgJ7ZpaeZPW9mZ1C1ts8GfmcC2cnTuV9PSa+h\narH/NvCPqYa3PhRZz8ZIejNw0Mz21INHyO0sLaE8w9LZ0i9mtj/8Pwh8leoheVLSWoDw/+AYPUeF\nrxsQPimz0GvYNRpjZk+GB/rvgc9Tpekkej4FHCvp6L7wFbLC8X8AHBqmk6QXUFXWXzSzr4y5187S\nc5CeOaZnDzP7BXAX1ZxEW9kx9W+q52YzO2AVzwL/ncnTM0a+vx54i6THqVZ/fyNVDybPtGwy7pjL\nRjVW+BjVpFNvgun0GVz3xcBLa/v/m8qD4o9YOfH28bD/L1g5uXdPCD+OystkTdh+BBwXjt0b4vYm\n9y5ood96Vs5dJNdr2DVa6rm2tv/vgZvD/umsnGB8jGpycWj+A3/BygnGPwj7l7FyEvOWEfoJuBH4\ndF94Vuk5Qs/c0vMVwLFh/0XAt4E3t5UdU/+Weq6tpfengWsyeY42cWTyPqu0PKxj08orl43KI+Ov\nqcZqPzKja74qJHTPHfEjIfzlVJNpe8P/XiES1YfLHqVyMdxYk/VuKne+fcC7auEbgQfDOX9K8wnm\nm6iGPf4fVavj0lnoNewaLfX8H0GP+6nWhqtXjB8J13yEmofcsPwPeXRP0P8vgGNC+K+H3/vC8VeN\n0PGfUHXz76fmsptbeo7QM7f0/EdUrrH3h3v+z5PKjqV/Sz3/KqTng8D/5IjnWGfPUThnE0cMS1Zp\n2dt8SRfHcRwnKqXNsTiO4ziZ44bFcRzHiYobFsdxHCcqblgcx3GcqLhhcRzHcaLihsVxHMeJihsW\nx3EcJyr/H0jSaZ9ge6ZbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1170a8438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAD8CAYAAABU4IIeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztnX/QJVWZ3z/fBWVdf4RBhBpnYAdT\n4+6ClUKYIFsm7sSNMBBLNKXJkCqZUiqjLlQ0MRVBtzIvklShq6tL1kLHdQIkLsiWukxREhxZKa2K\nyMwo4YfIzoAorzPFyI4KKbfYDPvkjz7vTL/3vT+67z3n9jn3Pp+qft++p08//fQ5p89zfjx9WmaG\n4ziO48Ti17pWwHEcx5kt3LA4juM4UXHD4jiO40TFDYvjOI4TFTcsjuM4TlTcsDiO4zhRGWlYJJ0m\n6ZuSHpH0sKT3h/CTJO2StC/8XxXCJel6SfslPSDpnJqsLSH+PklbauHnSnownHO9JKW4WcdxHCc9\nTXosR4APmtnvAOcDV0g6E7gKuNvM1gN3h98AFwHrw7YVuAEqQwRsA14HnAdsWzJGIc7W2nmbJr81\nx3EcpwtGGhYzO2hm3wv7zwKPAGuAS4CbQrSbgLeG/UuAm63iXuBESauBC4FdZnbYzH4O7AI2hWMv\nM7PvWPW25s01WY7jOE5hHN8msqR1wGuB7wKnmtlBqIyPpFNCtDXAk7XTFkPYsPDFPuH9rr+VqmeD\n9OJzzX67jfqO4zgOe582s1ekvEJjwyLpJcCXgQ+Y2TNDpkH6HbAxwlcGmm0Htlf6bDCzPaPUboUE\nMVe4iS0vB5lN4saKMyp+07Bx9Wgiv+3545Ii31PI77p8ppKZoky3JV4e6ceTSxlOI68wSS+gMipf\nNLOvhOCnwjAW4f+hEL4InFY7fS1wYET42j7hjuM4ToE08QoT8AXgETP749qhncCSZ9cW4PZa+GXB\nO+x84JdhyOwu4AJJq8Kk/QXAXeHYs5LOD9e6rCYrP6Rmm+M4zpzSZCjs9cA7gQcl3R/CPgxcB9wm\n6XLgJ8A7wrGvARcD+4FfAe8CMLPDkq4Fdod4HzWzw2H/fcCNwIuAO8PmOI4TH6kaa+87CO+rvcdA\npS6b39kcS9PeiFkWY8OxZfoci8+xTENOUpnDnuEBF5qxOZa9ZrZhckmDaeUV5jiOs4JaRb2sJ1Bo\no9WZHDcsbfGHxWlBzOm21FN348i3vuNJ48ubhqyUT3COedQFRRuWFIkcW2YJOjrpSNEOyWq4qf2o\nUvtLRB8Kiyirh9xfV1iSm5qiDUs2D9eU5KWU6TiOEwtf3dhxHMeJihsWx3EcJypFD4U58fDhMCcq\nOTu5JHoVwDmGGxYHaPceS8x4TuH0FByvsB1ww+IkosnLZPNAqvvM2dswZ90mkdkkfglepdPADYvT\nGSVVuuOSa+u9pJ7FPHtrurtxB+TQipm2vFQyu6Dt0itdy3XmjCaFoxTrPGWKNizz1orJpQXjlXFz\nck6rnHVzyqZow+J0g8+fNCfXBm1pQ2FOWfh7LI7jOE5UvMfiOAnJubWds269dDP32aBLV5Pj86nH\ncMPiOAnJdbiptKGwbBe0TCw3hXz3CptHhuT60W9dlFIjOI4zlxRtWGbRNbiJySi1ezwNDI1eFn2K\nhjnnvMpZN6dsRk7eS9oh6ZCkh2phX5J0f9iekHR/CF8n6W9rxz5bO+dcSQ9K2i/peqkq1pJOkrRL\n0r7wf1VT5c3ibbHlrZD5m+sw1H/7zXXLzhl53/1kTKink4bYZSrn8p5SV6csmniF3QhsqgeY2b82\ns7PN7Gzgy8BXaocfWzpmZu+thd8AbAXWh21J5lXA3Wa2Hrg7/HacFUjLt0FhjuN0y0jDYmbfAg73\nOxZ6Hf8KuGWYDEmrgZeZ2XfMzICbgbeGw5cAN4X9m2rhjrOMJi1vx3G6Z9I5ln8KPGVm+2phZ0j6\nPvAM8Idm9m1gDbBYi7MYwgBONbODAGZ2UNIpE+qUFcda0U8MjvRjjs4LNKkbFWIN+964k47entGw\nnlLOvaicdeslpq6lrlFXUn5NalguZXlv5SBwupn9jaRzgb+UdBb9p1Nbty8lbaUaTgNOb61sFzRu\nRfcrNWbLw23JoCydM4lmoy8/SdwuVoIdl7aGoV+e9nMHLcWltwQ9S9ARytAza3djSccD/xI4dynM\nzJ4Dngv7eyU9Bryaqoeytnb6WuBA2H9K0urQW1kNHBp0TTPbDmyvrr8h8+xLQMIS0eRhaPrQtHm4\n2siMRa+spoYith6OM6tM0mP558APzezoEJekVwCHzex5Sa+imqR/3MwOS3pW0vnAd4HLgP8WTtsJ\nbAGuC/9vn0AnZ84RNrLnYD3HkupTiCEqQc8SdAQaKar2AzZFMdKwSLoF2AicLGkR2GZmXwA2s3LS\n/g3ARyUdAZ4H3mtmSxP/76PyMHsRcGfYoDIot0m6HPgJ8I6mys/Seyz9ipnUfLywmIduzsh9WATK\nGb7JXUdo8bx3eC/TqCtkJeRWH6QNZrYngdzJM30sGSMG81M+WLGHuFINhTUZsmo719FU7jjhJVWG\nuetZgo6w1BjM+zsukvaa2YaU1yj6zfuZooSnxmlNKT3JEvQsQcemzNK99MMNi+MkpIT2Qgm9gRJ0\nhGAwGija5a1Mw6j591gcx3GcqLhhcRzHcaJS9FBYzm/Q5rBSstM9peRZCXqWoCOUo2dKijYsMcZc\nU4zdlvTRnyV5ThpKmRfIXc8SdIQy9Mz6zXtnPvGemOM4o3DD0hXDatOMmzwpll+Z9pIu06QUvUvQ\nswQdoRw9U+KGxXESknEbAch76CbFkG8Jw97TGEpPjXuFOY7jOFEpuscSy/J2MW8wrEFS0pLzTvl4\nWXJiU7RhybEL2lhGyymW1GuFlbCgZ4mUkA45D4U5ZVK0YSmazJ7m2BPobRahnGUyy+YVzHr6O93g\nhsVxHGcC3DivxA2L4ySkhEonZx1LGKJNMSSfc540wQ1LTBS+C9e0UPR+035Y1NpfpxxyHgrL2dUY\n4uqX6l5TGYDS3Y2LNiy5vQU+zre92pxTeivGcRozqpGWs0V0yjYsnXzpcajAdtHN2p2T0ivMSUPu\naZurfm42yqZow+I4MRhUuUqwjYWjv68J+8Pi95JzwzrrobBMDZ7TjJHfvJe0A3gzcMjMXhPCFoB/\nC/wsRPuwmX0tHLsauBx4Hvh3ZnZXCN8E/AlwHPBnZnZdCD8DuBU4Cfge8E4z+7uRikf45n1JyzHM\n9TfvR75tekxArG/e9722Wfx765Cs9WuR503FlfJcpl/SJY9v3t8I/Clwc0/4p8zsE/UASWcCm4Gz\ngFcC35D06nD4M8CbgEVgt6SdZvYD4GNB1q2SPktllG4Y834cpxv6VoSW7VDTEk31swZdCEUdwBoh\na4x0zfn7TdOQOU1GGhYz+5akdQ3lXQLcambPAT+StB84Lxzbb2aPA0i6FbhE0iPAG4F/E+LcBCzg\nhsWZEbpqzcbvNY6W06an2BZfkDKuzNRMMsdypaTLgD3AB83s58Aa4N5anMUQBvBkT/jrgJcDvzCz\nI33ir0DSVmBr9ev0ufvSY8oC0VR27HhN4rYdKesnr/QWoOOUxLiG5QbgWqpn/lrgk8C76d+uMfqv\nojzIJ2pgPWJm24HtsDTH0kzZrlp5JbVgBsoc47sx0dOm5XB723wcyu/9XjNBjuMcZSzDYmZPLe1L\n+jxwR/i5CJxWi7oWOBD2+4U/DZwo6fjQa6nHd6ZJzi5NAybNpzL5fM89iS/gOLPHWN9jkbS69vNt\nwENhfyewWdIJwdtrPXAfsBtYL+kMSS+kmuDfaVWf45vA28P5W4Dbx9HJcRzHyYORPRZJtwAbgZMl\nLQLbgI2SzqYatnoCeA+AmT0s6TbgB8AR4Aozez7IuRK4i8rdeIeZPRwu8SHgVkn/Bfg+8IVodzeK\nFa30TFrojjMAQ8uHBqfVq8yl9+oUwcj3WHKlzXssA4dMegyLsPmdY+k3kTFqLbNpzbEMiD9IRsp3\nTfrGHzFRE9MNd0U+TZAH9bglUIJXWBnk8R6LMw8U2sAogS6XDWpK7i8PlvIxulIakqlxw5KYUtyZ\nB8scUqqLaaE5s4D3WMrBDUtiSmnBxH5onTS0Ga51nK5ww+I4kxJhnsNxZgk3LM745Pzui+M4nTHf\nhmWFm1E3ajjd0nbozof6HGc4c2NY+jaueyyJf/53dmhT+af8HEGXzhsx12sbl5hyS/AKK0VmaubG\nsPR/j2XqaswNXT8ME73D0Xgt+dEX6cLRIsq7OROS82rEpXybJaWeqZkbw9KU2K28Ulow48gc9s5p\nk8qva+PjOE4a3LD0MLdv3o9z32O8rN8bd+TFlp/RTLDjOJ1StGHxFm/HzIsH2NJ9ZlDgVqwV1jfO\nsb+O0wVFG5ZZWwvJcYbiBdkphKINy8T0WKZqMUZnJmhcCTf8ApzjOI2Zb8PiOImJaZNSNnpycQiZ\nlrwS3KtTyJsWblicuaf34W38xeQGD30JKxtD/k4m7m4cV25q3LAkZtJM7PsSp+J+4yNbak9VymHK\nlQ9vsyu1dmpznDnBDUtiJm5xDKicfDXiCJiV9ZJsG19ux+kQNyzObDKjrtBNls13JxSna35tVARJ\nOyQdkvRQLeyPJP1Q0gOSvirpxBC+TtLfSro/bJ+tnXOupAcl7Zd0vVQ1vSSdJGmXpH3h/6oUNxqV\nhQW2sQALC11r4jiOkx0jv3kv6Q3A/wVuNrPXhLALgL8ysyOSPgZgZh+StA64Yylej5z7gPcD9wJf\nA643szslfRw4bGbXSboKWGVmHxqpeIxv3o8Zb9lwxJAToky+DRr6iNgij/2N+By+ed+GlItQ+lph\n+cnzyfsMvnlvZt8KBqMe9vXaz3uBtw+TIWk18DIz+074fTPwVuBO4BJgY4h6E3APMNKwdMq2bSxc\nAwvbulZkvuhnMJra3VhGqC1dufHm8CmAnF1v3d04LTHmWN4NfKn2+wxJ3weeAf7QzL4NrAEWa3EW\nQxjAqWZ2EMDMDko6ZdCFJG0Ftla/Tk+yLHizeAsAXHMNcE2c6w5iUIOl1AI3CYMMQz8j0oaUaTnt\nFmzs3mdTUn8ts6t774037LxS0jN7d2NJHwGOAF8MQQeB083sbySdC/ylpLPo73vTOtnMbDuwvbr2\nhlGjeGPRVeEYHHdlYIqHrItVnXMxjqkq4lzuz3GmzdiGRdIW4M3A71uo4s3sOeC5sL9X0mPAq6l6\nKGtrp68FDoT9pyStDr2V1cChcXVyxqfJ0NGkc1WTynQcpwxGeoX1Q9ImqnmQt5jZr2rhr5B0XNh/\nFbAeeDwMdT0r6fzgDXYZcHs4bSewJexvqYU7juM4BTKyxyLpFqrJ9ZMlLQLbgKuBE4BdwWv4XjN7\nL/AG4KOSjgDPA+81s8NB1PuAG4EXUU3a3xnCrwNuk3Q58BPgHVHuzHEcx+mEke7GudLG3bid3Nzm\nWMaLO6m8HIfCJpHZ5vy2OsSKH0OeT97Hk9cv3mxM3mfgbuw4KRn44apCGzy95Dw31OSjYcciz0Z+\npKYEl+hpULRhKcUX3Zlfjn3GeUihalhpRy2XLYX5M9GMMnosaeTWKdqwlOBu7Dgj6S0oGfYOmgxD\nJX0faFjvKsP0mneKNiyO4+RD0vrdG2lF4YbFcZzRjOiOWO2v47hhcZxp0Ota5DgzjBsWpxUT14kD\nBKz4IqbXvX2Jt1SOrfg6ac40+Q7NVClgXqxL3LA4rUi1/EqTyWEn7ns0UoPKOqOEj73+XO4rOmeU\n9K0p2rCU4G6cqqCXtMhj7qSsjPrFL6guz4q4RjVNIyn379osyUxN0YYl1VvQsTJz0rfKx43vq/C2\nY+pv3g9LbzP/VIJTPGMtQlk8S073A7aSxp5jsI2FY/efw+eWe/Nk1jBbsQnrZJx+xKNA9TQs38aW\nNcHmlEXRPRZnTpBwV9Y0jGfL+p9UX2SgxOGbVvhk/VDms8fiOI7jJGNuVjde1opq0vyJkC4lzbFk\ntbpxn/wZ5G46qe4p595iy+66nMS67jRkxpy872J147QyfXVjx3Fypk8jwGClg0KhDdi25O7CPC3c\nsAygxMx0nKniD8kKyuixxJXXDzcsA2g1bDbgZH/unNyZpIy2qe/8WZgv3LA4zhwzUWu4hbGYtNXt\nhqksGnmFSdoh6ZCkh2phJ0naJWlf+L8qhEvS9ZL2S3pA0jm1c7aE+PskbamFnyvpwXDO9VLiYtTn\nPYIV7xQ43dEvX5wkTOvdEn+PZb5o6m58I7CpJ+wq4G4zWw/cHX4DXASsD9tW4AaoDBGwDXgdcB6w\nbckYhThba+f1XisvvOJzZoQRbayhW5vIk1zHH6/yaGRYzOxbwOGe4EuAm8L+TcBba+E3W8W9wImS\nVgMXArvM7LCZ/RzYBWwKx15mZt+xyvf55posx3EcpzAmmWM51cwOApjZQUmnhPA1wJO1eIshbFj4\nYp/wFUjaStWzAU6fQPUJ8b6548wutee713X66G/vRg0lxeR9v1rXxghfGWi2HdgO1QuS44zzto07\nTvGpX6fpNWPfSwx5TcPqHJ2h6hNvKbx3DmvSNOobviLQ2p3fUocY8ZvE7bqcdMm07r3JMx+jLDWl\npDxaYhLD8pSk1aG3sho4FMIXgdNq8dYCB0L4xp7we0L42j7xR5LirfYV543vbRxNl1RykshLuKhB\nKz179Mj5bekU8UuROU785OU/QhmeJD1TrHrRGyc1k6wVthNY8uzaAtxeC78seIedD/wyDJndBVwg\naVWYtL8AuCsce1bS+cEb7LKarO7JeaZxYeGY20wOqxI7juPQsMci6Raq3sbJkhapvLuuA26TdDnw\nE+AdIfrXgIuB/cCvgHcBmNlhSdcCu0O8j5rZkkPA+6g8z14E3Bm2eDQx0T5mGp/eRZYcx5kLGhkW\nM7t0wKHf7xPXgCsGyNkB7OgTvgd4TRNdnIJwY+KUiDcyJ2bmVzeWaPbhrp50yGHsPbWcJPIS9g7b\nzbE0m2TJIZ9LmQ/J4b7GkTtteU1ldjfH4qsbd0rOHh6xdYslbxKPmiYc89ZrJuSoB9qQ6Nl1rIYp\nVGhDMHeyKwOF44ZlCF23ZFPJSCvTRsoc91LLZDasCFK2sIfJcdISO41TPFPzjBsWx3GKI/ZQmBOX\n+TAsS6WwXwnyoYXy6Tcg7ThOZ8yHYVnCjYhTAl5OncKZL8PizCa59lAa6mWA5MYEaJRmVvs7i1i1\n7lEtoLx7dcPiOBnQ1smgFErSNRpauV7evH3jqWjDkmqBx0nPSy0zV72mKXPc11lix58lUiyCOasG\nsy2D3ldpQ+r6LiZFG5ZUPcScXXpz1i2VvH7ym15gWKxUek7jwc7BCHS52GGu73FBukG6WGV1GuWz\naMPiOPNK7DfUc2jltiHrhlAXaZnZeoiTrG7sOM4ASquos6FB5Tdv8xUl4j0Wx0lAVT82qwCr9eyc\no4wcY5uOGl0hrERHsGW4YXGcWaSPZ9JASq/FcsMs+TzjqOsDnXabfSjMcRzHiUrRPZaSvE5iysxZ\nt1TyUsn3uZAyKa28TpUMbqZow7Kiq9nwGxyjKMGlN6bMEt2Nc03LJVlOWnIvr9Nu/KX+XEVbijYs\njuM0wBfpnAmGGr8W3sbTyP6x51gk/Zak+2vbM5I+IGlB0k9r4RfXzrla0n5Jj0q6sBa+KYTtl3TV\npDflOE4NafnmzB5mo7cpMnaPxcweBc4GkHQc8FPgq8C7gE+Z2Sfq8SWdCWwGzgJeCXxD0qvD4c8A\nbwIWgd2SdprZD8bVzXEcx+mOWENhvw88ZmY/1uAW0SXArWb2HPAjSfuB88Kx/Wb2OICkW0Pc9obF\n3SYdx3E6J5a78WbgltrvKyU9IGmHpFUhbA3wZC3OYggbFL4CSVsl7ZG0B34WSXXHmUHMqjfUMxsi\nceaDiXsskl4IvAW4OgTdAFxL5ahwLfBJ4N30n14y+hu3vqXdzLYD26vrbrB5dTdOIbM0982c03Ic\nWTmsdJtChxzuqwt5TWUOitMvvKTpsRhDYRcB3zOzpwCW/gNI+jxwR/i5CJxWO28tcCDsDwofyrLG\nVptU72mlleAeW4KOqeRO4y3m2O7Gpay8nXO+p5CVQm7u7u/9ZKcmxlDYpdSGwSStrh17G/BQ2N8J\nbJZ0gqQzgPXAfcBuYL2kM0LvZ3OI6ziO4xTIRD0WSb9B5c31nlrwxyWdTTWc9cTSMTN7WNJtVJPy\nR4ArzOz5IOdK4C7gOGCHmT08iV6O4zhOd8gKnbyTNpjZnnpA85N9KGxq8ku499jX8KGwMtJzfofC\ntNfMNqSRXuGLUDqO4zhR8SVdAiV4MZWgYyq5uXzutwtZqWXnnO8pZKWQm7uX4rQp2rAsT3jDGiyY\nI6yv43Npw0y5DomUcK+pZMaWXdoQYCqZuadnrvc7TH5qijYsKxN/dG70i1Fyy8BxHCc3fI7FmQ69\nCyGO2mLKd/Ji48bl+bNxY9caOZFxw+I4juNExQ2L4ziOE5Wi51gcxymQe+7pWgMnMUUblhLcGVPK\nz9XttJ+ctk4uvTJG6Tb043otFvqLRc5urNO4RkllM0e5pU8NFm1YcneTjCVzJtxOWz4o/XQZquMQ\n+U3ua9T9t0mfJnFjxZmEXMt7Kpk5pGdbHVKlZ2qKNixOQaS2jIMsUf3/cAFR1XGcecYNi+N0Ra/B\nK3TdPsfpxb3CHMdxnKh4j8VxRiFVA2X1Dob3LhxnIG5YnPliySCU7nbj5E1f18hIHghN6LjhU7Rh\nKcFNMpZMdzudT9wNPq7MEsujRKMFdnvP6ZKiDUuubo0uM2M3ya6fuJZ07R7rMtvJTIEZUdz1l3B3\nY8fJAbM070D4PI0zo0zsFSbpCUkPSrpf0p4QdpKkXZL2hf+rQrgkXS9pv6QHJJ1Tk7MlxN8nacuk\nejmO4zjdEMvd+J+Z2dm17yhfBdxtZuuBu8NvgIuA9WHbCtwAlSECtgGvA84Dti0ZI8cZG7Pqw27W\nYHOcnGlShjMqz6neY7kEuCns3wS8tRZ+s1XcC5woaTVwIbDLzA6b2c+BXcCmRLo5juOkJcPKfprE\nMCwGfF3SXklbQ9ipZnYQIPw/JYSvAZ6snbsYwgaFL0PSVkl7qiG3n0VQ3XEcx4lNjMn715vZAUmn\nALsk/XBI3H7+CIN8HlZ+eNhsO7AdQNpgObs1usx8dGx6zqh4ba7dJG6sOJOQSx7Nk8y21y3MkRGI\nYFjM7ED4f0jSV6nmSJ6StNrMDoahrkMh+iJwWu30tcCBEL6xJ/ye0deeVPvlpFz9NJXsmHKbyspl\nhdYc0jPmiraxV1BuQwkrCHeZjl3nYavrb1uAa67pH2HbNjTgUEwmGgqT9GJJL13aBy4AHgJ2Akue\nXVuA28P+TuCy4B12PvDLMFR2F3CBpFVh0v6CENY9w76l7t9WdxzHWcGkPZZTga+qqlSPB/7czP6X\npN3AbZIuB34CvCPE/xpwMbAf+BXwLgAzOyzpWmB3iPdRMzs8oW6O4zhOB8gK9VaQNpjZnsgy+3Q3\n2/REhqRlLkM3MWT5UFi7+F0Po0xTng+FZTAUNrI8am/t1ZAk+Jv3jtNLi8bE0VWPI3wmdMUKyssu\nVGYD0JlP3LCMwh9oJ1Ny9eLz6UanaMOSpAA3FKoxPmWb6oEr4UEuxV20JOZxKMwpg6INS5IC3LDw\ntr12LnMCo2SloqQ5FsdxJqNow+IUStvae9rDkS2ul7KF7jil4oZlGL01hjdnnVSkWprfcTog1SKU\njuM4zpziPZZheA9lIJ40juMMwg2LMxZtXyJMQS5yZ8XIluBuHGtB0ZQyU1w71fVTUbRhSZGATevL\nlKvs5iI3JnUd204jDLq/UuYj2iwll6pSakLunntdreDQ5Zv8beK3uXZqijYsaSqWZpOo41SOJTy4\nqQrdMj0jOIVNpGekZXraUEIF6zixKNqwOPGZSotnxEVSt7za1Mle4TpOe9ywOMUydqt9yh0WN07O\nvOGGxXGcchlmtUuZhJtB3LA4xTJ+T6DNrOm413BS0Jvnw3LSe4rdUbRhaVJwrEXNUF9YspRCWYqe\nKcjdGaIud15J7rk2YYdlnvMmJUUblkaVwBjj6am8bpK4R0f2NnLSkNsS9ysbXObl04lG0YalNNxF\ndH7Jzt24T9nx8unEYuy1wiSdJumbkh6R9LCk94fwBUk/lXR/2C6unXO1pP2SHpV0YS18UwjbL+mq\nyW7JcZyh9Knx2wwZZ4XZ4M3pjEl6LEeAD5rZ9yS9FNgraVc49ikz+0Q9sqQzgc3AWcArgW9IenU4\n/BngTcAisFvSTjP7wQS6HcMLmOM402DU56WXMdv10tiGxcwOAgfD/rOSHgHWDDnlEuBWM3sO+JGk\n/cB54dh+M3scQNKtIW4cw+JkwVCv0J4nsffhHOdrnY7jdEeUZfMlrQNeC3w3BF0p6QFJOyStCmFr\ngCdrpy2GsEHhTgcsORkM25rGq8efZMRi0DltdEhxX23kOnHoIs+bXtc5xsST95JeAnwZ+ICZPSPp\nBuBaqobntcAngXfTv4No9DdufasbSVuBrdWv05NmZgrZJciMvTDeUrwUxNS17YR4v/iDZKSoeCaR\nN+g25618NsnzVuWi5b2m6MXnYuAmMiySXkBlVL5oZl8BMLOnasc/D9wRfi4Cp9VOXwscCPuDwpdh\nZtuB7ZXsDRZ7FdE2FWWKCrgNqWQ6aYiVV3Hy3VZmtlnUAUcvn8MZmTYj7nWStJ1GOk7iFSbgC8Aj\nZvbHtfDVtWhvAx4K+zuBzZJOkHQGsB64D9gNrJd0hqQXUk3w7xxXL8dxGlAbV/Q5rEgspeWwcd85\n8VibpMfyeuCdwIOS7g9hHwYulXQ2Ve/tCeA9AGb2sKTbqCbljwBXmNnzAJKuBO4CjgN2mNnDE+jl\nzBopmr+O4yRDjcaTMqQaCtvTIJ4PhaWSGS29mvTNewTE1jX1HEteQ2HzIXPqcyxjxB8qaBgTXETS\nXjPbMLaABvib9073WPi4Wqkv6TlOZIQV3Ul3w1IwszSZ6TizQqznsuTn2w3LlEhRSNzrphxipm0J\nbsGpZJZArI/DpeqxTCNf3LBMCTcC843PscSR6ZSBGxYnfwoZbDbU9/0Dq/115oSFhf77c4J7hY0R\n173C2sXPwftmGl5hw1+ljpO+vdmmAAAHvElEQVRZOeR7KTI79Qqrl4UxEiLtUJh7hTmZ48MTjtOH\nbdu61qBT3LA4E+Gff3WcPszh8Fed2TYsGvB9hA6G/+bd66aJrm3vJ3V8Z/bJoVymltMFRRuWkS+n\njnFe08xsm+ldjzc3lZmKrua5lsUbdIOW5nvvvXrkKGvWZcZ8gT3FnEyMa44jMzVFG5ZxVwgddN44\nmdi0sJVCSbrGJgfD2pSulolJ4mTRcEmf2HrOc1lPTdGGxYlP20pjlpjne3ecmMy2YSnUldpxcmXQ\nuzrDznDmj9k2LI5TJ3FDY9DCgdUCm44zP7hhcRynW3xkYSglDru6YZkSOXjION0xs/k/qIs2+OfM\nE/t+S/QKG/vTxDkgxd1S0uRrpW22gTK3LVB9ILXPtm1hnr+W2ilTy/8xZZVw76UQO59LpOgeS4mW\n3HFmCn9onD4U3WNxHMdx8iOb1Y0lbQL+BDgO+DMzu254/GarG7fToX0vqIu3cXORmWI16KRv3keI\nNyz+IBk55FU0WS17KDE/sZviBcnY5b2MN+/Tr26cRY9F0nHAZ4CLgDOBSyWd2a1WjuP0InecdhqQ\nhWEBzgP2m9njZvZ3wK3AJR3r5DhOP+Zh9tmZiFwm79cAT9Z+LwKv640kaSuwNfx8TtJDsRUZZy5y\nxDknA08X4G46lp4pFu0cEvdk4OlU145x7yFsmZ7jyB73+i05qmek/IkSvw/J9Ix8L0me9QRl57ei\nS+whF8PSL+lWNHfMbDuwHUDSntTjhDFwPeNRgo7gesbG9YyLpLiT033IZShsETit9nstcKAjXRzH\ncZwJyMWw7AbWSzpD0guBzcDOjnVyHMdxxiCLoTAzOyLpSuAuKnfjHWb28IjTtqfXLAquZzxK0BFc\nz9i4nnFJrmc277E4juM4s0EuQ2GO4zjOjOCGxXEcx4mLmRW1AZuAR4H9wFVTvO4TwIPA/cCeEHYS\nsAvYF/6vCuECrg86PgCcU5OzJcTfB2yphZ8b5O8P56qhXjuAQ8BDtbDkeg26Rks9F4CfhjS9H7i4\nduzqcM1HgQtH5T9wBvDdoM+XgBeG8BPC7/3h+LohOp4GfBN4BHgYeH+O6TlEz9zS89eB+4D/E/S8\nZlzZsfRvqeeNwI9q6Xl2Bs/RccD3gTtyTMuj5zSpvHLZQqI+BrwKeGEoCGdO6dpPACf3hH18KQOA\nq4CPhf2LgTtDATwf+G6tED0e/q8K+0uV1H3A74Zz7gQuaqjXG4BzWF5hJ9dr0DVa6rkA/Mc+cc8M\neXtCKNSPhbwfmP/AbcDmsP9Z4H1h/w+Az4b9zcCXhui4mlBJAC8F/jroklV6DtEzt/QU8JKw/wKq\nyun8trJj6t9SzxuBt/eJ3+Vz9B+AP+eYYckqLY/q2aTyymULGXNX7ffVwNVTuvYTrDQsjwKrw/5q\n4NGw/zng0t54wKXA52rhnwthq4Ef1sKXxWug2zqWV9jJ9Rp0jZZ6LtC/IlyWr1Tegr87KP/Dw/o0\ncHxvOVk6N+wfH+I17Q3eDrwp1/Tso2e26Qn8BvA9qhU1WsmOqX9LPW+kv2HpJN+p3u+7G3gjcMc4\n+TSttCxtjqXf0i9rpnRtA74uaW9YWgbgVDM7CBD+nzJCz2Hhi33Cx2Uaeg26RluulPSApB2SVo2p\n58uBX5jZkT56Hj0nHP9liD8USeuA11K1XrNNzx49IbP0lHScpPuphkF3UbWK28qOqX8jPc1sKT3/\na0jPT0k6oVfPhvrEyvdPA/8J+Pvwe5x8Sp6WUN7kfaOlXxLxejM7h2oF5iskvWFI3EF6tg2PTW56\n3QD8Q+Bs4CDwyRAeU8/W9yDpJcCXgQ+Y2TPDokbUszV99MwuPc3seTM7m6q1fR7wO2PITp7OvXpK\neg1Vi/23gX9MNbz1och6NkbSm4FDZra3HjxEbmdpCeUZls6WfjGzA+H/IeCrVA/JU5JWA4T/h0bo\nOSx8bZ/wcZmGXoOu0Rgzeyo80H8PfJ4qTcfR82ngREnH94QvkxWO/wPg8CCdJL2AqrL+opl9ZcS9\ndpae/fTMMT2XMLNfAPdQzUm0lR1T/6Z6bjKzg1bxHPDfGT89Y+T764G3SHqCavX3N1L1YPJMyybj\njrlsVGOFj1NNOi1NMJ01heu+GHhpbf9/U3lQ/BHLJ94+Hvb/Bcsn9+4L4SdReZmsCtuPgJPCsd0h\n7tLk3sUt9FvH8rmL5HoNukZLPVfX9v89cGvYP4vlE4yPU00uDsx/4C9YPsH4B2H/CpZPYt42RD8B\nNwOf7gnPKj2H6Jlber4CODHsvwj4NvDmtrJj6t9Sz9W19P40cF0mz9FGjk3eZ5WWR3VsWnnlslF5\nZPw11VjtR6Z0zVeFhF5yR/xICH851WTavvB/qRCJ6sNlj1G5GG6oyXo3lTvffuBdtfANwEPhnD+l\n+QTzLVTDHv+PqtVx+TT0GnSNlnr+j6DHA1Rrw9Urxo+Eaz5KzUNuUP6HPLov6P8XwAkh/NfD7/3h\n+KuG6PhPqLr5D1Bz2c0tPYfomVt6/iMq19gHwj3/53Flx9K/pZ5/FdLzIeB/csxzrLPnKJyzkWOG\nJau0XNp8SRfHcRwnKqXNsTiO4ziZ44bFcRzHiYobFsdxHCcqblgcx3GcqLhhcRzHcaLihsVxHMeJ\nihsWx3EcJyr/H+s/XZQiueVoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a282b7dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# for learned kdtree，\n",
    "visualize_kdnodes_and_query(kdnodes, testing_set, domains_)\n",
    "visualize_kdnodes_and_query(fast_kdnodes, testing_set, domains_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 820,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEJCAYAAABVFBp5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGqpJREFUeJzt3X+sX3Wd5/Hny1bQVZGCYEiLC6Md\ntZodtHehs25cFAOFnVgmQVPiSHWb7ciCYeJsljIzWR11EzWrrCTobkdYipmxMIyzNATsdgXj7Ay/\nWkGgMMgdYOUOhB/bghizODDv/eP7KX4pt/d+77lfer/f5vlITu457/M557w/vbf3fc85n3O+qSok\nSeriVQudgCRpfFlEJEmdWUQkSZ1ZRCRJnVlEJEmdWUQkSZ0NXESSLEpyR5Lr2vLxSW5N8kCSq5Ic\n0uKHtuXJtv64vn1c1OL3JzmtL766xSaTbOyLT3sMSdJomMuZyAXAfX3LXwYurqrlwB5gfYuvB/ZU\n1duAi1s7kqwA1gLvAlYD32iFaRFwKXA6sAI4u7Wd6RiSpBEwUBFJsgz418C32nKADwLXtCabgTPb\n/Jq2TFt/Smu/BthSVc9V1UPAJHBimyar6sGq+iWwBVgzyzEkSSNg8YDt/gvwH4A3tOUjgaer6vm2\nPAUsbfNLgUcAqur5JM+09kuBW/r22b/NI/vET5rlGC+RZAOwAeB1sPIdA3ZKktSzE56qqqPmut2s\nRSTJbwFPVNXOJCfvDU/TtGZZt7/4dGdDM7V/ebBqE7AJYCKpHdM1kiTtV+D/dNlukDOR9wEfTnIG\n8BrgMHpnJocnWdzOFJYBj7b2U8CxwFSSxcAbgd198b36t5ku/tQMx5AkjYBZ74lU1UVVtayqjqN3\nY/zGqvoYcBNwVmu2Dri2zW9ty7T1N1bvLY9bgbVt9NbxwHLgNuB2YHkbiXVIO8bWts3+jiFJGgHz\neU7kQuAzSSbp3b+4rMUvA45s8c8AGwGqahdwNXAv8D3gvKp6oZ1lnA9sozf66+rWdqZjSJJGQA62\nV8F7T0SS5i6ws6om5rqdT6xLkjqziEiSOrOISJI6s4hIkjqziEiSOrOISJI6s4hIkjqziEiSOrOI\nSJI6s4hIkjqziEiSOrOISJI6s4hIkjqziEiSOrOISJI6s4hIkjqziEiSOrOISJI6m7WIJHlNktuS\n/DjJriR/3OJXJHkoyZ1tOqHFk+SSJJNJ7kry3r59rUvyQJvW9cVXJrm7bXNJkrT4EUm2t/bbkywZ\n/j+BJKmrWT9jvf1Cf11V/TzJq4H/DVwAfAq4rqqu2af9GcCngTOAk4CvV9VJSY4AdgATQAE7gZVV\ntSfJbW2ftwDXA5dU1Q1JvgLsrqovJdkILKmqC2fOd6JqDD5lPYFx+Hj7cciz9yfHeOQ56jmCeQ7b\n+OSZV+Yz1qvn523x1W2a6Z9kDXBl2+4W4PAkxwCnAdurandV7QG2A6vbusOq6ubqVbQrgTP79rW5\nzW/ui0uSRsDiQRolWUTvzOFtwKVVdWuSc4H/lOQ/At8HNlbVc8BS4JG+zadabKb41DRxgDdX1WMA\nVfVYkqMHy3eQVgvPPIdrXPLU8IzL93xc8uxioBvrVfVCVZ0ALANOTPJu4CLgHcA/B44A9l5mmu6f\nqzrEB5ZkQ5IdSXbAk3PZVDqgktGfxilPLbyBzkT2qqqnk/wAWF1V/7mFn0vy34F/35angGP7NlsG\nPNriJ+8T/0GLL5umPcDjSY5pZyHHAE/sJ69NwCbYe09kLr1aGONznXT08/SeyHCZ53CNU55dDDI6\n66gkh7f51wIfAv62/VLfe+P9TOCetslW4Jw2SmsV8Ey7JLUNODXJkjbK6lRgW1v3bJJVbV/nANf2\n7WvvKK51fXFJ0ggY5EzkGGBzuy/yKuDqqrouyY1JjqJ3OepOeqO1oDe66gxgEvgF8EmAqtqd5AvA\n7a3d56tqd5s/F7gCeC1wQ5sAvgRcnWQ98FPgI7Mlu5KdY3GuWzD9hbwRMw55vvhH3jjkOeI5gnkO\n27jk2dWsQ3zHzUQyBgN8JWm0BF6ZIb6SJO2PRUSS1JlFRJLUmUVEktSZRUSS1JlFRJLUmUVEktSZ\nRUSS1JlFRJLUmUVEktSZRUSS1JlFRJLUmUVEktSZRUSS1JlFRJLUmUVEktSZRUSS1JlFRJLUmUVE\nktTZrEUkyWuS3Jbkx0l2JfnjFj8+ya1JHkhyVZJDWvzQtjzZ1h/Xt6+LWvz+JKf1xVe32GSSjX3x\naY8hSRoNg5yJPAd8sKp+AzgBWJ1kFfBl4OKqWg7sAda39uuBPVX1NuDi1o4kK4C1wLuA1cA3kixK\nsgi4FDgdWAGc3doywzEkSSNg1iJSPT9vi69uUwEfBK5p8c3AmW1+TVumrT8lSVp8S1U9V1UPAZPA\niW2arKoHq+qXwBZgTdtmf8eQJI2Age6JtDOGO4EngO3A3wFPV9XzrckUsLTNLwUeAWjrnwGO7I/v\ns83+4kfOcIx989uQZEeSHU8O0iFJ0lAMVESq6oWqOgFYRu/M4Z3TNWtfs591w4pPl9+mqpqoqomj\npmsgSXpFzGl0VlU9DfwAWAUcnmRxW7UMeLTNTwHHArT1bwR298f32WZ/8admOIYkaQQMMjrrqCSH\nt/nXAh8C7gNuAs5qzdYB17b5rW2Ztv7GqqoWX9tGbx0PLAduA24HlreRWIfQu/m+tW2zv2NIkkbA\n4tmbcAywuY2iehVwdVVdl+ReYEuSLwJ3AJe19pcB304ySe8MZC1AVe1KcjVwL/A8cF5VvQCQ5Hxg\nG7AIuLyqdrV9XbifY0iSRkB6f/AfPCaS2rHQSUjSmAnsrKqJuW7nE+uSpM4sIpKkziwikqTOLCKS\npM4sIpKkziwikqTOLCKSpM4sIpKkziwikqTOLCKSpM4sIpKkziwikqTOLCKSpM4sIpKkziwikqTO\nLCKSpM4sIpKkziwikqTOZi0iSY5NclOS+5LsSnJBi38uyd8nubNNZ/Rtc1GSyST3JzmtL766xSaT\nbOyLH5/k1iQPJLkqySEtfmhbnmzrjxtm5yVJ8zPImcjzwO9X1TuBVcB5SVa0dRdX1Qltuh6grVsL\nvAtYDXwjyaIki4BLgdOBFcDZffv5ctvXcmAPsL7F1wN7quptwMWtnSRpRMxaRKrqsar6UZt/FrgP\nWDrDJmuALVX1XFU9BEwCJ7ZpsqoerKpfAluANUkCfBC4pm2/GTizb1+b2/w1wCmtvSRpBMzpnki7\nnPQe4NYWOj/JXUkuT7KkxZYCj/RtNtVi+4sfCTxdVc/vE3/Jvtr6Z1r7ffPakGRHkh1PzqVDkqR5\nGbiIJHk98BfA71XVz4BvAm8FTgAeA766t+k0m1eH+Ez7emmgalNVTVTVxFEz9kKSNEyLB2mU5NX0\nCsifVtV3Aarq8b71fwJc1xangGP7Nl8GPNrmp4s/BRyeZHE72+hvv3dfU0kWA28Eds+U605WQu0Y\npFsLKoF6WTkcPeOQ594LnOOQ56jnCOY5bOOSJx3vFAwyOivAZcB9VfW1vvgxfc1+G7inzW8F1raR\nVccDy4HbgNuB5W0k1iH0br5vraoCbgLOatuvA67t29e6Nn8WcGNrL0kaAYOcibwP+Dhwd5I7W+wP\n6I2uOoHe5aWHgd8FqKpdSa4G7qU3suu8qnoBIMn5wDZgEXB5Ve1q+7sQ2JLki8Ad9IoW7eu3k0zS\nOwNZO4++SpKGLAfbH/bJRJWXs4ZmHPL0ctZwmedwjU+e2VlVE3PdzifWJUmdWUQkSZ1ZRCRJnVlE\nJEmdWUQkSZ1ZRCRJnQ30xPo4WcnOzk9eHkgF07/UZcSMQ54vjp4chzxHPEcwz2Eblzy78kxEktSZ\nRUSS1JlFRJLUmUVEktSZRUSS1JlFRJLUmUVEktTZQfecyItG/N3L4/N66NHP01fBD5d5DtdY5DmP\nZ+sO3iIy4g8cjssDSOOQpw8bDpd5Dte45NmVl7MkSZ1ZRCRJnc1aRJIcm+SmJPcl2ZXkghY/Isn2\nJA+0r0taPEkuSTKZ5K4k7+3b17rW/oEk6/riK5Pc3ba5JOldi9rfMSRJo2GQM5Hngd+vqncCq4Dz\nkqwANgLfr6rlwPfbMsDpwPI2bQC+Cb2CAHwWOAk4EfhsX1H4Zmu7d7vVLb6/Y0iSRsCsRaSqHquq\nH7X5Z4H7gKXAGmBza7YZOLPNrwGurJ5bgMOTHAOcBmyvqt1VtQfYDqxu6w6rqpurqoAr99nXdMeQ\nJI2AOd0TSXIc8B7gVuDNVfUY9AoNcHRrthR4pG+zqRabKT41TZwZjrFvXhuS7Eiy48m5dEiSNC8D\nF5Ekrwf+Avi9qvrZTE2niVWH+MCqalNVTVTVxFFz2VCSNC8DFZEkr6ZXQP60qr7bwo+3S1G0r0+0\n+BRwbN/my4BHZ4kvmyY+0zEkSSNgkNFZAS4D7quqr/Wt2grsHWG1Dri2L35OG6W1CnimXYraBpya\nZEm7oX4qsK2tezbJqnasc/bZ13THkCSNgEGeWH8f8HHg7iR3ttgfAF8Crk6yHvgp8JG27nrgDGAS\n+AXwSYCq2p3kC8Dtrd3nq2p3mz8XuAJ4LXBDm5jhGJKkEZAa+Ze6zM1EUjsWOglJGjOBnVU1Mdft\nfGJdktSZRUSS1JlFRJLUmUVEktSZRUSS1JlFRJLUmUVEktSZRUSS1JlFRJLUmUVEktTZIO/OGis7\nWQlj8OKTBEb9jTNpL+k3z+EYh+85mOewjUueL/5HmiPPRCRJnR10ZyIr2dm5oh5IBdN/HNcIefGP\nJ/McinH4noN5Dtu45NmVZyKSpM4sIpKkziwikqTOLCKSpM4sIpKkzmYtIkkuT/JEknv6Yp9L8vdJ\n7mzTGX3rLkoymeT+JKf1xVe32GSSjX3x45PcmuSBJFclOaTFD23Lk239ccPqtCRpOAY5E7kCWD1N\n/OKqOqFN1wMkWQGsBd7VtvlGkkVJFgGXAqcDK4CzW1uAL7d9LQf2AOtbfD2wp6reBlzc2kmSRsis\nRaSqfgjsHnB/a4AtVfVcVT0ETAIntmmyqh6sql8CW4A1SQJ8ELimbb8ZOLNvX5vb/DXAKa29JGlE\nzOeeyPlJ7mqXu5a02FLgkb42Uy22v/iRwNNV9fw+8Zfsq61/prV/mSQbkuxIsuPJeXRIkjQ3XYvI\nN4G3AicAjwFfbfHpzhSqQ3ymfb08WLWpqiaqauKombKWJA1VpyJSVY9X1QtV9Y/An9C7XAW9M4lj\n+5ouAx6dIf4UcHiSxfvEX7Kvtv6NDH5ZTZJ0AHQqIkmO6Vv8bWDvyK2twNo2sup4YDlwG3A7sLyN\nxDqE3s33rVVVwE3AWW37dcC1ffta1+bPAm5s7SVJI2LWFzAm+Q5wMvCmJFPAZ4GTk5xA7/LSw8Dv\nAlTVriRXA/cCzwPnVdULbT/nA9uARcDlVbWrHeJCYEuSLwJ3AJe1+GXAt5NM0jsDWTunno14vRmH\n10OP0yvWYTzyHPUcwTyHbSzynMeYpRxsf9xPJGPwaSKSNFoCO6tqYq7b+cS6JKkzi4gkqTOLiCSp\nM4uIJKkzi4gkqbOD7jPWd7ISxmB81jgM+xunobMwHnmOeo5gnsM2Lnl2HeZ70BWRleyc15jnA6Vg\n+he7jJAXf+7NcyjG4XsO5jls45JnV17OkiR1ZhGRJHVmEZEkdWYRkSR1ZhGRJHVmEZEkdWYRkSR1\nZhGRJHVmEZEkdWYRkSR1ZhGRJHU2yGesXw78FvBEVb27xY4ArgKOo/cZ6x+tqj1JAnwdOAP4BfCJ\nqvpR22Yd8Edtt1+sqs0tvhK4AngtcD1wQVXV/o4xW76+gHF4xunFhjAeeY56jmCewzYueXZ+52BV\nzTgB7wfeC9zTF/sKsLHNbwS+3ObPAG6g97qxVcCtLX4E8GD7uqTNL2nrbgN+s21zA3D6TMeYbVrZ\n+345OTk5Oc1hAnYM8jt232nWy1lV9UNg9z7hNcDmNr8ZOLMvfmX13AIcnuQY4DRge1XtbmcT24HV\nbd1hVXVz9arFlfvsa7pjSJJGRNd7Im+uqscA2tejW3wp8Ehfu6kWmyk+NU18pmO8TJINSXYk2fFk\nxw5JkuZu2DfWp7uoVh3ic1JVm6pqoqomjprrxpKkzroWkcfbpSja1ydafAo4tq/dMuDRWeLLponP\ndAxJ0ojoWkS2Auva/Drg2r74OelZBTzTLkVtA05NsiTJEuBUYFtb92ySVW1k1zn77Gu6Y0iSRsQg\nQ3y/A5wMvCnJFPBZ4EvA1UnWAz8FPtKaX09vhNYkvSG+nwSoqt1JvgDc3tp9vqr23qw/l18N8b2h\nTcxwDEnSiEgbQnvQmEjG4CkRSRotgZ1VNTHX7XxiXZLUmUVEktSZRUSS1JlFRJLUmUVEktSZRUSS\n1JlFRJLUmUVEktSZRUSS1JlFRJLUmUVEktSZRUSS1JlFRJLUmUVEktSZRUSS1JlFRJLUmUVEktSZ\nRUSS1Nmsn7E+kyQPA88CLwDPV9VEkiOAq4DjgIeBj1bVniQBvk7vM9h/AXyiqn7U9rMO+KO22y9W\n1eYWX8mvPn/9euCCGvTzfEf8Y3+TkU+RpPd1HPIc9RzBPIfNPIdo73/2DuZVRJoPVNVTfcsbge9X\n1ZeSbGzLFwKnA8vbdBLwTeCkVnQ+C0wABexMsrWq9rQ2G4Bb6BWR1cANA2U1j3+UA6EARjtFXvy5\nH4c8RzxHMM9hM8/R8EpczloDbG7zm4Ez++JXVs8twOFJjgFOA7ZX1e5WOLYDq9u6w6rq5nb2cWXf\nviRJI2C+RaSA/5lkZ5INLfbmqnoMoH09usWXAo/0bTvVYjPFp6aJS5JGxHwvZ72vqh5NcjSwPcnf\nztB2uhO66hB/+Y57BWwDwFtmzleSNETzOhOpqkfb1yeAvwROBB5vl6JoX59ozaeAY/s2XwY8Okt8\n2TTx6fLYVFUTVTVx1Hw6JEmak85FJMnrkrxh7zxwKnAPsBVY15qtA65t81uBc9KzCnimXe7aBpya\nZEmSJW0/29q6Z5OsaiO7zunblyRpBMznctabgb/s/X5nMfBnVfW9JLcDVydZD/wU+Ehrfz294b2T\n9Ib4fhKgqnYn+QJwe2v3+ara3ebP5VdDfG9g0JFZkqQDIoM+djEuJpLasdBJSNKYCeysqom5bucT\n65KkziwikqTOLCKSpM4sIpKkzobx7qyRspOVMAa31sfhpWy+gHG4zHO4zHPIOr5v0DMRSVJnFhFJ\nUmcWEUlSZwfdPZGV7Bz5zxKB8fiMAT9PZLjMc7jMczR4JiJJ6swiIknqzCIiSerMIiJJ6uygu7Hu\nw4bD48OGw2Wew2WeQ+bDhpKkA80iIknqzCIiSersoLsn4sOGw+PDhsNlnsNlnqPBMxFJUmcjX0SS\nrE5yf5LJJBsXOh9J0q+MdBFJsgi4FDgdWAGcnWTFwmYlSdprpIsIcCIwWVUPVtUvgS3AmgXOSZLU\njPqN9aXAI33LU8BJ+zZKsgHY0BafC9xzAHJbKG8CnlroJF5BB3P/Dua+gf0bd2/vstGoF5HpxjS8\n7NnPqtoEbAJIsqOqJl7pxBaK/RtfB3PfwP6NuySdXvUx6pezpoBj+5aXAY8uUC6SpH2MehG5HVie\n5PgkhwBrga0LnJMkqRnpy1lV9XyS84FtwCLg8qraNctmm175zBaU/RtfB3PfwP6Nu079S43F6yUl\nSaNo1C9nSZJGmEVEktTZ2BaR2V6HkuTQJFe19bcmOe7AZ9nNAH37TJJ7k9yV5PtJ/ulC5NnVoK+y\nSXJWkkoyVsMqB+lfko+27+GuJH92oHOcjwF+Pt+S5KYkd7Sf0TMWIs8uklye5Ikk0z5rlp5LWt/v\nSvLeA53jfAzQv4+1ft2V5G+S/MasO62qsZvo3WT/O+DXgEOAHwMr9mnz74D/2ubXAlctdN5D7NsH\ngH/S5s8dl74N2r/W7g3AD4FbgImFznvI37/lwB3AkrZ89ELnPeT+bQLObfMrgIcXOu859O/9wHuB\ne/az/gzgBnrPsK0Cbl3onIfcv3/R93N5+iD9G9czkUFeh7IG2NzmrwFOScbgHfED9K2qbqqqX7TF\nW+g9PzMuBn2VzReArwD/70AmNwSD9O/fApdW1R6AqnriAOc4H4P0r4DD2vwbGaNnu6rqh8DuGZqs\nAa6snluAw5Mcc2Cym7/Z+ldVf7P355IBf7eMaxGZ7nUoS/fXpqqeB54Bjjwg2c3PIH3rt57eX0bj\nYtb+JXkPcGxVXXcgExuSQb5/vw78epK/TnJLktUHLLv5G6R/nwN+J8kUcD3w6QOT2gEx1/+f42yg\n3y0j/ZzIDAZ5HcpAr0wZQQPnneR3gAngX72iGQ3XjP1L8irgYuATByqhIRvk+7eY3iWtk+n9pfdX\nSd5dVU+/wrkNwyD9Oxu4oqq+muQ3gW+3/v3jK5/eK25cf6/MSZIP0Csi/3K2tuN6JjLI61BebJNk\nMb3T6plOU0fFQK96SfIh4A+BD1fVcwcot2GYrX9vAN4N/CDJw/SuO28do5vrg/5sXltV/1BVDwH3\n0ysq42CQ/q0HrgaoqpuB19B7eeHB4KB/FVOSfwZ8C1hTVf93tvbjWkQGeR3KVmBdmz8LuLHa3aIR\nN2vf2uWe/0avgIzT9XSYpX9V9UxVvamqjquq4+hdl/1wVXV6OdwCGORn83/QGxxBkjfRu7z14AHN\nsrtB+vdT4BSAJO+kV0SePKBZvnK2Aue0UVqrgGeq6rGFTmpYkrwF+C7w8ar6yUAbLfRogXmMMjgD\n+Am9kSJ/2GKfp/cLB3o/uH8OTAK3Ab+20DkPsW//C3gcuLNNWxc652H2b5+2P2CMRmcN+P0L8DXg\nXuBuYO1C5zzk/q0A/preyK07gVMXOuc59O07wGPAP9A761gPfAr4VN/37tLW97vH8Gdztv59C9jT\n97tlx2z79LUnkqTOxvVyliRpBFhEJEmdWUQkSZ1ZRCRJnVlEJGmMzfZSxX3aXpzkzjb9JMm8H3B1\ndJYkjbEk7wd+Tu+dXu+ew3afBt5TVf9mPsf3TESSxlhN81LFJG9N8r0kO5P8VZJ3TLPp2fSeG5mX\ncX13liRp/zbRe4DwgSQnAd8APrh3ZfsMouOBG+d7IIuIJB1Ekrye3ueC/Hnfp18cuk+ztcA1VfXC\nfI9nEZGkg8urgKer6oQZ2qwFzhvWwSRJB4mq+hnwUJKPwIsf6fvix9wmeTuwBLh5GMeziEjSGEvy\nHXoF4e1JppKsBz4GrE/yY2AXL/30ybOBLTWkobkO8ZUkdeaZiCSpM4uIJKkzi4gkqTOLiCSpM4uI\nJKkzi4gkqTOLiCSps/8PM+jmPJK2aXQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116edcc88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEJCAYAAABVFBp5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAHBhJREFUeJzt3X+sXOV95/H3JzYQtgnYgEGs7azd\nxE3iRK2D74KzWWUJRMb2rnKpBJXZNL5NrXVDoUqVSotpVyKFrBQiJd5FArpuYDFRE+PSdLlCOF5v\nAKVqseHe4ACGEt8CxTe2MKyNQ4QWavrdP85zyfEwP848d+bOzOXzko7mnO95zvNjZny/Pj9HEYGZ\nmVmO9/S6A2ZmNricRMzMLJuTiJmZZXMSMTOzbE4iZmaWzUnEzMyyVU4ikuZIelzS/Wl5qaS9kg5I\nukfSqSl+WlqeSOuXlOq4PsWflXRZKb4mxSYkbS7F67ZhZmb9oZ09kS8Dz5SWbwa2RMQy4BiwMcU3\nAsci4kPAllQOScuB9cDHgDXAbSkxzQFuBdYCy4GrUtlmbZiZWR+olEQkLQL+PfDttCzgEuDeVGQb\ncHmaH07LpPWXpvLDwPaIeCMingcmgAvTNBERz0XEm8B2YLhFG2Zm1gfmViz334D/DLw/LZ8NvBoR\nJ9LyJLAwzS8EDgJExAlJx1P5hcCeUp3lbQ7WxC9q0cZJJG0CNgH8Cqz8SMVBmZlZYRxeiYgF7W7X\nMolI+g/AkYgYl3TxVLhO0WixrlG83t5Qs/LvDEZsBbYCDEkxVq+QmZk1JPjHnO2q7Il8CvicpHXA\ne4EzKPZM5kmam/YUFgGHUvlJYDEwKWkucCZwtBSfUt6mXvyVJm2YmVkfaHlOJCKuj4hFEbGE4sT4\ngxHxeeAh4IpUbAS4L82PpmXS+gejeMrjKLA+Xb21FFgGPAo8BixLV2KdmtoYTds0asPMzPrAdO4T\nuQ74iqQJivMXd6T4HcDZKf4VYDNAROwHdgBPAz8AromIt9JexrXALoqrv3akss3aMDOzPqDZ9ih4\nnxMxM2ufYDwihtrdznesm5lZNicRMzPL5iRiZmbZnETMzCybk4iZmWVzEjEzs2xOImZmls1JxMzM\nsjmJmJlZNicRMzPL5iRiZmbZnETMzCybk4iZmWVzEjEzs2xOImZmls1JxMzMsjmJmJlZNicRMzPL\n1jKJSHqvpEcl/UTSfkl/muJ3SXpe0r40rUhxSbpF0oSkJyRdUKprRNKBNI2U4islPZm2uUWSUvws\nSbtT+d2S5nf+LTAzs2wR0XQCBLwvzZ8C7AVWAXcBV9Qpvw7YmbZbBexN8bOA59Lr/DQ/P617FPhk\n2mYnsDbFvwFsTvObgZtb93dlQHjy5MmTp7Ymxlr9fa03za2QZAL4RVo8JU3RZJNh4O603R5J8ySd\nD1wM7I6IowCSdgNrJD0MnBERj6T43cDlKZkMp+0AtgEPA9e16rNZv4pm/3L6hDQ4/YTO9rVc59T7\nUBubUrXdbryf5b5N1d2o7+X19ebLdeaodE5E0hxJ+4AjFIlgb1r1X9Mhqy2STkuxhcDB0uaTKdYs\nPlknDnBeRBwGSK/nVh6ZWR+S+n8apH52uq/lOsuvtfPttNuN97O2j836Xq8f9fqUq1ISiYi3ImIF\nsAi4UNLHgeuBjwD/muIQ1dQeQr3uREa8MkmbJI1JGoOX29nUzMymoa2rsyLiVYpDSmsi4nAU3gD+\nJ3BhKjYJLC5ttgg41CK+qE4c4KV0KIz0eqRBv7ZGxFBEDMGCdoZkZmbTUOXqrAWS5qX504HPAn9f\n+uMuinMYT6VNRoEN6SqtVcDxdChqF7Ba0vx0ldVqYFda95qkVamuDcB9pbpG0vxIKW5mZn2g5Yl1\n4Hxgm6Q5FElnR0TcL+lBSQsoDkftA76Uyj9AcYXWBPA68EWAiDgq6SbgsVTuxqmT7MDVFFd7nU5x\nQn1nin8d2CFpI/AicGWrzq5knLG6R8jMzKyR3L+aikG4DKMNQ1KM9boTZmYDRjBenBJoj+9YNzOz\nbE4iZmaWzUnEzMyyOYmYmVk2JxEzM8vmJGJmZtmcRMzMLJuTiJmZZXMSMTOzbE4iZmaWzUnEzMyy\nOYmYmVk2JxEzM8vmJGJmZtmcRMzMLJuTiJmZZXMSMTOzbE4iZmaWzUnEzMyytUwikt4r6VFJP5G0\nX9KfpvhSSXslHZB0j6RTU/y0tDyR1i8p1XV9ij8r6bJSfE2KTUjaXIrXbcPMzPpDlT2RN4BLIuI3\ngBXAGkmrgJuBLRGxDDgGbEzlNwLHIuJDwJZUDknLgfXAx4A1wG2S5kiaA9wKrAWWA1elsjRpw8zM\n+kDLJBKFX6TFU9IUwCXAvSm+Dbg8zQ+nZdL6SyUpxbdHxBsR8TwwAVyYpomIeC4i3gS2A8Npm0Zt\nmJlZH6h0TiTtMewDjgC7gX8AXo2IE6nIJLAwzS8EDgKk9ceBs8vxmm0axc9u0kZt/zZJGpM09nKV\nAZmZWUdUSiIR8VZErAAWUew5fLResfSqBus6Fa/Xv60RMRQRQwvqFTAzs65o6+qsiHgVeBhYBcyT\nNDetWgQcSvOTwGKAtP5M4Gg5XrNNo/grTdowM7M+UOXqrAWS5qX504HPAs8ADwFXpGIjwH1pfjQt\nk9Y/GBGR4uvT1VtLgWXAo8BjwLJ0JdapFCffR9M2jdowM7M+MLd1Ec4HtqWrqN4D7IiI+yU9DWyX\n9DXgceCOVP4O4DuSJij2QNYDRMR+STuAp4ETwDUR8RaApGuBXcAc4M6I2J/quq5BG2Zm1gdU/Id/\n9hiSYqzXnTAzGzCC8YgYanc737FuZmbZnETMzCybk4iZmWVzEjEzs2xOImZmls1JxMzMsjmJmJlZ\nNicRMzPL5iRiZmbZnETMzCybk4iZmWVzEjEzs2xOImZmls1JxMzMsjmJmJlZNicRMzPL5iRiZmbZ\nnETMzCxbyyQiabGkhyQ9I2m/pC+n+Fcl/UzSvjStK21zvaQJSc9KuqwUX5NiE5I2l+JLJe2VdEDS\nPZJOTfHT0vJEWr+kk4M3M7PpqbIncgL4o4j4KLAKuEbS8rRuS0SsSNMDAGndeuBjwBrgNklzJM0B\nbgXWAsuBq0r13JzqWgYcAzam+EbgWER8CNiSypmZWZ9omUQi4nBE/DjNvwY8AyxssskwsD0i3oiI\n54EJ4MI0TUTEcxHxJrAdGJYk4BLg3rT9NuDyUl3b0vy9wKWpvJmZ9YG2zomkw0mfAPam0LWSnpB0\np6T5KbYQOFjabDLFGsXPBl6NiBM18ZPqSuuPp/K1/dokaUzS2MvtDMjMzKalchKR9D7gr4A/jIif\nA7cDHwRWAIeBb04VrbN5ZMSb1XVyIGJrRAxFxNCCpqMwM7NOmlulkKRTKBLIX0TE9wEi4qXS+j8H\n7k+Lk8Di0uaLgENpvl78FWCepLlpb6NcfqquSUlzgTOBo836Os5KxFiVYZmZ2dvyzhS0TCLpHMQd\nwDMR8a1S/PyIOJwWfxN4Ks2PAt+V9C3gXwLLgEdTD5dJWgr8jOLk+3+MiJD0EHAFxXmSEeC+Ul0j\nwCNp/YMR8Y49kVqtS7RHqlZn1XLtlq2q0/2cOvv0bu1np3Wjn738HvVat97PqTqn3ofa2JR++Lde\nrrtR38vr682X68xRZU/kU8AXgCcl7UuxP6a4umoFxeGlF4DfKwYQ+yXtAJ6muLLrmoh4q+ikrgV2\nAXOAOyNif6rvOmC7pK8Bj1MkLdLrdyRNUOyBrM8bppmZdYMq/Md+oEhDEdHZw1neE5lev6bTvvdE\nWtfpPZHu1Pnu2xPReEQMtdsX37FuZmbZnETMzCybk4iZmWVzEjEzs2xOImZmls1JxMzMslW6Y32Q\nrGQ8/66ZBgIq3cxZtVy7ZavqdD/fvgLwXdrPTutGP3v5Peq1rr2fqc4ov9bOt9FuNz+jct0N+66T\nt+l0n7wnYmZm2ZxEzMwsm5OImZllcxIxM7NsTiJmZpbNScTMzLI5iZiZWbZZd5/IFL3zV3TNzKyO\nmMZNI7M2iUznTTEzs2p8OMvMzLI5iZiZWbaWSUTSYkkPSXpG0n5JX07xsyTtlnQgvc5PcUm6RdKE\npCckXVCqaySVPyBppBRfKenJtM0tUvHwq0ZtmJlZf6iyJ3IC+KOI+CiwCrhG0nJgM/DDiFgG/DAt\nA6wFlqVpE3A7FAkBuAG4CLgQuKGUFG5PZae2W5PijdowM7M+0DKJRMThiPhxmn8NeAZYCAwD21Kx\nbcDlaX4YuDsKe4B5ks4HLgN2R8TRiDgG7AbWpHVnRMQjERHA3TV11WvDzMz6QFvnRCQtAT4B7AXO\ni4jDUCQa4NxUbCFwsLTZZIo1i0/WidOkjdp+bZI0Jmns5XYGZGZm01I5iUh6H/BXwB9GxM+bFa0T\ni4x4ZRGxNSKGImJoQTsbmpnZtFRKIpJOoUggfxER30/hl9KhKNLrkRSfBBaXNl8EHGoRX1Qn3qwN\nMzPrA1WuzhJwB/BMRHyrtGoUmLrCagS4rxTfkK7SWgUcT4eidgGrJc1PJ9RXA7vSutckrUptbaip\nq14bZmbWB6rcsf4p4AvAk5L2pdgfA18HdkjaCLwIXJnWPQCsAyaA14EvAkTEUUk3AY+lcjdGxNE0\nfzVwF3A6sDNNNGnDzMz6gIoLomaPISnGet0JM7MBIxiPiKF2t/Md62Zmls1JxMzMsjmJmJlZNicR\nMzPL5iRiZmbZnETMzCybk4iZmWVzEjEzs2xOImZmls1JxMzMslV5dtZAGWcldPjBJxJ08ukwSg+/\n7/cnznSrn1Xfz3bKQe/6WbWuKYPwuffze9ntOqGod6r+2tiUem3X61OrWLmdquOpt02jvte206hP\nJw2uDd4TMTOzbLNuT2Ql49kZtZGA+j+dNZ36OlxnN3Srn1Xfz7bKVSzbjk5+7if9p28QPvc+fi+7\nXmeqN8qvtfMN2q7Xp1axk9qpOJ562zTse007jfqUy3siZmaWzUnEzMyyOYmYmVk2JxEzM8vmJGJm\nZtlaJhFJd0o6IumpUuyrkn4maV+a1pXWXS9pQtKzki4rxdek2ISkzaX4Ukl7JR2QdI+kU1P8tLQ8\nkdYv6dSgzcysM6rsidwFrKkT3xIRK9L0AICk5cB64GNpm9skzZE0B7gVWAssB65KZQFuTnUtA44B\nG1N8I3AsIj4EbEnlzMysj7RMIhHxI+BoxfqGge0R8UZEPA9MABemaSIinouIN4HtwLAkAZcA96bt\ntwGXl+ralubvBS5N5c3MrE9M55zItZKeSIe75qfYQuBgqcxkijWKnw28GhEnauIn1ZXWH0/l30HS\nJkljksZensaAzMysPblJ5Hbgg8AK4DDwzRSvt6cQGfFmdb0zGLE1IoYiYmhBs16bmVlHZSWRiHgp\nIt6KiH8G/pzicBUUexKLS0UXAYeaxF8B5kmaWxM/qa60/kyqH1YzM7MZkJVEJJ1fWvxNYOrKrVFg\nfbqyaimwDHgUeAxYlq7EOpXi5PtoRATwEHBF2n4EuK9U10iavwJ4MJU3M7M+0fIBjJK+B1wMnCNp\nErgBuFjSCorDSy8AvwcQEfsl7QCeBk4A10TEW6mea4FdwBzgzojYn5q4Dtgu6WvA48AdKX4H8B1J\nExR7IOvbGlkH840fBd/5et9tj4Kfqs86y4+C79Cj4Kfx5dRs+8/9kNThXxMxM5v9BOMRMdTudr5j\n3czMsjmJmJlZNicRMzPL5iRiZmbZnETMzCzbrPuN9XFWQoevz/Ilvp2v15f4Wif4Et8OXeJbrqBN\nsy6JrGS84/9aO/mj9m/X1+E6u6Fb/az6frZVrmLZdnTtc7fO6cZnnuqN8mvtfIO2631nWsVOaqfi\neOpt07DvNe006lMuH84yM7NsTiJmZpbNScTMzLI5iZiZWTYnETMzy+YkYmZm2ZxEzMwsm5OImZll\ncxIxM7NsTiJmZpbNScTMzLK1TCKS7pR0RNJTpdhZknZLOpBe56e4JN0iaULSE5IuKG0zksofkDRS\niq+U9GTa5hapePBVozZaKR7AGB2dROfr63Sd3Zi61c+qdbZTrpf97PX7+W6duvWZT9Vbfq2db9R2\nvXirWG177Yy9tp56fW/W5klTrohoOgGfBi4AnirFvgFsTvObgZvT/DpgJ8WjvVYBe1P8LOC59Do/\nzc9P6x4FPpm22QmsbdZGq6nzKcSTJ0+eZv8EjFX5G1s7tdwTiYgfAUdrwsPAtjS/Dbi8FL87CnuA\neZLOBy4DdkfE0Yg4BuwG1qR1Z0TEI1Fki7tr6qrXhpmZ9YnccyLnRcRhgPR6boovBA6Wyk2mWLP4\nZJ14szbeQdImSWOSxl7OHJCZmbWv0yfW6z2hPjLibYmIrRExFBFDC9rd2MzMsuUmkZfSoSjS65EU\nnwQWl8otAg61iC+qE2/WhpmZ9YncJDIKjKT5EeC+UnxDukprFXA8HYraBayWND9dZbUa2JXWvSZp\nVboqa0NNXfXaMDOzPtHy53ElfQ+4GDhH0iRwA/B1YIekjcCLwJWp+AMUV2hNAK8DXwSIiKOSbgIe\nS+VujIipk/VXA3cBp1NcnbUzxRu1YWZmfULpEtpZY0iKsV53wsxswAjGI2Ko3e18x7qZmWVzEjEz\ns2xOImZmls1JxMzMsjmJmJlZNicRMzPL5iRiZmbZnETMzCybk4iZmWVzEjEzs2xOImZmls1JxMzM\nsjmJmJlZNicRMzPL5iRiZmbZnETMzCybk4iZmWVzEjEzs2wtf2O9GUkvAK8BbwEnImJI0lnAPcAS\n4AXgtyLimCQB/53iN9hfB34nIn6c6hkB/kuq9msRsS3FV/LL319/APhyVP093w7+7K/U0eq6Vmen\nScVrr8beTjnobD8Hqc5efjc7/Vm2o1t1QlHvVP21sSn12q7Xp1axcjvtvu+19dTre207dftUHlib\nppVEks9ExCul5c3ADyPi65I2p+XrgLXAsjRdBNwOXJSSzg3AEBDAuKTRiDiWymwC9lAkkTXAzkq9\nmsabUisAOldd1+rstLe/Yz0ae1vlKpataqDq7OF3s9OfZTu6VmeqN8qvtfMN2q7Xp1axk9pp831/\nRz31+l7TTqM+5erG4axhYFua3wZcXorfHYU9wDxJ5wOXAbsj4mhKHLuBNWndGRHxSNr7uLtUl5mZ\n9YHpJpEA/rekcUmbUuy8iDgMkF7PTfGFwMHStpMp1iw+WSduZmZ9YrqHsz4VEYcknQvslvT3TcrW\n23mKjPg7Ky4S2CaADzTvr5mZddC09kQi4lB6PQL8NXAh8FI6FEV6PZKKTwKLS5svAg61iC+qE6/X\nj60RMRQRQwumMyAzM2tLdhKR9CuS3j81D6wGngJGgZFUbAS4L82PAhtUWAUcT4e7dgGrJc2XND/V\nsyute03SqnRl14ZSXWZm1gemczjrPOCvi7/vzAW+GxE/kPQYsEPSRuBF4MpU/gGKy3snKC7x/SJA\nRByVdBPwWCp3Y0QcTfNX88tLfHdS9cosMzObEap628WgGJJirNedMDMbMILxiBhqdzvfsW5mZtmc\nRMzMLJuTiJmZZXMSMTOzbJ14dlZfGWcldPjU+qA8PK7T/ADGwajTD2DsbJ3wLnwAY7mCNnlPxMzM\nsjmJmJlZtll3OAs6+hT4gauzG3o59nbaHpTPqNN19rqP3fgsO912br3l+mvbatR2vXirWL32cvpY\nL95oDJ1672ZdElnJOGP9/kMdZmZ9Jvevpg9nmZlZNicRMzPL5iRiZmbZnETMzCzbrDuxPs5KhJ/j\na2bWnrxT67MuicDg3BnrO9Y7Vw58d3kv6vQd6637NCh3rOde8jsrk4hZvxqUeyUG5T4R671ZmUR6\n/Y+ql3V2g282tH7lmw19s2HH+WZDM7P2+WZDMzObcX2fRCStkfSspAlJm3vdHzMz+6W+TiKS5gC3\nAmuB5cBVkpb3tldmZjalr5MIcCEwERHPRcSbwHZguMd9MjOzpN9PrC8EDpaWJ4GLagtJ2gRsSotv\nCJ6agb71yjnAK73uRBfN5vHN5rGBxzfoPpyzUb8nkXoXDLzjdpyI2ApsBZA0FhFD3e5Yr3h8g2s2\njw08vkEnKetRH/1+OGsSWFxaXgQc6lFfzMysRr8nkceAZZKWSjoVWA+M9rhPZmaW9PXhrIg4Iela\nYBcwB7gzIva32Gxr93vWUx7f4JrNYwOPb9BljU/R708BNDOzvtXvh7PMzKyPOYmYmVm2gU0irR6H\nIuk0Sfek9XslLZn5XuapMLavSHpa0hOSfijpX/Win7mqPspG0hWSQtJAXVZZZXySfit9hvslfXem\n+zgdFb6fH5D0kKTH03d0XS/6mUPSnZKOSKp7r5kKt6SxPyHpgpnu43RUGN/n07iekPR3kn6jZaUR\nMXATxUn2fwB+FTgV+AmwvKbM7wN/lubXA/f0ut8dHNtngH+R5q8elLFVHV8q937gR8AeYKjX/e7w\n57cMeByYn5bP7XW/Ozy+rcDVaX458EKv+93G+D4NXAA81WD9OmAnxT1sq4C9ve5zh8f3b0rfy7VV\nxjeoeyJVHocyDGxL8/cCl0oD8esQLccWEQ9FxOtpcQ/F/TODouqjbG4CvgH8v5nsXAdUGd9/Am6N\niGMAEXFkhvs4HVXGF8AZaf5MBujeroj4EXC0SZFh4O4o7AHmSTp/Zno3fa3GFxF/N/W9pOLflkFN\nIvUeh7KwUZmIOAEcB86ekd5NT5WxlW2k+J/RoGg5PkmfABZHxP0z2bEOqfL5/Rrwa5L+VtIeSWtm\nrHfTV2V8XwV+W9Ik8ADwBzPTtRnR7r/PQVbpb0tf3yfSRJXHoVR6ZEofqtxvSb8NDAH/rqs96qym\n45P0HmAL8Dsz1aEOq/L5zaU4pHUxxf/0/kbSxyPi1S73rROqjO8q4K6I+KakTwLfSeP75+53r+sG\n9e9KWyR9hiKJ/NtWZQd1T6TK41DeLiNpLsVudbPd1H5R6VEvkj4L/AnwuYh4Y4b61gmtxvd+4OPA\nw5JeoDjuPDpAJ9erfjfvi4h/iojngWcpksogqDK+jcAOgIh4BHgvxcMLZ4NZ/ygmSb8OfBsYjoj/\n26r8oCaRKo9DGQVG0vwVwIORzhb1uZZjS4d7/gdFAhmk4+nQYnwRcTwizomIJRGxhOK47OciIuvh\ncD1Q5bv5vygujkDSORSHt56b0V7mqzK+F4FLASR9lCKJvDyjveyeUWBDukprFXA8Ig73ulOdIukD\nwPeBL0TETytt1OurBaZxlcE64KcUV4r8SYrdSPEHB4ov7l8CE8CjwK/2us8dHNv/AV4C9qVptNd9\n7uT4aso+zABdnVXx8xPwLeBp4Elgfa/73OHxLQf+luLKrX3A6l73uY2xfQ84DPwTxV7HRuBLwJdK\nn92taexPDuB3s9X4vg0cK/1tGWtVpx97YmZm2Qb1cJaZmfUBJxEzM8vmJGJmZtmcRMzMLJuTiJnZ\nAGv1UMWaslsk7UvTTyVN+wZXX51lZjbAJH0a+AXFM70+3sZ2fwB8IiJ+dzrte0/EzGyARZ2HKkr6\noKQfSBqX9DeSPlJn06so7huZlkF9dpaZmTW2leIGwgOSLgJuAy6ZWpl+g2gp8OB0G3ISMTObRSS9\nj+J3Qf6y9OsXp9UUWw/cGxFvTbc9JxEzs9nlPcCrEbGiSZn1wDWdaszMzGaJiPg58LykK+Htn/R9\n+2duJX0YmA880on2nETMzAaYpO9RJIQPS5qUtBH4PLBR0k+A/Zz865NXAdujQ5fm+hJfMzPL5j0R\nMzPL5iRiZmbZnETMzCybk4iZmWVzEjEzs2xOImZmls1JxMzMsv1/qWaAzJPZfhAAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116edcac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# for traditional kdtree\n",
    "visualize_kdnodes_and_query(traditional_kdnodes, testing_set, domains_)\n",
    "visualize_kdnodes_and_query(traditional_variance_kdnodes, testing_set, domains_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 821,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# # === Baseline Performance 2D\n",
    "# #\n",
    "# # @dataset: the ordered dataset, contains only the columns to build, and ordered\n",
    "# # @threshold: the block/page size\n",
    "# #\n",
    "# # return @minmaxnodes: contains all the min max info of all the build dimensions (in the build_dim order)\n",
    "# #\n",
    "# def MinMaxOrdering(dataset, threshold):\n",
    "#     minmaxnodes = []\n",
    "    \n",
    "#     # currently sorted by A, B\n",
    "#     dataset = dataset[np.lexsort((dataset[:,1], dataset[:,0]))]\n",
    "    \n",
    "#     count = 0;\n",
    "#     total_size = len(dataset)\n",
    "#     block_amount = math.ceil(total_size / threshold)\n",
    "    \n",
    "#     lower = 0;\n",
    "#     upper = lower + threshold\n",
    "    \n",
    "#     # concat the results of each sorted nodes\n",
    "#     for i in range(block_amount):\n",
    "#         if upper > total_size:\n",
    "#             upper = total_size\n",
    "#         sub_dataset = dataset[lower:upper,:]\n",
    "        \n",
    "#         minmax_dim0 = [min(sub_dataset[:,0]),max(sub_dataset[:,0])]\n",
    "#         minmax_dim1 = [min(sub_dataset[:,1]),max(sub_dataset[:,1])]\n",
    "#         minmaxnodes.append([minmax_dim0, minmax_dim1])\n",
    "        \n",
    "#         lower += threshold\n",
    "#         upper += threshold\n",
    "        \n",
    "#     return minmaxnodes\n",
    "\n",
    "\n",
    "# def QueryMinMaxOrdering(queries, minmaxnodes):\n",
    "    \n",
    "#     counts = []\n",
    "#     count_single_query = 0;\n",
    "    \n",
    "#     # number of dimensions\n",
    "#     dims = int(len(queries[0]))\n",
    "    \n",
    "#     # for each query\n",
    "#     for i in range(len(queries)):\n",
    "        \n",
    "#         count_single_query = 0\n",
    "        \n",
    "#         # check for intersection for each kdnode\n",
    "#         for j in range(len(minmaxnodes)):\n",
    "            \n",
    "#             # for each dimension\n",
    "#             intersection_tag = True\n",
    "#             for k in range(dims):\n",
    "                \n",
    "#                 # an intersection holds if it intersecs in all dimensions\n",
    "#                 if queries[i][k][0] >= minmaxnodes[j][k][1] or queries[i][k][1] <= minmaxnodes[j][k][0]:\n",
    "#                     intersection_tag = False\n",
    "#                     break\n",
    "                \n",
    "#             # if the query intersect with this kdnode\n",
    "#             if intersection_tag:\n",
    "#                 count_single_query += 1\n",
    "            \n",
    "#         counts.append(count_single_query)\n",
    "    \n",
    "#     #print(\"blocks IO: \", counts)\n",
    "#     print(\"blocks IO(average): \", statistics.mean(counts))\n",
    "    \n",
    "# minmaxnodes = MinMaxOrdering(dataset, block_size)\n",
    "# print('min max nodes amount: ',len(minmaxnodes))\n",
    "# QueryMinMaxOrdering(testing_set, minmaxnodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 822,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "min max nodes amount:  375\n",
      "blocks IO(average):  337.84422110552765\n"
     ]
    }
   ],
   "source": [
    "# === Baseline Performance 3D\n",
    "#\n",
    "# @dataset: the ordered dataset, contains only the columns to build, and ordered\n",
    "# @threshold: the block/page size\n",
    "#\n",
    "# return @minmaxnodes: contains all the min max info of all the build dimensions (in the build_dim order)\n",
    "#\n",
    "def MinMaxOrdering(dataset, threshold):\n",
    "    minmaxnodes = []\n",
    "    \n",
    "    # currently sorted by A, B, C\n",
    "    dataset = dataset[np.lexsort((dataset[:,2], dataset[:,1], dataset[:,0]))]\n",
    "    \n",
    "    count = 0;\n",
    "    total_size = len(dataset)\n",
    "    block_amount = math.ceil(total_size / threshold)\n",
    "    \n",
    "    lower = 0;\n",
    "    upper = lower + threshold\n",
    "    \n",
    "    # concat the results of each sorted nodes\n",
    "    for i in range(block_amount):\n",
    "        if upper > total_size:\n",
    "            upper = total_size\n",
    "        sub_dataset = dataset[lower:upper,:]\n",
    "        \n",
    "        minmax_dim0 = [min(sub_dataset[:,0]),max(sub_dataset[:,0])]\n",
    "        minmax_dim1 = [min(sub_dataset[:,1]),max(sub_dataset[:,1])]\n",
    "        minmax_dim2 = [min(sub_dataset[:,2]),max(sub_dataset[:,2])]\n",
    "        minmaxnodes.append([minmax_dim0, minmax_dim1, minmax_dim2])\n",
    "        \n",
    "        lower += threshold\n",
    "        upper += threshold\n",
    "        \n",
    "    return minmaxnodes\n",
    "\n",
    "\n",
    "def QueryMinMaxOrdering(queries, minmaxnodes):\n",
    "    \n",
    "    counts = []\n",
    "    count_single_query = 0;\n",
    "    \n",
    "    # number of dimensions\n",
    "    dims = int(len(queries[0]))\n",
    "    \n",
    "    # for each query\n",
    "    for i in range(len(queries)):\n",
    "        \n",
    "        count_single_query = 0\n",
    "        \n",
    "        # check for intersection for each kdnode\n",
    "        for j in range(len(minmaxnodes)):\n",
    "            \n",
    "            # for each dimension\n",
    "            intersection_tag = True\n",
    "            for k in range(dims):\n",
    "                \n",
    "                # an intersection holds if it intersecs in all dimensions\n",
    "                if queries[i][k][0] >= minmaxnodes[j][k][1] or queries[i][k][1] <= minmaxnodes[j][k][0]:\n",
    "                    intersection_tag = False\n",
    "                    break\n",
    "                \n",
    "            # if the query intersect with this kdnode\n",
    "            if intersection_tag:\n",
    "                count_single_query += 1\n",
    "            \n",
    "        counts.append(count_single_query)\n",
    "    \n",
    "    #print(\"blocks IO: \", counts)\n",
    "    print(\"blocks IO(average): \", statistics.mean(counts))\n",
    "    \n",
    "minmaxnodes = MinMaxOrdering(dataset, block_size)\n",
    "print('min max nodes amount: ',len(minmaxnodes))\n",
    "QueryMinMaxOrdering(testing_set, minmaxnodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 823,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "blocks IO(average):  7.517587939698492\n",
      "blocks IO(average):  9.809045226130653\n",
      "blocks IO(average):  9.658291457286433\n",
      "blocks IO(average):  15.864321608040202\n",
      "blocks IO(average):  375\n"
     ]
    }
   ],
   "source": [
    "# === Ignore PK ===¶\n",
    "for i in range(len(testing_set)):\n",
    "    testing_set[i][0][0] = domains_[0][0]\n",
    "    testing_set[i][0][1] = domains_[0][1]\n",
    "Query(testing_set, kdnodes)\n",
    "QueryDense(testing_set, dense_kdnodes)\n",
    "Query(testing_set, traditional_kdnodes)\n",
    "QueryDense(testing_set, dense_traditional_kdnodes)\n",
    "QueryMinMaxOrdering(testing_set, minmaxnodes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}