Project 2 - Web Traffic Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project 2: Web Traffic Analysis\n",
    "**This is the second of three mandatory projects to be handed in as part of the assessment for the course 02807 Computational Tools for Data Science at Technical University of Denmark, autumn 2019.**\n",
    "\n",
    "#### Practical info\n",
    "- **The project is to be done in groups of at most 3 students**\n",
    "- **Each group has to hand in _one_ Jupyter notebook (this notebook) with their solution**\n",
    "- **The hand-in of the notebook is due 2019-11-10, 23:59 on DTU Inside**\n",
    "\n",
    "#### Your solution\n",
    "- **Your solution should be in Python**\n",
    "- **For each question you may use as many cells for your solution as you like**\n",
    "- **You should document your solution and explain the choices you've made (for example by using multiple cells and use Markdown to assist the reader of the notebook)**\n",
    "- **You should not remove the problem statements**\n",
    "- **Your notebook should be runnable, i.e., clicking [>>] in Jupyter should generate the result that you want to be assessed**\n",
    "- **You are not expected to use machine learning to solve any of the exercises**\n",
    "- **You will be assessed according to correctness and readability of your code, choice of solution, choice of tools and libraries, and documentation of your solution**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "In this project your task is to analyze a stream of log entries. A log entry consists of an [IP address](https://en.wikipedia.org/wiki/IP_address) and a [domain name](https://en.wikipedia.org/wiki/Domain_name). For example, a log line may look as follows:\n",
    "\n",
    "`192.168.0.1 somedomain.dk`\n",
    "\n",
    "One log line is the result of the event that the domain name was visited by someone having the corresponding IP address. Your task is to analyze the traffic on a number of domains. Counting the number of unique IPs seen on a domain doesn't correspond to the exact number of unique visitors, but it is a good estimate.\n",
    "\n",
    "Specifically, you should answer the following questions from the stream of log entries.\n",
    "\n",
    "- How many unique IPs are there in the stream?\n",
    "- How many unique IPs are there for each domain?\n",
    "- How many times was IP X seen on domain Y? (for some X and Y provided at run time)\n",
    "\n",
    "**The answers to these questions can be approximate!**\n",
    "\n",
    "You should also try to answer one or more of the following, more advanced, questions. The answers to these should also be approximate.\n",
    "\n",
    "- How many unique IPs are there for the domains $d_1, d_2, \\ldots$?\n",
    "- How many times was IP X seen on domains $d_1, d_2, \\ldots$?\n",
    "- What are the X most frequent IPs in the stream?\n",
    "\n",
    "You should use algorithms and data structures that you've learned about in the lectures, and you should provide your own implementations of these.\n",
    "\n",
    "Furthermore, you are expected to:\n",
    "\n",
    "- Document the accuracy of your answers when using algorithms that give approximate answers\n",
    "- Argue why you are using certain parameters for your data structures\n",
    "\n",
    "This notebook is in three parts. In the first part you are given an example of how to read from the stream (which for the purpose of this project is a remote file). In the second part you should implement the algorithms and data structures that you intend to use, and in the last part you should use these for analyzing the stream."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading the stream\n",
    "The following code reads a remote file line by line. It is wrapped in a generator to make it easier to extend. You may modify this if you want to, but your solution should remain parametrized, so that your notebook can be run without having to consume the entire file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from urllib.request import urlopen\n",
    "\n",
    "def stream(n):\n",
    "    i = 0\n",
    "    with urlopen('https://files.dtu.dk/fss/public/link/public/stream/read/traffic_2?linkToken=_DcyO-U3MjjuNzI-&itemName=traffic_2') as f:\n",
    "        for line in f:\n",
    "            element = line.rstrip().decode(\"utf-8\")\n",
    "            yield element\n",
    "            i += 1\n",
    "            if i == n:\n",
    "                break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "STREAM_SIZE = 10\n",
    "web_traffic_stream = stream(STREAM_SIZE)\n",
    "# list(web_traffic_stream)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data structures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our data structure is based on the Flajolet-Martin algorithm, which the Hyperloglog algorithm is also based on.\n",
    "To implement this algorithm, and estimate the number of unique IPs in the stream, we defined some helper functions: make_bucket, hash, and calc_prefix_length.\n",
    "Essentially, we are trying to estimate the number of unique elements in the stream by hashing the IP address into a bit string.\n",
    "We use a 64-bit hash function, which returns a tuple of two 32-bit hash values. The first is used to specify a bucket number, but taking modulus of BUCKET_SIZE. The second value is used to make the count estimate:\n",
    "We know that the random chance of having a 0 in the first position of the bit string is 1/2, thus we can assume that the probability will start with 'k' number of 0's is (1/2)^k. We then use 2^k as the estimate of the number of unique elements seen in the stream.\n",
    "Additionally, to reduce the variance of overestimating, we introduce X (BUCKET_SIZE) amount of different hash functions in the same hash familiy (murmurhash3).\n",
    "Each bucket represents the largest number of zeros in sequence seen from the start of any element in the stream for the corresponding bucket hash value.\n",
    "Finally, we take the harmonic mean of the bucket values (k) and use this to calculate our estimate (2^k). The harmonic mean is less influenced by large outliers, and is thus less prone to overestimation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimating the X most frequent IPs in the stream:\n",
    "The estimate was created by implementing the Count-min sketch. A 'd x w' matrix M was created, where d represents hash functions within the same family and w the number of possible hash positions.\n",
    "Thus implementing stacked bloom filters containing individual counts of hashed 'IP + domain'-strings. We chose to hash the IPs and domains together in order to look at each IP's occurance in the stream. Additionally, the total occurance of an IP in the stream was determined by hashing IP + the string 'total' and incrementing the appropriate position in the M matrix.\n",
    "\n",
    "The 'count_ip' function finds the estimated count of a specific IP in a particular domain or 'total', by looking up in the M matrix and taking the minimal value of the positions specified by the hash functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mmh3\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "# Function for making hashing more readable. Returning the absolute value of a 32-bit or 64-bit hashing.\n",
    "def hash(element, seed, bit_size = 32):\n",
    "    if bit_size == 64:\n",
    "        hash_val = mmh3.hash64(element, seed = seed)\n",
    "        return (abs(hash_val[0]), str(bin(abs(hash_val[1]))))\n",
    "    else:\n",
    "        return abs(mmh3.hash(element, seed = seed))\n",
    "        \n",
    "\n",
    "# Function to calculate the length of a bit string prefix of 0s:\n",
    "def calc_prefix_length(bit_string):\n",
    "    # The initial 3 values in a bit string are not directly related to the actual number.\n",
    "    bit_string = bit_string[3:]\n",
    "    prefix_length = bit_string.find('1')\n",
    "    return prefix_length\n",
    "\n",
    "\n",
    "# Function for creating or editing a bucket of prefix lengths, based on hashing in different seed values (0-BUCKET_SIZE).\n",
    "def make_bucket(element, bucket):\n",
    "    for i in range(BUCKET_SIZE):                      # Number of hash functions\n",
    "        bucket_hash, estimate = hash(element, seed = i, bit_size = 64)\n",
    "        bucket_hash = bucket_hash % BUCKET_SIZE       # Bucket location\n",
    "        prefix_len = calc_prefix_length(estimate)\n",
    "        if prefix_len > bucket[bucket_hash]:\n",
    "            bucket[bucket_hash] = prefix_len\n",
    "    return bucket\n",
    "\n",
    "\n",
    "# Function for estimating the count of a specific IP.\n",
    "def count_ip(IP, domain):\n",
    "    X_min = 10**10  # Some large initial value\n",
    "    for i in range(d):\n",
    "        # In order to distinguish IP and domain combinations we combine these in the hash function.\n",
    "        wi = hash(IP+domain, seed = i) % w            # w defines the size limit of returned hash values.\n",
    "        X_val = M[i, wi]\n",
    "        if X_val < X_min:\n",
    "            X_min = X_val\n",
    "    return X_min\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize constants\n",
    "BUCKET_SIZE = 32        # Number of hash functions and number of buckets the prefix lens are inserted in (HyperLogLog)\n",
    "STREAM_SIZE = 100000    # Number of elements to go through\n",
    "d = 10                  # Number of hash functions in Count-min sketch\n",
    "w = STREAM_SIZE         # Number of possible hash values in Count-min sketch\n",
    "M = np.zeros((d, w))    # Count matrix for Count-min sketch\n",
    "X = 10                  # X IPs with highest occurance (estimate)\n",
    "max_count = [0]*X       # Count of occurances of X IPs used for sorting the max_ip list\n",
    "max_ip = [0]*X          # List of X most frequent IPs\n",
    "\n",
    "\n",
    "domains = dict()        # Dictionary of domains with corresponding buckets\n",
    "domains['total'] = [0]*BUCKET_SIZE\n",
    "\n",
    "\n",
    "for element in stream(STREAM_SIZE):\n",
    "    element = element.split('\\t')\n",
    "    IP = element[0]\n",
    "    domain = element[1]\n",
    "    \n",
    "    # Add to total unique count\n",
    "    domains['total'] = make_bucket(IP, domains['total'])\n",
    "    # Add to domains dict according to domain \n",
    "    if domain in domains:\n",
    "        domains[domain] = make_bucket(IP, domains[domain])\n",
    "    else:\n",
    "        domains[domain] = [0]*BUCKET_SIZE\n",
    "    \n",
    "    \n",
    "    # Count-min sketch\n",
    "    for i in range(d):\n",
    "        wi = hash(IP+domain, i) % w     # Location of hashed IP with domain\n",
    "        wit = hash(IP+'total', i) % w   # Location of hashed IP + 'total'\n",
    "        M[i, wi] += 1\n",
    "        M[i, wit] += 1\n",
    "    \n",
    "    # Estimating the X most frequent IPs\n",
    "    if IP not in max_ip:\n",
    "        num_ip = count_ip(IP, 'total')\n",
    "        max_ip = [ip for count, ip in sorted(zip(max_count, max_ip))]\n",
    "        max_count.sort()\n",
    "        if num_ip > max_count[0]:\n",
    "            max_count[0] = num_ip\n",
    "            max_ip[0] = IP\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Short description of 'Estimating the X most frequent IPs':\n",
    "\n",
    "* We only consider a IP that is not already there.\n",
    "* 1) Estimate the count of a particular IP\n",
    "* 2) Sort max_ip and max_count in order to easily compare with lowest value and insert if higher.\n",
    "* 3) Add if it is higher than the total count of the current IP.\n",
    "\n",
    "The needed data structures (domains, M, and max_ip) are now ready for our analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### How many unique IPs are there in the stream?\n",
    "To answer this question, we define a function that exstracts the bucket of a specific domain from the domains dictionary.\n",
    "The function returns 2^k, where k is the harmonic mean of the bucket."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def harmonic_mean(alist):\n",
    "    n = len(alist)\n",
    "    return n/np.sum([1/val for val in alist if val != 0])\n",
    "\n",
    "def get_count(domain):\n",
    "    bucket = [val for val in domains[domain] if val != 0]\n",
    "    bucket.sort()\n",
    "    return int(round(2**harmonic_mean(bucket)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of unique IPs are: 87413\n"
     ]
    }
   ],
   "source": [
    "print('The number of unique IPs are:', get_count('total'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How many unique IPs are there for each domain?\n",
    "We loop through the domains the domains dictionary and run the get_count function for each domain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total has 87413 unique IPs\n",
      "python.org has 33523 unique IPs\n",
      "wikipedia.org has 59177 unique IPs\n",
      "pandas.pydata.org has 12100 unique IPs\n",
      "dtu.dk has 2645 unique IPs\n",
      "google.com has 2647 unique IPs\n",
      "databricks.com has 1282 unique IPs\n",
      "github.com has 1433 unique IPs\n",
      "spark.apache.org has 737 unique IPs\n",
      "datarobot.com has 203 unique IPs\n",
      "scala-lang.org has 3 unique IPs\n"
     ]
    }
   ],
   "source": [
    "# How many unique IPs are there for each domain?\n",
    "for domain in domains.keys():\n",
    "    print('{} has {} unique IPs'.format(domain, get_count(domain)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How many times was IP X seen on domain Y? (for some X and Y provided at run time)\n",
    "Here, the count_ip functions comes in handy, where we look up the IP + domain combination in the M Count-min matrix using the count_ip function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IP 72.187.84.158 was seen 26.0 times on scala-lang.org\n"
     ]
    }
   ],
   "source": [
    "# How many times was IP X seen on domain Y? (for some X and Y provided at run time)\n",
    "# Define input variables:\n",
    "IP_X = '72.187.84.158'\n",
    "domain_Y = 'scala-lang.org'\n",
    "\n",
    "count = count_ip(IP_X, domain_Y)\n",
    "print('IP {} was seen {} times on {}'.format(IP_X, count, domain_Y))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How many unique IPs are there for the domains d1,d2,…?\n",
    "This is done in the same way as before, however, a list of domain names can be defined if desired."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of unique IPs visiting python.org is: 33523\n",
      "The number of unique IPs visiting wikipedia.org is: 59177\n",
      "The number of unique IPs visiting pandas.pydata.org is: 12100\n",
      "The number of unique IPs visiting github.com is: 1433\n"
     ]
    }
   ],
   "source": [
    "# How many unique IPs are there for the domains d1,d2,…?\n",
    "domain_names = ['python.org', 'wikipedia.org', 'pandas.pydata.org', 'github.com'] # domains.keys()\n",
    "\n",
    "for domain in domain_names:\n",
    "    print('The number of unique IPs visiting {} is: {}'.format(domain, get_count(domain)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How many times was IP X seen on domains d1,d2,…?\n",
    "Again a list of domains of interest are defined and the tables are looped through, counting the number of times a specific IP occurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "python.org 0.0\n",
      "wikipedia.org 0.0\n",
      "pandas.pydata.org 0.0\n",
      "github.com 0.0\n",
      "scala-lang.org 26.0\n"
     ]
    }
   ],
   "source": [
    "# How many times was IP X seen on domains d1,d2,…?\n",
    "IP_X = '72.187.84.158'\n",
    "domains_to_look_at = ['python.org', 'wikipedia.org', 'pandas.pydata.org', 'github.com', 'scala-lang.org'] # domains.keys()\n",
    "for domain in domains_to_look_at:\n",
    "    print(domain, count_ip(IP_X, domain))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What are the X most frequent IPs in the stream?\n",
    "The X most frequent IPs were stored in the max_ip list. As the max_count was not reevaluated constantly, we need to use the count_ip function on the X IPs in the max_ip list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Top 10 IPs and their frequency:\n",
      "56.29.201.127 5.0\n",
      "56.30.200.129 5.0\n",
      "57.30.198.127 5.0\n",
      "204.141.72.187 7.0\n",
      "54.28.199.128 6.0\n",
      "55.29.200.127 6.0\n",
      "55.30.200.128 6.0\n",
      "56.29.201.129 6.0\n",
      "108.41.112.108 27.0\n",
      "72.187.84.158 26.0\n"
     ]
    }
   ],
   "source": [
    "print('Top', X, 'IPs and their frequency:')\n",
    "for ip in max_ip: # Not sorted\n",
    "    print(ip, count_ip(ip,'total'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Document the accuracy of your answers when using algorithms that give approximate answers\n",
    "The Count-min sketch error analysis is derived from this article: http://dimacs.rutgers.edu/~graham/pubs/papers/cmencyc.pdf\n",
    "\n",
    "HyperLogLog accuracy:\n",
    "The accuracy is determined as 1/sqrt(m), where m is the number of buckets.\n",
    "We assume the hash function will uniformly distribute the element hash values to the 32 buckets, and we therefore set the number of buckets, m, to BUCKET_SIZE.\n",
    "Reference: https://storage.googleapis.com/pub-tools-public-publication-data/pdf/40671.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error added with each element added to the M matrix: 2.7182818284590452e-05\n",
      "Probability of allowing a count estimate outside the error above: 4.5399929762484854e-05\n",
      "\n",
      "Error rate of HyperLogLog implementation: ±0.17677669529663687\n",
      "The variance is thus: 0.031249999999999993\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "# Count-min sketch error analysis:\n",
    "print('Error added with each element added to the M matrix:', math.e/w)\n",
    "print('Probability of allowing a count estimate outside the error above:', math.exp(-d))\n",
    "print('')\n",
    "\n",
    "# HyperLogLog accuracy:\n",
    "error = 1/math.sqrt(BUCKET_SIZE)\n",
    "print('Error rate of HyperLogLog implementation: ±{}'.format(error))\n",
    "print('The variance is thus: {}'.format(error**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can thus conclude that with the chosen parameters for d and w, we have a very high chance of hitting the true value, i.e. a high accuracy.\n",
    "\n",
    "Our HyperLogLog implementation has a theoretical variance of ~3%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Argue why you are using certain parameters for your data structures\n",
    "\n",
    "* **BUCKET_SIZE = 32**        -- We use 32 to minimize the risk of having a string hashing to a bit value with many initial 0s. With a higher bucket size, the chance of hitting critically wrong is diminished, as the harmonic mean of the buckets are used.\n",
    "* **STREAM_SIZE = 100000**    -- A stream size of 100000 is used to run through a significant number of elements in the stream. \n",
    "* **d = 10**         -- For the Count-min sketch, the chance of two strings hashing to the same position reduces with a higher d and w. In our case, a d of 10 gave very good estimates when compared to actual counts.\n",
    "* **w = STREAM_SIZE**         -- A w of 100000 means there are enough possible hash values for each single element in the stream. This combined with the d of 10 gives a total of 10 * 100000 positions in the M matrix and thus when using the Count-min sketch, a very good estimate of a count is achieved."
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}