willy_nilly.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $ \\mathcal{T}_{X} = \\{X_{i}, Y_{i}\\}_{i = 1}^{n_{1}} \\stackrel{iid}{\\sim} F_{X,Y} $ and $ \\mathcal{T}_{Z} = \\{Z_{j}, W_{j}\\}_{j=1}^{n_{2}} \\stackrel{iid}{\\sim} F_{Z,W} $. Suppose $ X_{i}, Z_{j} \\in \\mathcal{X} $ and $ Y_{i}, W_{j} \\in \\{0, 1\\} $. For the sake of simplicity let $ n_{1} = n_{2} = n $, and $ \\mathcal{X} = \\mathbb{R} $. We refer to $ \\mathcal{T}_{X} $ and $ \\mathcal{T}_{Z} $ as training data. The goal of any classification procedure is to predict the class label of an unlabeled observation in the input space, $ X_{0} \\in \\mathcal{X} $. Though we do not introduce the necessary notational machinery, it is assumed that any classifier knows what distribution, $ F_{X, Y} $ or $ F_{Z, W} $, a particular test observation came from. In this write up we make the distinction via different letters (X and Z). \n",
    "\n",
    "In this notebook we investigate two types of classifiers. The first type ignores potentially useful information in $ \\mathcal{T}_{Z} $ when classifying $ X_{0} \\sim F_{X,Y} $ (and the same is true when classifying $ Z_{0} \\sim F_{Z, W} $). The second type attempts to use information in both sets of training data. In particular, consider $ \\hat{\\eta}_{X}(\\cdot) $ and $ \\hat{\\eta}_{Z}(\\cdot) $, estimates of $ P_{Y=1|X}(Y=1|X) $ and $ P_{W=1|Z}(W=1|Z) $, respectively. \n",
    "\n",
    "The classifier we consider that is of the first type is simply Bayes', or, \n",
    "$$ \\hat{Y}_{0} = g_{X}(X_{0} | \\mathcal{T}_{X}, \\mathcal{T}_{Z}) = \\begin{cases} 1 \\mbox{ if $ \\hat{\\eta}_{X}(X_{0}) > \\frac{1}{2} $} \\\\ 0 \\mbox{ o/w} \\end{cases} $$ where $ g_{Z}(Z_{0} | \\mathcal{T}_{X}, \\mathcal{T}_{Z}) $ is defined analogously. \n",
    "\n",
    "We consider two classifiers of the second type. The first can be described as follows: Assume we are classifying a test observation from $ F_{X, Y} $. Estimate the posterior associated with $ X $, $ \\hat{\\eta}_{X} $. Then, using $ \\mathcal{T}_{Z} $, partition $ \\mathcal{X} $. Then use $ \\mathcal{T}_{X} $ to estimate the posterior associated with $ X $ -- simply use the propotion of points in $ \\mathcal{T}_{X} $ such that $ X_{i} $ is in the partition and $ Y_{i} = 1 $. Let $ \\hat{\\eta}_{Z,X} $ denote the posterior estimated from using the structure $ \\mathcal{T}_{Z} $ for observations from $ F_{X, Y} $. A natural classifier from this procedure is given by \n",
    "$$ \\hat{Y}_{0} = g_{X, L2M}(X_{0} | \\mathcal{T}_{X}, \\mathcal{T}_{Z}) = \\begin{cases} 1 \\mbox{ if $ \\hat{\\eta}_{X}(X_{0}) + \\hat{\\eta}_{X,Z}(X_{0}) > 1 $} \\\\ 0 \\mbox{ o/w} \\end{cases} $$ \n",
    "\n",
    "The second classifier uses $ \\mathcal{T}_{Z} $ a bit differently. Let $$ \\hat{\\alpha}_{X} = \\frac{\\sum_{i=1}^{n} \\left[ \\mathbb{1}\\{\\hat{\\eta}_{X}(X_{i}) > 0.5\\} \\mathbb{1}\\{\\hat{\\eta}_{Z}(X_{i}) > 0.5\\} + \\mathbb{1}\\{\\hat{\\eta}_{X}(X_{i}) < 0.5\\} \\mathbb{1}\\{\\hat{\\eta}_{Z}(X_{i}) < 0.5\\} \\right]}{n} $$\n",
    "\n",
    "Then consider the classifier $$ g_{X, cep}(X_{0} | \\mathcal{T}_{X}, \\mathcal{T}_{Z}) = \\begin{cases} 1 \\mbox{ if $ \\hat{\\eta}_{X}(X_{0}) + \\hat{\\alpha}_{X} \\hat{\\eta}_{Z}(X_{0}) > \\frac{(1 + \\hat{\\alpha}_{X})}{2} $} \\\\ 0 \\mbox{ otherwise} \\end{cases} $$\n",
    "\n",
    "Though i do not write it out, it is possible that $ F_{X,Y} $ and $ F_{W,Z} $ are such that $ \\eta_{X}(X) = 1 - \\eta_{Z}(X) $ for all $ X $. In such a case, there is a lot of information in $ \\mathcal{T}_{Z} $ about $ F_{X, Y} $ that $ \\hat{\\alpha}_{X} $ is not able to take advantage of. To this point, one should permute the labels, $ W $, such that the sum is maximized. This is simply a label alignment. Note that the other two classifiers do not require this check.\n",
    "\n",
    "The first experiment to compare the classifiers is as follows:\n",
    "Fix $ n = n_{1} = n_{2} $. Let $ f(X|Y=0) = f(Z|Y=1) = Uniform(0,1) $ and $ f(X|Y=1) = f(X|Y=0) = Uniform(1,2) $. Then $ L^{*}_{X} = L^{*}_{Z} = 0 $. \n",
    "\n",
    "To estimate posteriors, we use $ k = log(n) $ random partitions of $ \\mathcal{X} $. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_sample(n, pi, conditional_0, params0, conditional_1, params1, acorn=None):\n",
    "    if acorn is None:\n",
    "        acorn = np.random.seed(10**6)\n",
    "    np.random.seed(acorn)\n",
    "    \n",
    "    n0 = int(np.random.binomial(n, pi))\n",
    "    n1 = n - n0\n",
    "    ns = [n0, n1]\n",
    "    \n",
    "    X0 = conditional_0(*params0, size=n0)\n",
    "    X0 = np.sort(X0[::])\n",
    "    \n",
    "    X1 = conditional_1(*params1, size=n1)\n",
    "    X1 - np.sort(X1[::])\n",
    "    \n",
    "    labels = np.concatenate([i*np.ones(ns[i]) for i in range(len(ns))])\n",
    "    \n",
    "    return [[X0, X1], labels]\n",
    "\n",
    "class LinearClassifier:\n",
    "    def __init__(self, boundaries, boundary_left=-10**6, boundary_right=10**6):\n",
    "        boundaries = np.sort(boundaries)[::]\n",
    "        boundaries = np.concatenate(([boundary_left], boundaries, [boundary_right]))\n",
    "        self.boundaries = boundaries\n",
    "    \n",
    "    def fit(self, X, y):\n",
    "        posteriors = np.zeros(len(self.boundaries) - 1)\n",
    "        \n",
    "        unique_labels, label_counts = np.unique(y, return_counts=True)\n",
    "        \n",
    "        for i in range(len(self.boundaries) - 1):\n",
    "            temp_idx = ((X > self.boundaries[i]).astype(int) + (X < self.boundaries[i+1]).astype(int)) == 2\n",
    "            temp = X[temp_idx]\n",
    "            temp_labels = y[temp_idx]\n",
    "            \n",
    "            unique_, counts_ = np.unique(temp_labels, return_counts=True)\n",
    "            \n",
    "            if len(unique_) == 0:\n",
    "                posteriors[i] = label_counts[-1] / np.sum(label_counts)\n",
    "            elif len(unique_) == 2:\n",
    "                posteriors[i] = counts_[-1] / np.sum(counts_)\n",
    "            elif unique_[0] == 0:\n",
    "                posteriors[i] = 0\n",
    "            else:\n",
    "                posteriors[i] = 1\n",
    "\n",
    "        self.posteriors = posteriors\n",
    "        \n",
    "    def predict(self, X, return_posteriors=False):\n",
    "        posteriors = -1*np.ones(len(X))\n",
    "        partitions = np.sum([np.array(X > c_) for c_ in self.boundaries], axis=0)\n",
    "        for i, c_ in enumerate(self.boundaries[1:]):\n",
    "            posteriors[np.where(partitions == i+1)] = self.posteriors[i]\n",
    "\n",
    "        posteriors[posteriors == -1] = self.posteriors[-1]\n",
    "        \n",
    "        predictions = -1*np.ones(len(X))\n",
    "        predictions[posteriors < 0.5], predictions[posteriors >= 0.5] = 0, 1\n",
    "        predictions[posteriors == 0.5] = np.random.binomial(1, 0.5, size=np.sum(posteriors == 0.5))\n",
    "        \n",
    "        if return_posteriors:\n",
    "            return predictions, posteriors\n",
    "        else:\n",
    "            return predictions\n",
    "\n",
    "def estimate_alpha(train1, train2, labels1, labels2, splits1, splits2):\n",
    "    \n",
    "    class_11 = posterior1[posterior1 > 0.5].astype(int)\n",
    "    class_12 = posterior2[posterior2 > 0.5].astype(int)\n",
    "    \n",
    "    alpha = np.sum(class_11 + class_12 != 1)/1\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 15/15 [01:30<00:00, 10.13s/it]\n"
     ]
    }
   ],
   "source": [
    "mc_its = 100 # number of simulation repititions\n",
    "ns0 = 10*np.arange(1,6,step=1)\n",
    "ns1 = 100*np.arange(0.5, 5.5, step=0.5) # number of training samples \n",
    "ns = np.concatenate((ns0, ns1))\n",
    "m = 1000 # number of test samples each monte carlo iteration\n",
    "\n",
    "pi = 0.5\n",
    "\n",
    "cond_X_0, cond_X_1 = 2*[np.random.uniform]\n",
    "params_X_0, params_X_1 = [0, 0.5], [0.5, 1]\n",
    "\n",
    "cond_Z_0, cond_Z_1 = 2*[np.random.uniform]\n",
    "params_Z_0, params_Z_1 = [0.5, 1], [0, 0.5]\n",
    "\n",
    "mean_errors_X = np.zeros(len(ns))\n",
    "std_errors_X = np.zeros(len(ns))\n",
    "\n",
    "mean_errors_Z = np.zeros(len(ns))\n",
    "std_errors_Z = np.zeros(len(ns))\n",
    "\n",
    "mean_errors_HHX = np.zeros(len(ns))\n",
    "std_errors_HHX = np.zeros(len(ns))\n",
    "\n",
    "mean_errors_HHZ = np.zeros(len(ns))\n",
    "std_errors_HHZ = np.zeros(len(ns))\n",
    "\n",
    "mean_errors_joint = np.zeros(len(ns))\n",
    "std_errors_joint = np.zeros(len(ns))\n",
    "\n",
    "mean_errors_LLX = np.zeros(len(ns))\n",
    "std_errors_LLX = np.zeros(len(ns))\n",
    "\n",
    "mean_errors_LLZ = np.zeros(len(ns))\n",
    "std_errors_LLZ = np.zeros(len(ns))\n",
    "\n",
    "for j, n in enumerate(tqdm(ns)):\n",
    "    n = int(n)\n",
    "    k = int(np.floor(np.log(n)))\n",
    "    \n",
    "    T = int(np.floor(np.sqrt(n)))\n",
    "    \n",
    "    b = int(np.floor(0.627*n))\n",
    "    \n",
    "    error_X = []\n",
    "    error_Z = []\n",
    "    error_HHX = []\n",
    "    error_HHZ = []\n",
    "    error_joint = []\n",
    "    error_LLX = []\n",
    "    error_LLZ = []\n",
    "    \n",
    "    for i in range(mc_its):\n",
    "        in_bag_idx_X = np.random.choice(range(n), size=(b,T), replace=True)\n",
    "        in_bag_idx_Z = np.random.choice(range(n), size=(b,T), replace=True)\n",
    "        \n",
    "        #- Initialize predictions\n",
    "        predictionsXX = np.zeros(m)\n",
    "        predictionsZX = np.zeros(m)\n",
    "        \n",
    "        predictionsZZ = np.zeros(m)\n",
    "        predictionsXZ = np.zeros(m)\n",
    "        \n",
    "        predictionsHHX = np.zeros(m)\n",
    "        predictionsHHZ = np.zeros(m)\n",
    "        \n",
    "        predictionsjoint = np.zeros(2*m)\n",
    "        \n",
    "        predictionsLLX = np.zeros(m)\n",
    "        predictionsLLZ = np.zeros(m)\n",
    "        \n",
    "        #- Initialize posteriors\n",
    "        posteriorsXX = np.zeros(m)\n",
    "        posteriorsXZ = np.zeros(m)\n",
    "        \n",
    "        posteriorsZZ = np.zeros(m)\n",
    "        posteriorsZX = np.zeros(m)\n",
    "        \n",
    "        X, labelsX = generate_sample(n, pi, cond_X_0, params_X_0, cond_X_1, params_X_1)\n",
    "        X = np.concatenate(X)\n",
    "        testX, test_labelsX = generate_sample(m, pi, cond_X_0, params_X_0, cond_X_1, params_X_1)\n",
    "        testX = np.concatenate(testX)\n",
    "        \n",
    "        Z, labelsZ = generate_sample(n, pi, cond_Z_0, params_Z_0, cond_Z_1, params_Z_1)\n",
    "        Z = np.concatenate(Z)\n",
    "        testZ, test_labelsZ = generate_sample(m, pi, cond_Z_0, params_Z_0, cond_Z_1, params_Z_1)\n",
    "        testZ = np.concatenate(testZ)\n",
    "        \n",
    "        for _ in range(T):\n",
    "            temp_X, temp_labels_X = X[in_bag_idx_X[_]], labelsX[in_bag_idx_X[_]]\n",
    "            temp_Z, temp_labels_Z = Z[in_bag_idx_Z[_]], labelsZ[in_bag_idx_Z[_]]\n",
    "            \n",
    "            joint, joint_labels = np.concatenate((temp_X, temp_Z)), np.concatenate((temp_labels_X, temp_labels_Z))\n",
    "            joint_test = np.concatenate((testX, testZ), axis = 0)\n",
    "            \n",
    "            cX = np.random.uniform(min(temp_X), max(temp_X), size=k)\n",
    "            cZ = np.random.uniform(min(temp_Z), max(temp_Z), size=k)\n",
    "            \n",
    "            lc_X = LinearClassifier(cX)\n",
    "            lc_X.fit(temp_X, temp_labels_X)\n",
    "            Y_hat_XX, posteriors_XX = lc_X.predict(testX, return_posteriors=True)\n",
    "            Y_hat_XZ, posteriors_XZ = lc_X.predict(testZ, return_posteriors=True)\n",
    "            \n",
    "            lc_LLX = LinearClassifier(cZ)\n",
    "            lc_LLX.fit(temp_X, temp_labels_X)\n",
    "            Y_hat_LLX, posteriors_LLX = lc_LLX.predict(testX, return_posteriors=True)\n",
    "            \n",
    "            predictionsXX += Y_hat_XX\n",
    "            posteriorsXX += posteriors_XX\n",
    "            \n",
    "            predictionsLLX += Y_hat_LLX + Y_hat_XX\n",
    "            \n",
    "            lc_Z = LinearClassifier(cZ)\n",
    "            lc_Z.fit(temp_Z, temp_labels_Z)\n",
    "            Y_hat_ZZ, posteriors_ZZ = lc_Z.predict(testZ, return_posteriors=True)\n",
    "            Y_hat_ZX, posteriors_ZX = lc_Z.predict(testX, return_posteriors=True)\n",
    "            \n",
    "            lc_LLZ = LinearClassifier(cX)\n",
    "            lc_LLZ.fit(temp_Z, temp_labels_Z)\n",
    "            Y_hat_LLZ, posteriors_LLZ = lc_LLZ.predict(testZ, return_posteriors=True)\n",
    "            \n",
    "            predictionsZZ += Y_hat_ZZ\n",
    "            posteriorsZZ += posteriors_ZZ\n",
    "            \n",
    "            predictionsLLZ += Y_hat_LLZ + Y_hat_ZZ\n",
    "            \n",
    "            predictionsZX += Y_hat_ZX\n",
    "            predictionsXZ += Y_hat_XZ\n",
    "            \n",
    "            posteriorsZX += posteriors_ZX\n",
    "            posteriorsXZ += posteriors_XZ\n",
    "            \n",
    "            cjoint = np.random.uniform(min(joint), max(joint), size=2*k)\n",
    "            \n",
    "            lc_joint = LinearClassifier(cjoint)\n",
    "            lc_joint.fit(joint, joint_labels)\n",
    "            predictionsjoint += lc_joint.predict(joint_test)\n",
    "            \n",
    "        # average\n",
    "        predictionsXX, predictionsZZ = predictionsXX/T, predictionsZZ/T\n",
    "        predictionsZX, predictionsXZ = predictionsZX/T, predictionsXZ/T\n",
    "        predictionsjoint = predictionsjoint / T\n",
    "        \n",
    "        predictionsLLX = predictionsLLX / (2 * T)\n",
    "        predictionsLLZ = predictionsLLZ / (2 * T)\n",
    "        \n",
    "        posteriorsXX, posteriorsZZ = posteriorsXX/T, posteriorsZZ/T\n",
    "        posteriorsZX, posteriorsXZ = posteriorsZX/T, posteriorsXZ/T\n",
    "        \n",
    "        # threshold\n",
    "        predictionsXX[predictionsXX < 0.5], predictionsXX[predictionsXX >= 0.5] = 0, 1\n",
    "        predictionsZZ[predictionsZZ < 0.5], predictionsZZ[predictionsZZ >= 0.5] = 0, 1\n",
    "        \n",
    "        predictionsZX[predictionsZX < 0.5], predictionsZX[predictionsZX >= 0.5] = 0, 1\n",
    "        predictionsXZ[predictionsXZ < 0.5], predictionsXZ[predictionsXZ >= 0.5] = 0, 1\n",
    "        \n",
    "        predictionsLLX[predictionsLLX < 0.5], predictionsLLX[predictionsLLX >= 0.5] = 0, 1\n",
    "        predictionsLLZ[predictionsLLZ < 0.5], predictionsLLZ[predictionsLLZ >= 0.5] = 0, 1\n",
    "        \n",
    "        predictionsjoint[predictionsjoint < 0.5], predictionsjoint[predictionsjoint >= 0.5] = 0, 1\n",
    "        \n",
    "        alpha_X = max(np.sum(predictionsXX == predictionsZX), np.sum(predictionsXX == -1*(predictionsZX - 1))) / m\n",
    "        alpha_Z = max(np.sum(predictionsZZ == predictionsXZ), np.sum(predictionsZZ == -1*(predictionsXZ - 1))) / m\n",
    "        \n",
    "        if np.argmax([np.sum(predictionsXX == predictionsZX), np.sum(predictionsXX == -1*(predictionsZX - 1))]) == 1:\n",
    "            flip = True\n",
    "        else:\n",
    "            flip = False\n",
    "        \n",
    "        if flip:\n",
    "            posteriorsHHX = posteriorsXX + alpha_X*(1 - posteriorsZX)\n",
    "            posteriorsHHZ = posteriorsZZ + alpha_Z*(1 - posteriorsXZ)\n",
    "        else:\n",
    "            posteriorsHHX = posteriorsXX + alpha_X*posteriorsZX\n",
    "            posteriorsHHZ = posteriorsZZ + alpha_Z*posteriorsXZ\n",
    "        \n",
    "        predictionsHHX[posteriorsHHX < (alpha_X + 1)*0.5], predictionsHHX[posteriorsHHX >= (alpha_X + 1)*0.5] = 0, 1\n",
    "        predictionsHHZ[posteriorsHHZ < (alpha_Z + 1)*0.5], predictionsHHZ[posteriorsHHZ >= (alpha_Z + 1)*0.5] = 0, 1\n",
    "        \n",
    "        error_X.append(1 - np.sum(predictionsXX == test_labelsX)/m)\n",
    "        error_Z.append(1 - np.sum(predictionsZZ == test_labelsZ)/m)\n",
    "        error_HHX.append(1 - np.sum(predictionsHHX == test_labelsX)/m)\n",
    "        error_HHZ.append(1 - np.sum(predictionsHHZ == test_labelsZ)/m)\n",
    "        error_joint.append(1 - np.sum(predictionsjoint == np.concatenate((test_labelsX, test_labelsZ)))/(2*m))\n",
    "        error_LLX.append(1 - np.sum(predictionsLLX == test_labelsX)/m)\n",
    "        error_LLZ.append(1 - np.sum(predictionsLLZ == test_labelsZ)/m)\n",
    "        \n",
    "        \n",
    "    mean_errors_X[j] = np.mean(error_X)\n",
    "    std_errors_X[j] = np.std(error_X, ddof=1)\n",
    "    \n",
    "    mean_errors_Z[j] = np.mean(error_Z)\n",
    "    std_errors_Z[j] = np.std(error_Z, ddof=1)\n",
    "    \n",
    "    mean_errors_HHX[j] = np.mean(error_HHX)\n",
    "    std_errors_HHX[j] = np.std(error_HHX, ddof=1)\n",
    "    \n",
    "    mean_errors_HHZ[j] = np.mean(error_HHZ)\n",
    "    std_errors_HHZ[j] = np.std(error_HHZ, ddof=1)\n",
    "    \n",
    "    mean_errors_joint[j] = np.mean(error_joint)\n",
    "    std_errors_joint[j] = np.std(error_joint, ddof=1)\n",
    "    \n",
    "    mean_errors_LLX[j] = np.mean(error_LLX)\n",
    "    std_errors_LLX[j] = np.std(error_LLX, ddof=1)\n",
    "    \n",
    "    mean_errors_LLZ[j] = np.mean(error_LLZ)\n",
    "    std_errors_LLZ[j] = np.std(error_LLZ, ddof=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8FfW5+PHPk40AIaxBEhJZlEXZYqSARa27uAC1tbcgVEvt5VeXutwq2spFQG3x1rrdi721WkVblVZLBa1arSjqVTFAkEVACCJZkASSQEiALM/vj5mkh5jlkJyZyfK8X6+8cmY732eGcJ7z/c7MM6KqGGOMMU2JCjoAY4wxbYMlDGOMMWGxhGGMMSYsljCMMcaExRKGMcaYsFjCMMYYExZLGB2AiJSKyGD39dMicq/7+hwRyfGgvbNEZGvI9BcicoH7er6I/DHSbUZC3bjrWX6ieyyj/YyrNRCRd0Tkx0HHUZ/WHFt7YwmjA1DVBFXN9rG991R1mF/tNZeIqIicXDNdN+7QROcu/9I9llV+x+rG00lEnhSRXSJyUESyROSSRtb/oYi872eMx0NEHhCRz9192SIiVwcdk2mcJQzT4YhITNAxNFMMsBv4FtAdmAv8WUQGBhhTSxwCJuPsyzXAIyLyzWBDMo2xhNFGicgsEVkRMv25iPwlZHq3iKS7r4/5Jt3A+90uIi/VmfeoiDzSwPpfiMhtIvKpiJSIyFIRiXeXhTXUJSKvishP68z7VESuqGfdge5+zBaRPBHJF5HbQpaPE5EPRaTYXfY/IhIXslxF5AYR+Rz4XERWuYvWu8NM3w+NW0SeBU4EVrjL54TEEOOukyIiy0Vkv4hsF5F/D2lvvoj8WUSecb9BbxKRsSHL7xCRXHfZVhE5v6njpaqHVHW+qn6hqtWq+gqwEzi9nuN1CvC/wBlu/MXu/MtEZJ2IHHD/RuaHbBMvIn8UkX3ucfxERE6o572T3X+n293pH4pItrsvO0VkRlP74u7P3aq6xd2Xj4H3gDMaWl9Eprq9qgMiskNEJtWzzkki8ra7D4Ui8icR6RGyvN7j7v79ZLrv/ZWIPBjOPnQ4qmo/bfAHGAwU4yT9FGAXkBOyrAiIcqcVONl9/TRwr/v6nJBtknG+8fVwp2OAvcDpDbT/BbDabbsX8Bnwk7rvG7LuBe7r+cAf3df/Bnwcst4YYB8QV097A939eB7oCowCCkLe93Rgghv3QDeeW0K2V+BNN9bOdY9LU3HXiSHGnV4FPAbEA+luPOeF7Odh4FIgGvgV8JG7bBhOTyEl5H1PasbfwAluG8MbWP5D4P06885xj10UMBr4Cvi2u+z/ASuALm7MpwOJ7rJ3gB8Dg4BtwGx3flfgADAs5O9oRDP2pTOQD0xqYPk4oAS40I29f81+18Tmvj7ZXacTkOT+Gz3c1HEHPgR+4L5OACYE/X+8Nf5YD6ONUuecxEGcD6qzgTeAPBEZjjNk8Z6qVh/H++Xj/Of6njtrElCoqmsa2exRVc1T1f04HzTpx7kby4GhIjLEnf4BsFRVjzayzQJ1vmlvAJ4Cprvxr1HVj1S1UlW/AH6HcxxC/UpV96tq+XHG+TUikgZMBO5Q1cOqmgU8AYSOw7+vqn9X55zHszgJEaAK5wPtVBGJVafHsOM4248F/gQsUdUt4W6nqu+o6gZ1vtV/ipOAa45TBdAbJ4lWucf0QMjmpwIrgbtV9fGQ+dXASBHprKr5qrrpePbF9b/Aepy/4/pcC/xBVd90Y8+tb79Vdbu7zhFVLQAeDNm/xo57BXCyiPRR1VJV/agZ+9DuWcJo297F+cZ4tvv6HZz/HN9yp4/XEmCm+3omzodcY/aEvC7D+WYWNlU9DCwFZopIFM6Hf1Nt7g55vQunh4OIDBWRV0Rkj4gcAH4J9Glk25ZKAfar6sE68fQPma57fOJFJEZVtwO34PRC9orICyKSEm7D7rF6FjgK3Hg8QYvIeBFZKSIFIlIC/IR/HadncT6wX3CH/f7LTUw1ZgC5wIs1M1T1EPB9933y3WHG4ccZ06+BkcC/qfsVvx5pQJNJVUROcI9nrvt38Mea/WviuF8LDAW2uENxlx/PPnQUljDatpqEcZb7+l1aljD+BowWkZHA5TjfYL22BOeD6HygTFU/bGL9tJDXJwJ57uvfAluAIaqaCPwCkDrbHm9p5sbWzwN6iUi3OvHkhvXGqs+p6pnAALed+8PZTkQEeBJnOOq7qlrRWDP1zHsOp2eXpqrdcb7ZixtThaouUNVTgW/i/A2E9pjmA4XAcxJyabGqvqGqF+IMR20Bfh/Ovrj7swC4BLioTm+mrt3ASWG85S9x9nuU+3cwk5C/g4aOu6p+rqrTgb7uvBdFpGu4+9FRWMJo294FzsUZk8/BOWk4CWdYYd3xvpn7jf9FnA+V1ar6ZQRjbajND3GGNH5D070LgP8UkS4iMgKYhdNDAeiGM5Ze6n7DvS6M9/oK53zPcS9X1d3A/wG/ck8Wj8b5ltrkPSYiMkxEzhORTjjnIMpxjkHNBQONJarfAqcAk8MYWvsKSA09+Y9znPar6mERGQdcFRLXuSIyyk0GB3CGaUKHNStwhiy7As+ISJT7jX6q++F6BCgN2ZeaiwQGNnAcfu62f4Gq7mtiX54EZonI+W67/RvoyXRzYygRkf7A7SHtNXbcZ4pIkjuMW+xuEvaQbkdhCaMNU9VtOP853nOnDwDZwAfa/HsFluCcFA3nwztSnnHbDOeGvneB7cA/gQdU9R/u/NtwPnwO4nzDXVr/5seYDyxxrwj6t3qW/wqY6y6/rZ7l03FOnOYBy3DG9t8Ko91OwCKcb+t7cL7V/txdloaTiL5GRAbgnJhOB/a4Vz+VNnJV0tvAJnfdQnfe9cBCETkIzAP+HLJ+P5wvDAdwLhp4lzp/B+75pe/g9HD+gHORwX+4x2A/Tu+2Jlmn4QzTNdTr+iVOr2x7yL78or4VVXU1zheEh3BOfr+L00uoawGQ4a7zKvDXkGWNHfdJwCYRKQUeAaZF4lxXeyMNDxmajkhETsQZVujXxBBBJNu8GueqmzMbWWcgziWksapa6UdcQRCRJ4C/qGpDJ3/bDBGZCxSo6u+CjsVEhiUMU8s9mfogzqWUP/KpzS4434QfU9VnGllvIB0gYRjTmtmQlAHAHYM+gHMN+90+tXkxzr0LX+GcNzHGtGKe9zDcuzEfwbkR6AlVXVRn+X/g3BBUifPh8SNV3eUuqwI2uKt+qapTPA3WGGNMgzxNGO7VFttwvrXmAJ8A01V1c8g65+Lc7VsmItcB56jq991lpap6XNf2G2OM8YbXRdjGAdvdu5IRkReAqUBtwlDVlSHrf8S/bhw7bn369NGBAwc2d3NjjOmQ1qxZU6iqSU2t53XC6M+xd9fmAOMbWf9a4LWQ6XgRycQZrlqkqn+ru4GIzAZmA5x44olkZma2OGhjjOlIRGRXOOu1mjLPIjITGMux9X8GqGquOA//eVtENtStuePWtHkcYOzYsXbJlzHGeMTrq6RyObaUQyr13MQjzkNq7gKmqOqRmvmqmuv+zsapk3Sal8EaY4xpmNcJ4xNgiIgMcssTTMOpY1NLRE7DqSw6RVX3hszv6d7Cj4j0wakMuhljjDGB8HRISlUrReRGnAqY0TjliTeJyEIgU1WXA7/GqXL6F6euWu3ls6cAvxORapzEtij06ipjTPtSUVFBTk4Ohw8fDjqUdis+Pp7U1FRiY2ObXrke7epO77Fjx6qd9Dambdq5cyfdunWjd+/euF8eTQSpKvv27ePgwYMMGjTomGUiskZVxzawaS2709sY0yocPnzYkoWHRITevXu3qAdnCcMY02pYsvBWS4+vJQxjjDFhsYRhjDGuL774gpEjRwYdRqtlCcMYY0xYLGEYY0yIyspKZsyYwSmnnMKVV15JWVkZCxcu5Bvf+AYjR45k9uzZqCo7duwgIyOjdrvPP/+8dnrNmjV861vf4vTTT+fiiy8mPz8fgEcffZRTTz2V0aNHM23atED2ryVaTWkQY4ypccstt5CVlRXR90xPT+fhhx9ucr2tW7fy5JNPMnHiRH70ox/x2GOPceONNzJv3jwAfvCDH/DKK68wefJkunfvTlZWFunp6Tz11FPMmjWLiooKfvrTn/Lyyy+TlJTE0qVLueuuu/jDH/7AokWL2LlzJ506daK4uLiJSFof62GEyC0uZ8X6PHKL7VG+xnRUaWlpTJw4EYCZM2fy/vvvs3LlSsaPH8+oUaN4++232bRpEwA//vGPeeqpp6iqqmLp0qVcddVVbN26lY0bN3LhhReSnp7OvffeS05ODgCjR49mxowZ/PGPfyQmpu19X297EXskt7icix9ahaoiIrxx69n079E56LCM6ZDC6Ql4pe6lpyLC9ddfT2ZmJmlpacyfP7/2Xobvfve7LFiwgPPOO4/TTz+d3r17k5eXx4gRI/jwww+/9t6vvvoqq1atYsWKFdx3331s2LChTSUO62G41u4qQlU5dLQKVWXtrqKgQzLGBODLL7+s/bB/7rnnOPPMMwHo06cPpaWlvPjii7XrxsfHc/HFF3Pdddcxa9YsAIYNG0ZBQUHte1RUVLBp0yaqq6vZvXs35557Lvfffz8lJSWUlpb6vHct03ZSm8cyBvREROgaF42IkDGgZ9AhGWMCMGzYMBYvXsyPfvQjTj31VK677jqKiooYOXIk/fr14xvf+MYx68+YMYNly5Zx0UUXARAXF8eLL77ITTfdRElJCZWVldxyyy0MHTqUmTNnUlJSgqpy00030aNHjyB2sdmsllSIbdlfkvnxOsaOP42hg0+MYGTGmKZ89tlnnHLKKUGHcdweeOABSkpKuOeee4IOJSz1Hedwa0lZD8N1oHAvbyz8Garwxht/pN8Di0ns0zfosIwxrdgVV1zBjh07ePvtt4MOxReWMFx527agChWHy4mN70zeti2WMIwxjVq2bFnQIfjKTnq7UoYOR1Bi42IQlJShw4MOyRhjWhXrYbgSY45wzaA15B3qQkrXMhJjjjS9kTHGdCCWMGrkrCYx7iiJ7Ie4BMhZDT3Smt7OGGM6CBuSqpE6DhAnWSDutDHGmBqWMGr0SIPrP4Qpjzq/rXdhTIdSXFzMY4891uztBw4cSGFhYaPr3HXXXaSlpZGQkNDsdoJkCSNERVk0B77sTEVZdNChGGN81tKEEY7JkyezevVqT9vwkiUMV0VeHtmTp5A39z/JnjyFiry8oEMyxvjozjvvZMeOHaSnp3P77bdTWlrK+eefT0ZGBqNGjeLll18G4NChQ1x22WWMGTOGkSNHsnTp0mPep7y8nEsuuYTf//73X2tjwoQJJCcn+7I/XrCT3q7yrCxUFS0rgy5dKM/KIjYlJeiwjDGNyC0uZ+2uIjIG9GxxsdBFixaxcePG2rLqlZWVLFu2jMTERAoLC5kwYQJTpkzh9ddfJyUlhVdffRWAkpKS2vcoLS1l2rRpXH311Vx99dUtiqc1soTh6pye7lSp7NIFEaFzenrQIRljGuF1hWlV5Re/+AWrVq0iKiqK3NxcvvrqK0aNGsXPfvYz7rjjDi6//HLOOuus2m2mTp3KnDlzmDFjRsTiaE1sSMoVm5LC4BXLSbn3HgavWG69C2NaOa8rTP/pT3+ioKCANWvWkJWVxQknnMDhw4cZOnQoa9euZdSoUcydO5eFCxfWbjNx4kRef/112lONvlCWMELEpqSQeOmlliyMaQMiXWG6W7duHDx4sHa6pKSEvn37Ehsby8qVK9m1axcAeXl5dOnShZkzZ3L77bezdu3a2m0WLlxIz549ueGGG1oUS2tlCcMY0yb179GZN249m0XfHR2R4ajevXszceJERo4cye23386MGTPIzMxk1KhRPPPMMwwf7pQL2rBhA+PGjSM9PZ0FCxYwd+7cY97nkUceoby8nDlz5nytjTlz5pCamkpZWRmpqanMnz+/RTH7zcqbG2NahbZa3rytaUl5c+thGGOMCYslDGOMMWGxhBEivzSf13e+Tn5pftChGGNMq2P3YbjyS/O5YvkVtdd0L5uyjOSEtntHpjHGRJr1MFzrC9ajqpRVlqGqrC9YH3RIxhjTqljCcI1JGoOI0CXGudN7TNKYoEMyxphWxfOEISKTRGSriGwXkTvrWf4fIrJZRD4VkX+KyICQZdeIyOfuzzVexpmckMyyKctY8M0FNhxlTAfkdXnzsrIyLrvsMoYPH86IESO4886vfRy2ep4mDBGJBhYDlwCnAtNF5NQ6q60DxqrqaOBF4L/cbXsBdwPjgXHA3SLSsls5m5CckMykQZMsWRjTAflR3vy2225jy5YtrFu3jg8++IDXXnvN0/Yizesexjhgu6pmq+pR4AVgaugKqrpSVcvcyY+AVPf1xcCbqrpfVYuAN4FJHsdrjOmgvC5v3qVLF84991wA4uLiyMjIICcnx5+dixCvr5LqD+wOmc7B6TE05FqgJuXWt23/iEZnjGnbindDzmrnkcotfEqmn+XNi4uLWbFiBTfffHOLYvZbqznpLSIzgbHAr49zu9kikikimQUFBd4EZ4xpfYp3w2NnwPKbnN/Fu5ve5jjUlDcfPXo0F1xwwTHlzd98803uuOMO3nvvPbp37167zdSpU5k1a1ajyaKyspLp06dz0003MXjw4IjG7DWvE0YuEJr2U915xxCRC4C7gCmqeuR4tlXVx1V1rKqOTUpKiljgxphWLmc1oHC01PmdE9lHn3pV3nz27NkMGTKEW265JaLx+sHrhPEJMEREBolIHDANWB66goicBvwOJ1nsDVn0BnCRiPR0T3Zf5M4zxhhnGAqBuATnd+q4Fr2dH+XN586dS0lJCQ8//HCLYg2KpwlDVSuBG3E+6D8D/qyqm0RkoYhMcVf7NZAA/EVEskRkubvtfuAenKTzCbDQnWeMMc45i+s/hCmPOr9beA7D6/LmOTk53HfffWzevJmMjAzS09N54oknWhSz36y8eYiD+w+zJ7uEfoO7061XfAQjM8Y0xcqb+6Ml5c2tlpTr4P7DPL/wY1BAYPq88ZY0jDEmRKu5Sipoe7JLQKHiSBWoO22MMaaWJQxXv8HdQSC2UzSIO22MMaaWDUm5uvWKZ/q88XYOwxhjGmAJI0S3XvGWKIwxpgE2JGWMMSYsljCMMQbvy5sD3HXXXaSlpZGQkNDsdoJkCcMYY/CnvPnkyZNZvTqyJUz8ZAnDGGPwvrw5wIQJE0hObvx5O6WlpcyaNYtRo0YxevRoXnrpJQD+8Y9/cMYZZ5CRkcH3vvc9SktLAadnM2fOHEaNGsW4cePYvn17JA5HveyktzGmzcovzWd9wXrGJI1p8YPP/Cxv3ph77rmH7t27s2HDBgCKioooLCzk3nvv5a233qJr167cf//9PPjgg8ybNw+gdv1nnnmGW265hVdeeaUlh6JBljCMMW1Sfmk+Vyy/AlVFRCL+aOWa8uarVq0iKirqmPLmP/vZz7jjjju4/PLLOeuss2q3mTp1KnPmzGHGjBnNbvett97ihRdeqJ3u2bMnr7zyCps3b2bixIkAHD16lDPOOKN2nenTp9f+vvXWW5vddlNsSMoY0yatL1iPqlJWWYaqsr5gfUTf36vy5s2hqlx44YVkZWWRlZXF5s2befLJJ2uXi0i9ryPNEoYxpk0akzQGEaFLTBdEhDFJY1r0fn6UNw/HhRdeyOLFi2uni4qKmDBhAh988EHt+YlDhw6xbdu22nVqzqMsXbr0mJ5HpFnCMMa0SckJySybsowF31wQkeEor8ubA8yZM4fU1FTKyspITU1l/vz5X1tn7ty5FBUVMXLkSMaMGcPKlStJSkri6aefZvr06YwePZozzjiDLVu21G5TVFTE6NGjeeSRR3jooYdadBwaY+XNjTGtgpU3b56BAweSmZlJnz59wlq/JeXNrYcRhtziclaszyO3uDzoUIwxJjB2lVQTcovLufihVbVXYrxx69n079E56LCMMQaAL774wre2rIfRhLW7ilBVDh2tQlVZu6so6JCMMSYQ1sNoQsaAniTLPsbEbWW9DCNjQM+gQzLGmEBYwmhCfwr5R6c5VFVXEx0VRRRnAy172LwxxrRFNiTVlJzVRAGxlWXOwcppu4XDjDGmJSxhNCV1HCAQl+D8Th0XdETGGA94Xd68rKyMyy67jOHDhzNixAjuvPPOZrcVFEsYTemRBtd/CFMedX73sOEoY9ojP8qb33bbbWzZsoV169bxwQcf8Nprr3naXqRZwghHjzQY+V1LFsa0Y16XN+/SpQvnnnsuAHFxcWRkZJCTk/O1OFpzeXNUtd38nH766WqMaZs2b9583Nsczc3Vkldf1aO5uS1uf+fOnTpixIja6YqKCi0pKVFV1YKCAj3ppJO0urpaX3zxRf3xj39cu15xcbGqqg4YMEB37typ559/vi5ZsqTRtoqKinTQoEG6Y8eOry2bM2eO3nzzzbXT+/fv14KCAj3rrLO0tLRUVVUXLVqkCxYsqG333nvvVVXVJUuW6GWXXdZo2/UdZyBTw/iMtaukjDFtUkVeHtmTp9TeVDt4xXJiU1Ii9v7qUXnzyspKpk+fzk033cTgwYO/ttzKmxtjTISVZ2U533zLnPLm5e6DjyLFq/Lms2fPZsiQIdxyyy1hx6JW3twYY5qvc3o6IoJ0ccqbd05Pb9H7+VHefO7cuZSUlPDwww83GIeVNzfGmAiLTUlh8IrlpNx7T0SGo7wub56Tk8N9993H5s2bycjIID09nSeeeOJrcVh5c594Vd784P7D7Mkuod/g7nTrFR/x9zfGWHnz5vKzvLmd9G7Cwf2HeX7hx6CAwPR54y1pGGM6JEsYTdiTXQIKFUeqiO0UzZ7sEksYxphWw8qbtyL9BncHgdhO0SDutDHGdEDWw2hCt17xTJ833s5hGGM6PM97GCIySUS2ish2EflatS0ROVtE1opIpYhcWWdZlYhkuT/LvY61Id16xTNk7AmWLIwxHZqnPQwRiQYWAxcCOcAnIrJcVTeHrPYl8EPgtnreolxVW3ZxtTHGmIjwuocxDtiuqtmqehR4AZgauoKqfqGqnwLVHsdijDEN8rq8OcBdd91FWloaCQkJDa7z9NNPc+ONNx4zr7HS6PPnz0dEjik6+PDDDyMiRPo2A68TRn9gd8h0jjsvXPEikikiH4nIt+tbQURmu+tkFhQUtCRWY0wH5kd588mTJ7N6dfMewtZYafRRo0YdU3/qL3/5CyNGjGhxvHW19qukBrg3k1wFPCwiJ9VdQVUfV9Wxqjo2KSnJ/wiNMe2C1+XNASZMmEBycvJxx9ZUafRvf/vbtfHt2LGD7t27h30j3/Hw+iqpXI59AHaqOy8sqprr/s4WkXeA04AdkQzQGNN2RbIKw6JFi9i4cSNZbhHDyspKli1bRmJiIoWFhUyYMIEpU6bw+uuvk5KSwquvvgo4NadqlJaWMm3aNK6++mquvvrqFsXTkOLiYlasWMHNN99cOy8xMZG0tDQ2btzIyy+/zPe//32eeuqpiLftdQ/jE2CIiAwSkThgGhDW1U4i0lNEOrmv+wATgc2Nb+WNA4V72fJ/qzhQuDeI5o0x9aipwrDy2S08v/BjDu4/HNH3rylvPnr0aC644IJjypu/+eab3HHHHbz33nt07/6ve7OmTp3KrFmzPEsWjZVGnzZtGi+88AJ/+9vfuOKKKzxp39OEoaqVwI3AG8BnwJ9VdZOILBSRKQAi8g0RyQG+B/xORDa5m58CZIrIemAlsKjO1VW+OFC4lyW33cA/fvffLLntBksaxrQSoVUYUHc6grwqb94SjZVGv/zyy3n22Wc58cQTSUxM9KR9z2/cU9W/A3+vM29eyOtPcIaq6m73f8Aor+NrSt62LVRXV1N55AgxnTqRt20LiX36Bh2WMR1epKswHE958169ejFz5kx69OhxTMXZhQsXsnDhQm644YaIn0CvKY1eX4VbcM5z3H///QwdOjSi7YZq7Se9AxfTvxflVeVURFdTXlVOTP9eQYdkjOFfVRjO/cHwiBQF9bq8OcCcOXNITU2lrKyM1NRU5s+fX28sTz/9NKmpqbU/4ZZGnzZtGhkZGS06Do2x8uZNeH3n69z/1gISCqsp7RPFHRfczaRBkyLahjHGypv7xcqbe2hM0hjKuyplXUBEGZM0JuiQjDEmEJYwmpCckMyyKctYX7CeMUljSE44/muojTGmPbCEEYbkhGRLFMaYDs9OehtjjAlL2AlDRO4PZ54xxpj26Xh6GBfWM++SSAVijDGmdWsyYYjIdSKyARgmIp+G/OwEPvU+RGOM8Z7X5c0bK1Eeqq2XN38OmIxTA2pyyM/pqjozotEYY0xA/Chv3liJ8pZs22rKm6tqifuQo+mqugsoBxRIEJETIx6RMcYEwOvy5k2VKG9MmytvLiKTgQeBFGAvMACnoGDk05gxxoThQOFe8rZtIWXo8BbXePOzvHl9JcrD1VbKm98LTAC2qeog4Hzgo4hHZIwxYfC6krRX5c0bK1HelLZU3rxCVfcBUSISpaorgSZrj7QHFXl5HPj736nIyws6FGOMK2/bFlSh4nA5qs50JHlV3ryxEuVNaUvlzYtFJAFYBfxJRPYChzyJqhWpyMsje/IUVBURYfCK5cSmpAQdljEdXsrQ4YhAbHxnRJzplvCjvHlTJcob09bKm0/FOeF9K/A6zqNSJ3sRVGtSnpWFqqJlZagq5e74pjEmWIl9+nLNA4u56P/9lGseWNzicxhelzcPt0Q5WHlzX3hR3tx6GMb4w8qb+8OX8uYi8h3gfqAvIO6Pqqo3g2WtRGxKCoNXLKc8K4vO6emWLIwxHdbxnMP4L2Cyqn7mVTCtVWxKiiUKY0yHdzznML7qiMnCGOOf9jRE3hq19Pg22cNwh6IAMkVkKfA34EhIAH9tUQTGGAPEx8ezb98+evfujYgEHU67o6rs27eP+PjmP/s8nCGp0CuhyoCLQmMALGEYY1qs5mqggoKCoENpt+Lj40lNTW329k0mDFWdFc4bicjPVfVXzY7EGNOhxcbGMmjQoKDDMI2I5BP3vhfB9zLGGNPKRDJh2KCjMca0Y5FMGHZ5gzHGtGPWwzDGGBOWFiUMEQktmfiXFsacodNkAAAWO0lEQVRijDGmFWtpD+M/al6o6i9b+F7GGGNasZYmDBuGqk/xbtj4kvPbGGPaieOpJVUfO9FdV/FueOwMnEMjcP2H0CMt6KiMMabFwikNcpD6E4MAnSMeUVuXsxpQOFoKcQnOtCUMY0w7EM6d3t38CKTdSB0HiJMsEHfaGGPavpYOSZm6eqQ5w1A5q51kYb0LY0w7Ecn7MOolIpNEZKuIbBeRO+tZfraIrBWRShG5ss6ya0Tkc/fnGq9jjZgeaTDyu5YsjDHtiqcJQ0SigcXAJcCpwHQRObXOal8CPwSeq7NtL+BuYDwwDrhbRHp6Ga8xxpiGed3DGAdsV9VsVT0KvABMDV1BVb9Q1U+B6jrbXgy8qar7VbUIeBOY5HG8xhhjGuB1wugPhN6MkOPOi9i2IjJbRDJFJNPq6BtjjHc8P4fhNVV9XFXHqurYpKSkoMMBILe4nBXr88gtLg86FGOMiRivr5LKBULP/Ka688Ld9pw6274Tkag8lFtczsUPrUJVERHeuPVs+vew21WMMW2f1z2MT4AhIjJIROKAacDyMLd9A7hIRHq6J7svcue1amt3FaGqHDpahaqydldR0CEZY0xEeJowVLUSuBHng/4z4M+quklEForIFAAR+YaI5OA8se93IrLJ3XY/cA9O0vkEWOjOazUOFO5ly/+t4kDh3tp5GQN6IiJ0jYtGRMgYYBd2GWPaB1FtP+Wgxo4dq5mZmb60daBwL0t+dj1aXYlExXDNbx4jsU9fwBmWWruriIwBPW04yhjT6onIGlUd29R6bf6kd1Dy1n2AHj1ExdFK9Ogh8tZ9ULusf4/OTB6TYsnCGNOuWMJoppQuBxAgVioRd9oYY9ozqyXVTIkjzuMHw+4mt7QL/RPKSBxxXtAhGWOMpyxhNFMufbhG72NMl62s12EsoU/YdyQaY0xbZAmjmdbuKiJfe7P96Hi6xkWzdleRnbMwxrRrdg6jmezyWWNMR2M9jGbq36Mzb9x6tl0+a4zpMCxhtED/Hp0tURhjOgwbkvLAwf2H+TzzKw7uP/y1ZVaY0BjTVlkPowUO7j/MnuwS+g3uTrde8bXznl/4MSggMH3e+NplVpjQGNOWWcJopoYSw57sElCoOFJFbKdo9mSX1CaM0MKEdmWVMaatsSGpZgpNDKg7DfQb3B0EYjtFg7jTLruyyhjTllkPo5kaSgzdesUz+cbBbP0oi2ET0mt7F2BXVhlj2jarVtsC9Z3DOFC4lyW33YAqiMA1DyyurWJrjDGtUbjVaq2H0QLdesUf04MAyNu2BVWoOFxObHxn8rZtsYRhjGkX7BxGhKUMHY4IxMZ3RsSZNsaY9sB6GBGW2Kcv1zywmLxtW0gZOtx6F8aYdsMShgcS+/S1RGGMaXdsSMpvxbth40vOb2OMaUOsh+Gn4t3w2BnU3u13/YfQIy3oqIwxJizWw/BTzmpA4Wip8ztnddARGWNM2Cxh+Cl1HAcr+/B5xfkcrOwDqeOCjsgYY8JmQ1I+OlidxHMFD6FVFUh0LFdVJ9Et6KCMMSZM1sPw0ZbVW6k+spejh7OpPrKXLau3Bh2SMcaEzXoYPir96n3KSt4CoAIo/eoCYEygMRljTLgsYfioc88kYuQIlRpFjFTTuWdS0CEZY0zYbEjKRycOH0WlOoe8UqM4cfgo/xq3+z+MMS1kPQwfHSgsIKZTJyqPHCGmUycOFBaQMvQU7xu2+z+MMRFgPQwfpQwdTlRUFLHxnYmKivKvMKHd/2GMiQDrYfgosU9frpz7YO3DlXyrN5U6DhCIS3B+2/0fxphmsITho4P7D7Pif7JBE9mWmc30eUlfe56GJ3qkOcNQOaudZGHDUcaYZrCE4aPQ54DHdopmT3aJPwkDnCRhicIY0wJ2DsNHDT0H3Bhj2gLPexgiMgl4BIgGnlDVRXWWdwKeAU4H9gHfV9UvRGQg8BlQczv0R6r6E6/jjYSKvDzKs7LonJ5ObEpK7fxuveKZPm/8154DbowxbYGnCUNEooHFwIVADvCJiCxX1c0hq10LFKnqySIyDbgf+L67bIeqpnsZY0vUlxgq8vLYPnky1dVVREVFc/KKFV9LGpYojDFtkddDUuOA7aqarapHgReAqXXWmQoscV+/CJwvIuJxXC1WkZdH9uQp5M39T7InT6EiLw+AvI/eobyijKjyI5RXlJH30TvBBmqMMRHidcLoD4TeWpzjzqt3HVWtBEqA3u6yQSKyTkTeFZGz6mtARGaLSKaIZBYUFEQ2+kaUZ2WhqmhZGapKeVYWANv6O7muPJZjpo0xpq1rzVdJ5QMnquo+ETkd+JuIjFDVA6ErqerjwOMAY8eOVb+C65yejohAly6ICJ3TnZGzkSPO4drrHuKk3ZXsSIvhyRHn+BWSMcZ4yuuEkQuEXsuZ6s6rb50cEYkBugP7VFWBIwCqukZEdgBDgUyPYw5LbEoKCUufYMf7r3HSmZfUnqdITkjmyR+8zPqC9YxJGkNyQnLAkRpjTGR4nTA+AYaIyCCcxDANuKrOOsuBa4APgSuBt1VVRSQJ2K+qVSIyGBgCZHscb9jyS/P57uqfoLGKrH6ZZf2W1SaH5ITkVpcocovLWburiIwBPenfo3PQ4Rhj2iBPE4aqVorIjcAbOJfV/kFVN4nIQiBTVZcDTwLPish2YD9OUgE4G1goIhVANfATVd3vZbzHY33BelSVssoyusR0YX3B+laXJGrkFpdz8UOrUFVEhDduPduShjHmuHl+DkNV/w78vc68eSGvDwPfq2e7l4CXvI6vucYkjUFE6BLjnMMYk9R6H4S0dlcRqsqho1V0jYtm7a4iSxjGmOPWmk96t2rJCcksm7KsTZyryBjQk4TKUgYdzqcgKpmMAT2DDskY0wZZwmiB5pyryC/N9z3JdKs8yFU5S6mqqiI6OppulWcD1sMwxhwfSxg+yi/N54rlV9SeS1g2ZZkvSSN77QYqj1aCVlBZFUv22g2kX3S+5+0aY9oXKz7oo/UF6+lRXMmYT0vpUVzJ+oL1vrQrMcmAALGAuNP+yS0uZ8X6PHKLy31t1xgTWdbD8FFyaSy//O0hoquhKuoQ0WfF+tLu4PTBdE2aRdWRPKI7pTA4fbAv7YJdoWVMe2IJw0cHV66kV4XzXV+rYP/KlTDK+6Ghbr3imbHgwkCq5NoVWsa0H5YwfNQ/oT/ldab9ElSV3IwBPRERusZFIyK+X6FlNywaEzmWMHyUOmkq2//7caorK4mKiSF1Ut3Cve1P/x6deevaweRueJf+o75FPx8/tAMfDivebY/FNe2KJQwfxaakcPKrr9b7cKV2q3g3/f54Lv1QWC/Os8V9+vAMdDiseDc8dgaggL/7bYxXLGH4LDYlpWMkiho5qwGFo6UQl+BM+/TBGehwWID7bYxXLGEYb6WOA8T50ETcaX/079GZ5VeeyOcfZjHkjHR/h6NSx1ENVMV0IRqI8nG/jfGKJYwOoqHnjHuuRxoHr1rFnnUb6XfaSLr5+C374BfZvPVYNpDIrnXZTL+jE90G+nNJcS59uObIfzFGt7JehrGEPl97cphnbduJfuMRSxgdQFPPGffSwf2H+dNvNlF1JJ/ot4UZC1J8u1prz7qNQDQV2plYKWfPuo2+JYy1u4rI195sPzre1/MngZ/oD5AlSu/Znd4dQJDPGc/OyuZQwVMcLnmdQwVPkZ3l3yNN+p02EoBYKT9m2g8ZA3rSrVoYUx1Lt2r/zp+s3VVEPy3g3Mr36acFrN1V5Eu7NfZ8+TlrXn2CPV9+7mu7NYnyzpc+5eKHVllVAY9YD6MD2NZf6M2xzxkf4FPbWpmPc6VQBRDrTp/qS9vdBg5m+h38azjMp94FQGK1MOtgPFXV1UQfiSKx2p9nu3+j5yH+JrchsYoilPZ8z5d2wUkWCU+exXAUXS3sufY9+p04xJe2axLl8MotbIkZbjeIesQSRgcQ5HPGB2eMYtWfYqiuiiIqOorBGaN8axucpOFnoqixJ7sErTpAZVkOUV1S2ZNd4stQXL8Dn1ISlcTesjT6dtlNvwOf4jys0nu5G95lOEpXOcwhjSd3w7u+JYwgEyV0nOEwSxgdQJDPGU/s05cf/mYxedu2kDJ0OIl9+vrWdpASehzlUOFToEpFmZDQI92Xdg8mnMbSPYuciQMwPeEUuvnSMvQf9S10tXBI41GE/qO+5VPLTqKsjosiquIQ1bFdSfAzUXag80aWMDqIIJ8zntinb2CJ4kDh3kCSVcnencTERVF55DAxcfGU7N1J8sneXyG2p7Ar1VGVHC3PI65zCnsKu9JtoOfNAtDvxCHs+M5ytr33d4aedSkn+dS7ACB1HIeqkthTMYF+UTvo5uNlzB2pXpolDNNuHSjcy5LbbkAVROCaBxb7ljRShg4nKkqIje+MiDPth4QeRzm0b0lIzybDl3bBOd5/f+jXqML2zE2+Hu+D1Uk8v+9hqK6GqCimVyf51rMKul6anyVoLGGYditv2xZUoeJwObHxncnbtsW3D7DEPn255gH/h+KC6tmAc7yrq9Vpu1O8r8d7T3YJ1VWlHC3LIc7Hc0YQbL00indT/dgE5+KKqCiirv/I06RhCcN4LqibBlOGDkdRouJiUdS3b/k1ghiKC6pnA9C97yAqj1YDsVQeraZ730G+tR3UOSMg0Hpp+7e+T6cjlbUXGhzZ+j69xk/3rD1LGMZTQd40eCi+infScxi6M4ptg6r5t/gqEn1pOThB9WwASovj6NpnVu23/NLiON/aDrJnRc5q8suS2FqczLAe+ST7WDdsbfUQJiCUajwgrK0ewgUetmcJw3iq5qbBzkehPM6ZHvCdq3xpe+Omd/jPp0qc20DegY3ffIdkD7991ZW7PYsd77/GSWdeQv+T/fvGG9RFBv0GdycqJpH4bqeCONN+CbJnlV8xgOc+HwjAmoKBXFUxAL8uLznllBFMf/UhTjtayLq4PvzvKSM8bc8ShvFUkDcNDs1V9gGdK5xkNTRXfWrZSRZffWc6CcBXDz4Df33et6QR1BBgt17xTJ83PpAnOyb26cuM2+ex6523GXDOeb4mzK2b8t1Xzs2pWzflk3yaP20nVguXFidQcaSYSzsleH6DqCUM46kgbxpMmXAOZbG/oTq6is5R0aRM8K/tHe+/RgLU9qx2vP+aLwkjyCFAgPjD+zlhbxbxKemAf+1W5OXx1dU/JK66iq/+8AzdfNzvYRPSWfOKgMaCCMMm+NebzM7Kpqzm3I0I2VknMeY87yopWMIwngrypsHYlBROXrEikG/bJ515CV89+Azlcf+a9kOQQ4BBJqsg9zv55DSuuvcRtn6UxbAJ6f6dO8H/0juWMIzngrxpMKgHVvU/OR3++rzv5zCCHAIM8kM7yP0GoF8MhyfGQpK/H6l+l96xhGGMR/qfnO7ryW4IdggwyA/tIPc7vzSfa5+dykm7K3k0LYYnf/Cyb1+QEvv05dKFd7Nm3UpOP+1cz8/diKp/JwK9NnbsWM3MzAw6DGMClV+aH8gQYOgH5w6fPzhr2g9iv9/8+Hl6//vC2se37/v9PC706Wq8SB1zEVmjqmObWs96GMa0M0ENAQZ5vqqm/SD2O8ir8TZueod7f3uwNlltzPD20nFLGMaYiAnyfFVQgrwaz+9kZQnDGGNaIMir8fxOVpYwjDGmhYK6Gs/vZGUJwxhj2jA/k1WU1w2IyCQR2Soi20XkznqWdxKRpe7yj0VkYMiyn7vzt4rIxV7HaowxpmGeJgwRiQYWA5fg3H44XUTq3oZ4LVCkqicDDwH3u9ueCkwDRgCTgMfc9zPGGBMAr4ekxgHbVTUbQEReAKYCm0PWmQrMd1+/CPyPiIg7/wVVPQLsFJHt7vt92FBjW7du5Zxzzon0PhhjjMH7Ian+wO6Q6Rx3Xr3rqGolUAL0DnNbRGS2iGSKSGZFRUUEQzfGGBOqzZ/0VtXHgcfBudP7nXfeCTYgY4xpY5xBnaZ53cPIBUJLN6a68+pdR0RigO7AvjC3NcYY4xOvE8YnwBARGSQicTgnsZfXWWc5cI37+krgbXUKXC0HprlXUQ0ChgCrPY7XGGNMAzwdklLVShG5EXgDiAb+oKqbRGQhkKmqy4EngWfdk9r7cZIK7np/xjlBXgncoKpVjbW3Zs2aQhHZFWZ4fYDCZu1Y29UR9xk65n53xH0G2+/mCquwcLuqVns8RCQznOqM7UlH3GfomPvdEfcZbL+9bsfzG/eMMca0D5YwjDHGhKUjJ4zHgw4gAB1xn6Fj7ndH3Gew/fZUhz2HYYwx5vh05B6GMcaY42AJwxhjTFg6XMJoqtx6WyYifxCRvSKyMWReLxF5U0Q+d3/3dOeLiDzqHodPRSQjuMibT0TSRGSliGwWkU0icrM7v73vd7yIrBaR9e5+L3DnD3IfE7DdfWxAnDu/wccItDUiEi0i60TkFXe6I+zzFyKyQUSyRCTTnef733iHShhhlltvy57GKQUf6k7gn6o6BPinOw3OMRji/swGfutTjJFWCfxMVU8FJgA3uP+m7X2/jwDnqeoYIB2YJCITcB4P8JD7uIAinMcHQAOPEWijbgY+C5nuCPsMcK6qpofcb+H/37iqdpgf4AzgjZDpnwM/DzquCO/jQGBjyPRWINl9nQxsdV//Dphe33pt+Qd4GbiwI+030AVYC4zHuds3xp1f+/eOU23hDPd1jLueBB17M/Y1FefD8TzgFUDa+z678X8B9Kkzz/e/8Q7VwyDMkuntzAmqmu++3gOc4L5ud8fCHXI4DfiYDrDf7tBMFrAXeBPYARSr85gAOHbfGnqMQFvzMDAHqHane9P+9xlAgX+IyBoRme3O8/1vvM2XNzfhU1UVkXZ5HbWIJAAvAbeo6oHQcs3tdb/Vqa2WLiI9gGXA8IBD8pSIXA7sVdU1InJO0PH47ExVzRWRvsCbIrIldKFff+MdrYfREUumfyUiyQDu773u/HZzLEQkFidZ/ElV/+rObvf7XUNVi4GVOMMxPdzHBMCx+9bQYwTakonAFBH5AngBZ1jqEdr3PgOgqrnu7704Xw7GEcDfeEdLGOGUW29vQsvHX4Mzxl8z/2r3iooJQElI97bNEKcr8STwmao+GLKove93ktuzQEQ645y3+QwncVzprlZ3v+t7jECboao/V9VUVR2I83/3bVWdQTveZwAR6Soi3WpeAxcBGwnibzzokzkBnDy6FNiGM957V9DxRHjfngfygQqccctrccZs/wl8DrwF9HLXFZwrxnYAG4CxQcffzH0+E2d891Mgy/25tAPs92hgnbvfG4F57vzBOM+N2Q78Bejkzo93p7e7ywcHvQ8t3P9zgFc6wj67+7fe/dlU87kVxN+4lQYxxhgTlo42JGWMMaaZLGEYY4wJiyUMY4wxYbGEYYwxJiyWMIwxxoTFEoYxxpiwWMIwxhgTFksYxnhIRAaKyGci8nv3uRX/cO/MNqbNsYRhjPeGAItVdQRQDHw34HiMaRZLGMZ4b6eqZrmv1+A8s8SYNscShjHeOxLyugp7rIBpoyxhGGOMCYslDGOMMWGxarXGGGPCYj0MY4wxYbGEYYwxJiyWMIwxxoTFEoYxxpiwWMIwxhgTFksYxhhjwmIJwxhjTFj+P+6gbu+Ria/kAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1)\n",
    "ax.scatter(ns, mean_errors_X, s=7.5, label = 'task 1')\n",
    "ax.scatter(ns, mean_errors_Z, s=7.5, label = 'task 2')\n",
    "\n",
    "ax.scatter(ns, mean_errors_HHX, s=7.5, label = 'task 1 cep')\n",
    "ax.scatter(ns, mean_errors_HHZ, s=7.5, label = 'task 2 cep')\n",
    "\n",
    "# ax.scatter(ns, mean_errors_joint, s=7.5, label = 'silly willy nilly')\n",
    "\n",
    "ax.scatter(ns, mean_errors_LLX, s=7.5, label = 'task 1 L2M')\n",
    "ax.scatter(ns, mean_errors_LLZ, s=7.5, label = 'task 2 L2M')\n",
    "\n",
    "ax.set_xlabel('n')\n",
    "ax.set_ylabel('L_hat')\n",
    "ax.set_title('willy nilly partitions, 2 tasks, 2 class')\n",
    "ax.axhline(y = 0, c='k',label = 'bayes')\n",
    "ax.legend()\n",
    "plt.savefig('willy_nilly_no_swn.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#- Code graveyard\n",
    "\n",
    "def estimate_posteriors(c, data, labels, c_left=-10**6, c_right=10**6):\n",
    "    unique_labels, label_counts = np.unique(labels)\n",
    "    \n",
    "    temp_idx = (data > c_left).astype(int) + (data < c_right).astype(int) == 2\n",
    "    temp_data = data[temp_idx]\n",
    "    temp_label = labels[temp_idx]\n",
    "    \n",
    "    unique_temp_labels, temp_label_counts = np.unique(temp_label)\n",
    "    \n",
    "    if len(unique_temp_labels) == 0:\n",
    "        return label_counts[-1] / np.sum(label_counts)\n",
    "    elif len(unique_temp_labels) == 2:\n",
    "        return temp_label_counts[-1] / len(temp_label)\n",
    "    elif unique_temp_labels[0] == 0:\n",
    "        return 0\n",
    "    else:\n",
    "        return 1\n",
    "    \n",
    "def build_linear_classifier(c, train, labels, test, posteriors=False):\n",
    "    #Assumes 1d\n",
    "    #Assumes 2 classes\n",
    "    #Assumes type(c) is array-like \n",
    "    c = np.sort(c)\n",
    "    \n",
    "    if train.ndim == 1:\n",
    "        n = len(train)\n",
    "        m = len(test)\n",
    "        proj_train = train\n",
    "        proj_test = test\n",
    "#     elif train.ndim == 2:\n",
    "#         n, d = train.shape\n",
    "#         m, d = test.shape\n",
    "        \n",
    "#         proj_train = train @ u\n",
    "#         proj_test = test @ u\n",
    "#     else:\n",
    "#         raise ValueError('bad')\n",
    "        \n",
    "    \n",
    "    \n",
    "    temp = train[train < c[0]]\n",
    "    temp_labels = labels[train < c[0]]\n",
    "    \n",
    "    if len(temp) == 0:\n",
    "        if posteriors:\n",
    "            partition_posteriors[0] = 0.5\n",
    "        else:\n",
    "            partition_posteriors[0] = np.random.binomial(1, 0.5)\n",
    "    else:\n",
    "        unique_, counts_ = np.unique(temp_labels, return_counts=True)\n",
    "        if posteriors:\n",
    "            partition_posteriors[0] = counts_[np.where(unique_ == 1)[0][0]]/np.sum(counts_)\n",
    "        else:\n",
    "            partition_posteriors[0] = unique_[np.argmax(counts_)]\n",
    "    \n",
    "    for i in range(1, len(c)):\n",
    "        temp_idx = ((train < c[i]).astype(int) + (train > c[i-1]).astype(int)) == 2\n",
    "        temp = train[temp_idx]\n",
    "        temp_labels = labels[temp_idx]\n",
    "        \n",
    "        if len(temp) == 0:\n",
    "            if posteriors:\n",
    "                partition_posteriors[0] = 0.5\n",
    "            else:\n",
    "                partition_posteriors[i] = np.random.binomial(1, 0.5)\n",
    "            \n",
    "        else:\n",
    "            unique_, counts_ = np.unique(temp_labels, return_counts=True)\n",
    "            partition_posteriors[i] = unique_[np.argmax(counts_)]\n",
    "            \n",
    "            if posteriors:\n",
    "                partition_posteriors[i] = counts_[np.where(unique_ == 1)[0][0]]/np.sum(counts_)\n",
    "            else:\n",
    "                partition_posteriors[i] = unique_[np.argmax(counts_)]\n",
    "            \n",
    "    temp = train[train > c[-1]]\n",
    "    temp_labels = labels[train > c[-1]]\n",
    "    \n",
    "    if len(temp) == 0:\n",
    "        if posteriors:\n",
    "                partition_posteriors[0] = 0.5\n",
    "        else:\n",
    "            partition_posteriors[i] = np.random.binomial(1, 0.5)\n",
    "    else:\n",
    "        unique_, counts_ = np.unique(temp_labels, return_counts=True)\n",
    "        if posteriors:\n",
    "            partition_posteriors[-1] = counts_[np.where(unique_ == 1)[0][0]]/np.sum(counts_)\n",
    "        else:\n",
    "            partition_posteriors[-1] = unique_[np.argmax(counts_)]\n",
    "\n",
    "    # Now, classify test observations    \n",
    "    predictions = -1*np.ones(m)\n",
    "    \n",
    "    partitions = np.sum([np.array(test > c_) for c_ in c], axis=0)\n",
    "    for i, c_ in enumerate(c):\n",
    "        predictions[np.where(partitions == i)] = partition_posteriors[i]\n",
    "        \n",
    "    predictions[predictions == -1] = partition_posteriors[-1]\n",
    "    \n",
    "    return predictions"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "hh",
   "language": "python",
   "name": "hh"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}