cuped-for-starters/cuped-for-starters.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e66262fc-dfb4-473f-8e83-aeffbbf00a59",
   "metadata": {},
   "source": [
    "# CUPED for Starters: Enhancing Controlled Experiments with Pre-Experiment Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71376a2d-df18-4230-8d1f-943612297363",
   "metadata": {},
   "source": [
    "### Intro"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "833834d9-5d04-4d64-947e-c42d1be85249",
   "metadata": {},
   "source": [
    "In this article, I will briefly explain the randomized controlled experiments and why modern companies use them to make data-driven decisions. Then, I will introduce you to a CUPED procedure that improves the sensitivity of these experiments. After that, I will show you why it works in theory and practice. We will simulate an ordinary A/B test and compare results with the CUPED-adjusted procedure. We will see empirically that CUPED allows the detection of more minor effects on the same sample size."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c908fbd8-6a08-47bd-9204-9bbd18127b22",
   "metadata": {},
   "source": [
    "### What are Controlled Experiments, and why do we need them"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b923b119-7c0d-4f69-898e-6f679b92d20a",
   "metadata": {},
   "source": [
    "Controlled experiments, a cornerstone of data science, are systematic and controlled approaches used to investigate the impact of specific changes, interventions, or treatments on a target system. In these experiments, researchers manipulate one or more variables while keeping other factors constant, creating a controlled environment to accurately observe and measure the effects. The goal is to establish a causal relationship between the manipulated variables and the observed outcomes. Controlled experiments provide valuable insights into cause-and-effect relationships, enabling data scientists, researchers, and decision-makers to make informed choices based on empirical evidence rather than mere correlations or assumptions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11ac045c-9ac7-4241-8643-c3a1f9ef10e9",
   "metadata": {},
   "source": [
    "Usually, people encounter controlled experiments in the form of A/B tests. In an A/B test, a single variable (such as a webpage design, machine learning algorithm, or pricing plan) is changed between two versions: A (the control) and B (the variant). The goal of an A/B test is to compare the performance of these two versions to determine which one performs better in a specific metric (e.g., click-through rates, revenue). A/B tests are commonly used in marketing and web optimization to make data-driven decisions about changes to user interfaces or marketing material."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "783f1443-2b33-4cde-b081-9e89d2a24b07",
   "metadata": {},
   "source": [
    "To perform an A/B test, you have to gather some data. The amount of data you need is directly connected to the variance of the metric you want to compare between two groups. The more the chosen metric fluctuates naturally, the more data you need to spot the difference between the two groups. The more data you need, the longer you wait until you get the results. And time is money! That is why Microsoft researchers developed a [CUPED technique](https://dl.acm.org/doi/10.1145/2433396.2433413). This statistical trick helps you to decrease the variance of the tested metric. With lower variance, you need a smaller sample size to spot the same difference, or you can spot a smaller one with the same sample size."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19dc73a6-9264-4315-a53b-0c24c3e4be62",
   "metadata": {},
   "source": [
    "### What is CUPED and how it works?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e46ad9e-a5fd-4c94-a942-f445fef123dc",
   "metadata": {},
   "source": [
    "CUPED stands for \"Controlled-Experiment using Pre-Experiment Data.\" The idea is simple: modify the metric of interest in a way that its expectation does not change while variance decreases. The exciting part is how to do so! "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67eb3cd8-febe-4339-b9eb-e0cb2c4920d0",
   "metadata": {},
   "source": [
    "Suppose you want to test if there is a difference in metric $Y$. In a usual A/B test, you would calculate mean value of the metric in control group $\\bar{Y}_{A}$ and then compare it with mean value of the metric in test group $\\bar{Y}_{B}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af6f250c-f129-491b-972f-05adb04c00f7",
   "metadata": {},
   "source": [
    "Researchers suggested this transformation:\n",
    "$$\n",
    "\\hat{Y}_{CUPED} = \\bar{Y} - \\theta(\\bar{X} - E(X))\n",
    "$$\n",
    "Where $X$ is a random variable correlated with $Y$ but not affected by the A/B test, and $\\theta$ is some arbitrary (for now) constant."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd670862-9890-482d-9b0f-b5ac6c9b4b33",
   "metadata": {},
   "source": [
    "A great candidate for $X$ is the same metric $Y$ but before the start of the test. Hence, there is \"pre-experiment\" in the name of the method."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf5539f1-84f6-4252-a0dd-8d8cb12586f5",
   "metadata": {},
   "source": [
    "But how does this help to decrease variance and not mess up the results? The theoretical explanation in the paper needed to be more specific for me. So, let us derive results by ourselves."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31a0cab8-f15e-4cfb-bfe1-b0c587b09563",
   "metadata": {},
   "source": [
    "First of all, $\\hat{Y}_{CUPED}$ is an unbiased estimator of $E(Y)$:\n",
    "\n",
    "![CUPED Proof Part 1](./input/CUPED_proof_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3595730a-dc07-4676-8095-f8704bd5e166",
   "metadata": {},
   "source": [
    "Secondly, let us derive the $Var(\\hat{Y}_{CUPED})$:\n",
    "\n",
    "![CUPED Proof Part 1](./input/CUPED_proof_2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d3d34a4-4f35-4cb2-9af2-2260be9d587e",
   "metadata": {},
   "source": [
    "Now, let us find a value of $\\theta$ that minimizes this variance:\n",
    "\n",
    "![CUPED Proof Part 1](./input/CUPED_proof_3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d459ccfe-2e68-4950-b718-8ea7ed47467c",
   "metadata": {},
   "source": [
    "With such $\\theta$ we will have $Var(\\hat{Y}_{CUPED})$:\n",
    "\n",
    "![CUPED Proof Part 1](./input/CUPED_proof_4.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a756405e-8a27-40c7-aa28-dabf08d5e9c8",
   "metadata": {},
   "source": [
    "Cool! Variance is $(1 - \\rho^2)$ times smaller that the variance of our initial estimate $Var(\\bar{Y})$!\n",
    "The higher the correlation between $X$ and $Y$, the higher the variance reduction power."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccdceb30-dc3a-4f14-aaa7-182cbfcfd2dc",
   "metadata": {},
   "source": [
    "### How to use CUPED in practice?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9b9675e-2a8f-4524-8fe3-913676da973d",
   "metadata": {},
   "source": [
    "Enough theory for today! Let us apply this technique to a real-world problem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5770f6e0-7728-480e-a95e-0bd27c281368",
   "metadata": {},
   "source": [
    "Suppose we are working for some food delivery service (like Uber Eats, for example). Each restaurant has a display of meals users can order. There are so many meals that we want to recommend something to the user to speed up the choosing process. The faster he places an order, the quicker he gets his food, and the faster he becomes happy!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fcb88b2-83ca-4978-ba8a-6d75cfeab2c4",
   "metadata": {},
   "source": [
    "We already have a recommendation engine A, but the ML team has introduced another recommendation engine B. Now we want to understand which engine is better. We need an A/B test! Let us randomly expose each user to either engine A or engine B. Now, we need a metric to compare between two groups. Money is the best metric as each business tries to maximize its profits. So, let us measure the commission revenue from users within two weeks."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53743eec-4fcf-4f4a-b80a-a44db64c411f",
   "metadata": {},
   "source": [
    "I have prepared a simple dataset that represents this A/B test. Each row represents a user.\n",
    "Here is an explanation of what each column means:\n",
    "1. `group` - is either `test` or `control`; if it is `test`, then the user is in the group with new recommendation engine\n",
    "2. `user_id` - a unique identifier of a user\n",
    "3. `commission_revenue_pre` - commission revenue company got from orders of this user in a pre-experiment period\n",
    "4. `commission_revenue_exp` - commission revenue company for from order of this user during an experiment period"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0edfe023-0a6f-45ff-9671-37e488174909",
   "metadata": {},
   "source": [
    "The only thing different from real life here is that I know exactly how a new engine impacts the metric. If a user is exposed to engine B, his `commission_revenue_exp` will increase by 2%."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "397d8206-7f87-4cc1-91c3-396dfa055cac",
   "metadata": {},
   "source": [
    "Let us do the usual A/B test comparing average `commission_revenue_exp` for two engines using a t-test (assume that all assumptions hold, and we can use it here). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3180bdff-dab9-4626-8d20-463125f91eff",
   "metadata": {},
   "source": [
    "First of all, we need to import some useful libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c3f1b4f4-bfe3-4e12-8acb-0bee37eb5df3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e79b8df8-4a04-4c3d-839c-1d0d3f3aaf04",
   "metadata": {},
   "source": [
    "Now, let's load the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e0f6c8f5-b2e8-462d-9799-b1702988f8a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>group</th>\n",
       "      <th>user_id</th>\n",
       "      <th>commission_revenue_pre</th>\n",
       "      <th>commission_revenue_exp</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>277</th>\n",
       "      <td>test</td>\n",
       "      <td>d8236b3c-63b8-4c5a-b68a-aec573a443d3</td>\n",
       "      <td>90.304491</td>\n",
       "      <td>117.935468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1850</th>\n",
       "      <td>control</td>\n",
       "      <td>381587d3-0372-4208-87c9-ca9479e57e2a</td>\n",
       "      <td>74.032403</td>\n",
       "      <td>80.669517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1875</th>\n",
       "      <td>test</td>\n",
       "      <td>ac74ff20-ed82-4c54-9dd9-28f5989a39e0</td>\n",
       "      <td>139.042422</td>\n",
       "      <td>136.864395</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>145</th>\n",
       "      <td>control</td>\n",
       "      <td>d11da65b-194c-47b0-91c8-1e3082c9e9eb</td>\n",
       "      <td>138.068400</td>\n",
       "      <td>140.146841</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1951</th>\n",
       "      <td>test</td>\n",
       "      <td>60a40fd8-3b20-4867-9e03-ddb65ecbf119</td>\n",
       "      <td>83.290920</td>\n",
       "      <td>80.130643</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1871</th>\n",
       "      <td>test</td>\n",
       "      <td>be5c7f99-351c-49d9-9cb2-5893e7de036d</td>\n",
       "      <td>112.421797</td>\n",
       "      <td>93.468268</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1535</th>\n",
       "      <td>control</td>\n",
       "      <td>2d2731e8-88c9-431f-a4e0-fa2cc8fbf009</td>\n",
       "      <td>60.578493</td>\n",
       "      <td>53.202416</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1030</th>\n",
       "      <td>test</td>\n",
       "      <td>a0fa887a-8890-492d-ba30-248b103bd178</td>\n",
       "      <td>125.308029</td>\n",
       "      <td>131.100739</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1954</th>\n",
       "      <td>test</td>\n",
       "      <td>5f40cdcc-a9d9-4088-ad44-8cbd74e3468e</td>\n",
       "      <td>74.458162</td>\n",
       "      <td>70.715691</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1082</th>\n",
       "      <td>test</td>\n",
       "      <td>17f61cd0-6b32-4d30-ab3a-b96d84d056e5</td>\n",
       "      <td>68.946865</td>\n",
       "      <td>66.562188</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        group                               user_id  commission_revenue_pre  \\\n",
       "277      test  d8236b3c-63b8-4c5a-b68a-aec573a443d3               90.304491   \n",
       "1850  control  381587d3-0372-4208-87c9-ca9479e57e2a               74.032403   \n",
       "1875     test  ac74ff20-ed82-4c54-9dd9-28f5989a39e0              139.042422   \n",
       "145   control  d11da65b-194c-47b0-91c8-1e3082c9e9eb              138.068400   \n",
       "1951     test  60a40fd8-3b20-4867-9e03-ddb65ecbf119               83.290920   \n",
       "1871     test  be5c7f99-351c-49d9-9cb2-5893e7de036d              112.421797   \n",
       "1535  control  2d2731e8-88c9-431f-a4e0-fa2cc8fbf009               60.578493   \n",
       "1030     test  a0fa887a-8890-492d-ba30-248b103bd178              125.308029   \n",
       "1954     test  5f40cdcc-a9d9-4088-ad44-8cbd74e3468e               74.458162   \n",
       "1082     test  17f61cd0-6b32-4d30-ab3a-b96d84d056e5               68.946865   \n",
       "\n",
       "      commission_revenue_exp  \n",
       "277               117.935468  \n",
       "1850               80.669517  \n",
       "1875              136.864395  \n",
       "145               140.146841  \n",
       "1951               80.130643  \n",
       "1871               93.468268  \n",
       "1535               53.202416  \n",
       "1030              131.100739  \n",
       "1954               70.715691  \n",
       "1082               66.562188  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exp_data = pd.read_csv('./input/commission_revenue.csv')\n",
    "\n",
    "exp_data.sample(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "36b75334-7fe0-43bf-b8b8-4f276979b9f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 2000 entries, 0 to 1999\n",
      "Data columns (total 4 columns):\n",
      " #   Column                  Non-Null Count  Dtype  \n",
      "---  ------                  --------------  -----  \n",
      " 0   group                   2000 non-null   object \n",
      " 1   user_id                 2000 non-null   object \n",
      " 2   commission_revenue_pre  1963 non-null   float64\n",
      " 3   commission_revenue_exp  2000 non-null   float64\n",
      "dtypes: float64(2), object(2)\n",
      "memory usage: 62.6+ KB\n"
     ]
    }
   ],
   "source": [
    "exp_data.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38131e11-db21-44ff-ba30-9ce71cee9e3d",
   "metadata": {},
   "source": [
    "Let's have a look at the average commission revenue during the experiment in two groups:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1cb034e6-4714-42f1-88e8-c3db4a5663f3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In a test group (new engine B), the metric is on average equal to 101.73\n",
      "In a control group (old engine A), the metric is on average equal to 99.24\n",
      "Relative increase is equal to 2.5%\n"
     ]
    }
   ],
   "source": [
    "commission_revenue_exp_test_avg = exp_data[exp_data['group'] == 'test']['commission_revenue_exp'].mean()\n",
    "commission_revenue_exp_control_avg = exp_data[exp_data['group'] == 'control']['commission_revenue_exp'].mean()\n",
    "\n",
    "print(f'In a test group (new engine B), the metric is on average equal to {commission_revenue_exp_test_avg:0.2f}')\n",
    "print(f'In a control group (old engine A), the metric is on average equal to {commission_revenue_exp_control_avg:0.2f}')\n",
    "print(f'Relative increase is equal to {100.0 * (commission_revenue_exp_test_avg / commission_revenue_exp_control_avg - 1):0.1f}%')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f865679f-8b88-4440-9bbb-b0b770548941",
   "metadata": {},
   "source": [
    "We already know that engine B increases metric by two percent because data was synthetically created. But in the real world, we need to run some statistical procedures to understand whether this increase is statistically significant. Maybe we were just lucky to get this result."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fba67a7-317b-4673-b3b5-f58a7595fa12",
   "metadata": {},
   "source": [
    "Using [Student's t-tets](https://en.wikipedia.org/wiki/Student%27s_t-test), we can test a hypothesis of zero difference in the averages of two groups. Let's use $alpha = 0.05$ for this test, meaning that we will test this hypothesis at a 5% significance level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "13b7b3d2-c192-4930-9e35-1ffb7870c02f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=1.8021000484396446, pvalue=0.0716803352566017)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.ttest_ind(\n",
    "    exp_data[exp_data['group'] == 'test']['commission_revenue_exp'], \n",
    "    exp_data[exp_data['group'] == 'control']['commission_revenue_exp']\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55dba0da-2e5b-43ef-95e4-a9ec8da32abb",
   "metadata": {},
   "source": [
    "We got a p-value near $0.07$, which is definitely greater than $0.05 / 2 = 0.025$. Then, we cannot reject the null hypothesis of equal averages. This means we did not see any statistically significant improvement in the target metric from using the new recommendation engine B in production."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cfd4bc5-cdea-4cf2-95fa-518230205d6e",
   "metadata": {},
   "source": [
    "You may ask how this happened, if we know there must be an increase of 2%. When we have an infinite number of observations, sample means are equal to population means (because of [central limit theorem](https://en.wikipedia.org/wiki/Central_limit_theorem)). Unfortunately, we have to operate with finite samples in our experiments. This introduces some variation in our estimates of the parameters of distributions. So, we always have some margin of error in our estimates. And the smaller the sample size, the bigger the margin or error. When the sample size is not big enough, we cannot say that the difference between two samples is not present because of random associated with finite sample size."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02eefd3d-fda3-4a3e-ad0a-3760804b9431",
   "metadata": {},
   "source": [
    "As discussed earlier, CUPED can help us improve sensitivity by decreasing the metric variance we compare between two groups. Ultimately, this means we can detect smaller parameter changes for the same sample size. Let's make a CUPED transformation using `commission_revenue_pre`. We can do so because `commission_revenue_pre` is highly likely to be correlated with `commission_revenue_exp` because this is the same metric of the same users but in different time periods. But `commission_revenue_pre` is not affected by the new engine B because engine B was introduced after the metric was calculated. Hence, `commission_revenue_pre` is a fantastic candidate for `X`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2266ec70-d3ab-4b75-9aa1-ec1c81ad80c2",
   "metadata": {},
   "source": [
    "A closer look at the data shows that not all users have a `commission_revenue_pre` defined. This is absolutely normal in the real world - some users do not have pre-experiment data because they were not users before the experiment! We need to understand what to do with \"newbies\". Here, we can rely on the experience of Booking.com on how to handle missing data. They [suggest](https://booking.ai/how-booking-com-increases-the-power-of-online-experiments-with-cuped-995d186fff1d) to do nothing to the metric in case pre-experiment data is not available. The idea behind this trick is straightforward. When there is a missing value, you want to fill it with an average derived from the samples where this value is not missing. So, you will plug $\\bar{X}$ into the $X_i$ in CUPED transformation, and it will be canceled out with $E(X)$ with a big enough sample size!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "219e5b9c-d4b9-4876-9f61-e992695cb06a",
   "metadata": {},
   "source": [
    "Now, let's calculate $\\theta$, handle missing values, and make sure that the CUPED transformation yields an unbiased estimate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "eb8e3f3f-5f1c-41f1-8025-3d6fa2fae4a1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8818300441315146"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "not_null_pre = ~exp_data['commission_revenue_pre'].isna()\n",
    "\n",
    "cov_X_Y = np.cov(exp_data[not_null_pre]['commission_revenue_pre'], \n",
    "                 exp_data[not_null_pre]['commission_revenue_exp'],\n",
    "                 ddof=1)[0, 1]\n",
    "\n",
    "var_X = np.var(exp_data[not_null_pre]['commission_revenue_pre'], ddof=1)\n",
    "\n",
    "corr_X_Y = np.corrcoef([exp_data[not_null_pre]['commission_revenue_pre'], \n",
    "                        exp_data[not_null_pre]['commission_revenue_exp']])[0,1]\n",
    "\n",
    "theta = cov_X_Y / var_X\n",
    "\n",
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4f6d2db3-a997-4c17-ac7e-5b57d179c361",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_data['commission_revenue_exp_cuped'] = np.where(\n",
    "    not_null_pre,\n",
    "    exp_data['commission_revenue_exp'] - theta * (\n",
    "        exp_data['commission_revenue_pre'] - exp_data['commission_revenue_pre'].mean()\n",
    "    ),\n",
    "    exp_data['commission_revenue_exp']\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "473b8719-6109-40a9-9698-cb30a6c0c7c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In a test group (new engine B), the CUPED-metric is on average equal to 101.69\n",
      "In a control group (old engine A), the CUPED_metric is on average equal to 99.28\n",
      "Relative increase is equal to 2.4%\n"
     ]
    }
   ],
   "source": [
    "commission_revenue_exp_cuped_test_avg = exp_data[exp_data['group'] == 'test']['commission_revenue_exp_cuped'].mean()\n",
    "commission_revenue_exp_cuped_control_avg = exp_data[exp_data['group'] == 'control']['commission_revenue_exp_cuped'].mean()\n",
    "\n",
    "print(f'In a test group (new engine B), the CUPED-metric is on average equal to {commission_revenue_exp_cuped_test_avg:0.2f}')\n",
    "print(f'In a control group (old engine A), the CUPED_metric is on average equal to {commission_revenue_exp_cuped_control_avg:0.2f}')\n",
    "print(f'Relative increase is equal to {100.0 * (commission_revenue_exp_cuped_test_avg / commission_revenue_exp_cuped_control_avg - 1):0.1f}%')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4b31d79-3181-4f57-9e88-8adb8ac9f6da",
   "metadata": {},
   "source": [
    "In the test group, we got 101.73 before transformation and 101.69 after. A 0.04 absolute difference appears because we have missing pre-experiment data. If there were no missing values, the values would match exactly. Consider the proof of that fact as a take-home exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "852faa75-cf76-45a1-bca2-d96967350f38",
   "metadata": {},
   "source": [
    "Somewhere above, we have proved that the variance of sample average CUPED-metric is $(1-\\rho^2)$ times smaller, where $\\rho$ is a correlation between $Y$ and $X$. So, we must see this in the data also."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e712f0f5-9fb7-491d-a893-62f42891a125",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Var(Y_bar) = 1.00\n",
      "Var(Y_CUPED_bar) = 0.21\n",
      "(1 - rho^2) = 0.210\n",
      "Var(Y_CUPED_bar) / Var(Y_bar) = 0.210\n"
     ]
    }
   ],
   "source": [
    "var_Y = np.std(exp_data[exp_data['group'] == 'test']['commission_revenue_exp'])**2\n",
    "var_Y_bar = var_Y / len(exp_data[exp_data['group'] == 'test'])\n",
    "\n",
    "var_Y_CUPED = np.std(exp_data[exp_data['group'] == 'test']['commission_revenue_exp_cuped'])**2\n",
    "var_Y_CUPED_bar = var_Y_CUPED / len(exp_data[exp_data['group'] == 'test'])\n",
    "\n",
    "print(f'Var(Y_bar) = {var_Y_bar:0.2f}')\n",
    "print(f'Var(Y_CUPED_bar) = {var_Y_CUPED_bar:0.2f}')\n",
    "print(f'(1 - rho^2) = {1 - corr_X_Y**2:0.3f}')\n",
    "print(f'Var(Y_CUPED_bar) / Var(Y_bar) = {var_Y_CUPED_bar / var_Y_bar:0.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "926a8486-da75-4482-becd-5087966ca5f3",
   "metadata": {},
   "source": [
    "Amazing, empirics supports theory!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cdcfb18-0584-47e7-a732-374d4257cce6",
   "metadata": {},
   "source": [
    "A decrease in variance can also be visualized easily. We assume that the metric is normally distributed. Metric with lower variance would have a \"slimmer\" probability mass function with a \"taller\" peak:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "cb234be6-e9d0-4459-ab1c-fa9dd9a4c49b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1IAAAHWCAYAAAB9mLjgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABlb0lEQVR4nO3dd3gUVfv/8c+mJ4QkBFIhhN4kAQ0QI9IECUUEASmiFBF4MMCDgAULTTQKKljBRnkUVEABRQm9iFIURKoISBNSEEgiLSHJ/P7gl/2ypJCBJBvg/bquva7smbNn7pk9O7t3zswZi2EYhgAAAAAABeZg7wAAAAAA4GZDIgUAAAAAJpFIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACaRSAEAAACASSRSAAAAAGASiRQAAAAAmEQihRsybtw4WSyWYllX8+bN1bx5c+vztWvXymKxaMGCBcWy/r59+6pSpUrFsq7rdfbsWT3xxBMKDAyUxWLR8OHD86xbqVIl9e3b97rWc/V7kZ8b2W83EmNRKM7+XlwyMjL0zDPPKCQkRA4ODurUqZO9Qyo0Ja3/FLei2n4zn3/YX3F/T9etW7dY1gWUBCRSsJo1a5YsFov14ebmpuDgYEVHR+udd97Rv//+WyjrOXHihMaNG6ft27cXSnuFqSTHVhCvvvqqZs2apcGDB+uzzz7TY489Vizrvdn3W2GbO3eupk6dau8wCmTGjBmaPHmyunbtqtmzZ+upp56yd0h2cfHiRU2ZMkWRkZHy9vaWm5ubatSooSFDhujPP/+0d3jQ//3zLPvh6Ogof39/de3aVXv37rV3eLeNknC8L45jbEnYTpR8TvYOACXPhAkTVLlyZV26dEkJCQlau3athg8frrfeekvffvutwsPDrXVffPFFPffcc6baP3HihMaPH69KlSqpfv36BX7d8uXLTa3neuQX28cff6ysrKwij+FGrF69WnfffbfGjh17zbr79u2Tg8P1/S/l6vfiZt9vhW3u3LnatWtXviOCJcXq1atVvnx5TZkyxd6h2M0///yjNm3aaOvWrXrggQf0yCOPyNPTU/v27dOXX36pjz76SOnp6fYO07Qb+YznpziOxfkZNmyYGjZsqEuXLmnHjh2aPn261q5dq127dikwMNCusZVE1/M9nZ/r/Q4vTMVxjC0J24mSj0QKObRt21YNGjSwPh89erRWr16tBx54QA8++KD27t0rd3d3SZKTk5OcnIq2G50/f14eHh5ycXEp0vVci7Ozs13XXxBJSUmqU6dOgeq6urpe93rMvBc3w36zp4sXL8rFxaVIfvAWRFJSknx8fAqtvaysLKWnp8vNza3Q2ixqffv21W+//aYFCxaoS5cuNstefvllvfDCC3aK7MbcyGc8P/Y+Fjdp0kRdu3a1Pq9Zs6YGDx6s//3vf3rmmWfsGFnJVBzf07gxN+NxE5dxah8K5L777tNLL72kI0eO6PPPP7eW53bu9YoVK3TvvffKx8dHnp6eqlmzpp5//nlJl0/NaNiwoSSpX79+1lM0Zs2aJen/zq/eunWrmjZtKg8PD+tr8zovPzMzU88//7wCAwNVqlQpPfjggzp27JhNnbyuFbiyzWvFltu1PufOndPIkSMVEhIiV1dX1axZU2+88YYMw7CpZ7FYNGTIEC1atEh169aVq6ur7rjjDsXFxeW+w6+SlJSk/v37KyAgQG5ubqpXr55mz55tXZ59ysuhQ4f0/fffW2M/fPhwnm1evU+yT+386aefNGLECPn5+alUqVJ66KGHdPLkyULdb2+88YbuuecelS1bVu7u7oqIiLiha90K2p6Z92HDhg1q2LCh3NzcVLVqVX344YcFiqV58+b6/vvvdeTIEeu+yN7+7Pfpyy+/1Isvvqjy5cvLw8NDqampOn36tEaNGqWwsDB5enrKy8tLbdu21e+//27TfnYb8+bN0yuvvKIKFSrIzc1NLVu21IEDB2zq7t+/X126dFFgYKDc3NxUoUIF9ejRQykpKTp8+LAsFovWrFmj3bt3W2Ndu3atJPN9e86cObrjjjvk6uqquLg4a3/asGGDhg0bJj8/P/n4+GjQoEFKT09XcnKyevfurTJlyqhMmTJ65plncrSdlZWlqVOn6o477pCbm5sCAgI0aNAgnTlzxqaeYRiaOHGiKlSoIA8PD7Vo0UK7d+8u0Pu1efNmff/99+rfv3+OJEq6nIy88cYbNmWrV69WkyZNVKpUKfn4+Khjx445Ti3LPjb++eefevTRR+Xt7S0/Pz+99NJLMgxDx44dU8eOHeXl5aXAwEC9+eabNq+/8n0eP368ypcvr9KlS6tr165KSUlRWlqahg8fLn9/f3l6eqpfv35KS0uzaePqz/ilS5c0fvx4Va9eXW5ubipbtqzuvfderVixwlonISFB/fr1U4UKFeTq6qqgoCB17NjR5liS27H4WscoSdY+98Ybb+ijjz5S1apV5erqqoYNG+qXX37J8z26liZNmkiSDh48aFN+/PhxPf744woICLB+1mfMmGFdnpiYKCcnJ40fPz5Hm/v27ZPFYtF7771nLUtOTtbw4cOtn4lq1arp9ddftxlxN7ONeX2n5XbMLOhnITe5fU9f73fStY732fbs2aMWLVrIw8ND5cuX16RJk3K0lZaWprFjx6patWpydXVVSEiInnnmmRz9+Gr5HWPNtHsjv1Xy2z8NGjSw+d7Ib/9ffdyUpN9++01t27aVl5eXPD091bJlS23atMnm9Xld95Z93L3y81qpUiU98MADWr58uerXry83NzfVqVNH33zzTb7bgoLhXxQosMcee0zPP/+8li9frgEDBuRaZ/fu3XrggQcUHh6uCRMmyNXVVQcOHNBPP/0kSapdu7YmTJigMWPGaODAgdYvwHvuucfaxqlTp9S2bVv16NFDjz76qAICAvKN65VXXpHFYtGzzz6rpKQkTZ06Va1atdL27dutI2cFUZDYrmQYhh588EGtWbNG/fv3V/369bVs2TI9/fTTOn78eI5TpTZs2KBvvvlGTz75pEqXLq133nlHXbp00dGjR1W2bNk847pw4YKaN2+uAwcOaMiQIapcubLmz5+vvn37Kjk5Wf/9739Vu3ZtffbZZ3rqqadUoUIFjRw5UpLk5+dX4O3PNnToUJUpU0Zjx47V4cOHNXXqVA0ZMkRfffVVoew3SXr77bf14IMPqlevXkpPT9eXX36phx9+WEuWLFH79u1Nx2ymvYK8Dzt37lTr1q3l5+encePGKSMjQ2PHjr1mX5SkF154QSkpKfr777+tfcDT09OmzssvvywXFxeNGjVKaWlpcnFx0Z49e7Ro0SI9/PDDqly5shITE/Xhhx+qWbNm2rNnj4KDg23aeO211+Tg4KBRo0YpJSVFkyZNUq9evbR582ZJUnp6uqKjo5WWlqahQ4cqMDBQx48f15IlS5ScnCw/Pz999tlneuWVV3T27FnFxsZKuvx+mu3bq1ev1rx58zRkyBCVK1dOlSpVsl5XkL3u8ePHa9OmTfroo4/k4+Ojn3/+WRUrVtSrr76qH374QZMnT1bdunXVu3dva7uDBg3SrFmz1K9fPw0bNkyHDh3Se++9p99++00//fSTdbRzzJgxmjhxotq1a6d27dpp27Ztat26dYFOx/v2228lqcDXE65cuVJt27ZVlSpVNG7cOF24cEHvvvuuGjdurG3btuX4Ady9e3fVrl1br732mr7//ntNnDhRvr6++vDDD3Xffffp9ddf15w5czRq1Cg1bNhQTZs2tXl9bGys3N3d9dxzz+nAgQN699135ezsLAcHB505c0bjxo3Tpk2bNGvWLFWuXFljxozJM/Zx48YpNjZWTzzxhBo1aqTU1FT9+uuv2rZtm+6//35JUpcuXbR7924NHTpUlSpVUlJSklasWKGjR4/mOWlMQY5RV5o7d67+/fdfDRo0SBaLRZMmTVLnzp31119/XdcIdvaPxjJlyljLEhMTdffdd1t/sPr5+Wnp0qXq37+/UlNTNXz4cAUEBKhZs2aaN29ejtOhv/rqKzk6Ourhhx+WdPnMiGbNmun48eMaNGiQKlasqJ9//lmjR49WfHx8jut1CnsbC/pZMON6vpMKcrw/c+aM2rRpo86dO6tbt25asGCBnn32WYWFhalt27aSLieGDz74oDZs2KCBAweqdu3a2rlzp6ZMmaI///xTixYtyjPu/I6xBW23MH6rXO23335TmzZtFBQUpPHjxyszM1MTJkzI83s4t+Pm7t271aRJE3l5eemZZ56Rs7OzPvzwQzVv3lzr1q1TZGRknuvPz/79+9W9e3f95z//UZ8+fTRz5kw9/PDDiouLs372cZ0M4P+bOXOmIcn45Zdf8qzj7e1t3HnnndbnY8eONa7sRlOmTDEkGSdPnsyzjV9++cWQZMycOTPHsmbNmhmSjOnTp+e6rFmzZtbna9asMSQZ5cuXN1JTU63l8+bNMyQZb7/9trUsNDTU6NOnzzXbzC+2Pn36GKGhodbnixYtMiQZEydOtKnXtWtXw2KxGAcOHLCWSTJcXFxsyn7//XdDkvHuu+/mWNeVpk6dakgyPv/8c2tZenq6ERUVZXh6etpse2hoqNG+fft827uy7pX7JPv9b9WqlZGVlWUtf+qppwxHR0cjOTnZWnYj+80wDOP8+fM2z9PT0426desa9913X74x5qWg7RX0fejUqZPh5uZmHDlyxFq2Z88ew9HR0SjIYbN9+/Y5ttkw/q/PVqlSJUfMFy9eNDIzM23KDh06ZLi6uhoTJkzI0Ubt2rWNtLQ0a/nbb79tSDJ27txpGIZh/Pbbb4YkY/78+fnG2qxZM+OOO+6wKTPbtx0cHIzdu3fb1M3uT9HR0Tb9KSoqyrBYLMZ//vMfa1lGRoZRoUIFmz71448/GpKMOXPm2LQbFxdnU56UlGS4uLgY7du3t1nP888/b0i6Zv956KGHDEnGmTNn8q2XrX79+oa/v79x6tQpa9nvv/9uODg4GL1797aWZR8bBw4cmGM7LRaL8dprr1nLz5w5Y7i7u9vEmv0+161b10hPT7eW9+zZ07BYLEbbtm1t4oqKisrR567+/NSrVy/f48OZM2cMScbkyZPz3QdXf/4Leow6dOiQIckoW7ascfr0aWvdxYsXG5KM7777Lt/1Zu+TGTNmGCdPnjROnDhhxMXFGdWqVTMsFouxZcsWa93+/fsbQUFBxj///GPTRo8ePQxvb2/r5+/DDz+0+dxkq1Onjs3x4+WXXzZKlSpl/Pnnnzb1nnvuOcPR0dE4evSo6W28ej9mu/qYWdDPQl6u/p42jBv7TirId/j//vc/a1laWpoRGBhodOnSxVr22WefGQ4ODsaPP/5o8/rp06cbkoyffvop3xjyOsYWtN0b/a2Smw4dOhgeHh7G8ePHrWX79+83nJycct3/uR03O3XqZLi4uBgHDx60lp04ccIoXbq00bRpU2tZbu+pYfzfcffQoUPWstDQUEOS8fXXX1vLUlJSjKCgIJvfc7g+nNoHUzw9PfOdvS/7WovFixdf9wQDrq6u6tevX4Hr9+7dW6VLl7Y+79q1q4KCgvTDDz9c1/oL6ocffpCjo6OGDRtmUz5y5EgZhqGlS5falLdq1UpVq1a1Pg8PD5eXl5f++uuva64nMDBQPXv2tJY5Oztr2LBhOnv2rNatW1cIW/N/Bg4caHPKQJMmTZSZmakjR44U2jquHCk8c+aMUlJS1KRJE23btq3I27vW+5CZmally5apU6dOqlixorVe7dq1FR0dfV3xXa1Pnz45RktdXV2t10llZmbq1KlT1tNNctuOfv362Vyrkv0f0+zt8Pb2liQtW7ZM58+fNxWf2b7drFmzPK/N69+/v01/ioyMlGEY6t+/v7XM0dFRDRo0sPkszJ8/X97e3rr//vv1zz//WB8RERHy9PTUmjVrJF0eIUpPT9fQoUNt1lPQi9BTU1MlyeYYkpf4+Hht375dffv2la+vr7U8PDxc999/f67HnCeeeCLHdl69/T4+PqpZs2aux4LevXvbjDZk77/HH3/cpl5kZKSOHTumjIyMPOP38fHR7t27tX///lyXu7u7y8XFRWvXri3QKWPZzB6junfvbjN6dHXfvZbHH39cfn5+Cg4OVps2bZSSkqLPPvvMeiqWYRj6+uuv1aFDBxmGYdN/oqOjlZKSYv1Mde7cWU5OTjYj7rt27dKePXvUvXt3a9n8+fPVpEkTlSlTxqa9Vq1aKTMzU+vXry/UbbxSQT8LZl3vd9K1eHp66tFHH7U+d3FxUaNGjXJ8vmvXrq1atWrZbNN9990nSde9TQVttzB+q1wpMzNTK1euVKdOnWzOHqhWrZp1FO5qVx83MzMztXz5cnXq1ElVqlSxlgcFBemRRx7Rhg0brMcrs4KDg/XQQw9Zn3t5eal379767bfflJCQcF1t4jISKZhy9uzZfH9wdO/eXY0bN9YTTzyhgIAA9ejRQ/PmzTN1oCpfvrypi5mrV69u89xisahatWr5Xh9UGI4cOaLg4OAc+6N27drW5Ve68kd5tjJlylzzB8uRI0dUvXr1HJMR5LWeG3V1nNk/Bsz8sLqWJUuW6O6775abm5t8fX3l5+enadOmKSUlpcjbu9b7cPLkSV24cCFHv5IuX9ReGCpXrpyjLCsrS1OmTFH16tXl6uqqcuXKyc/PTzt27CjQdlz9PlWuXFkjRozQJ598onLlyik6Olrvv/9+gfax2b6d2/bkFWd2ghcSEpKj/Mo+tn//fqWkpMjf319+fn42j7NnzyopKckmlqvfLz8/P5sfsnnx8vKSpALd3iF7Xbn1g9q1a+uff/7RuXPnbMpz2343NzeVK1cuR3lunzEz+y8rKyvf93fChAlKTk5WjRo1FBYWpqefflo7duywLnd1ddXrr7+upUuXKiAgQE2bNtWkSZOu+UPL7DHqRo8xY8aM0YoVK7Rw4UL17t1bKSkpNus+efKkkpOT9dFHH+XoO9n/pMvuP+XKlVPLli01b9486+u/+uorOTk5qXPnztay/fv3Ky4uLkd7rVq1smmvsLbxSgX9LJh1vd9J11KhQoUc1+9c3e7+/fu1e/fuHNtTo0YNSTn3Z0EVtN3C+K1ypaSkJF24cEHVqlXLsSy3MinncfPkyZM6f/58nseXrKysHNd/F1S1atVyvCfZ+6Sofyvd6rhGCgX2999/KyUlJc+DgnT5P5rr16/XmjVr9P333ysuLk5fffWV7rvvPi1fvlyOjo7XXI+Z65oKKq+bEWZmZhYopsKQ13qMqy6wt7eijvPHH3/Ugw8+qKZNm+qDDz5QUFCQnJ2dNXPmTM2dO7fI2ysJ70NuffzVV1/VSy+9pMcff1wvv/yyfH195eDgoOHDh+f65V6Q7XjzzTfVt29fLV68WMuXL9ewYcMUGxurTZs2qUKFCkW6PdeKM7fyK2PPysqSv7+/5syZk+vrr+f6v9zUqlVL0uXr4rJHDQpTbttppg+a2X95tZGtadOmOnjwoLU/fPLJJ5oyZYqmT59uHTkbPny4OnTooEWLFmnZsmV66aWXFBsbq9WrV+vOO+/Ms20zbvQzGBYWZk1gOnXqpPPnz2vAgAG69957FRISYv28PProo+rTp0+ubVx5G48ePXqoX79+2r59u+rXr6958+apZcuWNsluVlaW7r///jxnBcz+UWpmGy0WS67bnJmZafO8qD4LRXUsLEi7WVlZCgsL01tvvZVr3av/UVBQBW23MH6r3Kgb+a2T328aFC8SKRTYZ599JknXPL3JwcFBLVu2VMuWLfXWW2/p1Vdf1QsvvKA1a9aoVatWhX6H9atPUzEMQwcOHLD5oixTpoySk5NzvPbIkSM2Q+hmYgsNDdXKlSv177//2vzn/o8//rAuLwyhoaHasWOHsrKybP7rWtjruRFm9tvXX38tNzc3LVu2zGZ65pkzZ17Xugu7PT8/P7m7u+d6+tO+ffsK1Mb19PEFCxaoRYsW+vTTT23Kk5OTc4xemBEWFqawsDC9+OKL+vnnn9W4cWNNnz5dEydOzPM1xdW381O1alWtXLlSjRs3zvcHR3Ys+/fvt/ksnzx5skD/We/QoYNiY2P1+eefXzORyl5Xbv3gjz/+ULly5VSqVKlrrtOefH191a9fP/Xr109nz55V06ZNNW7cOJtTEKtWraqRI0dq5MiR2r9/v+rXr68333zTZsbWK9n7GPXaa69p4cKFeuWVVzR9+nT5+fmpdOnSyszMtCZc+enUqZMGDRpkPb3vzz//1OjRo23qVK1aVWfPni1QewVVpkyZXE+ju3oEr6CfheJSGN/hVatW1e+//66WLVteV3t5vcZMu4X5W8Xf319ubm45Zk6VlGtZbvz8/OTh4ZHn8cXBwcGaCGaPcCYnJ9vcviKvM1QOHDggwzBstin7RuN5TSKDguHUPhTI6tWr9fLLL6ty5crq1atXnvVOnz6doyz7RnbZU49m/9DILbG5Hv/73/9sTstZsGCB4uPjbc5Lrlq1qjZt2mQzi9eSJUtyDJObia1du3bKzMy0mR5XkqZMmSKLxZLnedFmtWvXTgkJCTbn8GdkZOjdd9+Vp6enmjVrVijruRFm9pujo6MsFovNf84OHz6c7yxNxd1edHS0Fi1apKNHj1rL9+7dq2XLlhWojVKlSpk+TdHR0THHf4Lnz5+v48ePm2onW2pqao7rZcLCwuTg4HDN6YWLq2/np1u3bsrMzNTLL7+cY1lGRoa1r7Vq1UrOzs569913bfbf1bOo5SUqKkpt2rTRJ598kmufSU9P16hRoyRdvlahfv36mj17tk1f37Vrl5YvX6527doVePvs4dSpUzbPPT09Va1aNWt/OH/+vC5evGhTp2rVqipdunS+fcbex6iqVauqS5cumjVrlhISEuTo6KguXbro66+/1q5du3LUv/p2Dj4+PoqOjta8efP05ZdfysXFRZ06dbKp061bN23cuDHXY0BycnK+16blF/cff/xhE8/vv/9unTnuynUX5LNQXArjO7xbt246fvy4Pv744xzLLly4kOMU2dxiyO0YW9B2C/u3iqOjo1q1aqVFixbpxIkT1vIDBw7kuKY0vzZat26txYsX25xul5iYqLlz5+ree++1noqcfW3bldfmnTt3LsctB7KdOHFCCxcutD5PTU3V//73P9WvX996E+tLly7pjz/+UHx8fIHixWWMSCGHpUuX6o8//lBGRoYSExO1evVqrVixQqGhofr222/zvWHchAkTtH79erVv316hoaFKSkrSBx98oAoVKujee++VdPkA4OPjo+nTp6t06dIqVaqUIiMj873OIj++vr6699571a9fPyUmJmrq1KmqVq2azRTtTzzxhBYsWKA2bdqoW7duOnjwoD7//HObC23NxtahQwe1aNFCL7zwgg4fPqx69epp+fLlWrx4sYYPH56j7es1cOBAffjhh+rbt6+2bt2qSpUqacGCBfrpp580derUAl0kX9TM7Lf27dvrrbfeUps2bfTII48oKSlJ77//vqpVq2ZzvUZBFXZ7kjR+/HjFxcWpSZMmevLJJ60/Cu+4444CtRkREaGvvvpKI0aMUMOGDeXp6akOHTrk+5oHHnhAEyZMUL9+/XTPPfdo586dmjNnjs0oixmrV6/WkCFD9PDDD6tGjRrKyMjQZ599Zv2RmZ/i6tv5adasmQYNGqTY2Fht375drVu3lrOzs/bv36/58+fr7bffVteuXeXn56dRo0YpNjZWDzzwgNq1a6fffvtNS5cuLfBI3v/+9z+1bt1anTt3VocOHdSyZUuVKlVK+/fv15dffqn4+HjrvaQmT56stm3bKioqSv3797dOf+7t7a1x48YV4R65cXXq1FHz5s0VEREhX19f/frrr1qwYIGGDBki6fJ/qFu2bKlu3bqpTp06cnJy0sKFC5WYmKgePXrk2W5JOEY9/fTTmjdvnqZOnarXXntNr732mtasWaPIyEgNGDBAderU0enTp7Vt2zatXLkyxw/p7t2769FHH9UHH3yg6OjoHDepfvrpp/Xtt9/qgQceUN++fRUREaFz585p586dWrBggQ4fPmx65Pjxxx/XW2+9pejoaPXv319JSUmaPn267rjjDptJBQr6WSguhfEd/thjj2nevHn6z3/+ozVr1qhx48bKzMzUH3/8oXnz5mnZsmVq0KBBnq/P6xhb0HaL4rfKuHHjtHz5cjVu3FiDBw+2/jOqbt261ttBXMvEiROt97d68skn5eTkpA8//FBpaWk29+Jq3bq1KlasqP79++vpp5+Wo6OjZsyYIT8/P5t/AGarUaOG+vfvr19++UUBAQGaMWOGEhMTbc7cOH78uGrXrq0+ffpc835ZuEIxzhCIEi572szsh4uLixEYGGjcf//9xttvv20zzXa2q6fgXLVqldGxY0cjODjYcHFxMYKDg42ePXvmmDJ28eLFRp06dazTgmZPL5rbVMzZ8pr+/IsvvjBGjx5t+Pv7G+7u7kb79u1tpq3O9uabbxrly5c3XF1djcaNGxu//vprrtPP5hVbbtN4//vvv8ZTTz1lBAcHG87Ozkb16tWNyZMn20zDbBiXpzqNiYnJEVNBp/dOTEw0+vXrZ5QrV85wcXExwsLCcp2StTCmP796+vvs/bxmzRpr2Y3ut08//dSoXr264erqatSqVcuYOXNmrtO5FnT/FLQ9M+/DunXrjIiICMPFxcWoUqWKMX369DynnL3a2bNnjUceecTw8fExJFm3P3tf5jYl+cWLF42RI0caQUFBhru7u9G4cWNj48aNefb7q9vInnY5e7//9ddfxuOPP25UrVrVcHNzM3x9fY0WLVoYK1eutHldXp+5G+3befWn7H149bTDffr0MUqVKpWjnY8++siIiIgw3N3djdKlSxthYWHGM888Y5w4ccJaJzMz0xg/frx13zVv3tzYtWtXgfuPYVyeQv+NN94wGjZsaHh6ehouLi5G9erVjaFDh9pMEW0YhrFy5UqjcePGhru7u+Hl5WV06NDB2LNnzw1t59XvQ17vs5n9evX2T5w40WjUqJHh4+NjuLu7G7Vq1TJeeeUV6/Tq//zzjxETE2PUqlXLKFWqlOHt7W1ERkYa8+bNyxHr1Z//ghyjsvtobtOrSzLGjh2bo/xK+X1+DMMwmjdvbnh5eVlv1ZCYmGjExMQYISEhhrOzsxEYGGi0bNnS+Oijj3K8NjU11XB3d88xjfuV/v33X2P06NFGtWrVDBcXF6NcuXLGPffcY7zxxhvWfWh2Gz///HOjSpUqhouLi1G/fn1j2bJluR4zDaNgn4Xc3OixMDdmv8Nz26b09HTj9ddfN+644w7D1dXVKFOmjBEREWGMHz/eSElJyXf9eR1jC9rujf5WycuqVauMO++803BxcTGqVq1qfPLJJ8bIkSMNNzc3m3p57X/DMIxt27YZ0dHRhqenp+Hh4WG0aNHC+Pnnn3PU27p1qxEZGWm4uLgYFStWNN566608pz9v3769sWzZMiM8PNz6PZnXd0hBj5m4zGIYJexKdwAAAOAW0KlTp3xvO1DUKlWqpLp162rJkiV2Wf+tjmukAAAAgBt04cIFm+f79+/XDz/8oObNm9snIBQ5rpECAAAAblCVKlXUt29fValSRUeOHNG0adPk4uKS57T5uPmRSAEAAAA3qE2bNvriiy+UkJAgV1dXRUVF6dVXX831Bu+4NXCNFAAAAACYxDVSAAAAAGASiRQAAAAAmMQ1UpKysrJ04sQJlS5dWhaLxd7hAAAAALATwzD077//Kjg4WA4OeY87kUhJOnHihEJCQuwdBgAAAIAS4tixY6pQoUKey0mkJJUuXVrS5Z3l5eVl52gAAAAA2EtqaqpCQkKsOUJeSKQk6+l8Xl5eJFIAAAAArnnJD5NNAAAAAIBJJFIAAAAAYBKJFAAAAACYxDVSAAAAsDIMQxkZGcrMzLR3KECRcHR0lJOT0w3f9ohECgAAAJKk9PR0xcfH6/z58/YOBShSHh4eCgoKkouLy3W3QSIFAAAAZWVl6dChQ3J0dFRwcLBcXFxu+D/2QEljGIbS09N18uRJHTp0SNWrV8/3prv5IZECAACA0tPTlZWVpZCQEHl4eNg7HKDIuLu7y9nZWUeOHFF6errc3Nyuqx0mmwAAAIDV9f53HriZFEY/55MCAAAAACaRSAEAAACASSRSAAAAuK1YLBYtWrQo3zp9+/ZVp06dCtzm4cOHZbFYtH379huKDTcPu042MW3aNE2bNk2HDx+WJN1xxx0aM2aM2rZtK0m6ePGiRo4cqS+//FJpaWmKjo7WBx98oICAAGsbR48e1eDBg7VmzRp5enqqT58+io2NlZMT82gAAAAUhsXbjxfbujrWL2/6NX379lVycvI1k6Ns8fHxKlOmjKTLCVDlypX122+/qX79+tY6b7/9tgzDMB0Lbh92HZGqUKGCXnvtNW3dulW//vqr7rvvPnXs2FG7d++WJD311FP67rvvNH/+fK1bt04nTpxQ586dra/PzMxU+/btlZ6erp9//lmzZ8/WrFmzNGbMGHttEgAAAEq4wMBAubq65lvH29tbPj4+xRMQbkp2TaQ6dOigdu3aqXr16qpRo4ZeeeUVeXp6atOmTUpJSdGnn36qt956S/fdd58iIiI0c+ZM/fzzz9q0aZMkafny5dqzZ48+//xz1a9fX23bttXLL7+s999/X+np6fbcNAAAANhB8+bNNWzYMD3zzDPy9fVVYGCgxo0bZ1PnylP7KleuLEm68847ZbFY1Lx5c0k5T+2Li4vTvffeKx8fH5UtW1YPPPCADh48WAxbhJKqxFwjlZmZqS+//FLnzp1TVFSUtm7dqkuXLqlVq1bWOrVq1VLFihW1ceNGSdLGjRsVFhZmc6pfdHS0UlNTraNauUlLS1NqaqrNAwAAALeG2bNnq1SpUtq8ebMmTZqkCRMmaMWKFbnW3bJliyRp5cqVio+P1zfffJNrvXPnzmnEiBH69ddftWrVKjk4OOihhx5SVlZWkW0HSja7X0i0c+dORUVF6eLFi/L09NTChQtVp04dbd++XS4uLjmGVAMCApSQkCBJSkhIsEmispdnL8tLbGysxo8fX7gbAgDAlXYusHcEOYV1tXcEQLEIDw/X2LFjJUnVq1fXe++9p1WrVun+++/PUdfPz0+SVLZsWQUGBubZZpcuXWyez5gxQ35+ftqzZ4/q1q1biNHjZmH3EamaNWtq+/bt2rx5swYPHqw+ffpoz549RbrO0aNHKyUlxfo4duxYka4PAAAAxSc8PNzmeVBQkJKSkm6ozf3796tnz56qUqWKvLy8VKlSJUmXJz7D7cnuI1IuLi6qVq2aJCkiIkK//PKL3n77bXXv3l3p6elKTk62GZVKTEy0/rcgMDDQOhx75fLsZXlxdXW95gWGAAAAuDk5OzvbPLdYLDd8Cl6HDh0UGhqqjz/+WMHBwcrKylLdunW5Lv82ZvcRqatlZWUpLS1NERERcnZ21qpVq6zL9u3bp6NHjyoqKkqSFBUVpZ07d9r8h2HFihXy8vJSnTp1ij12AAAA3FxcXFwkXb5ePy+nTp3Svn379OKLL6ply5aqXbu2zpw5U1whooSy64jU6NGj1bZtW1WsWFH//vuv5s6dq7Vr12rZsmXy9vZW//79NWLECPn6+srLy0tDhw5VVFSU7r77bklS69atVadOHT322GOaNGmSEhIS9OKLLyomJoYRJwAAAFyTv7+/3N3dFRcXpwoVKsjNzU3e3t42dcqUKaOyZcvqo48+UlBQkI4eParnnnvOThGjpLBrIpWUlKTevXsrPj5e3t7eCg8P17Jly6wXAk6ZMkUODg7q0qWLzQ15szk6OmrJkiUaPHiwoqKiVKpUKfXp00cTJkyw1yYBAADccq7nJrk3CycnJ73zzjuaMGGCxowZoyZNmmjt2rU2dRwcHPTll19q2LBhqlu3rmrWrKl33nnHOlU6bk8Wg1s2KzU1Vd7e3kpJSZGXl5e9wwEA3AqYtQ83mYsXL+rQoUOqXLmy3Nzc7B0OUKTy6+8FzQ1K3DVSAAAAAFDSkUgBAAAAgEkkUgAAAABgEokUAAAAAJhEIgUAAAAAJpFIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACY52TsAAAAAlHA7FxTfusK6Ft+6cE3nz5/XY489phUrVujff//VmTNn5OPjY++wrqlSpUoaPny4hg8fXmTrYEQKAAAAN7WEhAQNHTpUVapUkaurq0JCQtShQwetWrXKWsdisWjRokU5Xtu3b1916tTJ+rx58+ayWCyyWCxyc3NTnTp19MEHH1iXz5o1y7r8yoebm5tNm9nlzs7OCggI0P33368ZM2YoKysr3225Oh57mz17tn788Uf9/PPPio+Pl7e3t71DKjFIpAAAAHDTOnz4sCIiIrR69WpNnjxZO3fuVFxcnFq0aKGYmJjranPAgAGKj4/Xnj171K1bN8XExOiLL76wLvfy8lJ8fLzN48iRIzZttGnTRvHx8Tp8+LCWLl2qFi1a6L///a8eeOABZWRk3NA2S9KlS5duuI2COHjwoGrXrq26desqMDBQFovFdBuZmZnXTCBvRiRSAAAAuGk9+eSTslgs2rJli7p06aIaNWrojjvu0IgRI7Rp06bratPDw0OBgYGqUqWKxo0bp+rVq+vbb7+1LrdYLAoMDLR5BAQE2LTh6uqqwMBAlS9fXnfddZeef/55LV68WEuXLtWsWbNyXe+4ceM0e/ZsLV682DqitXbtWh0+fFgWi0VfffWVmjVrJjc3N82ZM0enTp1Sz549Vb58eXl4eCgsLMwm4ZMuj7ANGzZMzzzzjHx9fRUYGKhx48ZZlxuGoXHjxqlixYpydXVVcHCwhg0bZn3tm2++qfXr18tisah58+aSpDNnzqh3794qU6aMPDw81LZtW+3fv9/a5qxZs+Tj46Nvv/1WderUkaurq44ePapKlSpp4sSJ6t27tzw9PRUaGqpvv/1WJ0+eVMeOHeXp6anw8HD9+uuvNtuwYcMGNWnSRO7u7goJCdGwYcN07tw56/KkpCR16NBB7u7uqly5subMmVPg9/pGkEgBAADgpnT69GnFxcUpJiZGpUqVyrG8sK7lcXd3V3p6+g23c99996levXr65ptvcl0+atQodevWzTqaFR8fr3vuuce6/LnnntN///tf7d27V9HR0bp48aIiIiL0/fffa9euXRo4cKAee+wxbdmyxabd2bNnq1SpUtq8ebMmTZqkCRMmaMWKFZKkr7/+WlOmTNGHH36o/fv3a9GiRQoLC5MkffPNNxowYICioqIUHx9vjbtv37769ddf9e2332rjxo0yDEPt2rWzGSU7f/68Xn/9dX3yySfavXu3/P39JUlTpkxR48aN9dtvv6l9+/Z67LHH1Lt3bz366KPatm2bqlatqt69e8swDEmXR8TatGmjLl26aMeOHfrqq6+0YcMGDRkyxLquvn376tixY1qzZo0WLFigDz74QElJSTf6dl0Tk00AAADgpnTgwAEZhqFatWoVSfuZmZn64osvtGPHDg0cONBanpKSIk9PT5u6TZo00dKlS6/ZZq1atbRjx45cl3l6esrd3V1paWkKDAzMsXz48OHq3LmzTdmoUaOsfw8dOlTLli3TvHnz1KhRI2t5eHi4xo4dK0mqXr263nvvPa1atUr333+/jh49qsDAQLVq1UrOzs6qWLGi9bW+vr7y8PCQi4uLNZ79+/fr22+/1U8//WRN8ubMmaOQkBAtWrRIDz/8sKTLpx5+8MEHqlevnk287dq106BBgyRJY8aM0bRp09SwYUPr65599llFRUUpMTFRgYGBio2NVa9evayTRlSvXl3vvPOOmjVrpmnTpuno0aNaunSptmzZooYNG0qSPv30U9WuXftab8UNI5ECAADATSl71KKwffDBB/rkk0+Unp4uR0dHPfXUUxo8eLB1eenSpbVt2zab17i7uxeobcMwrus6I0lq0KCBzfPMzEy9+uqrmjdvno4fP6709HSlpaXJw8PDpl54eLjN86CgIOuIzcMPP6ypU6eqSpUqatOmjdq1a6cOHTrIySn3NGHv3r1ycnJSZGSktaxs2bKqWbOm9u7day1zcXHJsd6rY8k+HTJ7BOzKsqSkJAUGBur333/Xjh07bE7XMwxDWVlZOnTokP788085OTkpIiLCurxWrVrFMrMgiRQAAABuStWrV5fFYtEff/xxzbqlS5dWSkpKjvLk5OQcM9H16tVLL7zwgtzd3RUUFCQHB9urYRwcHFStWrXrinnv3r2qXLnydb326tMXJ0+erLfffltTp05VWFiYSpUqpeHDh+c4DdHZ2dnmucVisU7+EBISon379mnlypVasWKFnnzySU2ePFnr1q3L8Toz3N3dc00Yr2wze3luZdnxnT17VoMGDbJet3WlihUr6s8//7zuGG8U10gBAADgpuTr66vo6Gi9//77NpMPZEtOTrb+XbNmTW3dutVmeWZmpn7//XfVqFHDptzb21vVqlVT+fLlcyRRN2L16tXauXOnunTpkmcdFxcXZWZmFqi9n376SR07dtSjjz6qevXqqUqVKteVWLi7u6tDhw565513tHbtWm3cuFE7d+7MtW7t2rWVkZGhzZs3W8tOnTqlffv2qU6dOqbXfS133XWX9uzZo2rVquV4uLi4qFatWsrIyLB5b/ft22fz3hcVRqQAAABw03r//ffVuHFjNWrUSBMmTFB4eLgyMjK0YsUKTZs2zXq62YgRI9S/f3/VqlVL999/v86dO6d3331XZ86c0RNPPGFqnYZhKCEhIUe5v7+/NfFKS0tTQkKCMjMzlZiYqLi4OMXGxuqBBx5Q796982y7UqVKWrZsmfbt26eyZcvme9+m6tWra8GCBfr5559VpkwZvfXWW0pMTDSV0MyaNUuZmZmKjIyUh4eHPv/8c7m7uys0NDTPdXbs2FEDBgzQhx9+qNKlS+u5555T+fLl1bFjxwKvt6CeffZZ3X333RoyZIieeOIJlSpVSnv27NGKFSv03nvvqWbNmmrTpo0GDRqkadOmycnJScOHDy/wqZY3gkQKAAAA+Qvrau8I8lSlShVt27ZNr7zyikaOHKn4+Hj5+fkpIiJC06ZNs9br2bOnDMPQW2+9peeee04eHh6KiIjQ+vXrc0xdfi2pqakKCgrKUR4fH2+dlCEuLk5BQUFycnJSmTJlVK9ePb3zzjvq06dPvqNcAwYM0Nq1a9WgQQOdPXtWa9asUaVKlXKt++KLL+qvv/5SdHS0PDw8NHDgQHXq1CnXUxjz4uPjo9dee00jRoxQZmamwsLC9N1336ls2bJ5vmbmzJnWe2Klp6eradOm+uGHH27oVMC8hIeHa926dXrhhRfUpEkTGYahqlWrqnv37jbxPPHEE2rWrJkCAgI0ceJEvfTSS4Uey9UsRlFdpXcTSU1Nlbe3t1JSUuTl5WXvcAAAt4KdC+wdQU4l+Mcw7O/ixYs6dOiQKleuLDc3N3uHAxSp/Pp7QXMDrpECAAAAAJNIpAAAAADAJBIpAAAAADCJRAoAAAAATCKRAgAAgBXzkOF2UBj9nEQKAAAA1qmrz58/b+dIgKKX3c9vZMp27iMFAAAAOTo6ysfHR0lJSZIkDw8PWSwWO0cFFC7DMHT+/HklJSXJx8dHjo6O190WiRQAAAAkyXoz2exkCrhV+fj4WPv79SKRAgAAgCTJYrEoKChI/v7+unTpkr3DAYqEs7PzDY1EZSORAgAAgA1HR8dC+aEJ3MqYbAIAAAAATCKRAgAAAACTSKQAAAAAwCQSKQAAAAAwiUQKAAAAAEwikQIAAAAAk0ikAAAAAMAkEikAAAAAMIlECgAAAABMIpECAAAAAJNIpAAAAADAJBIpAAAAADCJRAoAAAAATCKRAgAAAACTSKQAAAAAwCQSKQAAAAAwiUQKAAAAAEwikQIAAAAAk0ikAAAAAMAkuyZSsbGxatiwoUqXLi1/f3916tRJ+/bts6nTvHlzWSwWm8d//vMfmzpHjx5V+/bt5eHhIX9/fz399NPKyMgozk0BAAAAcBtxsufK161bp5iYGDVs2FAZGRl6/vnn1bp1a+3Zs0elSpWy1hswYIAmTJhgfe7h4WH9OzMzU+3bt1dgYKB+/vlnxcfHq3fv3nJ2dtarr75arNsDAAAA4PZg10QqLi7O5vmsWbPk7++vrVu3qmnTptZyDw8PBQYG5trG8uXLtWfPHq1cuVIBAQGqX7++Xn75ZT377LMaN26cXFxcinQbAAAAANx+StQ1UikpKZIkX19fm/I5c+aoXLlyqlu3rkaPHq3z589bl23cuFFhYWEKCAiwlkVHRys1NVW7d+/OdT1paWlKTU21eQAAAABAQdl1ROpKWVlZGj58uBo3bqy6detayx955BGFhoYqODhYO3bs0LPPPqt9+/bpm2++kSQlJCTYJFGSrM8TEhJyXVdsbKzGjx9fRFsCAAAA4FZXYhKpmJgY7dq1Sxs2bLApHzhwoPXvsLAwBQUFqWXLljp48KCqVq16XesaPXq0RowYYX2empqqkJCQ6wscAAAAwG2nRJzaN2TIEC1ZskRr1qxRhQoV8q0bGRkpSTpw4IAkKTAwUImJiTZ1sp/ndV2Vq6urvLy8bB4AAAAAUFB2TaQMw9CQIUO0cOFCrV69WpUrV77ma7Zv3y5JCgoKkiRFRUVp586dSkpKstZZsWKFvLy8VKdOnSKJGwAAAMDtza6n9sXExGju3LlavHixSpcubb2mydvbW+7u7jp48KDmzp2rdu3aqWzZstqxY4eeeuopNW3aVOHh4ZKk1q1bq06dOnrsscc0adIkJSQk6MUXX1RMTIxcXV3tuXkAAAAAblEWwzAMu63cYsm1fObMmerbt6+OHTumRx99VLt27dK5c+cUEhKihx56SC+++KLN6XhHjhzR4MGDtXbtWpUqVUp9+vTRa6+9JienguWJqamp8vb2VkpKCqf5AQAKx84F9o4gp7Cu9o4AAEq8guYGdh2RulYOFxISonXr1l2zndDQUP3www+FFRYAAAAA5KtETDYBAAAAADcTEikAAAAAMIlECgAAAABMIpECAAAAAJNIpAAAAADAJBIpAAAAADCJRAoAAAAATCKRAgAAAACTSKQAAAAAwCQSKQAAAAAwycneAQAAcMN2LrB3BACA2wwjUgAAAABgEokUAAAAAJhEIgUAAAAAJpFIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACaRSAEAAACASSRSAAAAAGASiRQAAAAAmEQiBQAAAAAmkUgBAAAAgEkkUgAAAABgEokUAAAAAJhEIgUAAAAAJpFIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACaRSAEAAACASSRSAAAAAGASiRQAAAAAmEQiBQAAAAAmkUgBAAAAgEkkUgAAAABgEokUAAAAAJhEIgUAAAAAJpFIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACaRSAEAAACASSRSAAAAAGASiRQAAAAAmEQiBQAAAAAmkUgBAAAAgEkkUgAAAABgEokUAAAAAJhEIgUAAAAAJtk1kYqNjVXDhg1VunRp+fv7q1OnTtq3b59NnYsXLyomJkZly5aVp6enunTposTERJs6R48eVfv27eXh4SF/f389/fTTysjIKM5NAQAAAHAbsWsitW7dOsXExGjTpk1asWKFLl26pNatW+vcuXPWOk899ZS+++47zZ8/X+vWrdOJEyfUuXNn6/LMzEy1b99e6enp+vnnnzV79mzNmjVLY8aMsccmAQAAALgNWAzDMOwdRLaTJ0/K399f69atU9OmTZWSkiI/Pz/NnTtXXbt2lST98ccfql27tjZu3Ki7775bS5cu1QMPPKATJ04oICBAkjR9+nQ9++yzOnnypFxcXK653tTUVHl7eyslJUVeXl5Fuo0AgCKwc4G9I7g5hHW1dwQAUOIVNDcoUddIpaSkSJJ8fX0lSVu3btWlS5fUqlUra51atWqpYsWK2rhxoyRp48aNCgsLsyZRkhQdHa3U1FTt3r071/WkpaUpNTXV5gEAAAAABVViEqmsrCwNHz5cjRs3Vt26dSVJCQkJcnFxkY+Pj03dgIAAJSQkWOtcmURlL89elpvY2Fh5e3tbHyEhIYW8NQAAAABuZSUmkYqJidGuXbv05ZdfFvm6Ro8erZSUFOvj2LFjRb5OAAAAALcOJ3sHIElDhgzRkiVLtH79elWoUMFaHhgYqPT0dCUnJ9uMSiUmJiowMNBaZ8uWLTbtZc/ql13naq6urnJ1dS3krQAAAABwu7DriJRhGBoyZIgWLlyo1atXq3LlyjbLIyIi5OzsrFWrVlnL9u3bp6NHjyoqKkqSFBUVpZ07dyopKclaZ8WKFfLy8lKdOnWKZ0MAAAAA3FbsOiIVExOjuXPnavHixSpdurT1miZvb2+5u7vL29tb/fv314gRI+Tr6ysvLy8NHTpUUVFRuvvuuyVJrVu3Vp06dfTYY49p0qRJSkhI0IsvvqiYmBhGnQAAAAAUCbsmUtOmTZMkNW/e3KZ85syZ6tu3ryRpypQpcnBwUJcuXZSWlqbo6Gh98MEH1rqOjo5asmSJBg8erKioKJUqVUp9+vTRhAkTimszAAAAANxmStR9pOyF+0gBwE2O+0gVDPeRAoBruinvIwUAAAAANwMSKQAAAAAwiUQKAAAAAEwikQIAAAAAk0ikAAAAAMAkEikAAAAAMIlECgAAAABMIpECAAAAAJNIpAAAAADAJBIpAAAAADCJRAoAAAAATCKRAgAAAACTSKQAAAAAwCQSKQAAAAAwiUQKAAAAAEwikQIAAAAAk0ikAAAAAMAkEikAAAAAMMnJ3gEAAIBisnOBvSPIKayrvSMAgOvCiBQAAAAAmEQiBQAAAAAmmU6ktm3bpp07d1qfL168WJ06ddLzzz+v9PT0Qg0OAAAAAEoi04nUoEGD9Oeff0qS/vrrL/Xo0UMeHh6aP3++nnnmmUIPEAAAAABKGtOJ1J9//qn69etLkubPn6+mTZtq7ty5mjVrlr7++uvCjg8AAAAAShzTiZRhGMrKypIkrVy5Uu3atZMkhYSE6J9//inc6AAAAACgBDKdSDVo0EATJ07UZ599pnXr1ql9+/aSpEOHDikgIKDQAwQAAACAksZ0IjV16lRt27ZNQ4YM0QsvvKBq1apJkhYsWKB77rmn0AMEAAAAgJLG1A15MzMzlZycrPXr16tMmTI2yyZPnixHR8dCDQ4AAAAASiJTI1KOjo5q3bq1kpOTcyxzc3OTs7NzYcUFAAAAACWW6VP76tatq7/++qsoYgEAAACAm4LpRGrixIkaNWqUlixZovj4eKWmpto8AAAAAOBWZ+oaKUnW6c4ffPBBWSwWa7lhGLJYLMrMzCy86AAAAACgBDKdSK1Zs6Yo4gAAAACAm4bpRKpZs2ZFEQcAAAAA3DRMXyMlST/++KMeffRR3XPPPTp+/Lgk6bPPPtOGDRsKNTgAAAAAKIlMJ1Jff/21oqOj5e7urm3btiktLU2SlJKSoldffbXQAwQAAACAkua6Zu2bPn26Pv74Y5v7RjVu3Fjbtm0r1OAAAAAAoCQynUjt27dPTZs2zVHu7e2d6416AQAAAOBWYzqRCgwM1IEDB3KUb9iwQVWqVCmUoAAAAACgJDOdSA0YMED//e9/tXnzZlksFp04cUJz5szRqFGjNHjw4KKIEQAAAABKFNPTnz/33HPKyspSy5Ytdf78eTVt2lSurq4aNWqUhg4dWhQxAgAAAECJYjEMw7ieF6anp+vAgQM6e/as6tSpI09Pz8KOrdikpqbK29tbKSkp8vLysnc4AACzdi6wdwS4XmFd7R0BANgoaG5wXfeRkiQXFxfVqVNHtWrV0sqVK7V3797rbQoAAAAAbiqmE6lu3brpvffekyRduHBBDRs2VLdu3RQeHq6vv/660AMEAAAAgJLGdCK1fv16NWnSRJK0cOFCZWVlKTk5We+8844mTpxY6AECAAAAQEljOpFKSUmRr6+vJCkuLk5dunSRh4eH2rdvr/379xd6gAAAAABQ0phOpEJCQrRx40adO3dOcXFxat26tSTpzJkzcnNzK/QAAQAAAKCkMT39+fDhw9WrVy95enoqNDRUzZs3l3T5lL+wsLDCjg8AAAAAShzTidSTTz6pRo0a6dixY7r//vvl4HB5UKtKlSpcIwUAAADgtmA6kZKkBg0aqEGDBjZl7du3L5SAAAAAAKCkM51IPf744/kunzFjRoHbWr9+vSZPnqytW7cqPj5eCxcuVKdOnazL+/btq9mzZ9u8Jjo6WnFxcdbnp0+f1tChQ/Xdd9/JwcFBXbp00dtvv31T3yAYAAAAQMlmOpE6c+aMzfNLly5p165dSk5O1n333WeqrXPnzqlevXp6/PHH1blz51zrtGnTRjNnzrQ+d3V1tVneq1cvxcfHa8WKFbp06ZL69eungQMHau7cuaZiAQAAAICCMp1ILVy4MEdZVlaWBg8erKpVq5pqq23btmrbtm2+dVxdXRUYGJjrsr179youLk6//PKL9VTDd999V+3atdMbb7yh4OBgU/EAAAAAQEGYnv4810YcHDRixAhNmTKlMJqzsXbtWvn7+6tmzZoaPHiwTp06ZV22ceNG+fj42Fyv1apVKzk4OGjz5s15tpmWlqbU1FSbBwAAAAAUVKEkUpJ08OBBZWRkFFZzki6f1ve///1Pq1at0uuvv65169apbdu2yszMlCQlJCTI39/f5jVOTk7y9fVVQkJCnu3GxsbK29vb+ggJCSnUuAEAAADc2kyf2jdixAib54ZhKD4+Xt9//7369OlTaIFJUo8ePax/h4WFKTw8XFWrVtXatWvVsmXL62539OjRNtuRmppKMgUAAACgwEwnUr/99pvNcwcHB/n5+enNN9+85ox+N6pKlSoqV66cDhw4oJYtWyowMFBJSUk2dTIyMnT69Ok8r6uSLl93dfWkFQAAAABQUKYTqTVr1hRFHAXy999/69SpUwoKCpIkRUVFKTk5WVu3blVERIQkafXq1crKylJkZKTd4gQAAABwa7uuG/IWlrNnz+rAgQPW54cOHdL27dvl6+srX19fjR8/Xl26dFFgYKAOHjyoZ555RtWqVVN0dLQkqXbt2mrTpo0GDBig6dOn69KlSxoyZIh69OjBjH0AAAAAikyhTTZxPX799VfdeeeduvPOOyVdvv7qzjvv1JgxY+To6KgdO3bowQcfVI0aNdS/f39FREToxx9/tDktb86cOapVq5Zatmypdu3a6d5779VHH31kr00CAAAAcBuwGIZh2DsIe0tNTZW3t7dSUlLk5eVl73AAAGbtXGDvCHC9wrraOwIAsFHQ3MCuI1IAAAAAcDMikQIAAAAAk64rkRoyZIhOnz5d2LEAAAAAwE2hwInU33//bf177ty5Onv2rKTLN8o9duxY4UcGAAAAACVUgac/r1WrlsqWLavGjRvr4sWLOnbsmCpWrKjDhw/r0qVLRRkjAAAAAJQoBR6RSk5O1vz58xUREaGsrCy1a9dONWrUUFpampYtW6bExMSijBMAAAAASowCJ1KXLl1So0aNNHLkSLm7u+u3337TzJkz5ejoqBkzZqhy5cqqWbNmUcYKAAAAACVCgU/t8/HxUf369dW4cWOlp6frwoULaty4sZycnPTVV1+pfPny+uWXX4oyVgAAAAAoEQo8InX8+HG9+OKLcnV1VUZGhiIiItSkSROlp6dr27Ztslgsuvfee4syVgAAAAAoESyGYRhmX1SmTBmtX79ee/fuVe/evRUYGKjExEQ1atRI69atK4o4i1RB714MACihdi6wdwS4XmFd7R0BANgoaG5w3Tfk9fb2Vrdu3eTs7KzVq1fr0KFDevLJJ6+3OQAAAAC4aRT4Gqkr7dixQ+XLl5ckhYaGytnZWYGBgerevXuhBgcAAAAAJdF1JVIhISHWv3ft2lVowQAAAADAzeC6T+0DAAAAgNsViRQAAAAAmEQiBQAAAAAmkUgBAAAAgEkkUgAAAABgEokUAAAAAJhEIgUAAAAAJpFIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACY52TsAAMBNZucCe0cAAIDdMSIFAAAAACaRSAEAAACASSRSAAAAAGASiRQAAAAAmEQiBQAAAAAmkUgBAAAAgEkkUgAAAABgEokUAAAAAJhEIgUAAAAAJpFIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACaRSAEAAACASSRSAAAAAGASiRQAAAAAmEQiBQAAAAAmkUgBAAAAgEkkUgAAAABgEokUAAAAAJhEIgUAAAAAJpFIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACaRSAEAAACASXZNpNavX68OHTooODhYFotFixYtslluGIbGjBmjoKAgubu7q1WrVtq/f79NndOnT6tXr17y8vKSj4+P+vfvr7NnzxbjVgAAAAC43TjZc+Xnzp1TvXr19Pjjj6tz5845lk+aNEnvvPOOZs+ercqVK+ull15SdHS09uzZIzc3N0lSr169FB8frxUrVujSpUvq16+fBg4cqLlz5xb35gBAoVq8/bi9Q8hVR0d7RwAAgP1ZDMMw7B2EJFksFi1cuFCdOnWSdHk0Kjg4WCNHjtSoUaMkSSkpKQoICNCsWbPUo0cP7d27V3Xq1NEvv/yiBg0aSJLi4uLUrl07/f333woODi7QulNTU+Xt7a2UlBR5eXkVyfYBgFklN5HaaO8QcCsJ62rvCADARkFzgxJ7jdShQ4eUkJCgVq1aWcu8vb0VGRmpjRsvf4lv3LhRPj4+1iRKklq1aiUHBwdt3rw5z7bT0tKUmppq8wAAAACAgiqxiVRCQoIkKSAgwKY8ICDAuiwhIUH+/v42y52cnOTr62utk5vY2Fh5e3tbHyEhIYUcPQAAAIBbmV2vkbKX0aNHa8SIEdbnqampJFMAgEL165HT9g4hhwahvvYOAQBuGSV2RCowMFCSlJiYaFOemJhoXRYYGKikpCSb5RkZGTp9+rS1Tm5cXV3l5eVl8wAAAACAgiqxI1KVK1dWYGCgVq1apfr160u6PHK0efNmDR48WJIUFRWl5ORkbd26VREREZKk1atXKysrS5GRkfYKHQAA3MRK4kQvHeuXt3cIAK5i10Tq7NmzOnDggPX5oUOHtH37dvn6+qpixYoaPny4Jk6cqOrVq1unPw8ODrbO7Fe7dm21adNGAwYM0PTp03Xp0iUNGTJEPXr0KPCMfQAAAABgll0TqV9//VUtWrSwPs++bqlPnz6aNWuWnnnmGZ07d04DBw5UcnKy7r33XsXFxVnvISVJc+bM0ZAhQ9SyZUs5ODioS5cueuedd4p9WwAAAADcPuyaSDVv3lz53cbKYrFowoQJmjBhQp51fH19ufkuAAAAgGJVYq+RAgCUTMxGBwAAiRQAALeNEpkEh9k7AgC4PiRSAICbXklMEAAAt7YSex8pAAAAACipSKQAAAAAwCQSKQAAAAAwiUQKAAAAAEwikQIAAAAAk5i1DwAkLd5+3N4hAACAmwgjUgAAAABgEokUAAAAAJjEqX0AAAAlXEk8/bhj/fL2DgGwK0akAAAAAMAkRqQAAIDdlMSRFgAoCEakAAAAAMAkEikAAAAAMIlECgAAAABMIpECAAAAAJNIpAAAAADAJBIpAAAAADCJRAoAAAAATCKRAgAAAACTSKQAAAAAwCQSKQAAAAAwiUQKAAAAAEwikQIAAAAAk0ikAAAAAMAkEikAAAAAMIlECgAAAABMIpECAAAAAJNIpAAAAADAJBIpAAAAADCJRAoAAAAATCKRAgAAAACTSKQAAAAAwCQSKQAAAAAwiUQKAAAAAExysncAAG4/i7cft3cIAAAAN4QRKQAAAAAwiUQKAAAAAEwikQIAAAAAk0ikAAAAAMAkEikAAAAAMIlECgAAAABMIpECAAAAAJNIpAAAAADAJBIpAAAAADCJRAoAAAAATCKRAgAAAACTSKQAAAAAwCQneweQn3Hjxmn8+PE2ZTVr1tQff/whSbp48aJGjhypL7/8UmlpaYqOjtYHH3yggIAAe4QLAABMKv/39/YOIYfjFdrbOwQAN4ESPyJ1xx13KD4+3vrYsGGDddlTTz2l7777TvPnz9e6det04sQJde7c2Y7RAgAAALgdlOgRKUlycnJSYGBgjvKUlBR9+umnmjt3ru677z5J0syZM1W7dm1t2rRJd999d3GHCgAAAOA2UeJHpPbv36/g4GBVqVJFvXr10tGjRyVJW7du1aVLl9SqVStr3Vq1aqlixYrauHFjvm2mpaUpNTXV5gEAAAAABVWiR6QiIyM1a9Ys1axZU/Hx8Ro/fryaNGmiXbt2KSEhQS4uLvLx8bF5TUBAgBISEvJtNzY2Nse1VwBQ0pTEa0cAAMBlJTqRatu2rfXv8PBwRUZGKjQ0VPPmzZO7u/t1tzt69GiNGDHC+jw1NVUhISE3FCsAAACA20eJP7XvSj4+PqpRo4YOHDigwMBApaenKzk52aZOYmJirtdUXcnV1VVeXl42DwAAAAAoqBI9InW1s2fP6uDBg3rssccUEREhZ2dnrVq1Sl26dJEk7du3T0ePHlVUVJSdIwUAALi1Ld5+3N4h5NCxfnl7h4DbSIlOpEaNGqUOHTooNDRUJ06c0NixY+Xo6KiePXvK29tb/fv314gRI+Tr6ysvLy8NHTpUUVFRzNgHAAAAoEiV6ETq77//Vs+ePXXq1Cn5+fnp3nvv1aZNm+Tn5ydJmjJlihwcHNSlSxebG/ICAAAAQFGyGIZh2DsIe0tNTZW3t7dSUlK4XgooBiXxdJCSiFn7APs4XqG9vUPAdeLUPhSGguYGN9VkEwAAAABQEpBIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACaRSAEAAACASSRSAAAAAGASiRQAAAAAmEQiBQAAAAAmkUgBAAAAgEkkUgAAAABgEokUAAAAAJjkZO8AABStxduP2zsEAACAWw4jUgAAAABgEokUAAAAAJhEIgUAAAAAJpFIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACaRSAEAAACASSRSAAAAAGASiRQAAAAAmEQiBQAAAAAmkUgBAAAAgEkkUgAAAABgEokUAAAAAJjkZO8AAAAASpLyf39v7xByOF6hvb1DAHAVRqQAAAAAwCQSKQAAAAAwiUQKAAAAAEwikQIAAAAAk0ikAAAAAMAkZu0DAJXMWboAAEDJxYgUAAAAAJhEIgUAAAAAJpFIAQAAAIBJJFIAAAAAYBKJFAAAAACYRCIFAAAAACaRSAEAAACASSRSAAAAAGASN+QFAADALWHx9uP2DiGHjvXL2zsEFBFGpAAAAADAJBIpAAAAADCJRAoAAAAATCKRAgAAAACTSKQAAAAAwCRm7QMKUUmcLQgAcPMr//f39g4hh+MV2ts7BMCuSKQAFLuS+IMAAADAjFvm1L73339flSpVkpubmyIjI7VlyxZ7hwQAAADgFnVLJFJfffWVRowYobFjx2rbtm2qV6+eoqOjlZSUZO/QAAAAANyCbolE6q233tKAAQPUr18/1alTR9OnT5eHh4dmzJhh79AAAAAA3IJu+muk0tPTtXXrVo0ePdpa5uDgoFatWmnjxo25viYtLU1paWnW5ykpKZKk1NTUog0Wher7HSfsHQKu09nzF+wdAgDgBnn/ucDeIdwUvjgbbe8QcmgfHmzvEEq07JzAMIx86930idQ///yjzMxMBQQE2JQHBATojz/+yPU1sbGxGj9+fI7ykJCQIokRAAAAwM3l33//lbe3d57Lb/pE6nqMHj1aI0aMsD7PysrS6dOnVbZsWVksFjtGdjkDDgkJ0bFjx+Tl5WXXWG4X7PPixf4ufuzz4sc+L17s7+LHPi9+7PPiYxiG/v33XwUH5z9yd9MnUuXKlZOjo6MSExNtyhMTExUYGJjra1xdXeXq6mpT5uPjU1QhXhcvLy8+JMWMfV682N/Fj31e/NjnxYv9XfzY58WPfV488huJynbTTzbh4uKiiIgIrVq1ylqWlZWlVatWKSoqyo6RAQAAALhV3fQjUpI0YsQI9enTRw0aNFCjRo00depUnTt3Tv369bN3aAAAAABuQbdEItW9e3edPHlSY8aMUUJCgurXr6+4uLgcE1DcDFxdXTV27Ngcpx6i6LDPixf7u/ixz4sf+7x4sb+LH/u8+LHPSx6Lca15/QAAAAAANm76a6QAAAAAoLiRSAEAAACASSRSAAAAAGASiRQAAAAAmEQiVYK8//77qlSpktzc3BQZGaktW7bYO6RbRmxsrBo2bKjSpUvL399fnTp10r59+2zqNG/eXBaLxebxn//8x04R3/zGjRuXY3/WqlXLuvzixYuKiYlR2bJl5enpqS5duuS4sTYKrlKlSjn2t8ViUUxMjCT6d2FYv369OnTooODgYFksFi1atMhmuWEYGjNmjIKCguTu7q5WrVpp//79NnVOnz6tXr16ycvLSz4+Purfv7/Onj1bjFtxc8lvn1+6dEnPPvuswsLCVKpUKQUHB6t37946ceKETRu5fTZee+21Yt6Sm8O1+njfvn1z7Ms2bdrY1KGPm3OtfZ7bcd1isWjy5MnWOvRx+yGRKiG++uorjRgxQmPHjtW2bdtUr149RUdHKykpyd6h3RLWrVunmJgYbdq0SStWrNClS5fUunVrnTt3zqbegAEDFB8fb31MmjTJThHfGu644w6b/blhwwbrsqeeekrfffed5s+fr3Xr1unEiRPq3LmzHaO9uf3yyy82+3rFihWSpIcffthah/59Y86dO6d69erp/fffz3X5pEmT9M4772j69OnavHmzSpUqpejoaF28eNFap1evXtq9e7dWrFihJUuWaP369Ro4cGBxbcJNJ799fv78eW3btk0vvfSStm3bpm+++Ub79u3Tgw8+mKPuhAkTbPr+0KFDiyP8m861+rgktWnTxmZffvHFFzbL6ePmXGufX7mv4+PjNWPGDFksFnXp0sWmHn3cTgyUCI0aNTJiYmKszzMzM43g4GAjNjbWjlHdupKSkgxJxrp166xlzZo1M/773//aL6hbzNixY4169erluiw5OdlwdnY25s+fby3bu3evIcnYuHFjMUV4a/vvf/9rVK1a1cjKyjIMg/5d2CQZCxcutD7PysoyAgMDjcmTJ1vLkpOTDVdXV+OLL74wDMMw9uzZY0gyfvnlF2udpUuXGhaLxTh+/HixxX6zunqf52bLli2GJOPIkSPWstDQUGPKlClFG9wtKLf93adPH6Njx455voY+fmMK0sc7duxo3HfffTZl9HH7YUSqBEhPT9fWrVvVqlUra5mDg4NatWqljRs32jGyW1dKSookydfX16Z8zpw5KleunOrWravRo0fr/Pnz9gjvlrF//34FBwerSpUq6tWrl44ePSpJ2rp1qy5dumTT52vVqqWKFSvS5wtBenq6Pv/8cz3++OOyWCzWcvp30Tl06JASEhJs+rS3t7ciIyOtfXrjxo3y8fFRgwYNrHVatWolBwcHbd68udhjvhWlpKTIYrHIx8fHpvy1115T2bJldeedd2ry5MnKyMiwT4C3gLVr18rf3181a9bU4MGDderUKesy+njRSkxM1Pfff6/+/fvnWEYftw8newcA6Z9//lFmZqYCAgJsygMCAvTHH3/YKapbV1ZWloYPH67GjRurbt261vJHHnlEoaGhCg4O1o4dO/Tss89q3759+uabb+wY7c0rMjJSs2bNUs2aNRUfH6/x48erSZMm2rVrlxISEuTi4pLjx05AQIASEhLsE/AtZNGiRUpOTlbfvn2tZfTvopXdb3M7jmcvS0hIkL+/v81yJycn+fr60u8LwcWLF/Xss8+qZ8+e8vLyspYPGzZMd911l3x9ffXzzz9r9OjRio+P11tvvWXHaG9Obdq0UefOnVW5cmUdPHhQzz//vNq2bauNGzfK0dGRPl7EZs+erdKlS+c4DZ4+bj8kUrjtxMTEaNeuXTbX60iyOYc7LCxMQUFBatmypQ4ePKiqVasWd5g3vbZt21r/Dg8PV2RkpEJDQzVv3jy5u7vbMbJb36effqq2bdsqODjYWkb/xq3s0qVL6tatmwzD0LRp02yWjRgxwvp3eHi4XFxcNGjQIMXGxsrV1bW4Q72p9ejRw/p3WFiYwsPDVbVqVa1du1YtW7a0Y2S3hxkzZqhXr15yc3OzKaeP2w+n9pUA5cqVk6OjY44ZyxITExUYGGinqG5NQ4YM0ZIlS7RmzRpVqFAh37qRkZGSpAMHDhRHaLc8Hx8f1ahRQwcOHFBgYKDS09OVnJxsU4c+f+OOHDmilStX6oknnsi3Hv27cGX32/yO44GBgTkmEMrIyNDp06fp9zcgO4k6cuSIVqxYYTMalZvIyEhlZGTo8OHDxRPgLaxKlSoqV66c9ThCHy86P/74o/bt23fNY7tEHy9OJFIlgIuLiyIiIrRq1SprWVZWllatWqWoqCg7RnbrMAxDQ4YM0cKFC7V69WpVrlz5mq/Zvn27JCkoKKiIo7s9nD17VgcPHlRQUJAiIiLk7Oxs0+f37duno0eP0udv0MyZM+Xv76/27dvnW4/+XbgqV66swMBAmz6dmpqqzZs3W/t0VFSUkpOTtXXrVmud1atXKysry5rYwpzsJGr//v1auXKlypYte83XbN++XQ4ODjlOQYN5f//9t06dOmU9jtDHi86nn36qiIgI1atX75p16ePFh1P7SogRI0aoT58+atCggRo1aqSpU6fq3Llz6tevn71DuyXExMRo7ty5Wrx4sUqXLm09V9vb21vu7u46ePCg5s6dq3bt2qls2bLasWOHnnrqKTVt2lTh4eF2jv7mNGrUKHXo0EGhoaE6ceKExo4dK0dHR/Xs2VPe3t7q37+/RowYIV9fX3l5eWno0KGKiorS3Xffbe/Qb1pZWVmaOXOm+vTpIyen/zu8078Lx9mzZ21G8A4dOqTt27fL19dXFStW1PDhwzVx4kRVr15dlStX1ksvvaTg4GB16tRJklS7dm21adNGAwYM0PTp03Xp0iUNGTJEPXr0sDkNE/8nv30eFBSkrl27atu2bVqyZIkyMzOtx3ZfX1+5uLho48aN2rx5s1q0aKHSpUtr48aNeuqpp/Too4+qTJky9tqsEiu//e3r66vx48erS5cuCgwM1MGDB/XMM8+oWrVqio6OlkQfvx7XOq5Il/8pM3/+fL355ps5Xk8ftzN7TxuI//Puu+8aFStWNFxcXIxGjRoZmzZtsndItwxJuT5mzpxpGIZhHD161GjatKnh6+truLq6GtWqVTOefvppIyUlxb6B38S6d+9uBAUFGS4uLkb58uWN7t27GwcOHLAuv3DhgvHkk08aZcqUMTw8PIyHHnrIiI+Pt2PEN79ly5YZkox9+/bZlNO/C8eaNWtyPY706dPHMIzLU6C/9NJLRkBAgOHq6mq0bNkyx3tx6tQpo2fPnoanp6fh5eVl9OvXz/j333/tsDU3h/z2+aFDh/I8tq9Zs8YwDMPYunWrERkZaXh7extubm5G7dq1jVdffdW4ePGifTeshMpvf58/f95o3bq14efnZzg7OxuhoaHGgAEDjISEBJs26OPmXOu4YhiG8eGHHxru7u5GcnJyjtfTx+3LYhiGUeTZGgAAAADcQrhGCgAAAABMIpECAAAAAJNIpAAAAADAJBIpAAAAADCJRAoAAAAATCKRAgAAAACTSKQAAAAAwCQSKQAAAAAwiUQKAGBXFotFixYtuuF2KlWqpKlTp95wOwAAFASJFABACQkJGjp0qKpUqSJXV1eFhISoQ4cOWrVqVZGvOz4+Xm3btr3hdn755RcNHDiwECLK29q1a2WxWKwPPz8/tWvXTjt37izS9QIASh4SKQC4zR0+fFgRERFavXq1Jk+erJ07dyouLk4tWrRQTExMka8/MDBQrq6uN9yOn5+fPDw8CiGia9u3b5/i4+O1bNkypaWlqX379kpPTy+WdQMASgYSKQC4zT355JOyWCzasmWLunTpoho1auiOO+7QiBEjtGnTJmu9o0ePqmPHjvL09JSXl5e6deumxMRE6/Jx48apfv36mjFjhipWrChPT089+eSTyszM1KRJkxQYGCh/f3+98sorNuu/8tS+9PR0DRkyREFBQXJzc1NoaKhiY2MlSYZhaNy4capYsaJcXV0VHBysYcOGWdu5+tS+gsb72WefqVKlSvL29laPHj3077//XnOf+fv7KzAwUHfddZeGDx+uY8eO6Y8//rAu37Bhg5o0aSJ3d3eFhIRo2LBhOnfunCTp+eefV2RkZI4269WrpwkTJliff/LJJ6pdu7bc3NxUq1YtffDBB9Zlhw8flsVi0TfffKMWLVrIw8ND9erV08aNG3Ns35WmTp2qSpUq2ZTltx4AQN5IpADgNnb69GnFxcUpJiZGpUqVyrHcx8dHkpSVlaWOHTvq9OnTWrdunVasWKG//vpL3bt3t6l/8OBBLV26VHFxcfriiy/06aefqn379vr777+1bt06vf7663rxxRe1efPmXON555139O2332revHnat2+f5syZY/3h//XXX2vKlCn68MMPtX//fi1atEhhYWG5tmMm3kWLFmnJkiVasmSJ1q1bp9dee63A+y8lJUVffvmlJMnFxcXaZps2bdSlSxft2LFDX331lTZs2KAhQ4ZIknr16qUtW7bo4MGD1nZ2796tHTt26JFHHpEkzZkzR2PGjNErr7yivXv36tVXX9VLL72k2bNn26z/hRde0KhRo7R9+3bVqFFDPXv2VEZGRoHjL+h6AAC5MAAAt63Nmzcbkoxvvvkm33rLly83HB0djaNHj1rLdu/ebUgytmzZYhiGYYwdO9bw8PAwUlNTrXWio6ONSpUqGZmZmdaymjVrGrGxsdbnkoyFCxcahmEYQ4cONe677z4jKysrRwxvvvmmUaNGDSM9PT3XGENDQ40pU6bcULxPP/20ERkZmed+WLNmjSHJKFWqlFGqVClDkiHJePDBB611+vfvbwwcONDmdT/++KPh4OBgXLhwwTAMw6hXr54xYcIE6/LRo0fbrLdq1arG3Llzbdp4+eWXjaioKMMwDOPQoUOGJOOTTz7JsX179+61bl+9evVs2pgyZYoRGhpa4PUAAPLGiBQA3MYMwyhQvb179yokJEQhISHWsjp16sjHx0d79+61llWqVEmlS5e2Pg8ICFCdOnXk4OBgU5aUlJTrevr27avt27erZs2aGjZsmJYvX25d9vDDD+vChQuqUqWKBgwYoIULF+Y5+nK98QYFBeUZ25V+/PFHbd26VbNmzVKNGjU0ffp067Lff/9ds2bNkqenp/URHR2trKwsHTp0SNLlUam5c+dKuvwefPHFF+rVq5ck6dy5czp48KD69+9v08bEiRNtRrEkKTw83CZ2SQWK3+x6AAA5Odk7AACA/VSvXl0Wi8Xm+p4b4ezsbPPcYrHkWpaVlZXr6++66y4dOnRIS5cu1cqVK9WtWze1atVKCxYsUEhIiPbt26eVK1dqxYoVevLJJzV58mStW7cuxzpuJN68YrtS5cqV5ePjo5o1ayopKUndu3fX+vXrJUlnz57VoEGDbK7fylaxYkVJUs+ePfXss89q27ZtunDhgo4dO2Y97fDs2bOSpI8//jjHtVSOjo55xm+xWCTJGr+Dg0OORPnSpUvWv82sBwCQE4kUANzGfH19FR0drffff1/Dhg3LcZ1UcnKyfHx8VLt2bR07dkzHjh2zjvLs2bNHycnJqlOnTqHG5OXlpe7du6t79+7q2rWr2rRpo9OnT8vX11fu7u7q0KGDOnTooJiYGNWqVUs7d+7UXXfdZdNGccYbExOj2NhYLVy4UA899JDuuusu7dmzR9WqVcvzNRUqVFCzZs00Z84cXbhwQffff7/8/f0lXR6xCw4O1l9//WUdpboefn5+SkhIkGEY1iRr+/bt1uWFtR4AuF2RSAHAbe79999X48aN1ahRI02YMEHh4eHKyMjQihUrNG3aNO3du1etWrVSWFiYevXqpalTpyojI0NPPvmkmjVrpgYNGhRaLG+99ZaCgoJ05513ysHBQfPnz1dgYKB8fHw0a9YsZWZmKjIyUh4eHvr888/l7u6u0NDQHO0UV7yS5OHhoQEDBmjs2LHq1KmTnn32Wd19990aMmSInnjiCZUqVUp79uzRihUr9N5771lf16tXL40dO1bp6emaMmWKTZvjx4/XsGHD5O3trTZt2igtLU2//vqrzpw5oxEjRhQorubNm+vkyZOaNGmSunbtqri4OC1dulReXl6Fuh4AuF1xjRQA3OaqVKmibdu2qUWLFho5cqTq1q2r+++/X6tWrdK0adMkXT5tbPHixSpTpoyaNm2qVq1aqUqVKvrqq68KNZbSpUtr0qRJatCggRo2bKjDhw/rhx9+kIODg3x8fPTxxx+rcePGCg8P18qVK/Xdd9+pbNmyOdoprnizDRkyRHv37tX8+fMVHh6udevW6c8//1STJk105513asyYMQoODrZ5TdeuXXXq1CmdP39enTp1sln2xBNP6JNPPtHMmTMVFhamZs2aadasWapcuXKBY6pdu7Y++OADvf/++6pXr562bNmiUaNGFfp6AOB2ZTEKeqUxAAAAAEASI1IAAAAAYBqJFAAAAACYRCIFAAAAACaRSAEAAACASSRSAAAAAGASiRQAAAAAmEQiBQAAAAAmkUgBAAAAgEkkUgAAAABgEokUAAAAAJhEIgUAAAAAJv0/hX16Zb7DTZMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "\n",
    "plt.hist(exp_data[exp_data['group'] == 'test']['commission_revenue_exp'], \n",
    "         bins=np.arange(0.0, 200.0, 10.0), \n",
    "         alpha=0.35,\n",
    "         label='Initial')\n",
    "\n",
    "plt.hist(exp_data[exp_data['group'] == 'test']['commission_revenue_exp_cuped'], \n",
    "         bins=np.arange(0.0, 200.0, 10.0), \n",
    "         alpha=0.35,\n",
    "         label='CUPED transformed')\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "plt.xlabel('Comission Revenue')\n",
    "plt.ylabel('# users')\n",
    "\n",
    "plt.title('Distribution of initial and transformed Commission Revenue in the test group.')\n",
    "\n",
    "plt.savefig('./output/distributions.png', bbox_inches='tight', dpi=100)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "237de3a1-de74-46d1-a883-95134f5cc714",
   "metadata": {},
   "source": [
    "We are sure the transformed metric is an unbiased estimator empirically and theoretically. We have seen that variance indeed decreases with such transformation. It is time to rerun our A/B test using the transformed metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b279c09c-42c1-43f0-9694-da2788bc8d21",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_indResult(statistic=3.7257956767943394, pvalue=0.00020010540422191647)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.ttest_ind(\n",
    "    exp_data[exp_data['group'] == 'test']['commission_revenue_exp_cuped'], \n",
    "    exp_data[exp_data['group'] == 'control']['commission_revenue_exp_cuped']\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b7b677e-0cec-4028-b8bf-46b6747b5f1b",
   "metadata": {},
   "source": [
    "The p-value is $0.0002$, lower than $0.025$ by at least one order! Finally, we've got a statistically significant result! We now have statistical and data-driven proof that the engine's B recommendations increase commission revenue. Based on the results, we can recommend rolling out new functionality to 100% of users to increase the metric by two 2% of the full scale of a company."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39568d7a-fc53-4f85-97b7-42dbf4dbd9cc",
   "metadata": {},
   "source": [
    "### Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f908be26-d6d7-47bb-8a44-07c6c08a600d",
   "metadata": {},
   "source": [
    "I hope now it is clear what CUPED is and how it can enhance your A/B experiments in production. In the article, we went through the theoretical foundations of CUPED and had hands-on experience applying it to a real-life A/B test.\n",
    "\n",
    "If you have any questions, do not hesitate to reach out on [LinkedIn](https://www.linkedin.com/in/igorkhomyanin/) :)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf7fa46f-c379-4d6b-9181-6f79638a4645",
   "metadata": {},
   "source": [
    "### References"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cf04ed2-87a2-4859-ad7b-6ddc319faf8b",
   "metadata": {},
   "source": [
    "1. Kohavi, Ron & Deng, Alex & Xu, Ya & Walker, Toby. (2013). Improving the Sensitivity of Online Controlled Experiments by Utilizing Pre-Experiment Data. 10.1145/2433396.2433413. \n",
    "2. https://en.wikipedia.org/wiki/Student%27s_t-test\n",
    "3. https://en.wikipedia.org/wiki/Central_limit_theorem\n",
    "4. \"How Booking.com increases the power of online experiments with CUPED\" by Simon Jackson at https://booking.ai/how-booking-com-increases-the-power-of-online-experiments-with-cuped-995d186fff1d "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.11.0 (ds_env)",
   "language": "python",
   "name": "ds_env"
  },
  "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}