Assignment2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os import getcwd\n",
    "TRAIN = os.path.join(getcwd(), 'MinHash/data/articles_1000.train')\n",
    "TRUTH = os.path.join(getcwd(), 'MinHash/data/articles_1000.truth')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mapping\n",
    "# cache the index of every docid for O(1) query later\n",
    "def docid_map(docids):\n",
    "    docid_col_map = dict()\n",
    "    for i in range(0, len(docids)):\n",
    "        docid_col_map[docids[i]] = i\n",
    "    return docid_col_map    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parse_data(filename):\n",
    "    if not filename:\n",
    "        raise ValueError(\"filename is empty.\")\n",
    "    \n",
    "    import string\n",
    "    \n",
    "    with open(filename) as f:\n",
    "        content = f.readlines()  # read lines from filename\n",
    "        content = [x.strip() for x in content]\n",
    "        exclude = set(string.punctuation)\n",
    "        tuples = list()\n",
    "        \n",
    "    for line in content:\n",
    "        id, article = line.split(\" \", 1)\n",
    "        article = ''.join(ch for ch in article if ch not in exclude).replace(' ', '').lower()\n",
    "        tuples.append((id, article))  # construct tuple of id and text\n",
    "    return tuples "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shingle_document(string, k):\n",
    "    if not string:\n",
    "        raise ValueError(\"String is empty.\")\n",
    "    if not k:\n",
    "        raise ValueError(\"k is empty.\")\n",
    "    \n",
    "    import binascii\n",
    "    \n",
    "    unique_hash = set() # initialize set data structure\n",
    "    for i in range(0, len(string) - k + 1):\n",
    "        unique_hash.add(binascii.crc32(string[i:i+k].encode(\"utf-8\"))) # has substring into 32-bit integer\n",
    "        \n",
    "    return unique_hash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################\n",
    "# test parse_data()\n",
    "dd = parse_data(TRAIN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shingled_list(document):\n",
    "    shingled_list = list()\n",
    "    for d in document:\n",
    "        shingled_list.append((d[0], shingle_document(d[1], 10))) # shingle_document()\n",
    "    return shingled_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################\n",
    "# test shingled_list()\n",
    "shingled_dd = shingled_list(dd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jaccard(a, b):\n",
    "    if not a:\n",
    "        raise ValueError(\"a is empty.\")\n",
    "    if not b:\n",
    "        raise ValueError(\"b is empty.\")\n",
    "        \n",
    "    a = set(a)\n",
    "    b = set(b)\n",
    "\n",
    "    return len(a.intersection(b))/len(a.union(b)) # return raio of union and intersection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def similarity(filename, k):\n",
    "    if not filename:\n",
    "        raise ValueError(\"filename is empty.\")\n",
    "    if not k:\n",
    "        raise ValueError(\"k is empty.\")\n",
    "        \n",
    "    tuples = parse_data(filename)  # parse a file with data\n",
    "    tuples_js = list()\n",
    "\n",
    "    n = len(tuples)\n",
    "\n",
    "    for i in range(0, n - 1):\n",
    "        id_i, article_i = tuples[i]\n",
    "        if k > len(article_i):\n",
    "            raise ValueError(\"k is out of bound.\")\n",
    "\n",
    "        for j in range(i + 1, n):\n",
    "            id_j, article_j = tuples[j]\n",
    "            if k > len(article_j):\n",
    "                raise ValueError(\"k is out of bound.\")\n",
    "\n",
    "            s = jaccard(shingle_document(article_i, k), shingle_document(article_j, k)) # shingle_document()\n",
    "            tuples_js.append((id_i, id_j, s))\n",
    "    \n",
    "    return tuples_js # return a list of tuples, each tuple contains: (id1, id2, s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################\n",
    "# test similarity()\n",
    "# similarity(TRAIN, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# inversion\n",
    "def _invert_shingles(shingled_documents): # input a list of tuples of form (id, [items])\n",
    "    if not shingled_documents:\n",
    "        raise ValueError(\"filename is empty\")\n",
    "        \n",
    "    tuples = list() # initialize list for tuples\n",
    "\n",
    "    for row in shingled_documents:\n",
    "        for item in row[1]:\n",
    "            tuples.append((item, row[0]))\n",
    "    return tuples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    " def getKey(item):\n",
    "        return item[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sorting\n",
    "def invert_shingles(shingled_documents):\n",
    "    if not shingled_documents:\n",
    "        raise ValueError(\"filename is empty\")\n",
    "    tuples = list()\n",
    "    ids = list() # initialize list for document ids\n",
    "    \n",
    "    for row in shingled_documents:\n",
    "        ids.append(row[0])\n",
    "        \n",
    "    tuples = sorted(_invert_shingles(shingled_documents), key = getKey) # invert_shingles_(), getKey()\n",
    "\n",
    "    return tuples, ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "############################\n",
    "# test invert_shingles()\n",
    "inverted = invert_shingles(shingled_dd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate hash functions\n",
    "import random        \n",
    "\n",
    "def make_random_hash_fn(p=2**33-355, m=4294967295):\n",
    "    a = random.randint(1,p-1)\n",
    "    b = random.randint(0, p-1)\n",
    "    return lambda x: ((a * x + b) % p) % m\n",
    "\n",
    "def make_hashes(num_hash):\n",
    "    hash_fns = list() # initialize list for hash functions\n",
    "    \n",
    "    for i in range(0, num_hash):\n",
    "        hash_fns.append(make_random_hash_fn()) # make_random_hash_fn()\n",
    "        \n",
    "    return hash_fns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_minhash_signature(shingled_data, num_hashes):\n",
    "    inv_index, docids = invert_shingles(shingled_data) # invert_shingles()\n",
    "    num_docs = len(docids)\n",
    "\n",
    "    # print(inv_index)\n",
    "    # print(type(inv_index))\n",
    "    # print(len(inv_index))\n",
    "    # print(len(inv_index[0]))\n",
    "    \n",
    "    docid_col_map = docid_map(docids) # docid_map\n",
    "    # print(docid_col_map)\n",
    "    \n",
    "    unique_items = {item[0] for item in inv_index}\n",
    "    sorted_unique_items = sorted(unique_items)\n",
    "  \n",
    "    # initialize the signature matrix with infinity in every entry\n",
    "    sigmatrix = np.full([num_hashes, num_docs], np.inf)\n",
    "    \n",
    "    # generate hash functions\n",
    "    hash_funcs = make_hashes(num_hashes) # make_hashes()\n",
    "    \n",
    "    # iterate over each non-zero entry of the characteristic matrix\n",
    "    # docid locates in both characteristics matrix and signature matrix\n",
    "    # row retains the row index value\n",
    "    # use the same hash functions calculate through rows with the same docid and put the minimum value into sigmatrix\n",
    "\n",
    "    for row, docid in inv_index:\n",
    "        # i = row_map[row] # index of rows\n",
    "        j = docid_col_map[docid] # index of columes\n",
    "        for n in range(0, num_hashes):\n",
    "            sigmatrix[n][j]= min(hash_funcs[n](row), sigmatrix[n][j]) # why don't we add % len(sorted_unique_items)?\n",
    "    \n",
    "    return sigmatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################\n",
    "# test make_minhash_signature()\n",
    "minhash_sigmat = make_minhash_signature(shingled_dd, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minhash_similarity(id1, id2, minhash_sigmat):\n",
    "\n",
    "    similarity = 0\n",
    "    \n",
    "    for i in range(0, len(minhash_sigmat)):\n",
    "        if  minhash_sigmat[i][id1] == minhash_sigmat[i][id2]:\n",
    "            similarity = similarity + 1\n",
    "    \n",
    "    return similarity/len(minhash_sigmat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################\n",
    "# test minhash_similarity()\n",
    "# minhash_similarity(0, 0, minhash_sigmat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minhash_similarity_all(shingled_doc, num_hashes):\n",
    "    inv_index, docids = invert_shingles(shingled_doc) # invert_shingles()\n",
    "    \n",
    "    MinHash_Matrix = np.zeros((len(docids), len(docids)))\n",
    "    tuples = list()\n",
    "    \n",
    "    minhash_sigmat = make_minhash_signature(shingled_dd, num_hashes) # make_minhash_signature()\n",
    "    \n",
    "    docid_col_map = docid_map(docids) # docid_map()\n",
    "    # print(docid_col_map)\n",
    "    \n",
    "    for d1 in docids:\n",
    "        i = docid_col_map[d1]\n",
    "        for d2 in docids:\n",
    "            j = docid_col_map[d2]\n",
    "            if i < j:\n",
    "                MinHash_Matrix[i][j] = minhash_similarity(i, j, minhash_sigmat) # minhash_similarity()\n",
    "                tuples.append((d1, d2, MinHash_Matrix[i][j]))\n",
    "            \n",
    "    return MinHash_Matrix, tuples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################\n",
    "# test minhash_similarity_all()\n",
    "# minhash_similarity_all = minhash_similarity_all(shingled_dd, 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# minhash_similarity_matrix = minhash_similarity_all[0] \n",
    "# minhash_similarity_tuples = minhash_similarity_all[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################\n",
    "# check for consistency\n",
    "# print(minhash_similarity_matrix)\n",
    "# print(minhash_similarity_tuples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "jaccard_similarity_tuples = similarity(TRAIN, 10)\n",
    "js = pd.DataFrame(jaccard_similarity_tuples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "js.columns = [\"id1\", \"id2\", \"js\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_hashes = 100\n",
    "\n",
    "dt = list()\n",
    "    \n",
    "for k in range(1, num_hashes + 1):\n",
    "    \n",
    "    # compare jaccard_similarity_tuples and minhash_similarity_tuples    \n",
    "    minhash_similarity_tuples = minhash_similarity_all(shingled_dd, k)[1]\n",
    "    mh = pd.DataFrame(minhash_similarity_tuples)\n",
    "    mh.columns = [\"id1\", \"id2\", \"mh\"]\n",
    "    \n",
    "    # all_tuples = jaccard_similarity_tuples + minhash_similarity_tuples \n",
    "    # data[k-1] = (len(all_tuples) - len(set(all_tuples)))/len(jaccard_similarity_tuples)\n",
    "       \n",
    "    # use min squared error\n",
    "    df = pd.merge(js, mh, on = [\"id1\", \"id2\"])\n",
    "    df[\"diff\"] = (df.js - df.mh) ** 2\n",
    "    mse = df[\"diff\"].sum()\n",
    "    rmse = math.sqrt(mse/len(df))\n",
    "    fraction = np.mean(js.js == mh.mh)\n",
    "    num = k\n",
    "    dt.append((num, rmse, fraction))\n",
    "\n",
    "dt = pd.DataFrame(dt, columns = [\"num\", \"rmse\", \"fraction\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xe4VNX1//H3h66IYMEGKgrEggELYiwx9hZ7iaLG3mJMYoyJmhijxhITjSbRmK8t2BKx/7ArauwFsKCARsACYkFEBKRdWL8/1pncc4cp517u3Dbr9Tz3mTl19pmZe9asvffZR2ZGCCGEUEq75i5ACCGEli+CRQghhLIiWIQQQigrgkUIIYSyIliEEEIoK4JFCCGEsiJYtAGSfiTpM0lzJK3S3OWpL0k7SJrawG1/LemGBm57hKTHU9MmqV8D97VO8v63b8j2lSRpmKSLlnEfLeL4luW7UmBf35X0bgO3rfN+SPqPpBOWoSzjJO3Q0O2bQgSLPJI+kDQv+SJ8mvyjrZBaPiw5qeybt91VyfxjkulOkq6QNDXZ1/uSrizyOrm/qxtQ3o7An4HdzGwFM5tRYJ0GnwRbAkn7SXpD0teSvpD0pKQ+AGZ2iZk16J/UzG43s90ao4xm9lHy/i9OyrxMJ4+mIumY5Pvx57z5+yfzh8HSx1dmn32SbTvkzV/moFVfkgZIelzSTElfSRojaS8AM3vOzDZoyH7r835k3N8AM/tPUubzJd3WGPttTBEsCtvHzFYANgU2A87JW/5f4OjcRPJPcQgwKbXOOcBgYAjQDdgReL3Q66T+TmtAWVcHugDjGrBtk8s/gWRYvx9wC/ALoDuwHvB3YEnjl65h6ntMLdAk4NC84zgK/563dg8AT+D/J6sBPwW+btYSpbSm704EixLM7FPgMTxopD0AbCtppWR6D2As8GlqnS2B+8xsmrkPzOyWhpRDUuckc5mW/F2VzPsWkEujv5L0VIZ9DZH0UvIr6xNJV0vqlFo+QNITkr5MqrZ+ncxvn1T5TJI0O/mFtnay7C+SpiS//MdI+m5qf+dLulvSbZK+Bo6RtFzyK3OmpPHJe1XMpsD7ZvZk8j7ONrN7zOyj1P5vS57nftEem5RnpqRTJG0paWxyzP/L3pJf1c8XeZ++L+n15JimSDo/tSz3OsdL+gh4Kv1rWtLFwHeBq3MZo6RrJF2R9xoPSDq9yOuXe0/vlHRL8lmMkzQ4tXwzSa8ly4bjPyZK+RR4C9g92X5lYBtgRIFj7pBM/0fS7yW9kLzO45JWLfM6+cd4lzx7nyXpWUkDUsv2kjQ+2ffHks7M2/YXkj5PvsPHFtn/qviPi+vNbGHy94KZPZ8sr1OlJc/2f5l8V+ZKulHS6pIeScoxMvc/n/9+5L1uX0lPSZohz4Rvl9Qj73XOkjQWmJt8Zz6QtIukPYBf48F7jqQ3JR0iaUyB47+/Pu/3sopgUYKk3sCewMS8RfPxf6TDkumj8F+/aS8DZ0g6VdK3JWkZivIb4Dv4iXMQnq2ca2b/BXL/YD3MbKcM+1oM/BxYFdga2Bk4FUBSN2Ak8CiwFtAPeDLZ7gxgKLAXsCJwHPBNsmxUUraVgX8Bd0lKn6D2A+4GegC3A78D+iZ/u5PK0gp4DdhQ0pWSdlSqSrCErYD+wKHAVfj7twv+Xv1A0vcy7GMu/rn2AL4P/EjS/nnrfA/YKDmG/zGz3wDPAaelMsabgaGS2sH/TmQ7A/8u8vrl3tN9gTuS8o0Ark722wm4H7g12fYu4KAMx3tLcrzg3+v/Bywos83hwLH4L/ZOwJmlV1/KI/jntBr+Od+eWnYjcLKZdQM2AdI/hNbAs8xewPHANar94ZY2A//fvU1erbZ6hjIdBOwKfAvYJynjr/H/l3Z4ZlKOgEvx/6GNgLWB8/PWGYp/r3qYWU1uppk9ClwCDE++O4Pwz3c9SRultj8S/4ybjpnFX+oP+ACYA8wGDD9Z9kgtHwZcBGwHvIR/aT8DlgOeB45J1msP/Bh4Af+nmwYcXeB1vkr9nVikTJOAvVLTuwMfJM/7JOXsUOKYDOhXZNnpeAYE/gV+vch67wL7ZXwPZwKDkufnA8/mLZ8M7JGaPgmYWmJ/3wHuBKbjgXoYsEJq/7flvRe9UtvOAA5NTd8DnJ48PwZ4PuP7dBVwZd7rrJ9aXudzAP4DnJC3jwnArsnz04CH6/G9zH9PR6aWbQzMS55vn3zXlFr+InBRkf0ek3xvl8O/x93xHzrb4t/zYSWO79zUfk4FHs1b96u8v4UlytEj2aZ7Mv0RcDKwYt56OwDzSH3fgc+B7xTZb288kE7Cqy6fBfqn9jU1te4HwBF535VrU9M/Ae7P+nmnttuf1P9V8jrH5a3zAbBL/nc6tfxa4OLk+YDk+9A56/enMf4isyhsf/NfNDsAG+K/KuowT2V7AucCD5rZvLzli83sGjPbFv9HuBi4Ke/Xwf5m1iP1d32R8qwFfJia/jCZV2+SviXpwST9/xr/FZM7vrWp2+6SVnRZkhJPSKoTvsJPOOn3bEreJmvlzfuQEszsZTP7gZn1xKt3tsezhWI+Sz2fV2C6bHYiaStJT0uaLmkWcApLfw/yj6ucm/FfhFDml2GG9zRd5fkN0CWpElkL+NiSs0qi5PsLkHx/H8K/z6ua2QsZjie/DPnv66rp7zeeIeWOr72kP8irNb/GT5ZQe4wH4Vnsh5KekbR1ar8zLPVrvMhr545rqpmdZmZ9gXXxjLFUdXBjfHdWk3RHUn32NXAbjfPdOTypofghcKeZlcv8GlUEixLM7Bn8V+zlRVa5DW94LdkWYWbzzOwa/NfAxg0oyjT8i56zTjKvIa4F3sF/Xa2Ip9i5KrIpeNVQIQWXJXXpZwE/AFZKTgqzUvsE/wWW9gkefHLWyVp4MxsF3ItXTVTSv/D0f20z6w78g7rHBEsfV7lltwH7SRqEV08UrHPO+J4W8wnQK6/aM+v7m+tI0BTVG4fj1ZO74IGwTzJf4J+zme2HV1Hdj2eWy8TMpgDXUPnvzqX45z8w+R87kmX87pjZy3hm9l38vWvaKigiWGRxFbCrpPxGboC/4vWbz+YvkHR60oC2XNKAdTTeKyq/R1QW/wbOldQzqes+Dz/xNEQ3vDfIHEkbAj9KLXsQWCMpe2dJ3SRtlSy7Afi9pP5yA+XXdHQDavAqog6SzsPbNEq5EzhH0kpJu9BPiq0oaTtJJ0paLZneEK+vf7neR14/3YAvzWy+pCH4P2h9fAasn55hZlPxtohbgXvys9G8167ve5rzUrLtT5Pv3YF4G1cWz+Df579lXH9ZdMOrZ2cAy+MZLvC/budHSOpuZovw72u9u6gm368LJPWT1C753zmOpvnuzME7nfQCflnP7T8D+uTat1JuwavUapKajSYVwaIMM5uOf0i/LbDsS0t66RTYdB5wBZ6qf4G3XxxkZpNT6zygutdZ3FekGBcBo/EeV2/hjYH17a+eK+OZ+IlvNnA9MDx1PLPxk8U+Sbnfw7v8gl/LcSfwOP7PeyNez/0Y3gj4X7y6Yz7lU+wLknXfT/ZX6lfSV3hweEvSHLzx/T7gj2VeY1mdClwoaTYenOv7y/YvwMHyHll/Tc2/Gfg2pY+5Ie8pAGa2EDgQb4uYiTfy35txW0u+z19mWX8Z3YIf28fAeJY+gf8Q+CCpxjmF2uq7+liIZywj8e/s23iAOqZBJc7uAmBzPBt8iIzvf8pdyeMMSa+l5t+KZ0VNnlVA0ggW2i5JK+Jf2pXM7KvmLk+1k7Q9nhX2MbMWc61IaPkkLYc35m9uZu819etHZtH2HQpMikDR/ORX2/8MuCECRWiAHwGjmiNQALSaqwdD/Ul6Ee+J1eKHnWjrkl5wo4E38WsTQshM0gd4I3n+tT5NV4aohgohhFBOVEOFEEIoq81UQ6266qrWp0+f5i5GCCG0KmPGjPkiueC1pDYTLPr06cPo0aObuxghhNCqSCp7hT9ENVQIIYQMIliEEEIoK4JFCCGEsiJYhBBCKCuCRQghhLIiWIQQQigrgkUIIYSyIlhk9MADMHVq+fVCCKEtimCRwZIlcOCBcN11zV2SEEJoHhEsMpg3D2pq/DGEEKpRBIsM5s71x0WLmrccIYTQXCJYZJALFgsXNm85QgihuUSwyCAyixBCtYtgkUFkFiGEahfBIoPILEII1S6CRQaRWYQQql0EiwwiswghVLsIFhlEZhFCqHYRLDKIzCKEUO0iWGQQwSKEUO0iWGQQ1VAhhGoXwSKDyCxCCNUugkUGkVmEEKpdBIsMIrMIIVS7CBYZRGYRQqh2ESwyiMwihFDtIlhkEJlFCKHaRbDIIDKLEEK1i2CRQWQWIYRqF8Eig8gsQgjVrqLBQtIekt6VNFHS2QWWd5Y0PFn+iqQ+ecvXkTRH0pmVLGc5kVmEEKpdxYKFpPbANcCewMbAUEkb5612PDDTzPoBVwKX5S2/EnikUmXMwqw2WJjB4sXNWZoQQmgelcwshgATzWyymS0E7gD2y1tnP+Dm5PndwM6SBCBpf2AyMK6CZSxrwQJYsgRWXNGnI7sIIVSjSgaLXsCU1PTUZF7BdcysBpgFrCKpK3AWcEGpF5B0kqTRkkZPnz690QqelssqevTwx2i3CCFUo0oGCxWYZxnXuQC40szmlHoBM7vOzAab2eCePXs2sJil5YLFSiv5Y2QWIYRq1KGC+54KrJ2a7g1MK7LOVEkdgO7Al8BWwMGS/gj0AJZImm9mV1ewvAVFZhFCCJUNFqOA/pLWAz4GDgMOz1tnBHA08BJwMPCUmRnw3dwKks4H5jRHoAD45ht/jMwihFDNKhYszKxG0mnAY0B74CYzGyfpQmC0mY0AbgRulTQRzygOq1R5GioyixBCqGxmgZk9DDycN++81PP5wCFl9nF+RQqXUX6bRQSLEEI1iiu4y8jPLKIaKoRQjSJYlBGZRQghRLAoKzKLEEKIYFFWNHCHEEIEi7JywaJ7d3+MzCKEUI0iWJQxdy4stxx06eLTkVmEEKpRBIsy5s6Frl2hY0efjswihFCNSgYLSe0l3dZUhWmJcsGiUyefjswihFCNSgYLM1sM9JTUqYnK0+JEZhFCCNmu4P4AeEHSCGBubqaZ/blShWpJIrMIIYRswWJa8tcO6FbZ4rQ8kVmEEEKGYGFmFwBI6uaTpe8x0dotXgzt29dOz50La64ZmUUIobqV7Q0laRNJrwNvA+MkjZE0oPJFa3offeQX3z3xRO28yCxCCCFb19nrgDPMbF0zWxf4BXB9ZYvVPO6+G+bMgddfr50XbRYhhJAtWHQ1s6dzE2b2H6BrxUrUjO65xx+nTq2dF5lFCCFka+CeLOm3wK3J9JHA+5UrUvP45BN48UV/XihYdEjeqcgsQgjVKEtmcRzQE7g3+VsVOLaShWoO993nj+usUxssamo8k+jaFSTPLiKzCCFUo5KZhaT2wK/N7KdNVJ5mc++98K1vwXbbwSOP+LzcIIJdk0q3Tp0iswghVKcsV3Bv0URlaTYzZsB//gMHHQS9e8Onn3pQyA8WHTtGsAghVKcsbRavJ1dv30XdK7jvrVipmtiIEX59xYEHek8oM2/DWLDAl6czi6iGCiFUoyzBYmVgBrBTap7h7Rdtwj33wLrrwhZbwPTpPu/jj31ocojMIoQQsrRZjDWzK5uoPE3u66/9Irwf/9gbsXv39vlTp8Jaa/nzdLCIzCKEUI2ytFns20RlaRZvv+3VSwce6NPpYBEN3CGE4LJUQ70o6WpgOHXbLF6rWKma0DbbeNVT7grtHj1g+eU9WPTp4/MiswghVLsswWKb5PHC1DyjbhtGq5a7ZSrUVkVFZhFCCLWyjDq7Y1MUpCXp1atwsIjMIoRQrYq2WUi6KvX8Z3nLhlWwTM0uMosQQqirVAP39qnnR+ctG1iBsrQYvXvDtGkwe7ZPR2YRQqh2pYKFijxv83r39nGh3n/fs4ncIIKRWYQQqlWpNot2klbCA0rueS5otC++WeuX6z777rveMyonMosQQrUqFSy6A2OoDRDprrJWsRK1AOlg0TV1547ILEII1aposDCzPk1YjhalVy9/nDkTevasnR+ZRQihWmW5n0XV6dmz9s54kVmEEEIEi4LatavNLtLBIjKLEEK1imBRRK7dIjKLEEIo0WYhaeVSG5rZl41fnJajULCIzCKEUK1K9YYag/d6ErAOMDN53gP4CFiv4qVrRpFZhBBCraLVUGa2npmtDzwG7GNmq5rZKsDetKEbHxVTLLOIYBFCqEZZ2iy2NLOHcxNm9gjwvcoVqWUo1MCdu62qtemrTEIIYWlZgsUXks6V1EfSupJ+g99mtSxJe0h6V9JESWcXWN5Z0vBk+SuS+iTzh0h6I/l7U9IB9TmoxlAsswC/X3cIIVSTLMFiKNATuC/565nMKym5Jes1wJ7AxsBQSRvnrXY8MNPM+gFXApcl898GBpvZpsAewP9JynLvjUZTKlhEI3cIodpkuZ/Fl8DPJK1gZnPqse8hwEQzmwwg6Q5gP2B8ap39gPOT53cDV0uSmX2TWqcLzTC8SK9ecOmlcMghtfNyd9OLdosQQrUpm1lI2kbSeJKTvKRBkv6eYd+9gCmp6anJvILrmFkNMAtYJXmdrSSNA94CTkmW55ftJEmjJY2ePn16hiJlJ8HZZ8N6qT5fkVmEEKpVlmqoK4HdSdopzOxN6t7rophCw5rnZwhF1zGzV8xsALAlcI6kLkutaHadmQ02s8E904M4VUhkFiGEapXpCm4zm5I3K0sT71Rg7dR0b2BasXWSNonuQJ2L/cxsAjAX2CRLWSspMosQQrXKEiymSNoGMEmdJJ0JTMiw3Sigv6T1JHUCDgNG5K0zgtq78B0MPGVmlmzTAUDSusAGwAcZXrOiIrMIIVSrLD2MTgH+grcvTAUeB35cbiMzq5F0Gn5RX3vgJjMbJ+lCYLSZjQBuBG6VNBHPKA5LNt8OOFvSImAJcKqZfVG/Q2t8kVmEEKpVyWCRdH/9oZkd0ZCdJxfzPZw377zU8/nAIQW2uxW4tSGvWUmRWYQQqlXJaigzW4x3bw1EZhFCqF5ZqqFekHQ1MBxvaAbAzF4rvknbFJlFCKFaZQkW2ySPF6bmGbBT4xenZYvMIoRQrbJcwb1jUxSkNYjMIoRQrTKNtyTp+8AAfOgNAMzswuJbtE2RWYQQqlWW4T7+ARwK/AS/4voQYN0Kl6tFiswihFCtslyUt42ZHYWPDnsBsDV1r8yuGrnMIoJFCKHaZAkW85LHbyStBSyijd9StZhcZhHVUCGEapOlzeJBST2APwGv4T2hbqhoqVqoyCxCCNUqS2+o3ydP75H0INDFzGZVtlgtU2QWIYRqVTZYSDqqwDzM7JbKFKnliswihFCtslRDbZl63gXYGa+OqtpgEZlFCKHaZKmG+kl6WlJ3WuAgf00hus6GEKpVppsf5fkG6N/YBWkNIrMIIVSrLG0WD1B7O9R2wMbAnZUsVEvVvr3fmzsyixBCtcnSZnF56nkN8KGZTa1QeVo0ybOLyCxCCNUmS5vFM01RkNaiU6fILEII1SdLNdRsaquh6iwCzMxWbPRStWCRWYQQqlGWaqgrgU/xHlACjgC6mdkfK1mwlioyixBCNcrSG2p3M/u7mc02s6/N7FrgoEoXrKWKzCKEUI2yBIvFko6Q1F5SO0lHAIsrXbCWKjKLEEI1yhIsDgd+AHwGfI7fz+LwShaqJYvMIoRQjbL0hvoA2K/yRWkdIrMIIVSjopmFpBMl9U+eS9JNkmZJGitp86YrYssSmUUIoRqVqob6GfBB8nwoMAhYHzgD+Etli9VyRWYRQqhGpYJFjZnlTot7A7eY2QwzGwl0rXzRWqaOHSNYhBCqT6lgsUTSmpJyw5KPTC1brrLFark6dYpqqBBC9SnVwH0eMBpoD4wws3EAkr4HTG6CsrVIHTvC7NnNXYoQQmhaRYOFmT0oaV38au2ZqUWjgUMrXrIWKjKLEEI1Ktl11sxqgJl58+ZWtEQtXLRZhBCqUUNuflTVoutsCKEaRbCop+g6G0KoRllGnUVSL2Dd9Ppm9mylCtWSRWYRQqhGWe5ncRneoD2e2gEEDajKYBGZRQihGmXJLPYHNjCzBZUuTGsQmUUIoRplabOYDHSsdEFai8gsQgjVKEtm8Q3whqQngf9lF2b204qVqgWLzCKEUI2yBIsRyV/AM4uaGjADqe6yadN83pprNk/ZQgihUrLcz+LmpihIa9ExqZBbtMgDR9rQodC1Kzz8cNOXK4QQKqlsm4Wk/pLuljRe0uTcX5adS9pD0ruSJko6u8DyzpKGJ8tfkdQnmb+rpDGS3koed6rvgVVKLkDkt1uYweuvw4cfNn2ZQgih0rI0cP8TuBaoAXYEbgFuLbeRpPbANcCewMbAUEkb5612PDDTzPoBVwKXJfO/APYxs28DR2d5vaaSyyzy2y2mTfMBBj//vOnLFEIIlZYlWCxnZk8CMrMPzex8IMsv/SHARDObbGYLgTtY+vas+wG5aq67gZ0lycxeN7NpyfxxQBdJnTO8ZsUVyyzGj/fHGTNg8WJCCKFNyRIs5ktqB7wn6TRJBwCrZdiuFzAlNT01mVdwnWTQwlnAKnnrHAS8Xug6D0knSRotafT06dMzFGnZFcssJkzwRzMPGCGE0JZkCRanA8sDPwW2AI7Eq4bKUYF5Vp91JA3Aq6ZOLvQCZnadmQ02s8E9e/bMUKRlVyyzyAULiKqoEELbk6U31CgArx2yY+ux76nA2qnp3sC0IutMldQB6A58mbxeb+A+4Cgzm1SP162odG+otPHjoUMH71YbwSKE0NZk6Q21taTxwIRkepCkv2fY9yigv6T1JHUCDmPp6zVGUJulHAw8ZWYmqQfwEHCOmb2Q8ViaRC6zKFQNteWW/ryJasRCCKHJZKmGugrYHZgBYGZvAtuX2yhpgzgNeAwPNHea2ThJF0raN1ntRmAVSROBM4Bc99rTgH7AbyW9kfxlaSepuEKZxRdfeIDYPnlXIrMIIbQ1mYYoN7Mpqnu5cqb+Pmb2MPBw3rzzUs/nA4cU2O4i4KIsr9HUCmUWufaK7baDP/0pgkUIoe3JEiymSNoGsKQ66ackVVLVqFBmkQsWm2wCq6wS1VAhhLYnSzXUKcCP8W6uU4FNk+mqVKjr7IQJsPzysM46sNpqkVmEENqeLL2hvgCOaIKytAqFus6OHw8bbgjt2kWwCCG0TUWDhaS/ltqwmocoh6Uzi1zjds+e8MYbTV+uEEKopFKZxSnA28Cd+PURhS6gqzr5mcXs2TBlCmy0kU9HZhFCaItKBYs18Z5Kh+KDCA4H7jGzmU1RsJYqP7N45x1/3DgZInG11eCrr3x5/hDmIYTQWhVt4DazGWb2DzPbETgG6AGMk/TDpipcS5SfWeR6QuUyi9yoI1980bTlCiGESspyBffm+PhQRwKPAGMqXaiWLD+zmDDBh/no29enV0suHYyqqBBCW1I0WEi6QNIY/MrqZ4DBZna8mY1vstK1QPmZxfjx0L9/bRApFCz+/W849VQfkTaEEFqjUpnFb/GB/QYBlwKvSRqb3L1ubJOUrgVKZxY1NfDqq/Dtb9cuz1VDpS/MGzYMrr0WHn20yYoZQgiNqlQD93pNVopWJJ1ZPPAAfPopHH547fJCmcXbb/vjOefA7rv79RghhNCaFA0WZhZ3ky4gnVlcey307g3f/37t8h49vA0jFyy+/NJvubrlljBqFAwfDkOHNn25QwhhWcRv3HrKBYvx4+GJJ+Ckkzw45EheFZWrhsplFeefDwMHwm9/u/S9MEIIoaWLYFFP7dt7NdKdd3qQOOGEpddJX5iXCxYDB8Ill8CkSXDjjU1X3hBCaAxZus7+LMu8atKpk1dD7b8/rLnm0svzg0X37tCrF+y1lw9jfvHF0TMqhNC6ZMksCt1v+5hGLkerkquK+tGPCi/Pr4baZBOvnpLgwANh6lRvywghhNai1ECCQ4HDgfUkpW+H2o3krnnVqlMn2GAD2HHHwstzmYWZB4tDD61dtv76/jhpkt/7IoQQWoNSXWdfBD4BVgWuSM2fDVTtdRYAZ54Jm27qmUIhq60Gc+Z4QJg50zOLnNyV3pMmwZAhxV/DrPj+QwihqZUaG+pDM/uPmW0NvINnFN2Aqcn9tavW2WfDHnsUX567MO/pp/0xHSxymcXkycW3f+UVWHllePnlZStnCCE0liwN3IcAr+Ij0P4AeEXSwZUuWGuWuzDvqaf8ccCA2mXLL++N4pMmFd/+pZd85Nojj/QMJYQQmluWBu5zgS3N7GgzOwoYgg8FEorIBYunn4Y11oBVV627fP31SweLSZO8EX3yZDjjjMqVM4QQssoSLNqZWXoM1RkZt6tauWqozz6rWwWV07dv6WAxcaKPN/XLX8L118OIEcXXDSGEppDlpP+opMckHSPpGOAh4OHKFqt1y2UWUHeQwZy+feHjj2H+/MLbT5rk61x4IQwa5Bf+zajq/mchhOZWNliY2S+B/wMG4iPQXmdmZ1W6YK3ZCitAly7+vFhmAfD++0svq6nx+f36QefOcPXVfs3Gk09WrrwhhFBO1uqkF4CngSeT56GE3PhQUDhYpK+1yPfRRx4wcgFliy18eJFx4ypT1hBCyCJLb6gf4L2hDiZ6Q2WWq4rK3Zs7LX2tRb7cvH79/HG55Ty4RLAIITSnUhfl5fwG7w31OYCknsBI4O5KFqy1W2MNP8mvsMLSy3r29PmFrrWYONEfc8ECvOvt+Kq+P2EIobllCRbRG6oBLr0UZs8uvEwq3iNq4kRv70gPUDhgADz0kA9emLv5UgghNKWG9oZ6pLLFav2+/W3YZpviy4tda5HrCZW+m96AAd6O8d//Znvt3/629oLAEEJoDA3tDfWrShesrevb13s9LVlSd/7EibVtGjm5K8CztFt8+SVcdBH85S+NU84QQoBs1VCY2b3AvQCS2ks6wsxur2jJ2ri+fWHBAr/lau/ePm/JEs8sdt+97robbJC9R9To0f74wgsxGGEIofEUzSwkrSjpHElXS9pN7jRgMt4rKiyDQt1nP/nEL9RLN26Dt2H07ZstWLz6qj/OmAHvvts4ZS3lnXdgo41KD4wYQmj9SlVD3QqBxBcMAAAZKElEQVRsALwFnAA8jg8muJ+Z7dcEZWvTCnWfzfWEyq+GAq+KyhIsXnkFVlzRn7/QBFfE3H23B4zbbqv8a4UQmk+pYLG+mR1jZv8HDAUGA3ub2RtNU7S2bZ11/H7e6WCRf41F2oABHkwWLCi+TzPPLPbf3wcvfP75xi1zISNH+uM991T+tUIIzadUsFiUe2Jmi4H3zaxIZ9BQXx07wrrr1q2+mTgROnTwQJJvwABYvLh0j6iPPvI79G21FWy7beUzi7lz4cUX/bqRsWNrM6MQQttTKlgMkvR18jcbGJh7LunrpipgW5bffXbiROjTxwNGviw9onLtFUOGeLB47z0f+bZSnnsOFi2CSy7x6cguQmi7St0pr72ZrZj8dTOzDqnnKzZlIduq/AvzJk0qXAUF3iOqffvywaJTJxg4ELbbzue9+GLjlTffyJE+2OERR8DgwREsQmjL4krsZrTppn5dxLXXentDoWsscjp39kBSLlhstpkHjM03920q2W4xcqRnMMstBwcdBKNGeVVYCKHtiWDRjE44AfbeG378Y7/J0ddfF88swAclzAWLOXPgiitg6lSfrqnxayyGDPHpzp39eaXaLT7/HN58E3bZxacPOsgf77uvMq8XQmheFQ0WkvaQ9K6kiZLOLrC8s6ThyfJXJPVJ5q8i6WlJcyRdXckyNqcOHWD4cG+QPvlkn1cqWOR6RL36qlf7nHkmHHKIB4oJE+Cbb2qDBfiv/jFjfH5jyw0nkgsW/fv7ECdRFRVC21SxYCGpPXANsCewMTBUUv6A3ccDM82sH3AlcFkyfz5+n+8zK1W+lmL55eHBB2HDDX26WDUUeLBYsgS23hpmzYKzzoKXX/YG5nTjds5223kgGTWq8cs9ciT06OHVXTkHHeTVXp9+2vivF0JoXpXMLIYAE81sspktBO4A8i/m2w+4OXl+N7CzJJnZXDN7Hg8abd4qq/jJ929/qw0ahQwe7I3cu+7qVUB/+AMceaTffvX66/3knc5Mtt7aHxu73cIMnngCdtrJy5Oz776+LAYxDKHtqWSw6AVMSU1PTeYVXMfMaoBZwCoVLFOL1asXnHZa6bGc+vXzNoqHH669udLVV/u2r7ziWUV6tNqVV/YAc9ttfo1G2s9/Dscf37CyTprkDdm5KqicTTbx60fGjm3YfkMILVclg0Wh0541YJ3iLyCdJGm0pNHTp0+vV+FaqzXWqBsQuneHW2/1IFNoSPRf/cqH47g7dauqV1+Fq66Cm25qWBbwSDJAfX6w6NjRG+EjWITQ9lQyWEwF1k5N9wamFVtHUgegO/Bl1hcws+vMbLCZDe6Zu+l1Fdp+e3jjDW/wznfQQX4C//3vvb3DzANIz56w9tpw9tk+rz7uvdfbT/r3X3rZwIERLEJoiyoZLEYB/SWtJ6kTcBgwIm+dEcDRyfODgafM6nvqCuAn6a5dl57frp3fDGncOD/JP/IIPPMM/O533tYxalT9ejBNnw7PPgsHHli8HB9/7KPetkY1NV51969/NXdJQmhZVMlzs6S9gKuA9sBNZnaxpAuB0WY2QlIXfHTbzfCM4jAzm5xs+wGwItAJ+ArYzcyK3ol68ODBNjp3M4dQx+LFngl06uRZxPz5fk/vdu1g0CC/Xeu4cV6NVM4NN8CJJ3omM2jQ0sufeAJ2282rt3bcsfGPpdLGjPFgscMO8PTTzV2aECpP0hgzG1xuvUw3P2ooM3sYeDhv3nmp5/PxYc8LbdunkmWrJu3bw7nnwg9/6NPDh9cGhksv9V5MN91Ue61HKffc42NaDRxYeHlu/tixrTNY5HqOPfcczJwJK63UvOUJoaWIK7irxGGHebfcrbbyC/ly9t7br8c4/3wfRbaUr76CJ5/0dpBivbZWX917alWq3eL++/3K9Up5/nnPwBYvhkcfrdzrhNDaRLCoEh06+NAfjz9e90QvwWWX+YV0V11Veh8PPuijzBZrr8gZONCvA2lsI0bAwQd7Q/777zf+/s08WBx0kAe8EfktbCFUsQgWVWTllWvvope2zTZ+w6TLLoMvvqidP2+et23k3HOPX9ORvkq8kEGDvA2kpsanzeC660rf72LRIr/A8Mgj4YMPll7+7LPwgx9424sEw4aVLkNDTJ7sQXP77T3jeuQRL1cIIYJFSFxyiVdDXXSRT3/6qVdPDRjgJ8977/VqmQMOqHudRyEDB3ojei44PP20t4fstFPtwIdpZnDqqXD77X49yIYbwjnneFXWs8/6/H32gfXW82qwXXeFf/5z6QsNC6mpgccey9Y9ONdesd12/nqzZnnbRQgBMLM28bfFFltYWDYnnGDWsaPZo4+arbee2fLLm511ltnaa5v56dbsqafK7+f1133d4cN9et99zVZayaxbN7NNNjGbObPu+pdd5uv/5jdmU6aYHXlk7evl/tZZx+yjj3z9O+7weY8/Xr4suX0/8US24+/Rw2zxYrM5c8w6dzY7/fTy24XQmuG9U8ueY5v9JN9YfxEslt3UqWbLLeffip49zV591ecvXGh2yy0eOGpqyu9n/nyzDh3Mfv1rs0mTzCQPBCNHejD63vfM3nzT7IUXzP76V3+9ww7zk3TOG294sBk50uy11/zknTNvngefww4rXY45c/w4wOy008qXe8MNzb7//drpvfYy69vXbMmS8tuG0FpFsAgNcvnlZptvbjZx4rLtZ5NNzPbe2+znP/fA8fHHPv9f/7KlsoZttvEAUB+nnea//L/80qfHjjW7/vq6AeeKK3z//fp5dlTqpP/5577upZfWzrv2Wp83blz9yhZCa5I1WESbRajjF7/wC9NKDZWexcCBfnX4jTd6V9211vL5Q4f6rV7vusvbQF580ds0unSp3/6POw4WLIB//ANOP93vEHjiid5Tyszv4fHHP8LOO3v7x5QpfiFh2uuv+wWJUHv72W23rV2+997+eMcdpcuSu6jxqad83ddfz34cb77pHQb+9CdvI2ku8+aV7oAQQrNnBI31F5lFy/KHP9RmDi+9VJnX2HRT379kdsop/pfLDq66yp8/84zZZ5/5Or/7Xe22Dz9cm9VMm2b2y1+adeq0dIZz4IFedfbii8XLsf/+dTOljh3Nnn8+2zGceKJZu3a+3QormP34x57d/f3vnoUtXFh82wULzH74Q7Obb872WqUcdZS3Uc2atez7Cq0LUQ0VmlPuZDxkSOVe4/77vY1h1CifXrzYbOjQ2hPvDjvUrrvttmabbVa73mabma2+up8g11rL2ya22Wbp1/jyS7P11zfr1curqvJ9+qlZ+/beKP/UU2Zjxpj17+9tJR98ULr8CxZ428sRR3i7zBFHeJVdOvD8/e/Ftz/zzNpgeeutpV+rlLff9n2Adx4I1SWCRWhW06ebrbyyn9Cb0oIFZrvvbkv13Mr1ivrwQ7O77vLnN9/sDel9+vj0r35VeJ+vvebtIzvvvHQD/1/+Yku1a0yYYNa9u9mgQWazZxcv6wMP+LYPPlg7b9Eis6+/9mzo2982Gzy48LaPPebbnnCC2U47ecC6997S700x++9vtuKKZquuanbooQ3bR2i9IliEZtdcvYjmzavNNnLeece/7VddZbbRRv6XO/F/8YX/Sn///eL7vPFG3/6CC+rO33LL2owl7dFHvXpp3309gBVy+OEeUIstz1Wlvflm3fmffeZZ0cYbm33zjQekrbf26q+RI4sfQyGvvOKv8fvfe+Dp1s17s4XqEcEihDwbbOC/+MGzi/o67DDPMCZN8ul33/V9XX554fWvucaX77KLZwtpc+eade1qdtJJxV9v+nQPAOlrPZYs8aq3zp3rBpGZM73r7wYbZOvenLPLLl5l9vXXZg895OV96KHs27cGS5aYvfVWdIEuJoJFCHl+9Sv/xm+6ad0utllNmeJtHAcc4NPnned1/bluwYUMG+ZVRJtv7u0bObkLC59+uvRrHnyw2Sqr1GYfuWzjr39det0777R6tTvkqrKuvNKn58/3zOKEE7JtX8js2d6RYOzYhu+jMc2fb3bMMX6c11/fNK+5cGHp70RLE8EihDxjxtReod5QF11k/7sifP31/Zd5OQ8/7EFm3XXNHnnE5+23nzesl8sCch0F7r7bbPRoL/+++xb+lbx4sVevbbJJ3WA4YoRfM/LVV7Xr/fnPnp2st17dHmCHHmq22mr1y05ypk3zoAh+JfzLL9d/H4UsWeL7KhdY833+udl223l5Vl7Z35tKZxcLF/p3okuX2gy0pYtgEUIB33yzbNvPm+cn2NVW8/+eYcOybffqq15FBN4dt1OnbEOJ1NR4T6wddvCLC3v39jaWYm691V/jvvt8+q67arvmdu1qdvLJZjvu6NP77FM32zEz+/e/fdlzz/mJ9c9/9kzsnXdKl3PsWL/wsWtXs+uu895lK6xQ/xN82uTJfuX/+uvb/3p9Pftstm0nTPDPqUsXP6ZbbvF9PPZYw8tTzpIl/v7muk8ffHC27b76yr8L+W1TTSWCRQgVcu+9/p/TpUv9rkuYP9/s4ot9O/DG5Sx+8xtfv1278ifLRYv8RL355p5BdezoXYKffdbs2GP9tVdYweyGGwr/yp41ywPZiSf6yQ68O2/fvt6Gkvbee2ZXX+1Z0vLLm625pmdvZl4Ns/HG/nr1bQOZMsXbcjp08GPedVcvb9++np3lMqRiXnrJM4nVVqt9j+fPN1tjDbM996xfWerjyiv9/TrrLO8IAdmCWy5b7dLFq8qaum0lgkUIFbJkiVfXNHSQwcmTvVop60lh0iQ/wV98cbb1cz23OnTw7rvpgRu/+qr8yXbPPWuD0x//6Bckdu7sVTrz5/sYYrnrWcC7Hp90Uu1AjznTp3vQ6tDB7LbbaucvWeL7/L//82M64wyz444zO+QQs91289fq2NEvUJw6tXa7l1/29p8jjihe9oce8vHN+vZdesiaCy/08k6YUDvvm288qP7852ZbbGF2/PGle8V98413Tsj34IOe+RxwgFfzzZ3rWeAWW5RuH8sFse9+16uvwHvJFXqNSolgEUIbUp+Tx4IFXnXTr9/S1UxZPPCAZwXpKqThw/1ssfXWXtXUubNnPO+9VzrozZpVW+115ZXeZpNrR8j95S6M3Ggjv4jzpJOKn7BzJ/xbbqltV5k716vbDj64cGeCnM8+83L/6Ef+Hl1xRW3vuM6d/cLNzp09uJ18sne/znUs+Pxzs3PP9baY/v19XznvvuvXqWy2Wd0BL2+7zcpWVd50k/1vBOWaGj8+yQNGsfd15kzPbhurTSSCRQhVbMaMxv91esklfsY44ADPjrKaN8/baXLBYe21vTfXlCn1H0By0SKvVsvtq1u32mq91Vc3++lPl+6mnHbccR6cvvUt32bPPT2A5d6rKVPMTj3VMxvwKrnNNqsdjXnvvf354MH+OrNmeZBbddWlr9hfvNhsq628Ouyee5Y++S9Z4p0RBg6su+z3v/fXuvbauvsaNswzr1zZunb1oLmsIliEEBrdJ580bLuaGg82N95Y/CLErKZP9yzld7/z6qMzzvCr9bP04Bo71qvXNtjAe5oVM22aZytnneXVQ8ceW3uV/ogRnsHsuqtf/d6+ffH7vIwdW9uxYYstPDDlAsOjj9r/RhJIW7zYbI89PFCNHu2jDuSys379fByzxx4z2357n3fssXUzmvqKYBFCCAW8/37pARqzyFUfgfcYK2XRIrN//rN2WJnNNjO7/XYfPmattQoHz+nTvc2jVy+v4irUKWHRIq8ak+reh6W+sgYL+bqt3+DBg2306NHNXYwQQpW4/nof+v6CC/y+8OUsXAi33gqXXw7vvOPzLr0Uzj678PovvQQ77ABbbeX3nF9//cLrjRwJK67oQ903hKQxZja47HoRLEIIoeksWQIPPQRPPOH3vF9xxeLrfvkl9OhR/r73yyJrsOhQuSKEEELI164d7LOP/5Wz8sqVL09Wcae8EEIIZUWwCCGEUFYEixBCCGVFsAghhFBWBIsQQghlRbAIIYRQVgSLEEIIZUWwCCGEUFabuYJb0nTgw3pssirwRYWK05JV43FX4zFDdR53NR4zLNtxr2tmPcut1GaCRX1JGp3lEve2phqPuxqPGarzuKvxmKFpjjuqoUIIIZQVwSKEEEJZ1RwsrmvuAjSTajzuajxmqM7jrsZjhiY47qptswghhJBdNWcWIYQQMopgEUIIoayqDBaS9pD0rqSJkorc1LB1k7S2pKclTZA0TtLPkvkrS3pC0nvJ40rNXdZKkNRe0uuSHkym15P0SnLcwyV1au4yNiZJPSTdLemd5DPfuho+a0k/T77fb0v6t6QubfGzlnSTpM8lvZ2aV/Dzlftrcn4bK2nzxihD1QULSe2Ba4A9gY2BoZI2bt5SVUQN8Asz2wj4DvDj5DjPBp40s/7Ak8l0W/QzYEJq+jLgyuS4ZwLHN0upKucvwKNmtiEwCD/2Nv1ZS+oF/BQYbGabAO2Bw2ibn/UwYI+8ecU+3z2B/snfScC1jVGAqgsWwBBgoplNNrOFwB3Afs1cpkZnZp+Y2WvJ89n4yaMXfqw3J6vdDOzfPCWsHEm9ge8DNyTTAnYC7k5WaVPHLWlFYHvgRgAzW2hmX1EFnzV+a+jlJHUAlgc+oQ1+1mb2LPBl3uxin+9+wC3mXgZ6SFpzWctQjcGiFzAlNT01mddmSeoDbAa8AqxuZp+ABxRgteYrWcVcBfwKWJJMrwJ8ZWY1yXRb+8zXB6YD/0yq3m6Q1JU2/lmb2cfA5cBHeJCYBYyhbX/WacU+34qc46oxWKjAvDbbf1jSCsA9wOlm9nVzl6fSJO0NfG5mY9KzC6zalj7zDsDmwLVmthkwlzZW5VRIUke/H7AesBbQFa+CydeWPussKvJ9r8ZgMRVYOzXdG5jWTGWpKEkd8UBxu5ndm8z+LJeSJo+fN1f5KmRbYF9JH+BVjDvhmUaPpKoC2t5nPhWYamavJNN348GjrX/WuwDvm9l0M1sE3AtsQ9v+rNOKfb4VOcdVY7AYBfRPekx0whvERjRzmRpdUk9/IzDBzP6cWjQCODp5fjTw/5q6bJVkZueYWW8z64N/tk+Z2RHA08DByWpt6rjN7FNgiqQNklk7A+Np4581Xv30HUnLJ9/33HG32c86T7HPdwRwVNIr6jvArFx11bKoyiu4Je2F/9psD9xkZhc3c5EanaTtgOeAt6itu/813m5xJ7AO/s92iJnlN5y1CZJ2AM40s70lrY9nGisDrwNHmtmC5ixfY5K0Kd6g3wmYDByL/xhs05+1pAuAQ/Hef68DJ+D1823qs5b0b2AHfCjyz4DfAfdT4PNNAufVeO+pb4BjzWz0MpehGoNFCCGE+qnGaqgQQgj1FMEihBBCWREsQgghlBXBIoQQQlkRLEIIIZQVwSK0KpJM0hWp6TMlnd9I+x4m6eDyay7z6xySjAz7dN78PulRRZN550s6swGvcYykq5e1rCHkRLAIrc0C4EBJqzZ3QdKS0YyzOh441cx2rFR5QmhsESxCa1OD32/45/kL8jMDSXOSxx0kPSPpTkn/lfQHSUdIelXSW5L6pnazi6TnkvX2TrZvL+lPkkYl9wc4ObXfpyX9C7/4Mb88Q5P9vy3psmTeecB2wD8k/ak+By7pxKQMb0q6R9LyyfxDktd4U9KzqU3WkvRocr+DP6b2s5uklyS9JumuZPwwkvdlfHKMl9enbKHt61B+lRBanGuAsekTYAaDgI3wYZ4nAzeY2RD5TaF+ApyerNcH+B7QF3haUj/gKHzIhC0ldQZekPR4sv4QYBMzez/9YpLWwu+rsAV+T4XHJe1vZhdK2gm/srzQVbV9Jb2Rml4DH1kV4F4zuz7Z/0V4hvI34DxgdzP7WFKP1Lab4qMNLwDelfQ3YB5wLrCLmc2VdBZwRlJldQCwoZlZ3n5CiGARWh8z+1rSLfiNb+Zl3GxUbnwcSZOA3Mn+LSBdHXSnmS0B3pM0GdgQ2A0YmMpauuM3llkIvJofKBJbAv8xs+nJa96O33Pi/jLlnGRmm+Ym8tpjNkmCRA9gBeCxZP4LwDBJd+KD6eU8aWazkv2MB9ZNtt0YD3jgw4O8BHwNzAdukPQQ8GCZcoYqE8EitFZXAa8B/0zNqyGpWk3Gx0nfTjM9NtCS1PQS6v4f5I9/Y/iQzz8xs8fSC5Kxp+YWKV+hYaKX1TBgfzN7U9Ix+FhBmNkpkrbCb/j0RjJOFNQ95sX4cQp4wsyGLlVgaQg+GN9hwGn4iL0hANFmEVqpZEC8O6l7y8wP8Gof8PscdGzArg+R1C5px1gfeBf/Bf8j+ZDvSPqW/OZCpbwCfE/Sqknj91DgmQaUJ60b8ElSjiNyMyX1NbNXzOw84AvqDk+d72Vg26R6jWTE1m8l7RbdzexhvEpu0xL7CFUoMovQml2B/wLOuR74f5Jexe9JXOxXfynv4if11YFTzGy+pBvwtozXkoxlOmVu1Wlmn0g6Bx8uW8DDZrasQ2X/Fg9CH+LVZ92S+X+S1D95nSeBNylysjez6UlW8u+k/QW8DWM2/t51SfazVAeCUN1i1NkQQghlRTVUCCGEsiJYhBBCKCuCRQghhLIiWIQQQigrgkUIIYSyIliEEEIoK4JFCCGEsv4/Dm4ps5NWgigAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(dt.num, dt.rmse, color = \"blue\")\n",
    "plt.title(\"RMSE of Jaccard Similarity and MinHash Similarity\")\n",
    "plt.xlabel(\"Number of Hash Functions\")\n",
    "plt.ylabel(\"Root Mean Squared Error\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What is the effect of the number of hash functions used to compute the Minhash signature on the accuracy of the Minhash estimate of Jaccard similarity. \n",
    "\n",
    "### My thoughts:\n",
    "#### 1. It is known to all that the probability of two min-hash values equalling is the Jaccard similarity for that two documents. However, the realization of the probability could be very different. Therefore, we are estimating Jaccard similarity instead of calculating it. Given the law of large numbers, we expect that with more experiements, which are hash functions in our case, we will have a better estimation which is closer to the expected value, namely a better accuracy of the Minhash estimate of Jaccard similarity.\n",
    "\n",
    "#### 2. After calculating the Root Mean Squared Error (RMSE) of Jaccard Similarity and MinHash Similarity, it is obvious that as the number of hash functions increases, the accurary of Minhash estimate of Jaccard similarity increases. We can see that the most dramatic decrease of RMSE happens around number of hash functions equaling to 5, and after the number of hash functions hits 60, RMSE levels out."
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}