12158054_11913754_12383201_code.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Silvan de Boer (12158054) - Oskar van der Wal (11913754) - Leon Lang (12383201)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import permutations\n",
    "from itertools import product\n",
    "import random\n",
    "from copy import deepcopy\n",
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "import math\n",
    "import statistics\n",
    "from matplotlib import pyplot as plt\n",
    "import pandas\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Offline Evaluation\n",
    "\n",
    "## Step 1: Simulate Ranking Pairs of E and P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_ranking_pairs():\n",
    "    \"\"\"\n",
    "    Outputs a list of tuples with all distinct ranking pairs using a pool of 12 document ids.\n",
    "    The first list in each pair corresponds to system E, the second to system P.\n",
    "    \"\"\"\n",
    "    pairs = []\n",
    "    \n",
    "    # add pairs without overlap\n",
    "    for rel_list in list(product([0, 1], repeat=6)): # all orders of relevance\n",
    "        doc = [-1, 5] # unused doc id with relevance 0 and 1\n",
    "        l_P = []\n",
    "        l_E = []\n",
    "        for r in rel_list[:3]:\n",
    "            doc[r] += 1\n",
    "            l_P.append(doc[r])\n",
    "        for r in rel_list[-3:]:\n",
    "            doc[r] += 1\n",
    "            l_E.append(doc[r])\n",
    "        pairs.append((l_P, l_E))\n",
    "    \n",
    "    # add pairs with one overlap\n",
    "    for overlap in list(product([0, 1, 2], repeat=2)): # all possible overlaps [pos in list P, pos in list E]\n",
    "        for rel_list in list(product([0, 1], repeat=5)): # all orders of relevance (excluding overlapping docs)\n",
    "            doc = [-1, 5] # unused doc id with relevance 0 and 1\n",
    "            l_P = [-1,-1,-1]\n",
    "            l_E = [-1,-1,-1]\n",
    "            i = 0\n",
    "            # construct list P, and overlapping doc of list E\n",
    "            for i in range(3):\n",
    "                r = rel_list[i]\n",
    "                doc[r] += 1\n",
    "                l_P[i] = doc[r]\n",
    "                if i == overlap[0]:\n",
    "                    l_E[overlap[1]] = doc[r]\n",
    "            # construct remainder of list E\n",
    "            j = 0\n",
    "            for i in range(3):\n",
    "                if l_E[i] == -1:\n",
    "                    r = rel_list[3+j]\n",
    "                    doc[r] += 1\n",
    "                    l_E[i] = doc[r]\n",
    "                    j += 1\n",
    "            pairs.append((l_P, l_E))\n",
    "                \n",
    "    # add pairs with two overlaps\n",
    "    overlaps = list(product([0, 1, 2], repeat=4)) # all possible double overlaps as before, concatenated\n",
    "    overlaps = [o for o in overlaps if o[0]!=o[2] and o[1]!=o[3] and o[0]<o[2]] # remove impossible overlaps\n",
    "    for overlap in overlaps:\n",
    "        for rel_list in list(product([0, 1], repeat=4)): # all orders of relevance (excluding overlapping docs)\n",
    "            doc = [-1, 5] # unused doc id with relevance 0 and 1\n",
    "            l_P = [-1,-1,-1]\n",
    "            l_E = [-1,-1,-1]\n",
    "            # construct list P, and overlapping doc of list E\n",
    "            for i in range(3):\n",
    "                r = rel_list[i]\n",
    "                doc[r] += 1\n",
    "                l_P[i] = doc[r]\n",
    "                if i == overlap[0]:\n",
    "                    l_E[overlap[1]] = doc[r]\n",
    "                if i == overlap[2]:\n",
    "                    l_E[overlap[3]] = doc[r]\n",
    "            # construct remainder of list E\n",
    "            j = 0\n",
    "            for i in range(3):\n",
    "                if l_E[i] == -1:\n",
    "                    r = rel_list[3+j]\n",
    "                    doc[r] += 1\n",
    "                    l_E[i] = doc[r]\n",
    "                    j += 1\n",
    "            pairs.append((l_P, l_E))\n",
    "\n",
    "    # add pairs with three overlaps\n",
    "    overlaps = list(product([0, 1, 2], repeat=6)) # all possible double overlaps as before, concatenated\n",
    "    overlaps = [o for o in overlaps if o[0]!=o[2]!=o[4]!=o[0] and o[1]!=o[3]!=o[5]!=o[1] and  o[0]<o[2]<o[4]] # remove impossible overlaps\n",
    "    for overlap in overlaps:\n",
    "        for rel_list in list(product([0, 1], repeat=3)): # all orders of relevance (excluding overlapping docs)\n",
    "            doc = [-1, 5] # unused doc id with relevance 0 and 1\n",
    "            l_P = [-1,-1,-1]\n",
    "            l_E = [-1,-1,-1]\n",
    "            # construct list P, and overlapping doc of list E (so all of them)\n",
    "            for i in range(3):\n",
    "                r = rel_list[i]\n",
    "                doc[r] += 1\n",
    "                l_P[i] = doc[r]\n",
    "                if i == overlap[0]:\n",
    "                    l_E[overlap[1]] = doc[r]\n",
    "                if i == overlap[2]:\n",
    "                    l_E[overlap[3]] = doc[r]\n",
    "                if i == overlap[4]:\n",
    "                    l_E[overlap[5]] = doc[r]\n",
    "            pairs.append((l_P, l_E))\n",
    "        \n",
    "    return pairs\n",
    "\n",
    "def get_relevance(doc):\n",
    "    \"\"\"\n",
    "    Outputs the integer relevance corresponding to a document id\n",
    "    \"\"\"\n",
    "    if (doc<0 or doc>11):\n",
    "        raise Exception('Document id should be from 0 to 12!')\n",
    "    if doc<6:\n",
    "        return 0\n",
    "    else:\n",
    "        return 1\n",
    "    \n",
    "def get_relevance_list(docs):\n",
    "    \"\"\"\n",
    "    Maps a list of document ids to integer relevance\n",
    "    \"\"\"\n",
    "    return [get_relevance(doc) for doc in docs]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that create_ranking_pairs() creates all distinct pairs, that is: all distinct combinations of documents shared between the two lists (including no shared documents) and relevance scores."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Compute the $\\Delta$ERR and form bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_click_probability(rel):\n",
    "    \"\"\"\n",
    "    Outputs the click-probability theta.\n",
    "    \"\"\"\n",
    "    rel_max = 1\n",
    "    return (2**rel - 1)/2**rel_max\n",
    "\n",
    "def compute_ERR(list):\n",
    "    \"\"\"\n",
    "    Computes Expected Reciprocal Rank for a list of form [doc_id1, doc_id2, doc_id3].\n",
    "    \"\"\"\n",
    "    ERR = 0\n",
    "    for r in range(len(list)):\n",
    "        theta_r = compute_click_probability(get_relevance(list[r]))\n",
    "        prod = 1\n",
    "        for i in range(r):\n",
    "            prod *= (1-compute_click_probability(get_relevance(list[i])))\n",
    "        ERR += prod*theta_r*1/(r+1)\n",
    "    return ERR\n",
    "\n",
    "def divide_pairs_over_bins(list_pairs):\n",
    "    \"\"\"\n",
    "    Divides a list of pairs over 10 bins according to the delta ERR.\n",
    "    The i'th element of dERRs contains all pairs with delta EER between i/10 and (i+1)/10\n",
    "    Note that the first bin only starts at 0.05 to make sure it does not happen too often that P outperforms E\n",
    "    \"\"\"\n",
    "    dERRs = [[] for i in range(10)]\n",
    "    for i in list_pairs:\n",
    "        deltaERR = compute_ERR(i[0]) - compute_ERR(i[1])\n",
    "        #Keep all the positive delta ERRs and put them in the respective bin\n",
    "        if deltaERR > 0.05:\n",
    "            if deltaERR < 0.1:\n",
    "                dERRs[0].append(i)\n",
    "            else:\n",
    "                bin = str(deltaERR*10)\n",
    "                dERRs[int(bin[0])].append(i)\n",
    "    return dERRs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Online Evaluation\n",
    "## Step 3: Implement Team-Draft and Probabilistic Interleaving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def coin_to_ranker(cointoss):\n",
    "    \"\"\"\n",
    "    Maps zero to E (experimental system) and one to P (production system)\n",
    "    \"\"\"\n",
    "    if (cointoss == 0):\n",
    "        return \"E\"\n",
    "    if (cointoss == 1):\n",
    "        return \"P\"\n",
    "    else:\n",
    "        raise Exception('This number should be either 0 or 1!')\n",
    "        \n",
    "def softmax(r, max_r):\n",
    "    \"\"\"\n",
    "    The softmax takes the rank of a document as well as the maximal rank and returns\n",
    "    the probability of choosing that document\n",
    "    \"\"\"\n",
    "    tau = 3\n",
    "    ranks = range(1,max_r+1)\n",
    "    normalizer = sum([rank**(-tau) for rank in ranks])\n",
    "    return r**(-tau)/normalizer\n",
    "\n",
    "def td_interleave(pair):\n",
    "    \"\"\"\n",
    "    Performs team draft interleaving to a ranking pair\n",
    "    Returns interleaved list as [(doc_id1, S1), (doc_id2, S2), (doc_id3, S3)] where\n",
    "    Si describes which list contributed the document: E or P\n",
    "    \"\"\"\n",
    "    new_pair = deepcopy(list([pair[0], pair[1]]))\n",
    "    result = []\n",
    "    for i in range(3):\n",
    "        cointoss = random.randint(0, 1)\n",
    "        result += [(new_pair[cointoss][0], coin_to_ranker(cointoss))]\n",
    "        if (new_pair[cointoss][0] in new_pair[1 - cointoss]):\n",
    "            new_pair[1 - cointoss].remove(new_pair[cointoss][0])\n",
    "        new_pair[cointoss] = new_pair[cointoss][1:]\n",
    "    return result\n",
    "\n",
    "def pr_interleave(pair):\n",
    "    \"\"\"\n",
    "    Performs probabilistic interleaving to a ranking pair\n",
    "    Returns interleaved list as [(doc_id1, S1), (doc_id2, S2), (doc_id3, S3)] where\n",
    "    Si describes which list contributed the document: E or P\n",
    "    \"\"\"\n",
    "    new_pair = deepcopy(list([pair[0], pair[1]]))\n",
    "    result = []\n",
    "    sm = [[],[]]\n",
    "    sm[0] = [softmax(r, 3) for r in range(1, 4)]\n",
    "    sm[1] = [softmax(r, 3) for r in range(1, 4)]\n",
    "    population = [[0, 1, 2], [0, 1, 2]]\n",
    "    for i in range(3):\n",
    "        cointoss = random.randint(0, 1)\n",
    "        index = random.choices(population[cointoss], sm[cointoss])[0]\n",
    "        sm[cointoss][index] = 0\n",
    "        result += [(new_pair[cointoss][index], coin_to_ranker(cointoss))]\n",
    "        for j in range(3):\n",
    "            if (new_pair[cointoss][index] == new_pair[1 - cointoss][j]):\n",
    "                sm[1 - cointoss][j] = 0       \n",
    "    return result\n",
    "\n",
    "def decide_winner(interl, clicks):\n",
    "    \"\"\"\n",
    "    Takes an interleaved list and a click pattern and returns the winner: E, P or NW (no winner)\n",
    "    \"\"\"\n",
    "    clicks_E = 0\n",
    "    clicks_P = 0\n",
    "    for i in range(3):\n",
    "        if (clicks[i] == 1):\n",
    "            if (interl[i][1] == \"E\"):\n",
    "                clicks_E += 1\n",
    "            elif (interl[i][1] == \"P\"):\n",
    "                clicks_P += 1\n",
    "    if (clicks_E > clicks_P):\n",
    "        return \"E\"\n",
    "    elif (clicks_E < clicks_P):\n",
    "        return \"P\"\n",
    "    else:\n",
    "        return \"NW\"\n",
    "    \n",
    "def decide_winner_marg(interl, ranking_pair, clicks):\n",
    "    \"\"\"\n",
    "    Returns winner of probabilistic interleaving according to the marginalised estimate\n",
    "    Note that the probability of an assignment given an interleaved list is proportional\n",
    "    to the probability of a list given an assignment, and thus only the latter is computed\n",
    "    \"\"\"\n",
    "    prob_win = [0,0] # unnormalised 'probability' that system E or P wins\n",
    "    \n",
    "    # loop over assignments\n",
    "    for assignment in [[a,b,c] \n",
    "                       for a in [0,1] \n",
    "                       for b in [0,1] \n",
    "                       for c in [0,1]]:\n",
    "        prob_list = 1 # probability of interleaved list given assignment\n",
    "        click_count = [0,0] # count how often each system's docs are clicked\n",
    "        sm = [[softmax(r, 3) for r in range(1, 4)], \n",
    "              [softmax(r, 3) for r in range(1, 4)]]\n",
    "        # loop over list positions\n",
    "        for i in range(3):\n",
    "            doc = interl[i][0]\n",
    "            ass = assignment[i]\n",
    "            cli = clicks[i]\n",
    "            # compute probability of list given assignment\n",
    "            if doc not in ranking_pair[ass]:\n",
    "                # assignment is not possible\n",
    "                prob_list = 0\n",
    "                break\n",
    "            else:\n",
    "                # count clicks\n",
    "                click_count[ass] += cli\n",
    "                # update probability\n",
    "                index = ranking_pair[ass].index(doc)\n",
    "                prob_list *= sm[ass][index]\n",
    "                # normalise softmax\n",
    "                sm[ass][index] = 0\n",
    "                sum_sm = sum(sm[ass])\n",
    "                if sum_sm != 0:\n",
    "                    sm[ass] = [s/sum_sm for s in sm[ass]]\n",
    "                # update softmax of other list if necessary\n",
    "                if doc in ranking_pair[1-ass]:\n",
    "                    index = ranking_pair[1-ass].index(doc)\n",
    "                    sm[1-ass][index] = 0\n",
    "                    sum_sm = sum(sm[1-ass])\n",
    "                    if sum_sm != 0:\n",
    "                        sm[1-ass] = [s/sum_sm for s in sm[1-ass]]\n",
    "        # assign win probability to appropriate system\n",
    "        if click_count[0] > click_count[1]:\n",
    "            prob_win[0] += prob_list\n",
    "        elif click_count[1] > click_count[0]:\n",
    "            prob_win[1] += prob_list\n",
    "            \n",
    "    # return winner\n",
    "    if (prob_win[1] > prob_win[0]):\n",
    "        return \"P\"\n",
    "    elif (prob_win[1] < prob_win[0]):\n",
    "        return \"E\"\n",
    "    else:\n",
    "        return \"NW\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Simulate User Clicks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def yandex_log_parser():\n",
    "    \"\"\"\n",
    "    Parses the Yandex Click Log File and returns a list with for each session a list with a 1\n",
    "    a list of the form [rank1_clicked, rank2_clicked, ..., rank10_clicked] where ranki_clicked\n",
    "    is either 1 if clicked or 0 if not clicked.\n",
    "    \"\"\"\n",
    "    sessionID = 0\n",
    "    links = [-1 for i in range(10)] #dummy list\n",
    "    clicks = [0 for i in range(10)]\n",
    "    list_clicks = []\n",
    "    with open(\"YandexRelPredChallenge.txt\") as f:\n",
    "        for line in f:\n",
    "            words = line.split()\n",
    "            sessionID_old = sessionID\n",
    "            sessionID = int(words[0])\n",
    "            if sessionID>sessionID_old:\n",
    "                list_clicks.append(clicks)\n",
    "                clicks = [0 for i in range(10)]\n",
    "            recordType = words[2]\n",
    "            if recordType==\"Q\":\n",
    "                links = [int(l) for l in words[5:]]\n",
    "            elif recordType==\"C\":\n",
    "                link_clicked = int(words[3])\n",
    "                if link_clicked in links:\n",
    "                    rank_clicked = links.index(link_clicked)\n",
    "                    clicks[rank_clicked] =1\n",
    "    return(list_clicks)\n",
    "\n",
    "def em():\n",
    "    \"\"\"\n",
    "    Expectation-maximization method for determining the parameters alpha and gamma, \n",
    "    using training data from the yandex_log_parser().\n",
    "    Using the update rules from [1].\n",
    "    \"\"\"\n",
    "    diffs = [1]*11 # differences in parameters (first alpha, then ten gammas)\n",
    "    diff_term = 1e-10 # minimal difference for termination\n",
    "    max_iter = 100\n",
    "    \n",
    "    gamma = [0.5 for x in range(10)] # initial probability for each rank whether observed\n",
    "    alpha = 0.2 # initial probability clicked if not relevant\n",
    "    epsilon = 0.05\n",
    "    tolerance = 0.01\n",
    "    click_log = yandex_log_parser()\n",
    "    i = 0\n",
    "    while i < max_iter and max(diffs) > diff_term:\n",
    "        i+=1\n",
    "        alpha_old = alpha\n",
    "        gamma_old = gamma.copy()\n",
    "        total = [2 for x in range(10)]\n",
    "        for j in click_log:\n",
    "            for rank in range(len(j)):\n",
    "                gamma_value = gamma[rank]/total[rank]\n",
    "                alpha_value = alpha/sum(total)\n",
    "                if j[rank]==1:\n",
    "                    gamma[rank] += 1\n",
    "                    alpha += 1\n",
    "                else:\n",
    "                    alpha += (1-gamma_value)*alpha_value/(1-gamma_value*alpha_value)\n",
    "                    gamma[rank] += (1-alpha_value)*gamma_value/(1-gamma_value*alpha_value)\n",
    "                total[rank] += 1\n",
    "        # New alpha and gamma\n",
    "        alpha = alpha/sum(total)\n",
    "        diffs[0] = abs(alpha-alpha_old)\n",
    "        for x in range(len(gamma)):\n",
    "            gamma[x] = gamma[x]/total[x]\n",
    "            diffs[x+1] = abs(gamma[x]-gamma_old[x])\n",
    "    return (alpha,gamma)\n",
    "\n",
    "def click_probabilities(l, alpha, gamma):\n",
    "    \"\"\"\n",
    "    Returns the probability of each position being clicked based on ranked relevance list l,\n",
    "    and parameters alpha and gamma\n",
    "    \"\"\"\n",
    "    return [abs(alpha - (1 - l[i]))*gamma[i] for i in range(3)]\n",
    "\n",
    "def produce_clicks(list, alpha, gamma):\n",
    "    \"\"\"\n",
    "    Performs the Position-Based Click Model (PBM)\n",
    "    Returns which positions are clicked based on the ranked relevance list and parameters alpha and gamma\n",
    "    \"\"\"\n",
    "    probabilities = click_probabilities(list, alpha, gamma)    \n",
    "    return [np.random.binomial(1, probabilities[i]) for i in range(3)]\n",
    "\n",
    "def produce_clicks_random(list, theta):\n",
    "    \"\"\"\n",
    "    Performs the Random Click Model\n",
    "    Returns which positions are clicked based on parameter theta (probability to click)\n",
    "    \"\"\"\n",
    "    return [np.random.binomial(1, theta) for i in range(3)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Simulate Interleaving Experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Simulation of Interleaving Experiment\n",
    "\n",
    "\n",
    "def estimate_win_proportion(ranking_pair, interleave_function, click_function, alpha, gamma, theta, k):\n",
    "    \"\"\"\n",
    "    Runs experiment to estimate the win proportion of E on a single pair\n",
    "    Interleaving and simulating clicks is done k times\n",
    "    The ranking pair, functions to interleave and click and their parameters are arguments to this function\n",
    "    Note that this function does not determine the winner by the marginalised estimation of probabilistic interleaving\n",
    "    \"\"\"\n",
    "    wins_E = 0\n",
    "    wins_P = 0\n",
    "    for i in range(k):\n",
    "        interl = interleave_function(ranking_pair)\n",
    "        relevance_list = get_relevance_list([elem[0] for elem in interl])\n",
    "        if (click_function == produce_clicks):\n",
    "            clicks = click_function(relevance_list, alpha, gamma)\n",
    "        elif (click_function == produce_clicks_random):\n",
    "            clicks = click_function(relevance_list, theta)\n",
    "        winner = decide_winner(interl, clicks)\n",
    "        if (winner == \"E\"):\n",
    "            wins_E += 1\n",
    "        elif (winner == \"P\"):\n",
    "            wins_P += 1\n",
    "            \n",
    "    if (wins_E + wins_P > 0):\n",
    "        return wins_E / (wins_E + wins_P)\n",
    "    else:\n",
    "        return 0.5 # if all simulations resulted in a tie, return 0.5\n",
    "    \n",
    "    \n",
    "# Determines win proportion of experimental system based on the marginalised estimator of probabilistic interleaving\n",
    "def estimate_win_proportion_marginalised(ranking_pair, interleave_function, click_function, alpha, gamma, theta, k):\n",
    "    \"\"\"\n",
    "    Runs experiment to estimate the win proportion of E on a single pair based on the marginalised estimation\n",
    "    of probabilistic interleaving\n",
    "    Interleaving and simulating clicks is done k times\n",
    "    The ranking pair, functions to interleave and click and their parameters are arguments to this function\n",
    "    \"\"\"\n",
    "    wins_E = 0\n",
    "    wins_P = 0\n",
    "    for i in range(k):\n",
    "        # interleave \n",
    "        interl = interleave_function(ranking_pair)\n",
    "        # simulate clicks\n",
    "        relevance_list = get_relevance_list([elem[0] for elem in interl])\n",
    "        if (click_function == produce_clicks):\n",
    "            clicks = produce_clicks(relevance_list, alpha, gamma)\n",
    "        elif (click_function == produce_clicks_random):\n",
    "            clicks = produce_clicks_random(relevance_list, theta)\n",
    "        # decide winning system\n",
    "        winner = decide_winner_marg(interl, ranking_pair, clicks)\n",
    "        if (winner == \"E\"):\n",
    "            wins_E += 1\n",
    "        elif (winner == \"P\"):\n",
    "            wins_P += 1\n",
    "    \n",
    "    if (wins_E + wins_P > 0):\n",
    "        return wins_E / (wins_E + wins_P)\n",
    "    else:\n",
    "        return 0.5 # if all simulations resulted in a tie, return 0.5\n",
    "    \n",
    "def run_interleaving_experiment(deltaERRs, estimator, interleaver, click_function, alpha, gamma, theta, k):\n",
    "    \"\"\"\n",
    "    Runs a full interleaving experiment on ranking pairs orginised in bins (deltaERRs), using an estimator\n",
    "    that runs the experiment for each pair based on an interleaver and click function\n",
    "    Returns for each bin a tuple containing the minimum, maximum and median number of the estimated minimum \n",
    "    number of impressions N, and the percentage of pairs that resulted in a win probability of E lower than 0.5\n",
    "    Note that in this case, N would be estimated as infinity\n",
    "    \"\"\"\n",
    "    result = []\n",
    "    counter = 0\n",
    "    total_pairs = sum([len(d) for d in deltaERRs])\n",
    "    print_at = np.ceil(np.linspace(total_pairs/5, total_pairs, 5))\n",
    "    start_time = time.time()\n",
    "        \n",
    "    for i in range(len(deltaERRs)): # loop over bins\n",
    "        if (len(deltaERRs[i])> 0):\n",
    "            Ns = []\n",
    "            count_Ninf = 0 # counts how often N is infinite (p1<=0.5)\n",
    "            for pair in deltaERRs[i]: # loop over pairs in bin\n",
    "                p1 = estimator(pair, interleaver, click_function, alpha, gamma, theta, k)\n",
    "                counter += 1\n",
    "                if counter in print_at:\n",
    "                    print(\"\\t{:.0f}%\\t({:.0f} secs to go)\".format(\n",
    "                        counter/total_pairs*100,\n",
    "                        math.ceil((time.time()-start_time)/counter*total_pairs-time.time()+start_time)))\n",
    "                N = compute_sample_size(p1)\n",
    "                Ns += [N]\n",
    "                if (N == np.inf):\n",
    "                    count_Ninf += 1\n",
    "                    \n",
    "            minimum = min(Ns)\n",
    "            maximum = max(Ns)\n",
    "            median = statistics.median(Ns)\n",
    "            result += [(minimum, median, maximum, count_Ninf/len(deltaERRs[i]))]\n",
    "        else:\n",
    "            result += [(None,None,None, None)]\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6: Compute Sample Size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_sample_size(p1):\n",
    "    \"\"\"\n",
    "    Returns the required minimal sample size according to the power analysis of the \n",
    "    one-sided one-sample proportion test with alpha 0.05, beta 0.1 and p0 0.5\n",
    "    \"\"\"\n",
    "    if p1<=0.5:\n",
    "        return np.inf\n",
    "    a = 0.05\n",
    "    b = 0.1\n",
    "    p0 = 0.5\n",
    "    delta = abs(p1-p0)\n",
    "    N_intermediate = ((norm.ppf(1 - a)*math.sqrt(p0*(1-p0)) + \n",
    "                      norm.ppf(1 - b)*math.sqrt(p1*(1 - p1)))\n",
    "                      /delta)**2\n",
    "    N = N_intermediate + 1/delta\n",
    "    return math.ceil(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tables and Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def results_table(results):\n",
    "    \"\"\"\n",
    "    Shows a table of min, median and max number of impressions in each bin, and the\n",
    "    percentage of experiments per bin that led to infinite N\n",
    "    \"\"\"\n",
    "    xlabels = [\"min\", \"median\", \"max\", \"infinite\"]\n",
    "    ylabels = []\n",
    "    table = []\n",
    "    i = 0\n",
    "    for row in results:\n",
    "        ylabels.append(\"{:.2f} - {:.2f}\".format(.05 if i==0 else i/10, (i+1)/10))\n",
    "        i += 1\n",
    "        if row[0] is None:\n",
    "            table.append([\"-\"]*4)\n",
    "        else:\n",
    "            table.append([\"{:.0f}\".format(row[0]),\n",
    "                          \"{:.1f}\".format(row[1]),\n",
    "                          \"{:.0f}\".format(row[2]),\n",
    "                          \"{:.1f}%\".format(row[3]*100)])\n",
    "            \n",
    "    return pandas.DataFrame(table, ylabels, xlabels)\n",
    "\n",
    "\n",
    "def plot_medians(td, pr, pr_mar):\n",
    "    \"\"\"\n",
    "    Displays plot of median number of impressions per interleaving method\n",
    "    \"\"\"\n",
    "    data_dict = {\n",
    "        \"team draft\": td,\n",
    "        \"probab\": pr,\n",
    "        \"probab, marg\": pr_mar}\n",
    "    \n",
    "    for name in data_dict.keys():\n",
    "        medians = [row[1] for row in data_dict[name]]\n",
    "        plt.plot([x+1 for x in range(10)],medians, \"o-\", label=name)\n",
    "    \n",
    "    plt.xlim((1,10))\n",
    "    plt.legend()\n",
    "    plt.ylabel(\"median\")\n",
    "    plt.xlabel(\"bin\")\n",
    "    plt.title(\"Median impressions per bin\")\n",
    "    plt.show()\n",
    "\n",
    "def plot_pairs_histogram(pairs_in_bins):\n",
    "    \"\"\"\n",
    "    Displays a histogram of the number of pairs in each bin\n",
    "    \"\"\"\n",
    "    height = [len(pair) for pair in pairs_in_bins]\n",
    "    bins = [x+1 for x in range(10)]\n",
    "    print(len(height), len(bins))\n",
    "    plt.bar(bins, height)\n",
    "    plt.xlabel(\"bin\")\n",
    "    plt.ylabel(\"number of pairs\")\n",
    "    plt.title(\"Number of ranking pairs in each bin\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Experiment execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Creating all ranking pairs\n",
      "\t688 pairs created\n",
      "\n",
      "Computing deltaERR for each pair and assigning pairs to bins\n",
      "\t7 non-empty bins created with in total 241 pairs\n",
      "\n",
      "Performing EM to estimate position based click model parameters\n",
      "\talpha: 0.950 (manually set after EM, result from EM: 0.597)\n",
      "\tgammas: 0.998, 0.781, 0.619\n",
      "\n",
      "Setting random click model parameter\n",
      "\ttheta: 0.500\n",
      "\n",
      "Starting online experiments with k=10000\n",
      "\n",
      "Running experiment 1: team draft interleaving & position-based click model\n",
      "\t20%\t(51 secs to go)\n",
      "\t40%\t(38 secs to go)\n",
      "\t60%\t(26 secs to go)\n",
      "\t80%\t(14 secs to go)\n",
      "\t100%\t(0 secs to go)\n",
      "\n",
      "Running experiment 2: team draft interleaving & random click model\n",
      "\t20%\t(49 secs to go)\n",
      "\t40%\t(36 secs to go)\n",
      "\t60%\t(25 secs to go)\n",
      "\t80%\t(13 secs to go)\n",
      "\t100%\t(0 secs to go)\n",
      "\n",
      "Running experiment 3: probabilistic interleaving & position-based click model\n",
      "\t20%\t(94 secs to go)\n",
      "\t40%\t(71 secs to go)\n",
      "\t60%\t(47 secs to go)\n",
      "\t80%\t(24 secs to go)\n",
      "\t100%\t(0 secs to go)\n",
      "\n",
      "Running experiment 4: probabilistic interleaving & random click model\n",
      "\t20%\t(92 secs to go)\n",
      "\t40%\t(69 secs to go)\n",
      "\t60%\t(46 secs to go)\n",
      "\t80%\t(23 secs to go)\n",
      "\t100%\t(0 secs to go)\n",
      "\n",
      "Running experiment 5: marginalised probabilistic interleaving & position-based click model\n",
      "\t20%\t(308 secs to go)\n",
      "\t40%\t(232 secs to go)\n",
      "\t60%\t(157 secs to go)\n",
      "\t80%\t(79 secs to go)\n",
      "\t100%\t(0 secs to go)\n",
      "\n",
      "Running experiment 6: marginalised probabilistic interleaving & random click model\n",
      "\t20%\t(305 secs to go)\n",
      "\t40%\t(228 secs to go)\n",
      "\t60%\t(155 secs to go)\n",
      "\t80%\t(78 secs to go)\n",
      "\t100%\t(0 secs to go)\n",
      "Showing results\n",
      "\n",
      "\n",
      "Table 1: team draft interleaving & position-based click model\n",
      "            min median  max infinite\n",
      "0.05 - 0.10  31  124.0  inf     4.4%\n",
      "0.10 - 0.20  14   65.0  457     0.0%\n",
      "0.20 - 0.30  10   32.0   66     0.0%\n",
      "0.30 - 0.40   9   28.0   37     0.0%\n",
      "0.40 - 0.50  12   13.0   32     0.0%\n",
      "0.50 - 0.60   8    8.0   13     0.0%\n",
      "0.60 - 0.70   7    8.0    8     0.0%\n",
      "0.70 - 0.80   -      -    -        -\n",
      "0.80 - 0.90   -      -    -        -\n",
      "0.90 - 1.00   -      -    -        -\n",
      "\n",
      "Table 2: team draft interleaving & random click model\n",
      "               min      median  max infinite\n",
      "0.05 - 0.10  17198         inf  inf    66.7%\n",
      "0.10 - 0.20  10564         inf  inf    54.5%\n",
      "0.20 - 0.30  19187  11327985.0  inf    45.5%\n",
      "0.30 - 0.40  10203         inf  inf    55.1%\n",
      "0.40 - 0.50   6561         inf  inf    54.5%\n",
      "0.50 - 0.60  34987         inf  inf    66.7%\n",
      "0.60 - 0.70  32973         inf  inf    60.0%\n",
      "0.70 - 0.80      -           -    -        -\n",
      "0.80 - 0.90      -           -    -        -\n",
      "0.90 - 1.00      -           -    -        -\n",
      "\n",
      "Table 3: probabilistic interleaving & position-based click model\n",
      "            min median   max infinite\n",
      "0.05 - 0.10  76  176.0  1710     0.0%\n",
      "0.10 - 0.20  13   69.5   253     0.0%\n",
      "0.20 - 0.30  10   42.0    81     0.0%\n",
      "0.30 - 0.40   9   26.0    40     0.0%\n",
      "0.40 - 0.50  14   17.0    25     0.0%\n",
      "0.50 - 0.60   8    8.0    15     0.0%\n",
      "0.60 - 0.70   7    8.0     8     0.0%\n",
      "0.70 - 0.80   -      -     -        -\n",
      "0.80 - 0.90   -      -     -        -\n",
      "0.90 - 1.00   -      -     -        -\n",
      "\n",
      "Table 4: probabilistic interleaving & random click model\n",
      "               min      median  max infinite\n",
      "0.05 - 0.10  13362         inf  inf    66.7%\n",
      "0.10 - 0.20  13215  56092747.0  inf    45.5%\n",
      "0.20 - 0.30  13070         inf  inf    58.2%\n",
      "0.30 - 0.40   7980         inf  inf    53.1%\n",
      "0.40 - 0.50  16094         inf  inf    63.6%\n",
      "0.50 - 0.60  31099         inf  inf    57.1%\n",
      "0.60 - 0.70  25084         inf  inf    60.0%\n",
      "0.70 - 0.80      -           -    -        -\n",
      "0.80 - 0.90      -           -    -        -\n",
      "0.90 - 1.00      -           -    -        -\n",
      "\n",
      "Table 5: marginalised probabilistic interleaving & position-based click model\n",
      "            min median   max infinite\n",
      "0.05 - 0.10   9   88.0   inf     6.7%\n",
      "0.10 - 0.20   9   41.0   inf     4.5%\n",
      "0.20 - 0.30   8   32.0  2362     0.0%\n",
      "0.30 - 0.40   8   21.0   128     0.0%\n",
      "0.40 - 0.50  10   17.0    26     0.0%\n",
      "0.50 - 0.60   8    8.0    17     0.0%\n",
      "0.60 - 0.70   7    8.0     8     0.0%\n",
      "0.70 - 0.80   -      -     -        -\n",
      "0.80 - 0.90   -      -     -        -\n",
      "0.90 - 1.00   -      -     -        -\n",
      "\n",
      "Table 6: marginalised probabilistic interleaving & random click model\n",
      "             min      median  max infinite\n",
      "0.05 - 0.10  216         inf  inf    51.1%\n",
      "0.10 - 0.20  181         inf  inf    61.4%\n",
      "0.20 - 0.30  204  11170942.0  inf    49.1%\n",
      "0.30 - 0.40  197  99997362.0  inf    49.0%\n",
      "0.40 - 0.50  364         inf  inf    59.1%\n",
      "0.50 - 0.60  182  29647135.0  inf    47.6%\n",
      "0.60 - 0.70  563         inf  inf    60.0%\n",
      "0.70 - 0.80    -           -    -        -\n",
      "0.80 - 0.90    -           -    -        -\n",
      "0.90 - 1.00    -           -    -        -\n",
      "\n",
      "Figure 1: Median impressions per bin, position-based click model\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd81dX5wPHPk0EmJIQdwggjzIRh2BERgQQHrmqdxeKs9VdrLa460NZW66pWW6uCUrVaRUEtIqCCCCgIgoAgeyUECCMBMiDj+f3xvQkh8ya5yc143q/Xfd17v+Oc5170Pvmec77niKpijDHGVMTH2wEYY4yp/yxZGGOMqZQlC2OMMZWyZGGMMaZSliyMMcZUypKFMcaYSlmyMHVCRFREerhevywiD9VCHWeLyGZPl1tXRORaEVng7Tg8qfi/exn7Gt3nbczE7rMwxYnILiASiFTVQ8W2rwUGANGquqsa5SrQU1W3eShU0wDYv3vjYVcWpiw7gasL34hILBDkvXDqHxHx83YM9Yl9H42fJQtTljeBXxR7Pxn4d/EDRCRARJ4WkT0icsDVtBRUbP9UEUkVkX0iMqXEuW+IyJ9cr1uKyP9EJE1EjrpeRxU7drGI/FFElonIcRFZICKtywpaRMaISHKx97tccawTkUwRmS4i7URknqusz0WkpevYrq4mk1tcMaeKyN3FypomIrNE5C0ROQbcICI+InKfiGwXkcMi8p6IRLiOD3Qde1hE0kXkOxFp59p3g4jscMWwU0SuLbZ9abE6R7rOy3A9j3Tne6mo7jK+s10icr+IbHR9/6+LSGCx/ReKyFpXOctFJK7EufeKyDogs4KEcb7r8x4SkadExKecz6sicpuIbHXF8pKISDllmjpmycKU5VughYj0ERFf4OfAWyWOeRKIAQYCPYCOwMMAIpIE/B4YD/QExlVQlw/wOtAF6AxkAy+WOOYa4JdAW6CZq2x3Xe6KIwa4CJgHPAC0dtX9mxLHn+uKeQJwn4gUj/1iYBYQDrztOvcS4BycprujwEuuYycDYUAnoBVwG5AtIiHAC8BEVW0OjATWlgzalXTmuo5tBTwLzBWRVsUOK+97KbPuCr6ja4FEoLvre3rQFcNgYAZwq6ucfwEfi0hAsXOvBi4AwlU1r5zyLwXigcE43+GUco4DuBAYgtPkeaUrLlMPWLIw5Sm8uhgP/ASkFO5w/bV3M3CXqh5R1ePAn4GrXIdcCbyuqhtUNROYVl4lqnpYVT9Q1SxXOY/j/PgW97qqblHVbOA9nATlrr+r6gFVTQG+Blao6hpVPQnMBgaVOP5RVc1U1fU4SezqYvu+UdU5qlrgiuVW4A+qmuwqbxrwM9df2Lk4P7A9VDVfVVer6jFXOQVAfxEJUtVUVf2xjLgvALaq6puqmqeq7+D8O1zkxvdSUd1leVFV96rqEZzvv/Az3wz8S1VXuMqZCZwEhhc79wXXuRUloydd/53sAf7Gmd9pSU+oarrr2EVU7d/a1CJLFqY8b+L85XoDJZqggDZAMLDa1TyRDnzm2g7OX9l7ix2/u7xKRCRYRP4lIrtdzTtLgHDXFU2h/cVeZwGhVfgcB4q9zi7jfcmySsYdWc4+cK6GZhf7DjYB+UA7nO9vPvCuq1nrryLi70qeP8f5az9VROaKSO8y4o6k9Pe2G+cKrlB530uZdZdRR2WfuQtwd+Hnc33GTlT8nVSl/LLU5N/a1CJLFqZMqrobp6P7fODDErsP4fzQ9lPVcNcjTFUL/8dOxflRKdS5gqruBnoBw1S1BTDatd1bbdUl495X7H3JoYN7cZqTwos9AlU1RVVzVfVRVe2L09R0Ia5+IFWdr6rjgQ44VwuvlhHHPpwf6+I6U+wKrzwV1V2O8j7zXuDxEp8v2HWVU1RdZfFUUL5pQCxZmIrcCIx1/TVcRFULcH7gnhORtgAi0lFECtuX38PpAO4rIsHAIxXU0Rwn8aS72ukrOrYuPOS62umH0x/w3wqOfRl4XES6AIhIGxG52PX6XBGJdV0hHcNpGsp3dbBPcvVdnARO4FyNlPQpECMi14iIn4j8HOgL/K+yD1Be3RWc8msRiXJ9/w8U+8yvAreJyDBxhIjIBSLSvLIYSpgqzkCGTsCdVPydmnrKkoUpl6puV9VV5ey+F9gGfOtqPvoc5woBVZ2H0zb9peuYLyuo5m84w3IP4XSsf+aZ6KvtK5yYvwCeVtWKbhp7HvgYWCAix3HiH+ba1x6nM/wYTvPUVziDBHxwrqb2AUdw+mduL1mwqh7GuSK4GzgM3ANcWPzelwqUV3d5/gMsAHa4Hn9yxbAKp9/iRZzO+204zZJV9RGwGqcjfy4wvRplGC+zm/KMwRk6i9Ps5l/BqJ5GR5ybMG9S1c+9HYup3+zKwhhjTKUsWRhjjKmUNUMZY4yplF1ZGGOMqVSDnvyrdevW2rVrV2+HYYwxDcrq1asPqWqbyo88rUEni65du7JqVXkjO40xxpRFRMqdVaE81gxljDGmUpYsjDHGVMqShTHGmEo16D4LY0z9kJubS3JyMjk5Od4OxRQTGBhIVFQU/v4VTTrsHksWxpgaS05Opnnz5nTt2hVb3K5+UFUOHz5McnIy0dHRNS6vYTdD7VsLz/WHde95OxJjmrScnBxatWpliaIeERFatWrlsau9WksWIjJDRA6KyIZi2/7rWs93rWv93rWu7V1FJLvYvpfdq0UhYy988htLGMZ4mSWK+seT/ya12Qz1Bs7UxkWrrKnqzwtfi8gzQEax47eravWWUMzNhi8eg7grqxepMcaYCtXalYWqLsGZr78U1xrOVwLvlLW/WjKSPVaUMaZhSU9P5x//+IdXY7jhhhuYNWtWlc5JS0tj2LBhDBo0iK+//po///nPtRRdzXmrz+Js4ICqbi22LVpE1ojIVyJydnknisgtIrJKRM68dTssqpZCNcZ42pw1KYx64kui75vLqCe+ZM6aSleLrVB9SBblyc8vf5HCL774gt69e7NmzRrOPvtsSxZluJozrypSgc6qOgj4HfAfEWlR1omq+oqqxqtqfNFG/yA47+HajNcY4yFz1qRw/4frSUnPRoGU9Gzu/3B9jRLGfffdx/bt2xk4cCBTp04F4KmnnmLIkCHExcXxyCOnV+u95JJLOOuss+jXrx+vvPJK0fbQ0FDuvfdezjrrLMaNG8fKlSsZM2YM3bp14+OPPy5Vp6pyxx130LdvXy644AIOHjxYtK9r16489thjJCQk8P777/Pqq68yZMgQBgwYwOWXX05WVhZr167lnnvu4dNPP2XgwIHce++9ZGdnM3DgQK699tpqfxe1pc6HzoqIH3AZcFbhNlU9ibMeMaq6WkS2AzGAexM/jX/c+iuMqSce/eRHNu47Vu7+NXvSOZVfcMa27Nx87pm1jndW7inznL6RLXjkon7llvnEE0+wYcMG1q5dC8CCBQvYunUrK1euRFWZNGkSS5YsYfTo0cyYMYOIiAiys7MZMmQIl19+Oa1atSIzM5MxY8bw5JNPcumll/Lggw+ycOFCNm7cyOTJk5k0adIZdc6ePZvNmzezfv16Dhw4QN++fZkyZUrR/sDAQJYuXQrA4cOHufnmmwF48MEHmT59Ov/3f//HY489xqpVq3jxxRcBeOmll4o+Q33jjfssxgE/qWpRJ4OItAGOqGq+iHQDeuKsBVyxNr2AFGgWVFuxGmM8rGSiqGx7dSxYsIAFCxYwaNAgAE6cOMHWrVsZPXo0L7zwArNnzwZg7969bN26lVatWtGsWTOSkpIAiI2NJSAgAH9/f2JjY9m1a1epOpYsWcLVV1+Nr68vkZGRjB079oz9P/950XgeNmzYwIMPPkh6ejonTpwgMTHRY5+1rtRashCRd4AxQGsRSQYeUdXpwFWU7tgeDTwmInlAPnCbqpbZOX4G/2Bo3gG2fAYDr/Fo/MaY6qnoCgBg1BNfkpKeXWp7x/Ag/nvrCI/EoKrcf//93HrrrWdsX7x4MZ9//jnffPMNwcHBjBkzpug+BH9//6Khpj4+PgQEBBS9zssre1n2ioamhoSEFL2+4YYbmDNnDgMGDOCNN95g8eLFNfl4XlGbo6GuVtUOquqvqlGuRIGq3qCqL5c49gNV7aeqA1R1sKp+4nZFPSfAti8h75SHP4ExpjZMTexFkL/vGduC/H2Zmtir2mU2b96c48ePF71PTExkxowZnDhxAoCUlBQOHjxIRkYGLVu2JDg4mJ9++olvv/222nWOHj2ad999l/z8fFJTU1m0aFG5xx4/fpwOHTqQm5vL22+/Xe5x/v7+5ObmVjum2tTwp/uISYLvZ8Ke5dBtjLejMcZU4pJBHQF4av5m9qVnExkexNTEXkXbq6NVq1aMGjWK/v37M3HiRJ566ik2bdrEiBHOlUpoaChvvfUWSUlJvPzyy8TFxdGrVy+GDx9e7TovvfRSvvzyS2JjY4mJieGcc84p99g//vGPDBs2jC5duhAbG3tGYivulltuIS4ujsGDB1eYVLyhQa/BHR8fr6uWfwVPRsOQGyHpL94OyZgmadOmTfTp08fbYZgylPVvIyKrzxhR6oaGPTcUQLMQiB4Nm+dBA058xhhTnzX8ZAEQkwhHd8Lhbd6OxBhjGqXGkyzAGRVljDHG4xpHsgjvDG37wZb53o7EGGMapcaRLMC5uti9HLLTvR2JMcY0Oo0oWSSB5sP2L7wdiTHGNDqNJ1lExUNQhDVFGWOqbMyYMaxa5d5UdABvvPEGd9xxRy1GVP80nmTh4+vczb11ARSUPyWwMaYeWPeesyTytPA6Wxq5oqnCTeUaT7IAp98i+ygkf+ftSIwx5Vn3nrMUcsZePLU08q5du+jduzeTJ08mLi6On/3sZ2RlZZWaKnzt2rUMHz6cuLg4Lr30Uo4ePVpUxltvvcXIkSPp378/K1euBGDlypWMHDmSQYMGMXLkSDZv3lx0/N69e0lKSqJXr148+uij1Y69oWj4030U130s+Pg5Q2g7V/82fmNMDcy7D/avL39/8neQf/LMbbnZ8NEdsHpm2ee0j4WJT1RY7ebNm5k+fTqjRo1iypQpRYshFZ8qPC4ujr///e+cc845PPzwwzz66KP87W9/AyAzM5Ply5ezZMkSpkyZwoYNG+jduzdLlizBz8+Pzz//nAceeIAPPvgAcBLJhg0bCA4OZsiQIVxwwQXEx1fppugGpXFdWQSFQ+cR1m9hTH1WMlFUtt1NnTp1YtSoUQBcd911RQmicKrwjIwM0tPTi+Zwmjx5MkuWLCk6/+qrrwacCQKPHTtGeno6GRkZXHHFFfTv35+77rqLH3/8sej48ePH06pVK4KCgrjsssuK6musGteVBTijohb8AY7uhpZdvB2NMU1PJVcAPNff1QRVQlgn+OXcaldbcrrwwvfFpwqv6vkPPfQQ5557LrNnz2bXrl2MGTOm0voaq8Z1ZQFOsgCno9sYU/+c97CzFHJxHlgaec+ePXzzzTcAvPPOOyQkJJyxPywsjJYtW/L1118D8Oabb54xU+x///tfAJYuXUpYWBhhYWFkZGTQsaMzG+4bb7xxRnkLFy7kyJEjZGdnM2fOnKKrmsaq8SWL1j0gortN/WFMfRV3JVz0gnMlgTjPF71Q46WR+/Tpw8yZM4mLi+PIkSP86le/KnXMzJkzmTp1KnFxcaxdu5aHHz6doFq2bMnIkSO57bbbmD59OgD33HMP999/P6NGjSo1miohIYHrr7+egQMHcvnllzfq/gpoDFOUlzU2+rMH4LtX4Z6dEBBa94EZ08R4e4ryXbt2ceGFF7JhwwavxVBf2RTlFYlJhPxTsPMrb0dijDGNQuNMFp1HQEALa4oypono2rWrXVXUssaZLPyaOfdcbFkABQXejsYYYxq8xpkswBkVdWI/7P/B25EYY0yDV2vJQkRmiMhBEdlQbNs0EUkRkbWux/nF9t0vIttEZLOIJNY4gJ7jAbEb9IwxxgNq88riDSCpjO3PqepA1+NTABHpC1wF9HOd8w8R8a1R7SGtIWqI9VsYY4wH1FqyUNUlwBE3D78YeFdVT6rqTmAbMLTGQcQkwr41cHx/jYsyxjRuNk15xbzRZ3GHiKxzNVO1dG3rCBS//z/Zta0UEblFRFaJyKq0tLSKa7K7uY2pl+bumMuEWROImxnHhFkTmLuj+tN8VEVjmaZcVSmo48E7dZ0s/gl0BwYCqcAzru1lTapS5t2CqvqKqsaranybNm0qrq1dP2gRZf0WxtQjc3fMZdryaaRmpqIoqZmpTFs+rUYJo7wpyoF6NU15165deeCBBxgxYgTx8fF8//33JCYm0r17d15++WUATpw4wXnnncfgwYOJjY3lo48+KvqMffr04fbbb2fw4MHs3buX6dOnExMTw5gxY7j55ptr9UqnTicSVNUDha9F5FXgf663yUCnYodGAftqXKGI0xT1w7uQmwP+gTUu0hhTsSdXPslPR34qd/+6tHWcKjh1xrac/BweXvYws7bMKvOc3hG9uXfovRXWW9YU5b///e+B+jVNeadOnfjmm2+46667uOGGG1i2bBk5OTn069eP2267jcDAQGbPnk2LFi04dOgQw4cPZ9KkSUWf8fXXX+cf//gH+/bt449//CPff/89zZs3Z+zYsQwYMKDCumuiTq8sRKRDsbeXAoUjpT4GrhKRABGJBnoCKz1SaUwS5GbC7sY9fbAxDUXJRFHZdneVN0U51K9pygt/+GNjYxk2bBjNmzenTZs2BAYGkp6ejqrywAMPEBcXx7hx40hJSeHAAefv7C5dujB8uLNWz8qVKznnnHOIiIjA39+fK664otrfnTtq7cpCRN4BxgCtRSQZeAQYIyIDcZqYdgG3AqjqjyLyHrARyAN+raqeaVyMPhv8gpymqB7jPFKkMaZ8lV0BTJg1gdTM1FLbO4R04PWk16tdb0VThtenacoDAgIA8PHxKXpd+D4vL4+3336btLQ0Vq9ejb+/P127diUnJ6fU56jref1qczTU1araQVX9VTVKVaer6vWqGquqcao6SVVTix3/uKp2V9VeqjrPY4H4B0G3Mc4Q2gY8aaIxjcWdg+8k0PfMJuFA30DuHHxnjcqtbIpyqLtpys877zxSUlKq9TkyMjJo27Yt/v7+LFq0iN27d5d53NChQ/nqq684evQoeXl5RU1jtaXxLX5UlphE2DIP0n6Ctt6bGdMYAxd0uwCA579/nv2Z+2kf0p47B99ZtL26Cqcov/XWW+nZs2eZU5SDM035bbfdRlZWFt26deP1109fzRROU37s2DFmzJgBONOUT548mWeffZaxY8eeUVbhNOXbtm3jmmuuIT4+noKCArZt20ZERES1Pse1117LRRddRHx8PAMHDqR3795lHtexY0ceeOABhg0bRmRkJH379iUsLKxadbpFVRvs46yzzlK3ZKSoPtJC9etn3TveGFMlGzdu9Gr9O3fu1H79+nk1hkLr16/Xu+66q07qOn78uKqq5ubm6oUXXqgffvhhqWPK+rcBVmkVf28b79xQxbWIhPZxNoTWGFPr+vfvz7PPPlsndU2bNo2BAwfSv39/oqOjueSSS2qtrqbRDAXOqKivn4asIxBcvctDY0z91FSnKH/66afrrK6mcWUBTrLQAtj2ubcjMaZRUhtAUu948t+k6SSLyEEQ0sYmFjSmFgQGBnL48GFLGPWIqnL48GECAz1zM3LTaYby8YGeifDTJ5CfC77+3o7ImEYjKiqK5ORkKp2vzdSpwMBAoqKiPFJW00kW4AyhXfsW7F0BXUuPwTbGVI+/vz/R0dHeDsPUoqbTDAXQ/Vzw8bemKGOMqaKmlSwCmjtXFDaE1hhjqqRpJQtwRkUd2gKHt3s7EmOMaTCaYLKY4DzbgkjGGOO2ppcsIrpB617Wb2GMMVXQ9JIFOKOidi2DnGPejsQYYxqEJposkqAgF3Ys8nYkxhjTIDTNZNFpGASG2agoY4xxU9NMFr5+0GO8kywKCrwdjTHG1HtNM1mA0xSVdQj2fe/tSIwxpt5rusmix3kgPjYqyhhj3NB0k0VwBHQabsnCGGPcUGvJQkRmiMhBEdlQbNtTIvKTiKwTkdkiEu7a3lVEskVkrevxcm3FdYaYRNi/HjKqt7C6McY0FbV5ZfEGkFRi20Kgv6rGAVuA+4vt266qA12P29ypYH1KBqOe+JI5a6r5Yx/jCm+rjYoyxpiK1FqyUNUlwJES2xaoap7r7bdAjSdaT0nP5v4P11cvYbTpBeFdbAitMcZUwpt9FlOAecXeR4vIGhH5SkTOLu8kEblFRFaJyKrCbdm5+Tw1f3PVIxBxri52LIZTWVU/3xhjmgivJAsR+QOQB7zt2pQKdFbVQcDvgP+ISIuyzlXVV1Q1XlXji2/fl55dvWBiEiEvB3Z9Xb3zjTGmCajzZCEik4ELgWvVtWCvqp5U1cOu16uB7UBMVcqNDA+qXkBdE8A/xEZFGWNMBeo0WYhIEnAvMElVs4ptbyMivq7X3YCewA53y23m68PUxF7VC8ovwFlBb8t8sMXmjTGmTLU5dPYd4Bugl4gki8iNwItAc2BhiSGyo4F1IvIDMAu4TVWPlFlwCb4+QqsQfyYNiKx+sDFJcCwFDmyo/FhjjGmC/GqrYFW9uozN08s59gPgg6rWEdsxjIeuHMCd767lox9SuHRQNQdX9XQtiLTlM2gfW70yjDGmEWvwd3BfFBdJv8gWPLNgCyfz8qtXSPN2EDnYhtAaY0w5Gnyy8PER7knqTfLRbP6zYk/1C4pJguRVcCLNc8EZY0wj0eCTBcDonq0Z0a0VL365jRMn8yo/oSwxiYDCtoUejc0YYxqDRpEsRIT7JvbmcOYpXl3i9iCqM3UYAM072BBaY4wpQ6NIFgADOoVzfmx7Xvt6B2nHT1a9ABGno3vbl5B3yvMBGmNMA9ZokgXA7yf0IievgBe/3Fq9AmKS4NRx2LPcs4EZY0wD16iSRbc2oVwZ34n/rNzDnsPVmOup2zngG2CjoowxpoRGlSwAfjuuJ74+wjMLqzGxYLMQiB4Nm+fZ3dzGGFNMo0sW7VoEMmVUNB+t3ceP+zKqXkBMIhzdCYe3eT44Y4xpoBpdsgC49ZzuhAX589fPqnF1EZPoPNuoKGOMKdIok0VYkD+/Prc7X21JY/n2Q1U7ObwztO1n/RbGGFNMo0wWAL8Y0ZUOYYE8+dlmtKr9DzGJsHs5ZKfXTnDGGNPANNpkEejvy13jYvhhbzqfbdhftZNjkkDzYfsXtROcMcY0MI02WQBcNrgjPdqG8tSCzeTlF7h/YlQ8BEVYU5Qxxrg06mTh5+vDPYm92JGWyfurk90/0ccXeo6HrQugoJoz2RpjTCPSqJMFwPi+7TirS0v+9vkWsk9V4Yc/JhGyj0Lyd7UXnDHGNBCNPlmICPcm9ebAsZO8vnyn+yd2Pw/E14bQGmMMTSBZAAyNjmBs77b8c/F20rPcnCQwKBy6jLR+C2OMwc1kISJtROQBEXlFRGYUPmo7OE+6J6kXJ07m8c/F290/KSYRDm6Eo7trLzBjjGkA3L2y+AgIAz4H5hZ7NBi927fg0kEdeWP5LlIzst07KSbJed66oPYCM8aYBsDdZBGsqveq6nuq+kHho1YjqwW/Gx+DKvxtoZtTmLfqARHdrN/CGNPkuZss/ici51e1cFdz1UER2VBsW4SILBSRra7nlq7tIiIviMg2EVknIoOrWl9loloGc93wLry/ei9bDxx35wM4Vxc7l8DJE54OxxhjGgx3k8WdOAkjW0SOichxETnmxnlvAEkltt0HfKGqPYEvXO8BJgI9XY9bgH+6GVuV3DG2B8HN/HhqvpuTDMYkQv4p2PlVbYRjjDENglvJQlWbq6qPqgapagvX+xZunLcEOFJi88XATNfrmcAlxbb/Wx3fAuEi0sG9j+G+iJBm3DK6Gws2HmD17qOVn9B5JDRrbk1Rxpgmze2hsyLSUkSGisjowkc162ynqqkArue2ru0dgb3Fjkt2bSsZxy0iskpEVqWlpVUrgBsTomkdGsCTn/1U+SSDfs2gx1jYssAWRDLGNFnuDp29CVgCzAcedT1P83AsUsa2Ur/OqvqKqsaranybNm2qVVFIgB93nteDlTuPsHizGwknJglO7IfUH6pVnzHGNHRV6bMYAuxW1XOBQUD1/qyHA4XNS67ng67tyUCnYsdFAfuqWUelrhramS6tgnnys5/IL6jkiqHHeEDsBj1jTJPlbrLIUdUcABEJUNWfgF7VrPNjYLLr9WScezgKt//CNSpqOJBR2FxVG/x9fbh7Qi9+2n+cj9amVHxwaBtnJlrrtzDGNFHuJotkEQkH5gALReQj3PirX0TeAb4BeolIsojcCDwBjBeRrcB413uAT4EdwDbgVeD2Kn2SargwtgP9IlvwzIItnMyrZJLBmETY9z0cP1DbYRljTL0jVV1FTkTOwbmb+zNVdXOipdoRHx+vq1atqlEZS7ak8YsZK3n4wr5MSYgu/8D96+HlBJj0Igy+vkZ1GmOMN4nIalWNr8o5FV5ZiEgL13NE4QNYDywFQqsdaT1yds/WjOrRihcXbeN4Tm75B7brDy06WlOUMaZJqqwZ6j+u59XAKtfz6mLvG7zCKcyPZJ7i1a8rmMJcxGmK2r4I8k7WXYDGGFMPVJgsVPVC13O0qnZzPRc+utVNiLUvLiqcC2I78NrXO0g7XkEiiEmC3EzYtbTugjPGmHqgsmaowRU96irIunD3hBhO5hXw9y8rmGQwejT4BdkQWmNMk1NZM9QzrsdLwArgFZyRSiuAF2o3tLrVrU0oPx/Sif+s2MPuw5llH+QfBN3Ocfot7G5uY0wTUlkz1Lmum/B2A4Ndd06fhXNT3ra6CLAiGw9vZMKsCczd4ZmlNe48ryd+vsIzC7aUf1BMIqTvhjQ3JyI0xphGwN37LHqr6vrCN6q6ARhYOyG5T1FSM1OZtnyaRxJGuxaB3JgQzcc/7GNDSkbZB/VMdJ5tVJQxpglxN1lsEpHXRGSMiJwjIq8Cm2ozsKrIyc/h+e+f90hZt57TnfBgf/5a3hTmYR2hfaz1WxhjmhR3k8UvgR9x5oj6LbDRta3e2J+53yPltAj059djerBkSxrLtx0q+6CYJNjWMd6HAAAgAElEQVT7LWSVnH3dGGMaJ3fXs8gBXgbuU9VLVfW5wrmi6ov2Ie09Vtb1I7oQGRZY/hTmMUmgBbDtC4/VaYwx9Zm7U5RPAtYCn7neDxSRj2szsKqa3G9y5Qe5KdDfl9+Oj+GH5AzmbSjjiiVyMAS3tn4LY0yT4W4z1CPAUCAdQFXXAl1rKSa3CULboLb4iz9LkpdUvpBRFVw+OIqYdqE8PX8zefkFZ+708XFGRW1bCPl5HqvTGGPqK3eTRZ6qljM8yHv6turLF1d+wb1D72X5vuXM2jrLY2X7+ghTE3uz41Am761KLn1ATCLkZMDeFR6r0xhj6it3k8UGEbkG8BWRniLyd2B5LcZVJVf0uoJhHYbx9HdPk3KikrUpqmBcn7bEd2nJ3z7fQvapElOYdzsXfPytKcoY0yS4myz+D+gHnMSZXDADZ2RUveAjPjw28jEAHln2CAVaUMkZ7hER7p3Ym4PHTzJjWYlJBgNbQNdRNoTWGNMkuJss+roefkAgcDHwXW0FVR2RoZFMHTKVFftX8N7m9zxW7pCuEZzXuy0vf7Wd9KwSy3fEJMGhzXBkh8fqM8aY+sjdZPE2MAO4DLjQ9biotoKqrst7Xs6oyFE8u/pZ9h7f67Fypyb14sTJPP6xePuZO2IK7+Ze4LG6jDGmPnI3WaSp6iequlNVdxc+ajWyahARpo2chp/48dCyhzzWHNW7fQsuGxTFG8t3sS89+/SOiG7QOsb6LYwxjZ7bQ2dd031cLSKXFT5qNbJqah/SnnuG3sPqA6v5z6b/VH6Cm+4a3xMU/vZ5iUkGYxKd9S1OHvdYXcYYU99UZbqPgUASTvPTRThNUfXSxd0vZnTUaJ7//nl2ZezySJlRLYO5fkQXZq1OZuuBYokhJgkKcp0V9IwxppFyN1kMcE1PPllVf+l6TKlOhSLSS0TWFnscE5Hfisg0EUkptv386pTvqoNHRjyCv68/Dy17iPyC/MpPcsOvz+1BSDO/MycZ7DQMAsNsVJQxplFzN1l8KyJ9PVGhqm5W1YGqOhA4C8gCZrt2P1e4T1U/rUk9bYPbcv/Q+1mbtpY3N75Z07ABiAhpxi2ju7Fw4wFW73ZNIujrDz3Gwdb5UOCZPhJjjKlv3E0WCcBaEdksIutEZL2IrPNA/ecB22urs/zCbhcyttNY/r7m7+xI98zw1hvPjqZ1aABPztt8enqRmCTITIN9azxShzHG1DfuJoskoCcwgdP9FZ4YOnsV8E6x93e4ktEMEWlZ1gkicouIrBKRVWlpaRUWLiI8NOIhgv2DeXDZg+QV1Hwep+Bmftw5ricrdx1h0eaDzsYe40B8bFSUMabRcneK8t1lPWpSsYg0AyYB77s2/RPojtORnoqz9ndZsbzi6j+Jb9OmTaX1tA5qzR+G/YH1h9bzxo9v1CTkIlcN6UTXVsE8OW8z+QUKwRFO34UlC2NMI+XulUVtmAh8r6oHAFT1gKrmq2oB8CrOLLcekdg1kfFdxvOPtf9g69GtNS7P39eHuyf0YvOB48xZ45qLKiYR9q+DY/tqXL4xxtQ33kwWV1OsCUpEOhTbdymwwVMViQgPDn+Q5s2a84elfyC3ILfGZV4Q24H+HVvw7MItnMzLd/otwEZFGWMaJa8kCxEJBsYDHxbb/NdiHefnAnd5ss6IwAgeHP4gm45sYvr66TUuz8dHuDepNynp2bz17R5o0xuCIuCz+2BaODzXH9Z5bo4qY4zxJj9vVKqqWUCrEtuur+16x3cZz8Toifzrh38xptMYekf0rlF5Z/dsQ0KP1ry0aBvXBH1L0MljUNiJnrEXPvmN8zruyhpGbowx3uXNZiiveGDoA4QFhPHg0gfJza95c9S9Sb05knmK3AXTTieKQrnZ8MVjNa7DGGO8rckli/DAcB4Z8Qibj27mX+v+VePyYqPCuCCuA6E5B8o+IKOMVfaMMaaBaXLJAuDczucyqfskXlv/Gj8e/rHG5f1+Qi9Sz2xVO615h7K3G2NMA9IkkwXAPUPuoVVgKx5c+iCn8k9VfkIFoluHsKzz7WRpszO2KzhNUQdqnpCMMcabmmyyCAsI45GRj7AtfRv//OGfNS4vp89l3Jd7E8kFrSlQIbmgNX8r+DnZ6gczkmxWWmNMg+aV0VD1xeio0Vza41JmbJjB2E5jiW0TW+2y/rVkJykFCXx8KuGM7UuCxjE77Dl4+2dw0Qsw6Nqahm2MMXWuyV5ZFJo6ZCptg9vyh2V/ICcvp9rlnLGCXjFrM0JgymfQNQE+uh0W/RkKJyA0xpgGoskni+bNmvPoyEfZmbGTl9a+VO1yIsODyt8eGAbXzoKB18JXT8Ls2yCvZv0kxhhTl5p8sgAYGTmSK2KuYOaPM1l7cG21ypia2Isgf99S268Z1tl54esPF78E5/4B1r0Lb10G2ek1CdsYY+qMJQuXu+PvpkNIBx5c9iDZeWU3KVXkkkEd+ctlsXQMD0KAdi0CaBHox+vLdrLt4AnnIBE45x649BXY8y3MSISjtbKUhzHGeJRoA24/j4+P11WrVnmsvBWpK7hpwU1c1+c67h16b43L2552gp//6xv8fHx4/7YRdIoIPr1z5xJ49zrwC4Br/gsdB9e4PmOMcYeIrFbV+KqcY1cWxQzrMIyrel3FW5ve4rv939W4vO5tQnnrpmHk5OVz9avfkppR7IolejTcuAD8AuGNC2CzrYVhjKm/LFmUcNdZdxEVGsXDyx4mKzerxuX1bt+Cf08ZSkZWLte+uoK04ydP72zbG276HFrHwLtXw8pXa1yfMcbUBksWJQT7B/OnhD+RciKF51Y/55Ey46LCef2XQ0jNyOG611ZwNLPYSKjm7eCXn0LPRPj09zD/D1BQ4JF6jTHGUyxZlOGsdmdxbZ9reXfzu6xIXeGRMuO7RvDa5Hh2Hs7kFzNWciyn2Iy3zULgqrdhyM3wzYsw6wZnmhBjjKknLFmU4zeDf0OXFl14eNnDnDh1wiNljurRmpevG8xP+4/xy9e/I/NksSnNfXzh/KdgwuOw8WOYOQkyD3mkXmOMqSlLFuUI8gviT6P+xP6s/Tyz+hmPlTu2dzuev2oQa/Yc5eZ/ryInN//0ThEYeQdcOdNZz3v6eDi83WN1G2NMdVmyqMDAtgOZ3Hcys7bMYnnKco+Ve35sB565cgDf7DjMr95azam8En0UfS+GyZ9ATga8Ns65J8MYY7zIkkUlfj3o10SHRfPw8oc5fuq4x8q9dFAUj18Sy6LNadz57hry8kskjE5D4caFENTSaZL6cbbH6jbGmKqyZFGJAN8AHh/1OGnZafz1u796tOxrhnXmoQv7Mm/DfqbOWkdBQYkbJFt1d4bWRg6C92+AZc/bJITGGK/wWrIQkV0isl5E1orIKte2CBFZKCJbXc8tvRVfcbFtYpnSfwpzts1hSfISj5Z9Y0I0v58Qw+w1KfxhzgZK3VEfHAG/+Aj6XQoLH4a5v4P8vLILM8aYWuLtK4tzVXVgsdvO7wO+UNWewBeu9/XCrwb8ih7hPZi2fBoZJzM8WvYdY3vy63O7887KPfzxf5tKJwz/QLh8Boy6E1bNgHevgZOeGaFljDHu8HayKOliYKbr9UzgEi/GcoZmvs14POFxjuQc4cmVT3q8/N9P6MUvR3VlxrKdPLNgS+kDfHxg/GNwwbOwbSG8PhGOpXo8DmOMKYs3k4UCC0RktYjc4trWTlVTAVzPbUueJCK3iMgqEVmVlpZWh+FC31Z9uTnuZj7Z8Qlf7vnSo2WLCA9f2Jerh3bixUXbeGnRtrIPHHIjXP1fZ0jta+PgwEaPxmGMMWXxZrIYpaqDgYnAr0VktDsnqeorqhqvqvFt2rSp3QjLcEvsLfRq2YvHvnmM9BzPrkchIvzpklguGRjJU/M3M33pzrIPjJngTBFSkOdMc75jsUfjMMaYkryWLFR1n+v5IDAbGAocEJEOAK7ng96Krzz+vv48nvA4GScz+POKP3u8fF8f4ekrBjCxf3v++L+NvLNyT9kHRg50RkqFRcFbl8Pa/3g8FmOMKeSVZCEiISLSvPA1MAHYAHwMTHYdNhn4yBvxVaZXRC9uHXAr83bNY8GuBR4v38/Xh+evGsS5vdrwwOz1zF6TXPaB4Z1Or+8951ew6C82tNYYUyu8dWXRDlgqIj8AK4G5qvoZ8AQwXkS2AuNd7+ulG2NvpG+rvvzp2z9xOPuwx8tv5ufDP687ixHdWnH3ez8wb305ndmBYXDN+671vZ+AObfb+t7GGI+zlfJqYOvRrfz8fz9nTKcxPHPOM4iIx+vIPJnHL2asZF1yOq9cH8+5vUv1+TtUYclTsOhxZ2GlK9+EoHCPx2OMafhspbw61rNlT24feDsLdy/ks121s9JdSIAfr/9yCL3bt+DWt1azfFs5M9EWre/9L9j9DcxIgvRy+juMMaaKLFnU0A39biC2dSyPr3icQ9m1M6V4i0B//j1lKNGtQrjp36tYtetI+QcPuAqu+wCO7XOG1n71FDzXH6aFO8/r3quVGI0xjZs1Q3nAjvQdXPHJFXQP7076yXT2Z+6nfUh77hx8Jxd0u8Bj9Rw8nsNV//qWtOMnefvmYcRFVdDMdPAn58a97BKJxT8ILnoB4q70WFzGmIbFmqG8pFt4N8Z1GcemI5tIzUxFUVIzU5m2fBpzd8z1WD1tmwfy1k3DaBHkzy9mrOSn/ccqOLg3+AWW3p6bDV885rGYjDFNgyULD1lzcE2pbTn5OTz//fMerScyPIh3bh5OgJ8P1722kh1pFcwRdbycEVQZ5QzFNcaYcliy8JD9mfvL3J6amcrSlKXk5OV4rK7OrYJ5+6bhqCrXvraCvUeyyj4wLKr8Qr74I2RV0PdhjDHFWLLwkPYh7cvd96vPf8XZ757Nrz7/FW9veps9x2o+SqlH21DeumkYWafyuea1b0nNyC590HkPO30UxfkFQuRg+PppeH4ALPozZHt22hJjTONjHdweMnfHXKYtn0ZO/ukriEDfQB4Y9gCtg1qzbN8ylqYsZfex3QB0bt6ZhI4JJHRMIL59PEF+QeUVXaEf9qZz7WsraNsigP/eMoI2zQPOPGDde04fRUayc6Vx3sNO5/b+Dc5NfJs+gYAwGPFrGH6bc5OfMaZRq04HtyULD5q7Yy7Pf/98haOh9h7by9cpX7M0ZSnf7f+OnPwcAnwDiG8XX5Q8urToUqUb/L7bdYRfTF9Jl1bBvHvLcMKDm7kfdOo6WPwEbJ4LgeEw8v9g2K0Q0Nz9MowxDYoliwYmJy+H1QdWszRlKUtTlrLr2C4AokKjSOiYwNlRZzOk/RC3rjqWbj3ElJnf0bt9c2fEVKB/1YLZt8ZJGls+g6AIGPUbGHIzBIRW45MZY+ozSxYN3N7je1mW4jRXrdy/kuy8bJr5NOOsdmc5Vx1RCUS3iC73quOLTQe49c3VDOwUzr9vHEpwM7+qB5G8Ghb/xVlgKbgVjPotDLkJmgXX8NMZY+oLSxaNyMn8k2dcdezMcNa26Bjasai5amj7oQT7n/kjPnddKv/3zveM7N6a1ybHE+jvW70A9q50ksb2LyGkDSTcBfFTSneYG2MaHEsWjVjKiRSWJjuJY8X+FWTnZePv43/6qqNjAt3CuiEifLA6mbvf/4Gxvdvy8nVn0cyvBoPe9nzrjJja+RWEtoOE38FZNzjrghtjGiRLFk3EqfxTfH/w+6LksT1jOwCRIZGM6jiKhI4J7EqO5LGPt9Ot608c8v+IAt+j+OS35PLom3lk7PVVr3TXMidp7F4KzTvA2XfD4F+AX0Dl5xpj6hVLFk3UvhP7ipqrVqSuICsvCz8fP5oVtCWzYD/iU1B0rBb4c0WXu6qXMAB2LnGSxp5voEVHJ2kMuh78qjACyxjjVZYsDLn5uaw5uIalKUt5fcNMkIJSx2huOLd1m0GniGCiWgbRqWUwrUMD8PFxc7iuqrPu96I/Q/JKCOsEo3/vLMDkW8VRWMaYOmfJwpyh/xuxlDVwShVOHpxI3rEBaJ4zc20zPx+iWgYR1fJ0AolqGVSUUFqFNCs9CksVtn/hJI2U1RDexVlTI+4q8K3GSCxjTJ2wZGHOEDd9NOp3tPQO9QXJB6Bb8/50C0qgeX48hzL82Hskm+SjWRzNyj3jlCB/X1cyOZ1AoloGO0klPJDwfYuRRX+G1LXQMhrOuRdirygzacxZk8JT8zezLz2byPAgpib24pJBHWvjKzDGlMGShTnDo1++yfu7n0N8Tv/wF/ZZTIkfw7xd8/h0x6dsz9iOr/gyvMNwJkZP5LzO56EFAaSkZ5N8JJu9R7NIPuokkcJkciwn74y6QgP8iAoP5ILAH/jZsX/TIXsrmaFdOTrkLprHX0VYiDN6as6aFO7/cD3ZuflF5wb5+/KXy2ItYRhTRyxZmFIe/fJNPtj5armjoVSVLUe3MG/nPD7b9RkpJ1Jo5tOM0VGjmRg9kdFRowksY12MjOxckl1JZO+R08nEeZ9JQt4Kfuv3AX189rCtIJJXfK5gY/hYth3OJie3dD9Kx/Aglt03tla/C2OMw5KFqRFV5Ye0H5i3cx7zd83ncM5hQvxDGNtpLBOjJzI8cjj+PpV3YKsq6Vm5JB/JJHfDHLpu+DsRmdtJ8e/C45kX408uU/3eJ1IOsU9b89e8K/mkIIGdT3huVUFjTPkaRLIQkU7Av4H2QAHwiqo+LyLTgJuBNNehD6jqpxWVZcmi9uQV5PHd/u+Yt3Men+/+nOO5x2kZ0JIJXScwMXoig9oOwkfcvNmvoAA2zobFT8KhzRSo4COn/7vL0mY8mH8LF//it4zu2bpKkygaY6quoSSLDkAHVf1eRJoDq4FLgCuBE6r6tLtlWbKoG6fyT7E0ZSnzds5j8d7F5OTn0D6kPUldk5gYPZE+EX3c+4EvyOfkE90IOFV6/YxUWjEi5+/079iC28f0ILFfe3zdHcprjKmSBpEsSgUg8hHwIjAKSxb1XlZuFov2LmLeznksS1lGnubRtUVXJkZPZGL0RKLDoisuYFo4UPq/OQUOthrKnON9+PB4H3Jb9ea2MT24ZGDHmk1XYowppcElCxHpCiwB+gO/A24AjgGrgLtVtdS4TxG5BbgFoHPnzmft3r27jqI1JaXnpLNwz0Lm7ZzHqv2rUJQ+EX2KEkeZqwc+1x8y9pbeHtAcwjrDwR8BSJPWfJ4by7rAIfQ7exKXj+hLULNqTopojDlDg0oWIhIKfAU8rqofikg74BDOH5l/xGmqmlJRGXZlUX8cyDzA/F3zmbdzHhsObwBgcNvBnB99PuO7jiciMMI5cN178MlvILfYMrD+QXDRC84Kfsf2wbbP0a0Lyd/2JX65J8hVX36QXpyKHkvsmJ/RvPNAyrzb0BjjlgaTLETEH/gfMF9Vny1jf1fgf6rav6JyLFnUT3uO7WHeznnM2znv9D0ckcM5P/p8xnYay1fLn+T5HbPZ7wPtC+DObpdywZg/li4oPxf2riR19Sfk/jSfzrk7ADju3wa/XuMI6pME3c+1pWCNqaIGkSzE6QmdCRxR1d8W295BVVNdr+8ChqnqVRWVZcmifivrHg5f8UVRCvT0vRaBvoFMGzmt1BK0JW3eupnvPp9Fy31fcbbPelpIFiq+SKdh0HMc9BgP7WPtqsOYSjSUZJEAfA2sxxk6C/AAcDUwEKcZahdwa2HyKI8li4aj8B6OWxfeSlZeVqn9Ab4BXBFzBR1DO9IxtCORoZF0DO1IaLPSy7ruPJTJq4s3s2PtYkbLD0wK/pGok1udnaHtoMc459H9XAhqWdsfzZgGp0EkC0+yZNHwxM2MQ8sYDQUQ5BdEdl72GdvCAsKIDIkkqnkUkSGRRIaefu1TEMHb3+znPyv2EJp7mNs77eSS0I20TF0KORkgvhA1pNhVRxz42MgqYyxZmHpvwqwJpGaWvmDsENKB+ZfPJ/1kOiknUooe+07sO+P1yfyTZ5wXERhB26AOnMwOY+eBAHKywujTqhO3xwRxTu42grYvdiY3BAhp67rqOA+6j4XgiKJy5i5+yL1+FGMaAUsWpt6bu2Mu05ZPIyc/p2ibu30WqsrhnMNlJpHC59yCM2fLbe4fQXRoOzqq0DHzKJGHdtIx6xgd8wvo0G4AzXpOYO6JHUw7uJScYjcBBhYo06ItYZjGyZKFaRDm7pjL898/z/7M/bQPac+dg++sNFG4o0ALSMtKY1fGXj7a8CMLt27iWN4BQkOOERJ8jON5h8jT07PlikKb/DyO+viQW0bzVId8ZcGUDTWOy5j6xpKFMcXk5Rcwd30q/1y8nZ/2HycyPIBrRrZgcHc4lLOflOMppKTv4KNd88oeQaXKf4Y8Qr/el+JjizmZRsSShTFlUFW+/OkgLy3axvd70mkd2owpCdFcN7wLLQL9mTCjP6m+5SwpKELLAmVkszYkdExgVNwNtIzoXvcfwhgPsmRhTAVUlRU7j/CPxdtZsiWN5oF+/GJEF6J93uTJ5I9K9VlMjRxHSEALlqZ8zfJTaRzxEUSV/urPqPDeJPS8mP69L8PXr5kXP5UxVWfJwhg3rU/O4J9fbWPehv34CpzV4r8car2aND+hTZ7S7kg8V47/c9HqfQX5eWzcPJuvt8xh2dFNrJdTFIgQXqCM8G/F2R0TGBk7mVatY7z8yYypnCULY6po28ETTHpxKVmn8kvta9M8gK+mjiG4Wen+ivSjO/lm3UyWpixh6cmDHHFdlfQr8GVUeC/O7jGJ2D5X2FWHqZcsWRhTDdH3zS3nNkFHx/AgurcNpXubEHq0DaV7m1B6tA2lVUgzRISC/Dw2bf2YpZtns+zoRn7gJAUitChQRvpHkBA5klFxN9C6de86+0zGVKQ6ycKGeJgmLzI8iJT07FLbI4L9+eWoaLannWBb2gm+23mE7NzTVyBhQf6u5BFCj7YD6BEziqS2obTwOcyK9W+wNHkJy3L289neubB3Ln0KfEkIi3FddfwMP//Sa5sbU1/ZlYVp8uasSeH+D9efkQiC/H35y2WxRX0WAAUFSuqxHLYfPMG2gyecJHLwBNvTMjl04vSd5c38fOjWOoTubULp1jqQtrKCQxkLWZu5mR84Sb4IzQuUEf4tSegwgoTYybRp269OP7Np2qwZyphqmrMmhafmb2ZfejaR4UFMTex1RqKoTEZWLtvSnASyvVgi2XMki4Ji/4v1DM+kf/gSTvmuY6PvEQ65huz2LvAhoUVPEnpcxIC+P8fPP5AX3v8dH2XML+p0vzgskd9cUWpG/zpXX+My7rNkYUw9czIvn12HsopdhZx+zsnNo3vgGiJDV5Admsz2wHzyRAgtKCAyz5ed/gXkSvHhvAVcH+rdH+YX3v8db56YT06xO97rQ1ymaixZGNNAFBQo+zKy2Z6WWZQ89uzfhZz4HwSsZ12LYxSUcVe5nyodcssosI6k+kNeGXG1zS3gi5t+9EJEpjqsg9uYBsLHR4hqGUxUy2DOiWnj2hoLXER61ilGvze4zPPygPYaUldhlrKXzDK3p/nZglONnSULY+qZ8OBmtMlTDvqX8Rd8njLjlhVeiMpx3mv9yoyrTV7DbaEw7rGVYIyphy4OSySwoOCMbYEFBVwcluiliBz1NS5T+yxZGFMP/eaKZ7k+NJG2uQWIKm1z60cncn2Ny9Q+6+A2xpgmpjod3HZlYYwxplL1LlmISJKIbBaRbSJyn7fjMcYYU8+ShYj4Ai8BE4G+wNUi0te7URljjKlXyQIYCmxT1R2qegp4F7jYyzEZY0yTV9+SRUdgb7H3ya5tRUTkFhFZJSKr0tLS6jQ4Y4xpqurbTXll3QZ6xnAtVX0FeAVARI6LyOa6CKyKWgOHvB1ECRaTeywm99XHuCwm9/Sq6gn1LVkkA52KvY8C9lVw/OaqDv+qCyKyqr7FZTG5x2JyX32My2Jyj4hU+Z6D+tYM9R3QU0SiRaQZcBXwsZdjMsaYJq9eXVmoap6I3AHMB3yBGapqU1kaY4yX1atkAaCqnwKfunn4K7UZSw3Ux7gsJvdYTO6rj3FZTO6pckwNeroPY4wxdaO+9VkYY4yphyxZGGOMqVSDTBYiMkNEDorIBm/HUkhEOonIIhHZJCI/isid9SCmQBFZKSI/uGJ61NsxFRIRXxFZIyL/83YshURkl4isF5G11RlaWBtEJFxEZonIT67/tkZ4OZ5eru+n8HFMRH7rzZhccd3l+m98g4i8IyKB9SCmO13x/OjN76is30sRiRCRhSKy1fXcsrJyGmSyAN4AkrwdRAl5wN2q2gcYDvy6HsxrdRIYq6oDgIFAkogM93JMhe4ENnk7iDKcq6oD69G4+OeBz1S1NzAAL39nqrrZ9f0MBM4CsoDZ3oxJRDoCvwHiVbU/zkjKq7wcU3/gZpwpjAYAF4pITy+F8walfy/vA75Q1Z7AF673FWqQyUJVlwBHvB1Hcaqaqqrfu14fx/mfumPFZ9V6TKqqJ1xv/V0Pr49oEJEo4ALgNW/HUp+JSAtgNDAdQFVPqWq6d6M6w3nAdlXd7e1AcEZ2BomIHxBMxTfz1oU+wLeqmqWqecBXwKXeCKSc38uLgZmu1zOBSyorp0Emi/pORLoCgwDvLZbs4mruWQscBBaqqtdjAv4G3AMUVHZgHVNggYisFpFbvB0M0A1IA153Ndm9JiIh3g6qmKuAd7wdhKqmAE8De4BUIENVF3g3KjYAo0WklYgEA+dz5uwU3tZOVVPB+UMXaFvZCZYsPExEQoEPgN+q6jFvx6Oq+a4mgyhgqOvy2GtE5ELgoKqu9mYc5RilqoNxpsj/tYiM9nI8fsBg4J+qOgjIxI3mgrrgmmFhEvB+PYilJc5fytFAJBAiItd5MyZV3QQ8CSwEPgN+wGmqbrAsWXiQiPjjJIq3VfVDb8dTnKv5YjHe7+sZBUwSkV04U9CPFZG3vBuSQ1X3uZ4P4rTDD/VuRCQDycWuBmfhJI/6YCLwvaoe8HYgwDhgp6qmqWou8AjpchIAAAJnSURBVCEw0ssxoarTVXWwqo7GaQba6u2YijkgIh0AXM8HKzvBkoWHiIjgtC1v0v9v7+5ZowijMAzfj6QxQRA/ESSks9PGQtBCCPoDBIsUEkRBUPwB2lgJ/gNBsAj4AUEROysRUdBChQhqZWEs7ASxMsixmAmGEBhlozOB+2r2XYZdDgszz847M+etGsTq9Ul2JtnajjfT7FQf+qypqi5V1d6qmqKZxnhcVb3+CwRIMpFky/IYOE4zldCbqvoCLCZZ7hA6DbzrsaSVZhjAFFTrE3AoyXi7H04zgJsnkuxqXyeBEwzn94Km595sO54FHnZ9YHDtPv5EkrvAUWBHks/Alaq62W9VHAZOAW/bawQAl9v2JX3ZA8y1KxBuAuarajC3qg7MbuBBc6xhDLhTVY/6LQmAi8DtdtrnI3C653po5+CPAef6rgWgql4muQe8ppnqecMwWmzcT7IdWAIuVNXXPopY63gJXAPmk5yhCduTnd9juw9JUhenoSRJnQwLSVInw0KS1MmwkCR1MiwkSZ0MC2kESabW6n7ctubou5GktG425HMW0tBV1dm+a5DWk2cW0ujGkswlWWjXnhhP8iTJQYAk35NcbdcVeZFkd98FS3/LsJBGtw+4UVX7gW/A+VXbJ2jaVR8AntKscyBtKIaFNLrFqnrejm8BR1Zt/wEst1l5BUz9p7qkdWNYSKNb3TNn9ful+t1X5ydeK9QGZFhIo5tcsTb2DPCsz2Kkf8GwkEb3HphNsgBsA673XI+07uw6K0nq5JmFJKmTYSFJ6mRYSJI6GRaSpE6GhSSpk2EhSepkWEiSOv0CakExC0R2w4oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1518291710>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Figure 2: Median impressions per bin, random click model\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW9//HXhxAIOwFxYZGAsgoRMCqCRdyxKOq1trW2P7z2utxeW663df1ZRXsXe+2vrd7aWntRvOrVWqrWFtu6F9cqIFWUGUFACCDbJGGNJOHz++NM4hAmySRkOLO8n49HHpk5c+acz0xg3nO+33O+X3N3REREGusQdgEiIpKZFBAiIpKUAkJERJJSQIiISFIKCBERSUoBISIiSSkgJG3MzM3s6Pjt+8zs+2nYxxfMLNre2z1YzOxSM3su7DraU+LfPcljOfd6c5npOggxs9VAf6C/u29JWL4EOBYY4u6r27BdB4a5+4p2KlWygP7uuUNHEFJvFXBJ/R0zGwt0Ca+czGNmHcOuIZPo/ch9Cgip9zDwfxLuzwT+J3EFM+tsZj8yszVmtjHebNQl4fHrzGyDma03s8sbPXeumf1r/Haxmf3BzDabWUX89sCEdV8xsx+Y2etmtt3MnjOzQ5IVbWZTzaw84f7qeB3vmdlOM5tjZoeZ2R/j23rBzIrj65bEm0OujNe8wcy+m7Ct2WY2z8weMbNtwGVm1sHMbjSzj81sq5k9YWZ94usXxdfdamaVZvaOmR0Wf+wyM1sZr2GVmV2asPy1hH1Oij+vKv57UirvS3P7TvKerTazm8zsw/j7/6CZFSU8fq6ZLYlv5w0zK2303BvM7D1gZzMh8cX4691iZneZWYcmXq+b2dVmtjxey71mZk1sUw4yBYTUewvoaWajzKwA+ArwSKN1fggMB8YBRwMDgFsBzGwa8D3gTGAYcEYz++oAPAgMBo4EdgM/a7TO14C/Bw4FOsW3naqL4nUMB84D/gjcDBwS3/d3Gq1/arzms4AbzSyx9vOBeUBv4NH4cy8ATiFolqsA7o2vOxPoBQwC+gJXA7vNrBtwD3COu/cAJgFLGhcdD5r58XX7Aj8G5ptZ34TVmnpfku67mffoUuBs4Kj4+3RLvIYJwAPAVfHt/BJ4xsw6Jzz3EmA60Nvda5vY/oVAGTCB4D28vIn1AM4FjidozvxyvC7JAFkZEGb2gJltMrOlKax7pJm9bGbvxr9VfvFg1Jil6o8izgQiwLr6B+Lf6q4ArnX3mLtvB/4d+Gp8lS8DD7r7UnffCcxuaifuvtXdf+vuu+Lb+TeCD9xED7r7R+6+G3iCIJRS9V/uvtHd1wGvAn9193fd/TPgKWB8o/Vvd/ed7v4+QXBdkvDYm+7+tLvvjddyFfB/3b08vr3ZwJfi36RrCD5Uj3b3Ondf5O7b4tvZC4wxsy7uvsHdP0hS93Rgubs/7O617v4Ywd/hvBTel+b2nczP3H2tu8cI3v/613wF8Et3/2t8Ow8BnwETE557T/y5zQXQD+P/TtYAP2Xf97SxO929Mr7uy7Tuby1plJUBAcwFpqW47i3AE+4+nuDD7OfpKioHPEzwDfUyGjUvAf2ArsCieNNDJfCn+HIIvk2vTVj/k6Z2YmZdzeyXZvZJvOlmAdA7fuRS79OE27uA7q14HRsTbu9Ocr/xthrX3b+JxyA46nkq4T1YBtQBhxG8f38GHo83Wf2nmRXGA/MrBN/qN5jZfDMbmaTu/uz/vn1CcKRWr6n3Jem+k+yjpdc8GPhu/euLv8ZBNP+etGb7yRzI31rSKCsDwt0XALHEZWZ2lJn9ycwWmdmrCf8BHegZv90LWH8QS80q7v4JQWf1F4EnGz28heDD9Rh37x3/6eXu9f+ZNxB8kNQ7spldfRcYAZzo7j2BKfHlYbU9N6478d9I49P81hI0FfVO+Cly93XuXuPut7v7aIJmpHOJ9+u4+5/d/UzgCIKjgl8lqWM9wQd0oiNJOJJrSnP7bkJTr3kt8G+NXl/X+NFMw+5aqqeZ7UsWycqAaML9wLfd/TiCdtn6I4XZwNct6Mh8Fvh2OOVljW8Cp8W/9TZw970EH2o/MbNDAcxsgJnVtxc/QdCJO9rMugK3NbOPHgRhUxlvd29u3YPh+/GjmmMI2vd/3cy69wH/ZmaDAcysn5mdH799qpmNjR8JbSNo9qmzoJN8Rrwv4jNgB8FRR2PPAsPN7Gtm1tHMvgKMBv7Q0gtoat/NPOWfzGxg/P2/OeE1/wq42sxOtEA3M5tuZj1aqqGR6yw4GWEQMIvm31PJUDkREGbWneBb028sOHf/lwTf1CBo+5zr7gMJvhk/XH9GhezP3T9294VNPHwDsAJ4K9409ALBkQDu/keCtuaX4uu81MxufkpwCu0Wgs7xP7VP9W32F4KaXwR+5O7NXch1N/AM8JyZbSeo/8T4Y4cTdGhvI2h6+gtBR38HgqOm9QRHvqcA32q8YXffSvDN/7vAVuB64NzEa1Oa0dS+m/K/wHPAyvjPv8ZrWEjQD/Ezgg74FQRNjq31O2ARQWf8fGBOG7YhIcvaC+XMrAT4g7uPMbOeQNTdj0iy3gfANHdfG7+/Epjo7psOZr2SeeL/hlYBhc2cjZNzLLgw8h/c/YWwa5HMlhPfpONna6wys4shOOPGzI6NP7wGOD2+fBRQBGwOpVARkSySlQFhZo8BbwIjzKzczL5JcF73N83sb8AHBOdeQ3C4fkV8+WPAZZ6th00iIgdR1jYxiYhIemXlEYSIiKRf1g22dcghh3hJSUnYZYiIZJVFixZtcfd+La/5uawLiJKSEhYubOosTBERScbMmhzdoClqYhIRkaQUECIikpQCQkREksq6PggRCV9NTQ3l5eVUV1eHXYo0UlRUxMCBAyksbG4w39QoIESk1crLy+nRowclJSVoArjM4e5s3bqV8vJyhgwZcsDbS1sTU0uT+sSHw7jHzFbEJ/KZkK5aJHzzV87nrHlnUfpQKWfNO4v5K+eHXZIcgOrqavr27atwyDBmRt++fdvtyC6dfRBzaX5Sn3MIpnkcBlwJ/CKNtUiI5q+cz+w3ZrNh5wYcZ8PODcx+Y7ZCIsspHDJTe/5d0hYQySb1aeR84H888BbBjGL7jcYq2e/uxXdTXbfvN5rqumruXnx3SBWJSCrCPItpAPtOS1jOvlMrNjCzK81soZkt3LxZA7Fmm093ftqq5SItqays5Oc/D3f24Msuu4x58+a16jmbN2/mxBNPZPz48bz66qv8+7//e5qqax9hBkSy46CkIwe6+/3uXubuZf36tepKcckAh3c7vFXLJfc8/e46Jt/5EkNunM/kO1/i6XdbnEW1WZkQEE2pq2t6Ir8XX3yRkSNH8u677/KFL3xBAdGMcvadt3Ygmrc2J82aMIuOtu8Jc0UFRcyaMCukiuRgevrdddz05Pusq9yNA+sqd3PTk+8fUEjceOONfPzxx4wbN47rrrsOgLvuuovjjz+e0tJSbrvt81lsL7jgAo477jiOOeYY7r///obl3bt354YbbuC4447jjDPO4O2332bq1KkMHTqUZ555Zr99ujvXXHMNo0ePZvr06Wza9PmcYyUlJdxxxx2cfPLJ/OY3v+FXv/oVxx9/PMceeywXXXQRu3btYsmSJVx//fU8++yzjBs3jhtuuIHdu3czbtw4Lr300ja/F+kU5mmuzwDXmNnjBFM2Vrn7hhDrkTSZPnQ6Dy59kI8rP6bO6zi82+HMmjCL6UOnh12atIPbf/8BH67f1uTj766pZE/d3n2W7a6p4/p57/HY22uSPmd0/57cdt4xTW7zzjvvZOnSpSxZsgSA5557juXLl/P222/j7syYMYMFCxYwZcoUHnjgAfr06cPu3bs5/vjjueiii+jbty87d+5k6tSp/PCHP+TCCy/klltu4fnnn+fDDz9k5syZzJgxY599PvXUU0SjUd5//302btzI6NGjufzyyxseLyoq4rXXXgNg69atXHHFFQDccsstzJkzh29/+9vccccdLFy4kJ/97GcA3HvvvQ2vIROlLSDik/pMBQ4xs3KCiekLAdz9PoIJ2r9IMOftLoLJ4iVHVXxWwbQh0/iPL/xH2KXIQdY4HFpa3hbPPfcczz33HOPHjwdgx44dLF++nClTpnDPPffw1FNPAbB27VqWL19O37596dSpE9OmBSdajh07ls6dO1NYWMjYsWNZvXr1fvtYsGABl1xyCQUFBfTv35/TTjttn8e/8pWvNNxeunQpt9xyC5WVlezYsYOzzz673V7rwZS2gHD3S1p43IF/Stf+JXPEqmNs2rWJkX1Ghl2KpEFz3/QBJt/5Eusqd++3fEDvLvz6qpPapQZ356abbuKqq67aZ/krr7zCCy+8wJtvvknXrl2ZOnVqwzUChYWFDaeEdujQgc6dOzfcrq1NPkV5c6eQduvWreH2ZZddxtNPP82xxx7L3LlzeeWVVw7k5YVGYzFJ2kVjUQAFRJ667uwRdCks2GdZl8ICrjt7RJu32aNHD7Zv395w/+yzz+aBBx5gx44dAKxbt45NmzZRVVVFcXExXbt2JRKJ8NZbb7V5n1OmTOHxxx+nrq6ODRs28PLLLze57vbt2zniiCOoqanh0UcfbXK9wsJCampq2lxTummoDUm7+oAYUdz2DwTJXheMD85ev+vPUdZX7qZ/7y5cd/aIhuVt0bdvXyZPnsyYMWM455xzuOuuu1i2bBknnRQckXTv3p1HHnmEadOmcd9991FaWsqIESOYOHFim/d54YUX8tJLLzF27FiGDx/OKaec0uS6P/jBDzjxxBMZPHgwY8eO3SfMEl155ZWUlpYyYcKEZoMkLFk3J3VZWZlrwqDscsOCG1i8aTHPf+n5sEuRdrJs2TJGjRoVdhnShGR/HzNb5O5lrdmOmpgk7aKxKCOL1bwkkm0UEJJW1bXVrNq2ihF91Lwkkm0UEJJWKypXsNf3MqqPmiNEso0CQtIqEosA6AhCJAspICStIrEI3Qu7M6B7289YEZFwKCAkrSKxCCP6jNDcASJZSAEhaVO3t46PKj7SBXKSkaZOnUprTpmfO3cu11xzTRoryjwKCEmbtdvXsrt2ty6QE3jvCfjJGJjdO/j93hMHZbfNDb0tLVNASNpEKoIO6lF9dQZTXnvvCfj9d6BqLeDB799/54BDYvXq1YwcOZKZM2dSWlrKl770JXbt2rXf0NtLlixh4sSJlJaWcuGFF1JRUdGwjUceeYRJkyYxZswY3n77bQDefvttJk2axPjx45k0aRLRaLRh/bVr1zJt2jRGjBjB7bfffkD1ZwMNtSFpE41F6dihI0f1OirsUiSd/ngjfPp+04+XvwN1n+27rGY3/O4aWPRQ8uccPhbOubPFXUejUebMmcPkyZO5/PLLGyYRShx6u7S0lP/6r//ilFNO4dZbb+X222/npz/9KQA7d+7kjTfeYMGCBVx++eUsXbqUkSNHsmDBAjp27MgLL7zAzTffzG9/+1sgCI+lS5fStWtXjj/+eKZPn05ZWasuTs4qCghJm2WxZRzV6ygKCwrDLkXC1DgcWlreCoMGDWLy5MkAfP3rX+eee+4BPh96u6qqisrKyoZxk2bOnMnFF1/c8PxLLgkGnZ4yZQrbtm2jsrKS7du3M3PmTJYvX46Z7TOY3plnnknfvn0B+Lu/+ztee+01BYRIW0RjUSb1nxR2GZJuLX3T/8mYePNSI70Gwd/PP6BdNz47rv5+4tDbrX3+97//fU499VSeeuopVq9ezdSpU1vcX65SH4SkxZbdW9iye4vOYBI4/VYo7LLvssIuwfIDtGbNGt58800AHnvsMU4++eR9Hu/VqxfFxcW8+uqrADz88MP7jML661//GoDXXnuNXr160atXL6qqqhgwILhuZ+7cufts7/nnnycWi7F7926efvrphqOXXKWAkLTQHBDSoPTLcN49wREDFvw+755g+QEaNWoUDz30EKWlpcRiMf7xH/9xv3UeeughrrvuOkpLS1myZAm33vp5MBUXFzNp0iSuvvpq5syZA8D111/PTTfdxOTJk/c7C+rkk0/mG9/4BuPGjeOiiy7K6eYl0HDfkiZz3p/DTxf/lNcveZ2enXqGXY60s0wY7nv16tWce+65LF26NNQ6MpGG+5aMFo1FGdB9gMJBJIspICQtlsWW6QI5SauSkhIdPaSZAkLa3a6aXXyy7RP1P4hkOQWEtLvllctxXEN8i2Q5BYS0u/ozmDRJkEh2U0BIu4vEIvTs1JPDux0edikicgAUENLuIrEII/uMzPmrTCW7abjvlikgpF3V7q3lo4qP1P8g+5i/cj5nzTuL0odKOWveWcxfeWBDbKQqV4b7dnf27t170PergJB2tWbbGj6r+0xnMEmD+SvnM/uN2WzYuQHH2bBzA7PfmH1AIdHUUN9ARg33XVJSws0338xJJ51EWVkZixcv5uyzz+aoo47ivvvuA2DHjh2cfvrpTJgwgbFjx/K73/2u4TWOGjWKb33rW0yYMIG1a9cyZ84chg8fztSpU7niiivSfkSjwfqkXUViwRwQugYif/zw7R82/N2TeW/ze+zZu2efZdV11dz6+q3M+2he0ueM7DOSG064odn9Jhvq+3vf+x6QWcN9Dxo0iDfffJNrr72Wyy67jNdff53q6mqOOeYYrr76aoqKinjqqafo2bMnW7ZsYeLEicyYMaPhNT744IP8/Oc/Z/369fzgBz9g8eLF9OjRg9NOO41jjz222X0fKB1BSLuKVEQo7FDI0N5Dwy5FMkTjcGhpeaoaD/VdHwjQ/HDfCxYsaFgv2XDfVVVVXHzxxYwZM4Zrr72WDz74oGH9+uG+u3Tp0jDcd0vqP+zHjh3LiSeeSI8ePejXrx9FRUVUVlbi7tx8882UlpZyxhlnsG7dOjZu3AjA4MGDmThxIhCE0ymnnEKfPn0oLCzcZ9jydNERhLSraCzK0b2PprCD5oDIFy190z9r3lls2Llhv+VHdDuCB6c92Ob9Njf0diYN9925c2cAOnTo0HC7/n5tbS2PPvoomzdvZtGiRRQWFlJSUkJ1dfV+ryOMcfN0BCHtxt0bzmASqTdrwiyKCor2WVZUUMSsCbMOaLstDfUNB2+479NPP51169a16XVUVVVx6KGHUlhYyMsvv8wnn3ySdL0TTjiBv/zlL1RUVFBbW9vQ7JVOOoKQdrN592Zi1TGdwST7mD50OgB3L76bT3d+yuHdDmfWhFkNy9uqfqjvq666imHDhiUd6huC4b6vvvpqdu3axdChQ3nwwc+PWuqH+962bRsPPPAAEAz3PXPmTH784x9z2mmn7bOt+uG+V6xYwde+9jXKysrYu3cvK1asoE+fPm16HZdeeinnnXceZWVljBs3jpEjk3/BGjBgADfffDMnnngi/fv3Z/To0fTq1atN+0yZu2fVz3HHHeeSmf6y9i8+Zu4YX/jpwrBLkTT78MMPQ93/qlWr/Jhjjgm1hnrvv/++X3vttQdlX9u3b3d395qaGj/33HP9ySefTLpesr8PsNBb+XmrJiZpN/VDbOgMJsknY8aM4cc//vFB2dfs2bMZN24cY8aMYciQIVxwwQVp3Z+amKTdRGIRBvUYRPdO3cMuRXJcvg71/aMf/eig7k9HENJu1EGdXzzLZqPMF+35d0lrQJjZNDOLmtkKM7sxyeNHmtnLZvaumb1nZl9MZz2SPjtrdrJm+xo1L+WJoqIitm7dqpDIMO7O1q1bKSoqannlFKSticnMCoB7gTOBcuAdM3vG3T9MWO0W4Al3/4WZjQaeBUrSVZOkz0cVHwHoCCJPDBw4kPLycjZv3hx2KdJIUVERAwcObJdtpbMP4gRghbuvBDCzx4HzgcSAcKB+0uJewPo01iNp1DDEhk5xzQuFhYUMGTIk7DIkzdLZxDQAWJtwvzy+LNFs4OtmVk5w9PDtZBsysyvNbKGZLdQ3lswUjUUp7lzMYV0PC7sUEWkn6QyIZNegN26wvASY6+4DgS8CD5vZfjW5+/3uXubuZf369UtDqXKgIrEII/qM0BwQIjkknQFRDgxKuD+Q/ZuQvgk8AeDubwJFwCFprEnSoGZvDcsrlqv/QSTHpDMg3gGGmdkQM+sEfBV4ptE6a4DTAcxsFEFAqA0py6yuWs2evXvU/yCSY9IWEO5eC1wD/BlYRnC20gdmdoeZzYiv9l3gCjP7G/AYcJnrvLmsU99BPbJYRxAiuSStV1K7+7MEnc+Jy25NuP0hMDmdNUj6RWNROhd0pqRXSdiliEg70pXUcsAiFRGG9R5Gxw4auUUklygg5IB4fA4I9T+I5B4FhByQjbs2UvVZlc5gEslBCgg5IA0d1AoIkZyjgJADEolFMIxhxcPCLkVE2pkCQg5INBblyJ5H0q0wtUniRSR7KCDkgGgOCJHcpYCQNtu+ZzvlO8oVECI5SgEhbaY5qEVymwJC2ixaEQSEjiBEcpMCQtosEovQp6gPh3TRALwiuUgBIW0WjUUZ2Wek5oAQyVEKCGmTmroalldqDgiRXKaAkDZZWbWS2r21CgiRHKaAkDapH2JDg/SJ5C4FhLRJJBahqKCIwT0Gh12KiKSJAkLaJFoRZXjxcAo6FIRdioikiQJCWq1+Dgj1P4jkNgWEtNr6nevZvme7+h9EcpwCQlpNc0CI5AcFhLRaNBalg3XQHBAiOU4BIa0WiUUY3HMwXTp2CbsUEUkjBYS0WjQWZWSxmpdEcp0CQlql6rMq1u9cz8i+CgiRXKeAkFapnwNCRxAiuU8BIa1SfwbT8D7DQ65ERNJNASGtEq2I0q9LP80BIZIHFBDSKpFYRBfIieQJBYSkbE/dHlZWrtQFciJ5QgEhKVtRuYJa1xwQIvlCASEpaziDSQEhkhcUEJKySCxCl45dGNRjUNiliMhBoICQlEViEUYUj6CD6Z+NSD7omMpKZtYPuAIoSXyOu1+enrIk0+z1vUQropw79NywSxGRgySlgAB+B7wKvADUpa8cyVTrdqxjZ81ORvUZFXYpInKQpBoQXd39htZu3MymAXcDBcB/u/udSdb5MjAbcOBv7v611u5H0k9zQIjkn1QD4g9m9kV3fzbVDZtZAXAvcCZQDrxjZs+4+4cJ6wwDbgImu3uFmR3aitrlIIrEIhRYAUf1PirsUkTkIEm1t3EWQUjsNrNtZrbdzLa18JwTgBXuvtLd9wCPA+c3WucK4F53rwBw902tKV4OnmgsypBeQyjqWBR2KSJykKQUEO7ew907uHsXd+8Zv9+zhacNANYm3C+PL0s0HBhuZq+b2VvxJqn9mNmVZrbQzBZu3rw5lZKlnWmIDZH8k2oTE2ZWDAwDGr5CuvuC5p6SZJkn2f8wYCowEHjVzMa4e+U+T3K/H7gfoKysrPE2JM0qqivYuGujhvgWyTOpnub6DwTNTAOBJcBE4E3gtGaeVg4kXlE1EFifZJ233L0GWGVmUYLAeCel6uWgaOig1iRBInmlNX0QxwOfuPupwHigpbaed4BhZjbEzDoBXwWeabTO08CpAGZ2CEGT08oUa5KDpH6IjRHFamISySepBkS1u1cDmFlnd48AzX5auHstcA3wZ2AZ8IS7f2Bmd5jZjPhqfwa2mtmHwMvAde6+tS0vRNInUhHhsK6HUVxUHHYpInIQpdoHUW5mvQm+8T9vZhXs31y0n/hpsc82WnZrwm0H/iX+IxkqGovq+geRPJRSQLj7hfGbs83sZaAX8Ke0VZUt3nsCXrwDqsqh10A4/VYo/XLYVbWr6tpqVlWt4rQjm+tuEpFc1GxAmFlPd99mZn0SFr8f/90diKWtskz33hPw++9Aze7gftXa4D7kVEh8XPkxdV6nIwiRPNTSEcT/AucCiwhOUU08ddWBoWmqK/O9eMfn4VCvZnewPIcCYllsGaAhNkTyUbMB4e7nxn8POTjlZJGq8tYtz1KRWITuhd0Z0L3xNY4ikutaamKa0Nzj7r64fcvJIr0GBs1KyZbnkGgsyvDi4ZoDQiQPtdTE9P/iv4uAMuBvBM1MpcBfgZPTV1qGO/3WffsgAAq7BMtzRP0cEBcefWHLK4tIzmn2a6G7nxq/MO4TYIK7l7n7cQQXyq04GAVmrNIvw3n3QK9BgAW/z7snp/of1m5fy+7a3ep/EMlTqV4HMdLd689ewt2Xmtm4NNWUPUq/nFOB0Fh9B7UG6RPJT6kGxDIz+2/gEYKzl75OcHW05LBoLEpH68jRvY8OuxQRCUGqAfH3wD8SjMkEsAD4RVoqkowRiUUY2nsonQo6hV2KiIQg1Supq83sPuBZd4+muSbJENFYlJP6nxR2GSISkpTOXYwPrreE+PAaZjbOzBqPzCo5ZMvuLWzevVkjuIrksVRPbr+NYArRSgB3XwKUpKkmyQAfxT4CdAW1SD5LNSBq3b0qrZVIRtEZTCKSaif1UjP7GlBgZsOA7wBvpK8sCVs0FqV/t/706twr7FJEJCSpHkF8GzgG+IxgAL8qPj+jSXJQpCKioweRPJdqQIyO/3QkGHbjfDRvdM7aVbOL1VWr1f8gkudSbWJ6FPgesBTYm75yJBOsqFyB4zqCEMlzqQbEZnf/fVorkYwRiUUAncEkku9SDYjb4kNtvEjQDwGAuz+ZlqokVJFYhB6detC/W/+wSxGRELVmqI2RQCGfNzE5oIDIQdFYlJF9RmJmLa8sIjkr1YA41t3HprUSyQh1e+v4qOIjvjT8S2GXIiIhS/UsprfMbHRaK5GM8Mn2T6iuq1b/g4ikfARxMjDTzFYR9EEY4O5emrbKJBTRWDAWowJCRFINiGlprUIyxrLYMjp26MjQXkPDLkVEQpbqcN+fpLsQyQzRWJRhvYdRWFAYdikiErJU+yAkD7g7kZiG2BCRgAJCGmzZvYVYdUz9DyICKCAkQf0V1JokSERAASEJGgJCTUwiggJCEkRiEQZ2H0iPTj3CLkVEMoACQhpEK6LqfxCRBgoIAWBnzU7WbFuj5iURaaCAEACWVyzHcR1BiEgDBYQAmgNCRPangBAgCIhenXtxWNfDwi5FRDJEWgPCzKaZWdTMVpjZjc2s9yUzczMrS2c90rRILKI5IERkH2kLCDMrAO4FzgFGA5ckGzLczHoA3wH+mq5apHm1e2uDGzasAAAMFElEQVRZXrGckcVqXhKRz6XzCOIEYIW7r3T3PcDjwPlJ1vsB8J9AdRprkWasrlrNnr17dAaTiOwjnQExAFibcL88vqyBmY0HBrn7H5rbkJldaWYLzWzh5s2b27/SPBepUAe1iOwvnQGRrDHbGx406wD8BPhuSxty9/vdvczdy/r169eOJQoEQ3x36tCJkl4lYZciIhkknQFRDgxKuD8QWJ9wvwcwBnjFzFYDE4Fn1FF98C2LLePo4qMp7KA5IETkc+kMiHeAYWY2xMw6AV8Fnql/0N2r3P0Qdy9x9xLgLWCGuy9MY03SiLsTjUUZ1WdU2KWISIZJW0C4ey1wDfBnYBnwhLt/YGZ3mNmMdO1XWmfjro1UflapDmoR2U+qc1K3ibs/CzzbaNmtTaw7NZ21SHLRWBRQB7WI7E9XUue5+iE2hhcPD7kSEck0Cog8F4lFOLLHkXQr7BZ2KSKSYRQQeS4Si6j/QUSSUkDkse17tlO+o1xnMIlIUgqIPPZRxUeA5qAWkeQUEHlMc0CISHMUEHksGovSp6gP/bpo+BIR2Z8CIo9FYhFGFI/QHBAikpQCIk/V1NWwonIFI/uqeUlEklNA5KmVVSup2VujSYJEpEkKiDwVrdAQGyLSPAVEnorEIhQVFDG45+CwSxGRDKWAyFORWIRhxcMo6FAQdikikqEUEHnI3TXEhoi0SAGRhzbs3MD2Pds1xEYGmr9yPmfNO4vSh0o5a95ZzF85P+ySJI+ldT4IyUz1V1DrCCKzzF85n9lvzKa6rhoIgnz2G7MBmD50eoiVSb7SEUQeisaiGMaw3sPCLkUS3L347oZwqFddV83di+8OqSLJdwqIPBSJRRjcczBdC7uGXYok+HTnp61aLpJuCog8FIlFdP1DBjq82+GtWi6SbgqIPFP1WRXrd65X/0MGmjVhFkUFRfssKyooYtaEWSFVJPlOndR5pn4OCJ3BlHnqO6LvXnw3n+78lMO7Hc6sCbPUQS2hUUDkGZ3BlNmmD52uQJCMoSamPBOJRTikyyEc0uWQsEsRkQyngMgzuoJaRFKlgMgje+r2sLJypYb4FpGUKCDyyMeVH1PrtZokSERSooDII/Ud1DqCEJFUKCDySLQiSpeOXRjUY1DYpYhIFlBA5JFILMLw4uGaA0JEUqKAyBN7fS/RWFRDbIhIyhQQeWLdjnXsqNmhU1xFJGUKiDwRjUUBDbEhIqlTQOSJSCxCB+vA0b2PDrsUEckSCog8EY1FGdJzCEUdi1peWUQEBUTeWBZbpv4HEWmVtAaEmU0zs6iZrTCzG5M8/i9m9qGZvWdmL5rZ4HTWk68qqivYuGujzmASkVZJW0CYWQFwL3AOMBq4xMxGN1rtXaDM3UuBecB/pquefBatCDqodQQhIq2RziOIE4AV7r7S3fcAjwPnJ67g7i+7+6743beAgWmsJ2/Vn8GkIwgRaY10BsQAYG3C/fL4sqZ8E/hjsgfM7EozW2hmCzdv3tyOJeaHSCzCoV0PpU9Rn7BLEZEsks6AsCTLPOmKZl8HyoC7kj3u7ve7e5m7l/Xr168dS8wPkVhERw8i0mrpDIhyIHFUuIHA+sYrmdkZwP8FZrj7Z2msJy9V11azqmoVI4rV/yAirZPOgHgHGGZmQ8ysE/BV4JnEFcxsPPBLgnDYlMZa8tbHlR9T53U6ghCRVktbQLh7LXAN8GdgGfCEu39gZneY2Yz4ancB3YHfmNkSM3umic1JGzXMAaGAEJFW6pjOjbv7s8CzjZbdmnD7jHTuX4KA6FbYjYE9dIKYiLSOrqTOcdGKKCOKR9DB9KcWkdbRp0YOq58DQhfIiUhbKCBy2Nrta9lVu0v9DyLSJgqIHFbfQa0jCBFpCwVEDovGonS0jpoDQkTaRAGRwyKxCEN6D6FzQeewSxGRLKSAyGHRWJSRxep/EJG2UUDkqK27t7Jp9yb1P4hImykgcpSG+BaRA6WAyFGRCg2xISIHRgGRoyKxCEd0O4JenXuFXYqIZCkFRI65/aWHKZ0zhWdX/pEN27dw+0sPh12SiGQpBUQOuf2lh/nNJz/BO1ZgBnSo4Tef/EQhISJtooDIIb9d9SusQ80+y6xDDb9d9auQKhKRbKaAyCF7CypatVxEpDkKiBzSoa64VctFRJqjgMghFw25At9buM8y31vIRUOuCKkiEclmCogccttp3+DiwdditcW4g9UWc/Hga7nttG+EXZqIZCFz97BraJWysjJfuHBh2GWIiGQVM1vk7mWteY6OIEREJCkFhIiIJKWAEBGRpBQQIiKSlAJCRESSyrqzmMxsOxANu45GDgG2hF1EI5lYE2RmXaopNaopdZlY1wh379GaJ3RMVyVpFG3tqVrpZmYLVVNqMrEu1ZQa1ZS6TKzLzFp9fYCamEREJCkFhIiIJJWNAXF/2AUkoZpSl4l1qabUqKbUZWJdra4p6zqpRUTk4MjGIwgRETkIFBAiIpJU1gSEmT1gZpvMbGnYtdQzs0Fm9rKZLTOzD8xsVgbUVGRmb5vZ3+I13R52TfXMrMDM3jWzP4RdC4CZrTaz981sSVtOAUwXM+ttZvPMLBL/t3VSyPWMiL9H9T/bzOyfw6wpXte18X/jS83sMTMryoCaZsXr+SDM9yjZ56WZ9TGz581sefx3izOJZU1AAHOBaWEX0Ugt8F13HwVMBP7JzEaHXNNnwGnufiwwDphmZhNDrqneLGBZ2EU0cqq7j8uwc9bvBv7k7iOBYwn5PXP3aPw9GgccB+wCngqzJjMbAHwHKHP3MUAB8NWQaxoDXAGcQPB3O9fMhoVUzlz2/7y8EXjR3YcBL8bvNytrAsLdFwCxsOtI5O4b3H1x/PZ2gv/IA0Kuyd19R/xuYfwn9DMRzGwgMB3477BryWRm1hOYAswBcPc97l4ZblX7OB342N0/CbsQggt9u5hZR6ArsD7kekYBb7n7LnevBf4CXBhGIU18Xp4PPBS//RBwQUvbyZqAyHRmVgKMB/4abiUNTTlLgE3A8+4eek3AT4Hrgb1hF5LAgefMbJGZXRl2MXFDgc3Ag/HmuP82s25hF5Xgq8BjYRfh7uuAHwFrgA1Albs/F25VLAWmmFlfM+sKfBEYFHJNiQ5z9w0QfLkFDm3pCQqIdmBm3YHfAv/s7tvCrsfd6+LNAQOBE+KHvqExs3OBTe6+KMw6kpjs7hOAcwiaB6eEXRDBt+IJwC/cfTywkxSaAg4GM+sEzAB+kwG1FBN8Ix4C9Ae6mdnXw6zJ3ZcBPwSeB/4E/I2gGTprKSAOkJkVEoTDo+7+ZNj1JIo3TbxC+H03k4EZZrYaeBw4zcweCbckcPf18d+bCNrUTwi3IgDKgfKEo755BIGRCc4BFrv7xrALAc4AVrn7ZnevAZ4EJoVcE+4+x90nuPsUgiae5WHXlGCjmR0BEP+9qaUnKCAOgJkZQVvxMnf/cdj1AJhZPzPrHb/dheA/UiTMmtz9Jncf6O4lBE0UL7l7qN/2zKybmfWovw2cRdBEECp3/xRYa2Yj4otOBz4MsaREl5ABzUtxa4CJZtY1/v/wdDLgBAgzOzT++0jg78ic9wvgGWBm/PZM4HctPSFrRnM1s8eAqcAhZlYO3Obuc8KtisnAN4D3423+ADe7+7Mh1nQE8JCZFRB8AXjC3TPitNIMcxjwVPDZQkfgf939T+GW1ODbwKPxJp2VwN+HXA/xNvUzgavCrgXA3f9qZvOAxQTNOO+SGcNb/NbM+gI1wD+5e0UYRST7vATuBJ4ws28SBOzFLW5HQ22IiEgyamISEZGkFBAiIpKUAkJERJJSQIiISFIKCBERSUoBIdJKZlaSbFTh+LAYYQ/WKNJusuY6CJFM5+7/EHYNIu1JRxAibdPRzB4ys/ficzd0NbNXzKwMwMx2mNm/xefleMvMDgu7YJHWUkCItM0I4H53LwW2Ad9q9Hg3gqGfjwUWEMwTIJJVFBAibbPW3V+P334EOLnR43uA+iFOFgElB6kukXajgBBpm8Zj1DS+X+Ofj2NTh/r7JAspIETa5siEuaIvAV4LsxiRdFBAiLTNMmCmmb0H9AF+EXI9Iu1Oo7mKiEhSOoIQEZGkFBAiIpKUAkJERJJSQIiISFIKCBERSUoBISIiSSkgREQkqf8P5d/ojjb1wREAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f150bca4978>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Figure 3: Number of pairs per bin\n",
      "10 10\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAGQBJREFUeJzt3Xm0ZWV9p/HnKwgogoAUiBRFQSCIdjNZKgmYiKhBwUBWlBbFlAiiphVsR6Lt2BrLZVqNaWOaxWCpqBAcIEIMBEGcgAAigmiLWMgoJbOgQMGv/9i7wuFS995zq2qfc6v281nrrHvOPnv4neGe737fPaWqkCT112PGXYAkabwMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQCT5bJIPjWnZSXJiktuTXDSC5T03yfWTPPecJD/ruoaZSvJPSd6zmub12yTbr455rQ5JzktyxJDjTvrZtc+vtvepb9YddwF6tCRLgMcB21fVPe2wI4BDq+q5YyytC3sDLwDmLn+t41JV3wF2GmcNK1JVr1+N83rC6prXbLM636e+sUUwe60LHD3uImYqyToznGRbYMmwIZDElZcBvh9aHQyC2etjwNuSbDLxiSTzk9Tgj8BgEzvJq5N8L8knktyR5Jokf9wOvy7JLUkWTpjt5knOTnJ3km8n2XZg3k9tn7styc+SHDzw3GeTfCbJmUnuAfZZQb1PSXJ6O/3VSV7bDj8cOA74o7bL4gMrmHbwtdwGvD/JHyT5VpJbk/wmyUmD71OSJUneluTyJHcmOTnJBit6k5McleQnSeZO7HqYbj5J3pHkpiQ3Jjmi/Ux2mGQ55yX5SJKL2nmdlmSzgef/OcnN7XPnJ3n6hPf4Q+395ya5Psk7k9wMnJhk8yTfaD/r25J8J8kK/7cHa2zn++kkZ7Sf+4VJ/mBF07Xj75nk++1yfpTkuQPPHZbkqnY+1yR53YRpD0xyWZK7kvwiyX4DT2/bfsZ3JzkryeaT1dDO613t574kySuneZ/e2n7fb0py2FTz7TODYPa6GDgPeNtKTv9s4HLgScAXgS8DzwR2AA4F/k+SwW6CVwL/C9gcuAw4CSDJhsDZ7Ty2AA4B/nHwhwp4BfBhYCPguyuo5UvA9cBTgJcCf5tk36o6Hng98IOqekJVvW+K13JNu/wPAwE+0s5vZ2Ab4P0TpjkY2A/YDtgFePXEmabpT3418KdVNVnf8wrn0/6QvQV4Ps17+qeTTD/or4DXtHUvAz418Ny/Aju2r/FS2vd/Ek8GNqNpTR0JvJXm/Z0DbAm8Cxj23DGHAB8ANgWupnl/HyXJ1sAZwIfaZb8N+EqSOe0otwAHABsDhwGfSLJHO+2zgM8Bbwc2Af4EWDIw+1e002wBrMfU3/kn03xHtwYWAscmmaw778nAE9txDwc+nWTTKebdWwbB7PZe4E0D/2wz8cuqOrGqHgROpvmx/GBV3VdVZwH30/yALXdGVZ1fVfcB76ZZS9+G5p97STuvZVV1KfAVmh/05U6rqu9V1UNV9fvBItp57A28s6p+X1WX0bQCXjWD13JjVf1Du/zfVdXVVXV2+1qWAh/n0T/En6qqG6vqNuBfgN0eWVY+DvwZsE87j8lMNp+DgROr6sqqupfmx3Q6n6+qK9pusPcAB6ftSquqE6rq7vb9fz+wa5InTjKfh4D3ta//d8ADwFbAtlX1QFV9p4Y/idhXq+qiqlpGEz67TTLeocCZVXVm+zmfTbOy8uK2/jOq6hfV+DZwFvCcdtrDgRPaz+yhqrqhqn46MO8Tq+r/ta/llClqWO497Wv/Nk04HTzJeA/QfOcfqKozgd8yC7cBzQYGwSxWVVcA3wCOWYnJfz1w/3ft/CYOG2wRXDew3N8Ct9GsuW4LPLvtDrgjyR00rYcnr2jaFXgKcFtV3T0w7FqatbRhPWL+SbZI8uUkNyS5C/gCzVrioJsH7t/LI1/rJjRr0h+pqjunWfZk83nKhLqmeg9WNM61wGNpuuTWSbKo7TK5i4fXlifrIlk6IXA/RrM2f1bbLTOT78tU79OgbYGXTfge7E0TQCR5UZIL2q6pO2gCYnn92wC/WA01ANw+YXvStTSfxYrc2gbcsPPuLYNg9nsf8Foe+cO5/B/h8QPDBn+YV8Y2y++0XUabATfS/Hh9u6o2Gbg9oareMDDtVGufNwKbJdloYNg84IYZ1DZx/h9ph+1SVRvTrK1mBvO7naalc2KSvWYw3aCbgLkDj7eZbMRJxplHs8b6G5qukQNpupmeCMxvx5nsNT3i/WhbEm+tqu2BlwBvSbLvEPXMxHU0LZrB78GGVbUoyfo0rcS/A7asqk2AMwfqvw6YdNvDDG3adlcuN4/mO6ZVYBDMclV1NU3XzlEDw5bS/JAe2q5NvoZV/0d7cZK9k6xHs63gwqq6jqZF8odJXpXkse3tmUl2HrL+64DvAx9JskGSXWi6CqbqA5/ORjTN/Dvavuu3z3QGVXUeTcvma0mevRI1nAIclmTnJI+n6cabzqFJntaO/0Hg1LbrbiPgPuBWmnD/25kUkuSAJDskCXAX8GB7W52+ALwkyZ+137kN2g2yc2n69dcHlgLLkrwIeOHAtMfTvFf7JnlMkq2TPHUVavlAkvWSPIcm0P95FeYlDII1xQeBDScMey3ND+CtwNNpfmxXxRdpWh+3Ac+g+ZGk7dJ5IfBymjWvm4GP0vzjD+sQmrXcG4Gv0fRvn70KtX4A2AO4k6aP+KsrM5O2hsOA05M8Y4bT/ivNxt5zabplftA+dd8Uk30e+CzNe7gBD4f752i6OG4AfgJcMJNaaDYy/ztNOP4A+Mc26FabNtAPpNkQvZRmLf/twGPa78hRNOF4O00L5/SBaS+i3YBM85l9m6araWXc3C7jRpqViddP2N6glRAvTCOturaFdAWw/oR+6eXPnwd8oaqOG3Vt0nRsEUgrKclftF0Um9K0kv5lRSEgzXYGgbTyXkfTTfILmj75N0w9ujQ72TUkST1ni0CSem6NOGHV5ptvXvPnzx93GZK0Rrnkkkt+U1XTnplgjQiC+fPnc/HFF4+7DElaoyS5dpjx7BqSpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknlsjjizWzM0/5ozOl7Fk0f6dL0NS92wRSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HOdnoY6yRLgbuBBYFlVLUiyGXAyMB9YAhxcVbd3WYckaXKjaBHsU1W7VdWC9vExwDlVtSNwTvtYkjQm4+gaOhBY3N5fDBw0hhokSa2ug6CAs5JckuTIdtiWVXUTQPt3i45rkCRNoetLVe5VVTcm2QI4O8lPh52wDY4jAebNm9dVfZLUe522CKrqxvbvLcDXgGcBv06yFUD795ZJpj22qhZU1YI5c+Z0WaYk9VpnQZBkwyQbLb8PvBC4AjgdWNiOthA4rasaJEnT67JraEvga0mWL+eLVfXNJP8BnJLkcOBXwMs6rEFjMP+YMzpfxpJF+3e+DKkvOguCqroG2HUFw28F9u1quZKkmfHIYknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqua6vUDZ2nhJZkqZmi0CSes4gkKSeMwgkqefW+m0E4+T2CUlrAlsEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST3XeRAkWSfJD5N8o328XZILk/w8yclJ1uu6BknS5EbRIjgauGrg8UeBT1TVjsDtwOEjqEGSNIlOgyDJXGB/4Lj2cYDnAae2oywGDuqyBknS1LpuEXwSeAfwUPv4ScAdVbWsfXw9sPWKJkxyZJKLk1y8dOnSjsuUpP7qLAiSHADcUlWXDA5ewai1oumr6tiqWlBVC+bMmdNJjZKkbi9Msxfw50leDGwAbEzTQtgkybptq2AucGOHNUiSptFZi6Cq/qaq5lbVfODlwLeq6pXAucBL29EWAqd1VYMkaXrjOI7gncBbklxNs83g+DHUIElqTRsESY5OsnEaxye5NMkLZ7KQqjqvqg5o719TVc+qqh2q6mVVdd/KFi9JWnXDtAheU1V3AS8E5gCHAYs6rUqSNDLDBMHyPX1eDJxYVT9ixXv/SJLWQMMEwSVJzqIJgn9LshEPHxcgSVrDTbn7aHsk8HtpuoSuqap7kzyJpntIkrQWmDIIqqqSfL2qnjEw7Fbg1s4rkySNxDBdQxckeWbnlUiSxmKYI4v3AV6X5FrgHpoNxVVVu3RamSRpJIYJghd1XoUkaWwmDYIkG7fHD9w9wnokSSM2VYvgi8ABwCU0ZwgdPHaggO07rEuSNCKTBsHAKSG2G105kqRRG+o01Ek2BXakOZ00AFV1fldFSZJGZ9ogSHIEzXWH5wKXAXsCP6C55KQkaQ03zHEERwPPBK6tqn2A3QGvHSlJa4lhguD3VfV7gCTrV9VPgZ26LUuSNCrDbCO4PskmwNeBs5PcjpeXlKS1xrRBUFV/0d59f5JzgScC3+y0KknSyAy719AewN40xw98r6ru77QqSdLIDHOpyvcCi2muL7w5cGKS/9l1YZKk0RimRXAIsPvABuNFwKXAh7osTJI0GsPsNbSEgQPJgPWBX3RSjSRp5IZpEdwHXJnkbJptBC8AvpvkUwBVdVSH9UmSOjZMEHytvS13XjelSJLGYZjdRxePohBJ0ngMs41AkrQWMwgkqecmDYIkn2//Hj26ciRJozZVi+AZSbYFXpNk0ySbDd5GVaAkqVtTbSz+J5pzCm1Pc7lKL1UpSWuhSVsEVfWpqtoZOKGqtq+q7QZu04ZAkg2SXJTkR0muTPKBdvh2SS5M8vMkJydZbzW+HknSDE27sbiq3pBk1yRvbG+7DDnv+4DnVdWuwG7Afkn2BD4KfKKqdgRuBw5f2eIlSatumJPOHQWcBGzR3k5K8qbppqvGb9uHj21vRXOJy1Pb4YuBg1aibknSajLMkcVHAM+uqnsAknyU5prF/zDdhEnWodm+sAPwaZpzFN1RVcvaUa4Htp5k2iOBIwHmzZs3RJmSpJUxzHEEAR4cePwgj9xwPKmqerCqdqO58P2zgJ1XNNok0x5bVQuqasGcOXOGWZwkaSUM0yI4EbgwyfLzDR0EHD+ThVTVHUnOA/YENkmybtsqmIuXvZSksRpmY/HHgcOA22g27h5WVZ+cbrokc9prHZPkccDzgauAc4GXtqMtBE5budIlSavDUJeqrKpLaS5GMxNbAYvb7QSPAU6pqm8k+Qnw5SQfAn7IDFsXkqTVa6ggWBlVdTmw+wqGX0OzvUCSNAt40jlJ6rkpgyDJOkn+fVTFSJJGb8ogqKoHgXuTPHFE9UiSRmyYbQS/B37cXrP4nuUDvVaxJK0dhgmCM9qbJGktNNQ1i9vjAOZV1c9GUJMkaYSGOencS4DLaK5NQJLdkpzedWGSpNEYZvfR99Ps938HQFVdBmzXYU2SpBEaJgiWVdWdE4at8ERxkqQ1zzAbi69I8gpgnSQ7AkcB3++2LEnSqAzTIngT8HSaK459CbgLeHOXRUmSRmeYvYbuBd7dXpCmquru7suSJI3KMHsNPTPJj4HLaQ4s+1GSZ3RfmiRpFIbZRnA88NdV9R2AJHvTXKxm2IvYS5JmsWG2Edy9PAQAquq7gN1DkrSWmLRFkGSP9u5FSf4vzYbiAv4bcF73pUmSRmGqrqH/PeHx+wbuexyBJK0lJg2CqtpnlIVIksZj2o3F7QXo/wqYPzi+p6GWpLXDMHsNnQlcAPwYeKjbciRJozZMEGxQVW/pvBJJ0lgMs/vo55O8NslWSTZbfuu8MknSSAzTIrgf+Bjwbh7eW6iA7bsqSpI0OsMEwVuAHarqN10XI0kavWGC4Erg3q4LkVaH+cd0f3ntJYv273wZ0igNEwQPApclOZfmVNSAu49K0tpimCD4enuTJK2FhrkeweJRFCJJGo9hjiz+JSs4t1BVTbnXUJJtgM8BT6Y5EO3Yqvr7dtfTk2mOVF4CHFxVt8+4cmmWcfuE1lTDdA0tGLi/AfAyYJjjCJYBb62qS5NsBFyS5Gzg1cA5VbUoyTHAMcA7Z1a2JGl1mfaAsqq6deB2Q1V9EnjeENPdVFWXtvfvBq4CtgYOBJZ3Ny0GDlrp6iVJq2yYrqE9Bh4+hqaFsNFMFpJkPrA7cCGwZVXdBE1YJNlikmmOBI4EmDdv3kwWJ0magWG6hgavS7CMtl9/2AUkeQLwFeDNVXVXkqGmq6pjgWMBFixY4PUPJKkjw+w1tNLXJUjyWJoQOKmqvtoO/nWSrdrWwFbALSs7f0nSqhuma2h94C959PUIPjjNdKG58P1VVfXxgadOBxYCi9q/p824aknSajNM19BpwJ3AJQwcWTyEvYBXAT9Oclk77F00AXBKksOBX9HshSRJGpNhgmBuVe030xlX1XeByTYI7DvT+UmSujHM9Qi+n+S/dl6JJGkshmkR7A28uj3C+D6atfyqql06rUySNBLDBMGLOq9CkjQ2w+w+eu0oCpEkjccw2wgkSWsxg0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeq5zoIgyQlJbklyxcCwzZKcneTn7d9Nu1q+JGk4XbYIPgvsN2HYMcA5VbUjcE77WJI0Rp0FQVWdD9w2YfCBwOL2/mLgoK6WL0kazqi3EWxZVTcBtH+3mGzEJEcmuTjJxUuXLh1ZgZLUN7N2Y3FVHVtVC6pqwZw5c8ZdjiSttUYdBL9OshVA+/eWES9fkjTBqIPgdGBhe38hcNqIly9JmqDL3Ue/BPwA2CnJ9UkOBxYBL0jyc+AF7WNJ0hit29WMq+qQSZ7at6tlSpJmbtZuLJYkjYZBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9VxnZx+VNDrzjzmj82UsWbR/58vQeNgikKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqubEEQZL9kvwsydVJjhlHDZKkxsiDIMk6wKeBFwFPAw5J8rRR1yFJaoyjRfAs4Oqquqaq7ge+DBw4hjokSUCqarQLTF4K7FdVR7SPXwU8u6reOGG8I4Ej24c7AT8baaHjsznwm3EXMQa+7n7xdY/GtlU1Z7qRxnHN4qxg2KPSqKqOBY7tvpzZJcnFVbVg3HWMmq+7X3zds8s4uoauB7YZeDwXuHEMdUiSGE8Q/AewY5LtkqwHvBw4fQx1SJIYQ9dQVS1L8kbg34B1gBOq6spR1zGL9a47rOXr7hdf9ywy8o3FkqTZxSOLJannDAJJ6jmDYBZIsk2Sc5NcleTKJEePu6ZRSrJOkh8m+ca4axmVJJskOTXJT9vP/Y/GXdMoJPkf7Xf8iiRfSrLBuGvqSpITktyS5IqBYZslOTvJz9u/m46zxuUMgtlhGfDWqtoZ2BP47z077cbRwFXjLmLE/h74ZlU9FdiVHrz+JFsDRwELquq/0Ows8vLxVtWpzwL7TRh2DHBOVe0InNM+HjuDYBaoqpuq6tL2/t00Pwpbj7eq0UgyF9gfOG7ctYxKko2BPwGOB6iq+6vqjvFWNTLrAo9Lsi7weNbiY4iq6nzgtgmDDwQWt/cXAweNtKhJGASzTJL5wO7AheOtZGQ+CbwDeGjchYzQ9sBS4MS2S+y4JBuOu6iuVdUNwN8BvwJuAu6sqrPGW9XIbVlVN0GzAghsMeZ6AINgVknyBOArwJur6q5x19O1JAcAt1TVJeOuZcTWBfYAPlNVuwP3MEu6CLrU9ocfCGwHPAXYMMmh461KYBDMGkkeSxMCJ1XVV8ddz4jsBfx5kiU0Z6F9XpIvjLekkbgeuL6qlrf6TqUJhrXd84FfVtXSqnoA+Crwx2OuadR+nWQrgPbvLWOuBzAIZoUkoekvvqqqPj7uekalqv6mquZW1XyajYbfqqq1fg2xqm4GrkuyUztoX+AnYyxpVH4F7Jnk8e13fl96sJF8gtOBhe39hcBpY6zlP43j7KN6tL2AVwE/TnJZO+xdVXXmGGtSt94EnNSeb+sa4LAx19O5qrowyanApTR7yv2QWXrKhdUhyZeA5wKbJ7keeB+wCDglyeE0wfiy8VX4ME8xIUk9Z9eQJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgTSPJ/MEzSA4MP65nJwfUWsrjCKSVVFVHjLsGaXWwRSANZ90ki5Nc3l5H4PFJzkuyACDJb5N8OMmPklyQZMtxFywNyyCQhrMTcGxV7QLcBfz1hOc3BC6oql2B84HXjrg+aaUZBNJwrquq77X3vwDsPeH5+4HlV1i7BJg/orqkVWYQSMOZeC6WiY8fqIfP1/Igbn/TGsQgkIYzb+C6wocA3x1nMdLqZBBIw7kKWJjkcmAz4DNjrkdabTz7qCT1nC0CSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknvv/Ni5HeLf7/TkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f150bea6710>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "random.seed(42)\n",
    "np.random.seed(42)\n",
    "print(\"Creating all ranking pairs\")\n",
    "pairs = create_ranking_pairs() #list of all distinct ranking pairs\n",
    "print(\"\\t{} pairs created\\n\".format(len(pairs)))\n",
    "\n",
    "print(\"Computing deltaERR for each pair and assigning pairs to bins\")\n",
    "deltaERRs = divide_pairs_over_bins(pairs) #list of all pairs (tuples of two lists) divided over the 10 bins\n",
    "print(\"\\t{} non-empty bins created with in total {} pairs\\n\".format(\n",
    "    len([bin for bin in deltaERRs if len(bin) > 0]),\n",
    "    sum([len(p) for p in deltaERRs])))\n",
    "\n",
    "print(\"Performing EM to estimate position based click model parameters\")\n",
    "(alpha, gamma) = em()\n",
    "epsilon = 0.05\n",
    "alpha_em = alpha\n",
    "alpha = 1 - epsilon    # alpha is manually determined\n",
    "gamma = gamma[:3]      # only first 3 gammas are used\n",
    "print(\"\\talpha: {:.3f} (manually set after EM, result from EM: {:.3f})\".format(alpha, alpha_em))\n",
    "print(\"\\tgammas: {}\\n\".format(\", \".join([\"{:.3f}\".format(g) for g in gamma])))\n",
    "\n",
    "print(\"Setting random click model parameter\")\n",
    "theta = 0.5\n",
    "print(\"\\ttheta: {:.3f}\\n\".format(theta))\n",
    "\n",
    "k = 10000\n",
    "print(\"Starting online experiments with k={}\".format(k))\n",
    "\n",
    "# Experiment 1: team draft interleaving, position-based click model\n",
    "estimator = estimate_win_proportion\n",
    "interleaver = td_interleave\n",
    "clicker = produce_clicks\n",
    "print(\"\\nRunning experiment 1: team draft interleaving & position-based click model\")\n",
    "results_td_pb = run_interleaving_experiment(deltaERRs, \n",
    "                                    estimator, \n",
    "                                    interleaver, \n",
    "                                    clicker, \n",
    "                                    alpha, gamma, theta, k)\n",
    "\n",
    "# Experiment 2: team draft interleaving, random click model\n",
    "estimator = estimate_win_proportion\n",
    "interleaver = td_interleave\n",
    "clicker = produce_clicks_random\n",
    "print(\"\\nRunning experiment 2: team draft interleaving & random click model\")\n",
    "results_td_rn = run_interleaving_experiment(deltaERRs, \n",
    "                                    estimator, \n",
    "                                    interleaver, \n",
    "                                    clicker, \n",
    "                                    alpha, gamma, theta, k)\n",
    "\n",
    "# Experiment 3: probabilistic interleaving, position-based click model\n",
    "estimator = estimate_win_proportion\n",
    "interleaver = pr_interleave\n",
    "clicker = produce_clicks\n",
    "print(\"\\nRunning experiment 3: probabilistic interleaving & position-based click model\")\n",
    "results_pr_pb = run_interleaving_experiment(deltaERRs, \n",
    "                                    estimator, \n",
    "                                    interleaver, \n",
    "                                    clicker, \n",
    "                                    alpha, gamma, theta, k)\n",
    "\n",
    "# Experiment 4: probabilistic interleaving, random click model\n",
    "estimator = estimate_win_proportion\n",
    "interleaver = pr_interleave\n",
    "clicker = produce_clicks_random\n",
    "print(\"\\nRunning experiment 4: probabilistic interleaving & random click model\")\n",
    "results_pr_rn = run_interleaving_experiment(deltaERRs, \n",
    "                                    estimator, \n",
    "                                    interleaver, \n",
    "                                    clicker, \n",
    "                                    alpha, gamma, theta, k)\n",
    "\n",
    "# Experiment 5: marginalised probabilistic interleaving, position-based click model\n",
    "estimator = estimate_win_proportion_marginalised\n",
    "interleaver = pr_interleave\n",
    "clicker = produce_clicks\n",
    "print(\"\\nRunning experiment 5: marginalised probabilistic interleaving & position-based click model\")\n",
    "results_pr_pb_mar = run_interleaving_experiment(deltaERRs, \n",
    "                                    estimator, \n",
    "                                    interleaver, \n",
    "                                    clicker, \n",
    "                                    alpha, gamma, theta, k)\n",
    "\n",
    "# Experiment 6: marginalised probabilistic interleaving, random click model\n",
    "estimator = estimate_win_proportion_marginalised\n",
    "interleaver = pr_interleave\n",
    "clicker = produce_clicks_random\n",
    "print(\"\\nRunning experiment 6: marginalised probabilistic interleaving & random click model\")\n",
    "results_pr_rn_mar = run_interleaving_experiment(deltaERRs, \n",
    "                                    estimator, \n",
    "                                    interleaver, \n",
    "                                    clicker, \n",
    "                                    alpha, gamma, theta, k)\n",
    "\n",
    "print(\"Showing results\\n\\n\")\n",
    "\n",
    "print(\"Table 1: team draft interleaving & position-based click model\")\n",
    "print(results_table(results_td_pb))\n",
    "\n",
    "print(\"\\nTable 2: team draft interleaving & random click model\")\n",
    "print(results_table(results_td_rn))\n",
    "    \n",
    "print(\"\\nTable 3: probabilistic interleaving & position-based click model\")\n",
    "print(results_table(results_pr_pb))\n",
    "\n",
    "print(\"\\nTable 4: probabilistic interleaving & random click model\")\n",
    "print(results_table(results_pr_rn)) \n",
    "\n",
    "print(\"\\nTable 5: marginalised probabilistic interleaving & position-based click model\")\n",
    "print(results_table(results_pr_pb_mar))\n",
    "\n",
    "print(\"\\nTable 6: marginalised probabilistic interleaving & random click model\")\n",
    "print(results_table(results_pr_rn_mar))\n",
    "\n",
    "print(\"\\nFigure 1: Median impressions per bin, position-based click model\")\n",
    "plot_medians(results_td_pb, results_pr_pb, results_pr_pb_mar)\n",
    "      \n",
    "print(\"\\nFigure 2: Median impressions per bin, random click model\")\n",
    "plot_medians(results_td_rn, results_pr_rn, results_pr_rn_mar)\n",
    "\n",
    "print(\"\\nFigure 3: Number of pairs per bin\")\n",
    "plot_pairs_histogram(deltaERRs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 7: Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results\n",
    "We refer to the tables and figures in the previous output cell.\n",
    "\n",
    "## Analysis and Conclusion\n",
    "The main conclusion is that the offline evaluation can inform the design of the online comparison of two ranking algorithms in a meaningful way, given the assumptions made in the experiment. The strongest assumptions are that PBM accurately describes user click behaviour; that the same set of queries are used in the offline and online evaluation; and that there is no query variance. Figure 1 shows a correlation between the result of the offline evaluation ($\\Delta$ERR) and the median number of impressions needed to significantly compare two ranking algorithms. As is expected, an experiment with an algorithm pair with a larger $\\Delta$ERR requires less impressions to make a significant conclusion about the relative performance.\n",
    "\n",
    "### Click models\n",
    "Experiments using the random click model are used as sanity check. About half of the pairs in each bin result in an infinite N. This confirms the check, because systems E and P have equal chances of winning. \n",
    "\n",
    "With the results for the position-based click model, there are significant differences between the interleaving methods, as can be seen in Table 1, 3 and 5.\n",
    "Probabilistic marginalised interleaving seems to require less impressions to show a significant difference in the lower bins 1 and 2 (see Figure 1). With higher $\\Delta$ERR values the difference is small.\n",
    "\n",
    "We first analyze the results of the random click-model. In this model, there is a fixed probability of, in our case, 50% that a document is clicked. This probability is independent from the rank and the relevance of a document. It therefore makes sense that neither the $\\Delta$ERR, nor the interleaving method has any significant influence on the results. This explains why we consistently see approximately the same numbers in the tables.\n",
    "\n",
    "With the position-based click model, the click-probability of a document depends on both its rank and its relevance. In particular, the dependence on the relevance effectively means that a ranker will receive more clicks if it more often produces relevant results. If those relevant results are higher in the ranking, they receive even more clicks. In our case, the experimental algorithm E more often produces relevant, high-ranked results in the cases where the $\\Delta$ERR is big, and that means that differences in the algorithms are more easily detectable in a statistically significant way. Therefore, we could expect a priori that the results for N drop with increasing Delta-ERR. That’s also what we observe in the tables.\n",
    "\n",
    "### Team draft interleaving\n",
    "About team-draft interleaving: This follows the trend already described in the click-models section. The required number of impressions drops almost exponentially with an increasing ΔERR. Furthermore, we observe that in the first bin, roughly 4.4% of pairs lead to a value of infinity. This can potentially be explained by the insufficient sensitivity of team draft interleaving that moves the average win-proportion of E closer to 0.5.\n",
    "\n",
    "### Probabilistic interleaving\n",
    "Based on Figure 1, team draft interleaving seems to perform slightly better than probabilistic interleaving. This is against our expectations. Probabilistic interleaving was developed to be more sensitive which means that is it better at distinguishing between systems that are very close (and are not distinguished by team draft)[3]. The source of the differences could be further investigated. \n",
    "\n",
    "Probabilistic marginalised interleaving seems to perform better than probabilistic interleaving. This is what we expected. In probabilistic interleaving, luck plays a larger role in winning a simulation. Marginalised interleaving looks at the interleaved list and the click pattern, and then decides which system most probably contributed the clicked documents. Normal probabilistic interleaving decides based on the accidental assignment per simulation. Thus, by marginalising we increase the probability of the better system winning, and therefore need a smaller amount of impressions. This effect is larger for systems with smaller differences, as is confirmed by Figure 1. If the differences are already large, the worse system can already hardly win and the number of required impressions is generally small.\n",
    "\n",
    "Nevertheless, we find the following difference unexpected: Probabilistic marginalised interleaving produced more max-values of infinity than usual probabilistic interleaving: 6.7% and 4.5% in the first two bins, compared to no cases at all. While it can be that $\\Delta$ERR favors one ranker that is actually not better on the online-evaluation, the difference between the two versions of probabilistic interleaving requires an explanation. Further experiments are needed to explain this phenomenon.\n",
    "\n",
    "## Possible improvements\n",
    "\n",
    "### Different click-models\n",
    "\n",
    "A different click model could be implemented in order to potentially better simulate actual user behaviour. Examples of widely used click models are the Dependent Click Model (DCM), Click Chain Model (CCM), and a Dynamic Bayesian Network Model (DBN)[1].\n",
    "\n",
    "### Better parameter-estimation\n",
    "\n",
    "To find the parameters of the click model, EM is used to infer the examination (gamma) and the attractiveness (alpha) probabilities. Currently, we use a fixed constant epsilon to define how probable it is that a user clicks on a non-relevant document when he sees it, because it is difficult to infer the value for alpha for each possible document-query pair and we have no information about which documents are relevant and which are not. It could be tried to acquire click logs where multiple sessions occur with the same document-query pairs to compute meaningful alpha values for these document-query pairs. The alpha inferred for documents which are clicked a significant number of times can inform a sensible choice of alpha for relevant documents. This could be done by for instance calculating the average over all these document-query pairs with enough clicks. These click logs could be acquired by setting up an experiment where users have to choose the best documents given a query.\n",
    "\n",
    "### Closed-form solutions\n",
    "\n",
    "A closed form solution can be found, which gives a more accurate and computationally efficient method for finding the number of impressions. This has also been suggested in [2] and our own description of the experiment (see the appendix of part A of the programming assignment). With a closed-form solution the variance is reduced and we therefore get more precise values about the number of impressions to choose. This can only be approached in the current experiment by choosing a very high value for *k*, which has computational disadvantages compared to the closed-form solution.\n",
    "\n",
    "### Beyond the no-query-variability assumption\n",
    "\n",
    "Currently, the assumption is made that there is no query variability, which means that all queries would lead to the same document ranking. It is also assumed that the queries in the offline and online experiment are the same. The assumptions make the experiment design easier, but at the expense of less realistically modelling real search algorithms. Finding a way to take the query variability into account, as is also suggested by [2], would be an improvement for the experimental design."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sources\n",
    "-----------------------------\n",
    "[1] Chuklin, Aleksandr, Ilya Markov, and Maarten de Rijke. \"Click models for web search.\" Synthesis Lectures on Information Concepts, Retrieval, and Services 7.3 (2015): 1-115.\n",
    "\n",
    "[2] Azarbonyad, Hosein, and Evangelos Kanoulas. \"Power Analysis for Interleaving Experiments by means of Offline Evaluation.\" Proceedings of the 2016 ACM International Conference on the Theory of Information Retrieval. ACM, 2016.\n",
    "\n",
    "[3] Hofmann, Katja, Shimon Whiteson, and Maarten De Rijke. \"A probabilistic method for inferring preferences from clicks.\" Proceedings of the 20th ACM international conference on Information and knowledge management. ACM, 2011."
   ]
  }
 ],
 "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}