experiment_2b.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Section 7.2 part 2\n",
    "\n",
    "In this part of the experiment we consider fixed (p, q) in case where they "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "Wd-5j4FHr9LZ"
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%reload_ext autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "Wd-5j4FHr9LZ"
   },
   "outputs": [],
   "source": [
    "import torch\n",
    "from torch import nn\n",
    "from torch.nn import functional as F\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import cf.gadget_2_jax as JP_Z\n",
    "from cf import utils\n",
    "from cf import models\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "suRH4mqhr9Le"
   },
   "outputs": [],
   "source": [
    "def train_gsm_with_gadget1(num_epochs, model, logits_p, logits_q, R, lr=0.001, tau=1.0, verbose=True):\n",
    "    optimizer = torch.optim.Adam(model.parameters(), lr=lr)\n",
    "    lr_scheduler = torch.optim.lr_scheduler.ExponentialLR(optimizer=optimizer, gamma=schedualer_decay)\n",
    "    bs = 64\n",
    "\n",
    "    losses = []\n",
    "    l_avg = []\n",
    "\n",
    "    k = logits_p.size(-1)\n",
    "\n",
    "    RiRj = torch.tensor([(R(i) - R(j)) ** 2 for i in range(k) for j in range(k)])  # /normalizer\n",
    "    RiRj = RiRj.view(1, k, k).repeat([bs, 1, 1])\n",
    "\n",
    "    ma_l, ma_c = 0, 0\n",
    "    logits_p = logits_p.repeat([bs, 1])\n",
    "    logits_q = logits_q.repeat([bs, 1])\n",
    "\n",
    "    for epoch in range(num_epochs):\n",
    "        p, q = add_noise_if_generalized_pq(logits_p, logits_q, GENERALIZED)\n",
    "        \n",
    "        bs = p.size(0)\n",
    "        u = torch.rand(bs, p.size(-1), q.size(-1))\n",
    "        joint_pq = model(p, is_pq=True)\n",
    "        joint_qp = model(q, is_pq=False)\n",
    "\n",
    "\n",
    "        gumbel_noise = - torch.log(1e-8 - torch.log(u + 1e-8))\n",
    "        logprobs_p = utils.marginal_gumbel_log_softmax(joint_pq, gumbel_noise, rows=True, tau=tau, hard=False)\n",
    "        logprobs_q = utils.marginal_gumbel_log_softmax(joint_qp, gumbel_noise.transpose(-2, -1), rows=True, tau=tau,\n",
    "                                                    hard=False)\n",
    "\n",
    "        covariance = torch.exp(logprobs_p[:, :, None] + logprobs_q[:, None, :]) * RiRj.view(bs, k, k)\n",
    "\n",
    "        loss = torch.sum((logprobs_p[:, :, None] + logprobs_q[:, None, :]) * covariance.detach())\n",
    "\n",
    "        l_avg.append((loss.item(), covariance.sum().item()))\n",
    "\n",
    "        optimizer.zero_grad()\n",
    "        loss.backward()\n",
    "\n",
    "        optimizer.step()\n",
    "        if epoch % 100 == 0:\n",
    "            lr_scheduler.step() if schedualer_decay != 1.0 else None\n",
    "            tau = tau * tau_decay\n",
    "\n",
    "        if epoch % 50 == 0 or np.remainder(np.log2(epoch), 1) == 0:\n",
    "            l, c = zip(*l_avg)\n",
    "            l, c = np.mean(l), np.mean(c)\n",
    "            losses.append((l, c))\n",
    "            l_avg = []\n",
    "            ma_l = .99 * ma_l + .01 * l\n",
    "            ma_c = .99 * ma_c + .01 * c\n",
    "            if verbose:\n",
    "                print(f'{epoch} {tau} {l} {c}')\n",
    "\n",
    "    return losses\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def train_gsm_with_gadget2(num_epochs, model, logits_p, logits_q, reward_f, s_prime_obs=None, cf=False,\n",
    "                            verbose=True, seed=0):\n",
    "    if cf:\n",
    "        assert s_prime_obs is not None, 's_prime_obs must be given in the CF method'\n",
    "    else:\n",
    "        #set dummy s_prime\n",
    "        s_prime_obs = 0\n",
    "\n",
    "    logits_p = logits_p.numpy().clip(min=-80.0)\n",
    "    logits_q = logits_q.numpy().clip(min=-80.0)\n",
    "    k = logits_p.shape[-1]\n",
    "\n",
    "    rewards = np.array([reward_f(i) for i in range(k)])\n",
    "\n",
    "    model.train(logits_p, logits_q,\n",
    "                np.array(s_prime_obs), rewards, num_iter=num_epochs,\n",
    "                cf=cf, noise_scale=GENERALIZED, verbose=verbose, seed=seed)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "yn7tnhmqr9Ld"
   },
   "outputs": [],
   "source": [
    "def add_noise_if_generalized_pq(p, q, std):\n",
    "    # add noise to logits    \n",
    "    if std != 0:\n",
    "        return torch.clamp(p + std*torch.randn_like(p), max=0.0), torch.clamp(q + std*torch.randn_like(q), max=0.0)\n",
    "    else:\n",
    "        return p, q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "gOLSQo-er9Ld"
   },
   "outputs": [],
   "source": [
    "\n",
    "# Non-monotonic reward function\n",
    "def R_sin(a, s=None):\n",
    "    return np.sin(15 * a)\n",
    "\n",
    "\n",
    "def R_sin30(a, s=None):\n",
    "    return np.sin(30 * a)\n",
    "\n",
    "\n",
    "def R_random_sin30(a, s=None):\n",
    "    return RANDOM_NON_MONOTONE_R[a]\n",
    "\n",
    "\n",
    "# Monotonic (in the region)\n",
    "def R_random(a, s=None):\n",
    "    return RANDOM_MONOTONE_R[a] \n",
    "\n",
    "\n",
    "def R_sqr(a, s=None):\n",
    "    return a ** 2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "R_monotone = R_random\n",
    "R_non_monotone = R_random_sin30\n",
    "\n",
    "#figpath = \"./figs/toy_example\"\n",
    "figpath = \"figs/\"\n",
    "coupler_save_path = R_monotone.__name__ + '-' + R_non_monotone.__name__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_random_monotone():\n",
    "    return np.cumsum(np.random.rand(K))\n",
    "\n",
    "def create_random_non_monotone():\n",
    "    return R_sin30(np.random.normal(loc=np.zeros(K), scale=1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "id": "fJnT8Lifr9Le"
   },
   "outputs": [],
   "source": [
    "K = 10\n",
    "\n",
    "schedualer_decay = 1.0\n",
    "tau_decay = 1.0\n",
    "learning_rate = 0.00001\n",
    "\n",
    "torch.manual_seed(100)\n",
    "np.random.seed(100)\n",
    "\n",
    "    \n",
    "RANDOM_MONOTONE_R = create_random_monotone()\n",
    "RANDOM_NON_MONOTONE_R = create_random_non_monotone()\n",
    "\n",
    "GENERALIZED = 0\n",
    "p_logits = (torch.arange(K, 0, step=-1).float()/4).unsqueeze(0)\n",
    "q_logits = (torch.arange(K).float()/4).unsqueeze(0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "8zKQu6uLr9Lf",
    "outputId": "35dd0aad-3b26-421b-c079-7454163e1f8d",
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1.0 -167.9389398643856 52.324004913466865\n",
      "1 1.0 -155.5490832190536 58.144065714950116\n",
      "2 1.0 -175.50655121120153 59.31430416194602\n",
      "4 1.0 -174.3931358170284 57.931231495763356\n",
      "8 1.0 -171.2437877525198 56.5110724927502\n",
      "16 1.0 -166.5353348913696 55.54475491220046\n",
      "32 1.0 -168.85153891476392 57.14526575097915\n",
      "50 1.0 -169.4767376565317 55.98763226258352\n",
      "64 1.0 -169.31475324757358 55.113062645211464\n",
      "100 1.0 -164.41959151759173 54.019143479658695\n",
      "128 1.0 -165.26732003016073 53.82468731915931\n",
      "150 1.0 -165.04225421833783 51.671158613726725\n",
      "200 1.0 -163.61047298397025 51.035025649660696\n",
      "250 1.0 -161.26902005378014 50.115249191829996\n",
      "256 1.0 -159.01353535161397 48.389123530846824\n",
      "300 1.0 -159.76142380584668 48.23573796995071\n",
      "350 1.0 -160.69120360056476 47.56453717837677\n",
      "400 1.0 -158.93637818423574 46.71928901230769\n",
      "450 1.0 -157.81276452439903 46.53609680013909\n",
      "500 1.0 -158.0909916578559 45.89234425061499\n",
      "512 1.0 -157.43416324083364 45.26634913084891\n",
      "550 1.0 -158.1204285643671 45.74081610857735\n",
      "600 1.0 -156.9141194732875 45.23666733931128\n",
      "650 1.0 -156.8586592254624 45.25010335831754\n",
      "700 1.0 -156.92556481956836 45.01469659620476\n",
      "750 1.0 -155.03927477408624 44.55721698000532\n",
      "800 1.0 -157.15975105399386 44.91858328926049\n",
      "850 1.0 -158.09796597677655 44.97354693793965\n",
      "900 1.0 -155.5617853333851 44.699870913116584\n",
      "950 1.0 -157.22951817334027 45.145065348889084\n",
      "1000 1.0 -156.6875545501713 44.82030193832495\n",
      "1024 1.0 -157.59837031599463 44.85810511915518\n",
      "1050 1.0 -154.76051555127432 44.11968055966038\n",
      "1100 1.0 -156.87297367211121 44.67180876476678\n",
      "1150 1.0 -155.94493520518998 44.24772393407728\n",
      "1200 1.0 -153.88557666541502 43.69531531172071\n",
      "1250 1.0 -154.95565139279594 43.95512604270458\n",
      "1300 1.0 -157.92625303334577 44.55745594025657\n",
      "1350 1.0 -157.88450224850615 44.55992621542202\n",
      "1400 1.0 -157.38159467204937 44.2828321332011\n",
      "1450 1.0 -156.75150242435487 44.143808089869715\n",
      "1500 1.0 -155.88769704200877 43.825779094911994\n",
      "1550 1.0 -156.18577143683825 44.01710350872093\n",
      "1600 1.0 -157.463260993944 44.30975921668558\n",
      "1650 1.0 -155.862410261882 43.55259504986849\n",
      "1700 1.0 -157.71052758792263 44.209162766792176\n",
      "1750 1.0 -157.00496263952482 44.04210757934873\n",
      "1800 1.0 -157.7075515945058 44.158097065844125\n",
      "1850 1.0 -155.3105103883583 43.40403581769887\n",
      "1900 1.0 -155.1838096844505 43.32602595684878\n",
      "1950 1.0 -159.68145683367354 44.778886925639206\n",
      "1: {'loss': 0.6844744086265564}\n",
      "2: {'loss': 0.6906892657279968}\n",
      "4: {'loss': 0.745762050151825}\n",
      "8: {'loss': 0.5971688032150269}\n",
      "16: {'loss': 0.6934273838996887}\n",
      "32: {'loss': 0.6648648381233215}\n",
      "64: {'loss': 0.6260109543800354}\n",
      "128: {'loss': 0.6422197222709656}\n",
      "200: {'loss': 0.7151163220405579}\n",
      "256: {'loss': 0.6568819284439087}\n",
      "400: {'loss': 0.42272570729255676}\n",
      "512: {'loss': 0.3662804663181305}\n",
      "600: {'loss': 0.3243226110935211}\n",
      "800: {'loss': 0.4315519630908966}\n",
      "1000: {'loss': 0.3905391991138458}\n",
      "1024: {'loss': 0.252959281206131}\n",
      "1200: {'loss': 0.2858180105686188}\n",
      "1400: {'loss': 0.3500309884548187}\n",
      "1600: {'loss': 0.26849132776260376}\n",
      "1800: {'loss': 0.3126137852668762}\n",
      "2000: {'loss': 0.26454639434814453}\n"
     ]
    }
   ],
   "source": [
    "gd1_non_monotone = models.JointPredictorGivenLogits(K, 1024)\n",
    "gd2_non_monotone = JP_Z.Gadget2Jax(s_dim=K, z_dim=20, hidden_features=[1024,1024], tmp=1.0)\n",
    "\n",
    "#\n",
    "train_gsm_with_gadget1(2000, gd1_non_monotone, p_logits, q_logits, R_non_monotone, tau=1.0, lr=learning_rate)\n",
    "train_gsm_with_gadget2(2000, gd2_non_monotone, p_logits, q_logits, R_non_monotone, s_prime_obs=None)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "fKyqlS7nr9Lf",
    "outputId": "1885958f-91c3-4fa1-d7e5-ad96c95d145a",
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1.0 -708.3702203004975 257.28724122019713\n",
      "1 1.0 -714.2217068482007 290.8004172156075\n",
      "2 1.0 -753.3472956603676 299.76778768428403\n",
      "4 1.0 -708.1431515712784 258.0230642630415\n",
      "8 1.0 -730.0606884518941 272.48550828098683\n",
      "16 1.0 -746.6774640116587 265.24287909825273\n",
      "32 1.0 -727.4646328265603 268.18615157844954\n",
      "50 1.0 -728.1139133561257 268.38696042390285\n",
      "64 1.0 -721.2240044477591 265.44342458868744\n",
      "100 1.0 -742.7275416705138 276.1595332724956\n",
      "128 1.0 -738.357437253195 266.94377085422633\n",
      "150 1.0 -738.2012007359391 263.8609905896301\n",
      "200 1.0 -734.4423230492101 262.83799644136997\n",
      "250 1.0 -726.378975315848 257.9427140494324\n",
      "256 1.0 -717.7958512069263 255.59375785658526\n",
      "300 1.0 -718.2819254039861 252.68372884570383\n",
      "350 1.0 -712.2109944213058 250.38463151719677\n",
      "400 1.0 -711.2767435224341 250.20797012678813\n",
      "450 1.0 -714.559915158114 245.85613786038687\n",
      "500 1.0 -717.3158270961245 249.63337041872532\n",
      "512 1.0 -701.2417934659508 238.98988057754005\n",
      "550 1.0 -711.2260216975068 245.21129744098857\n",
      "600 1.0 -713.6530013371101 244.69355784121268\n",
      "650 1.0 -700.996323438405 239.144140988994\n",
      "700 1.0 -700.2268783181659 239.67689095378498\n",
      "750 1.0 -705.9910369460968 238.73798245140722\n",
      "800 1.0 -709.3956521477359 239.11058218801165\n",
      "850 1.0 -705.7881713671451 238.91753868789715\n",
      "900 1.0 -707.1674707328118 239.59287733938245\n",
      "950 1.0 -704.3141523760751 237.01176089281674\n",
      "1000 1.0 -706.6391898860466 238.29438875074487\n",
      "1024 1.0 -697.1903763094655 233.902683456727\n",
      "1050 1.0 -699.8460045452167 235.96670074362888\n",
      "1100 1.0 -707.0642030895808 240.36679100859467\n",
      "1150 1.0 -704.6669506711667 237.26963233913125\n",
      "1200 1.0 -700.247289106904 235.51633611016615\n",
      "1250 1.0 -700.9067773806497 236.3480851519661\n",
      "1300 1.0 -706.32892674509 237.89509205159786\n",
      "1350 1.0 -705.7739661041653 235.57342286999196\n",
      "1400 1.0 -699.51440516983 235.54461154725658\n",
      "1450 1.0 -702.1153791800501 234.3490190463634\n",
      "1500 1.0 -694.3719105069067 234.64385626327658\n",
      "1550 1.0 -698.4868298025486 234.85343254608233\n",
      "1600 1.0 -698.9776719438714 233.00517784427208\n",
      "1650 1.0 -697.9982864840376 233.5488430840647\n",
      "1700 1.0 -699.5397315327506 234.1597721981606\n",
      "1750 1.0 -698.8666657339589 232.91221500158414\n",
      "1800 1.0 -699.6735973734739 235.22347697513868\n",
      "1850 1.0 -699.2420766216197 233.4381610399753\n",
      "1900 1.0 -695.6983501751738 231.72481298467258\n",
      "1950 1.0 -706.7312199445464 235.0940361329063\n",
      "1: {'loss': 4.113441467285156}\n",
      "2: {'loss': 3.7668910026550293}\n",
      "4: {'loss': 3.6756746768951416}\n",
      "8: {'loss': 3.3424761295318604}\n",
      "16: {'loss': 4.172622203826904}\n",
      "32: {'loss': 3.753202438354492}\n",
      "64: {'loss': 3.8878116607666016}\n",
      "128: {'loss': 4.503621578216553}\n",
      "200: {'loss': 3.5473196506500244}\n",
      "256: {'loss': 4.609161853790283}\n",
      "400: {'loss': 3.6037137508392334}\n",
      "512: {'loss': 3.6930065155029297}\n",
      "600: {'loss': 2.837578535079956}\n",
      "800: {'loss': 4.042028427124023}\n",
      "1000: {'loss': 3.341987371444702}\n",
      "1024: {'loss': 2.1955745220184326}\n",
      "1200: {'loss': 3.5798184871673584}\n",
      "1400: {'loss': 3.021186113357544}\n",
      "1600: {'loss': 3.020780086517334}\n",
      "1800: {'loss': 3.120607852935791}\n",
      "2000: {'loss': 2.3093464374542236}\n"
     ]
    }
   ],
   "source": [
    "\n",
    "gd1_monotone = models.JointPredictorGivenLogits(K, 1024)\n",
    "gd2_monotone = JP_Z.Gadget2Jax(s_dim=K, z_dim=20, hidden_features=[1024,1024], tmp=1.0)\n",
    "\n",
    "train_gsm_with_gadget1(2000, gd1_monotone, p_logits, q_logits, R_monotone, tau=1.0, lr=learning_rate)\n",
    "train_gsm_with_gadget2(2000, gd2_monotone, p_logits, q_logits, R_monotone, s_prime_obs=None)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "I90imHEpr9Lg"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Df6y9oiTr9Lg"
   },
   "source": [
    "Let's compare our joint predictor to Gumbel-Max and inverseCDF coupling:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "id": "ThnqYt8_r9Lg"
   },
   "outputs": [],
   "source": [
    "def f_gumbel(u, logits):\n",
    "    # u is a vector of size size(logits)\n",
    "    return torch.argmax(-torch.log(-torch.log(u)) + logits)\n",
    "\n",
    "\n",
    "def f_cdf_inv(u, logits):\n",
    "    u = u[0]\n",
    "    cum = 0\n",
    "    for i, p in zip(torch.arange(logits.size(-1)), torch.softmax(logits, dim=-1)):\n",
    "        cum += p\n",
    "        if cum > u:\n",
    "            return i\n",
    "\n",
    "\n",
    "def get_couple(f, logits_p, logits_q, r, antithetic=False, indep=False):\n",
    "    u = torch.rand_like(logits_p)\n",
    "    ap = f(u, logits_p)\n",
    "    \n",
    "    if not indep:\n",
    "        aq = f(u, logits_q) if not antithetic else f(1 - u, logits_q)\n",
    "    else:\n",
    "        aq = f(torch.rand_like(logits_p), logits_q)\n",
    "\n",
    "    rp = r(ap)\n",
    "    rq = r(aq)\n",
    "\n",
    "    return (ap.item(), aq.item()), (rp - rq).item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "id": "4tmaTseMr9Lg"
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import linprog\n",
    "\n",
    "\n",
    "class MaxCovSolver:\n",
    "    def __init__(self, logits, r):\n",
    "        p_logits, q_logits = logits\n",
    "        probs_p = torch.softmax(p_logits, dim=-1).detach().numpy()\n",
    "        probs_q = torch.softmax(q_logits, dim=-1).detach().numpy()\n",
    "\n",
    "        self.obj = lambda r1, r2: (r1 - r2)**2\n",
    "\n",
    "        self.K = len(probs_p)\n",
    "        assert self.K == len(probs_q)\n",
    "\n",
    "        self.ij_to_idx = {(i, j): i * self.K + j for i in range(self.K) for j in range(self.K)}\n",
    "        self.idx_to_ij = {i * self.K + j: (i, j) for i in range(self.K) for j in range(self.K)}\n",
    "\n",
    "        self.bnd = [(0.0, 1.0) for i in range(self.K) for j in range(self.K)]\n",
    "        self.lhs, self.rhs = self.create_constrains(probs_p, probs_q)\n",
    "        self.psi_obj = self.create_obj(r)\n",
    "\n",
    "        np.set_printoptions(threshold=np.inf)\n",
    "\n",
    "\n",
    "    def __call__(self, method=\"simplex\", tol=1.4e-07):\n",
    "        self.opt = linprog(c=self.psi_obj,\n",
    "                           A_eq=self.lhs,\n",
    "                           b_eq=self.rhs,\n",
    "                           bounds=self.bnd,\n",
    "                           method=method,\n",
    "                           options={'ipm_optimality_tolerance': 1e-12,\n",
    "                                    'dual_feasibility_tolerance':1e-10,\n",
    "                                    'primal_feasibility_tolerance':1e-10}\n",
    "                           #options={'rr': False, 'tol': tol}  # , 'autoscale':True\n",
    "                           )\n",
    "        return self.opt.x\n",
    "\n",
    "    def get_opt_coupling(self):\n",
    "        \"\"\"returns (a_p, a_q, prob)\"\"\"\n",
    "        return [(self.idx_to_ij[m][0], self.idx_to_ij[m][1], aij) for m, aij in enumerate(self.opt.x)]\n",
    "\n",
    "    def create_obj(self, r):\n",
    "        obj = [0 for _ in range(self.K) for _ in range(self.K)]\n",
    "        for i in range(self.K):\n",
    "            for j in range(self.K):\n",
    "                obj[self.ij_to_idx[(i, j)]] = self.obj(r(i), r(j))  #(r(i)-r(j))**2 #-r(i)*r(j)\n",
    "        return obj\n",
    "\n",
    "    def create_constrains(self, p, q):\n",
    "        K = self.K\n",
    "        lhs = [[0 for _ in range(K * K)] for __ in range(2 * K)]\n",
    "        rhs = []\n",
    "        for i, eq in enumerate(lhs):\n",
    "            for k1 in range(K):\n",
    "                for k2 in range(K):\n",
    "                    if i == k1:\n",
    "                        lhs[i][self.ij_to_idx[(k1, k2)]] = 1\n",
    "\n",
    "                    elif i - K == k1:\n",
    "                        lhs[i][self.ij_to_idx[(k2, k1)]] = 1\n",
    "                if i == k1:\n",
    "                    rhs.append(p[k1])\n",
    "                elif i - K == k1:\n",
    "                    rhs.append(q[k1])\n",
    "\n",
    "        return lhs, rhs\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "id": "38sJk24gr9Lh"
   },
   "outputs": [],
   "source": [
    "def sample_from_joint(joint_pq, joint_qp, r):\n",
    "    u = torch.rand_like(joint_pq)  # bs, k, k\n",
    "    gumbel_noise = - torch.log(1e-8 - torch.log(u + 1e-8))\n",
    "    ap = utils.marginal_gumbel_log_softmax(joint_pq, gumbel_noise, rows=True, tau=1.0).argmax(-1)\n",
    "    aq = utils.marginal_gumbel_log_softmax(joint_qp, gumbel_noise.transpose(-2, -1), rows=True, tau=1.0).argmax(-1)\n",
    "\n",
    "    rp = r(ap).squeeze()\n",
    "    rq = r(aq).squeeze()\n",
    "\n",
    "    return (ap.item(), aq.item()), (rp - rq).item()\n",
    "\n",
    "\n",
    "def sample_from_policy(model, logits_p, logits_q, r):\n",
    "    joint_pq = model(logits_p.view(1, logits_p.size(-1)), is_pq=True)\n",
    "    joint_qp = model(logits_q.view(1, logits_p.size(-1)), is_pq=False)\n",
    "\n",
    "    return sample_from_joint(joint_pq, joint_qp, r)\n",
    "\n",
    "\n",
    "def get_opt_joint(logits_p, logits_q, r):\n",
    "    m = MaxCovSolver(logits=(logits_p, logits_q), r=r)\n",
    "    sol = m(method='highs-ipm') #, tol=1.4e-6)  #\"interior-point\", revised simplex,highs-ipm\n",
    "    #print(m.opt)\n",
    "    _, _, probs_opt = zip(*m.get_opt_coupling())\n",
    "    -m.opt.fun\n",
    "    k = logits_p.size(-1)\n",
    "    return torch.tensor(probs_opt).log().reshape(1, k, k)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "id": "cSvE-eG-r9Lh"
   },
   "outputs": [],
   "source": [
    "\n",
    "def compare_coupling_strategies(n, p_logits, q_logits, R, model, model_z, params_z, monotone=False):\n",
    "    logits_pq_opt = get_opt_joint(p_logits, q_logits, R) if not monotone else None\n",
    "    ates = []\n",
    "    couples = []\n",
    "    r = np.array([R(i) for i in range(q_logits.shape[-1])])\n",
    "    for i in range(n):\n",
    "        # Independant \n",
    "        (ap_i, aq_i), ate_i = get_couple(f_gumbel, p_logits, q_logits, R, indep=True)\n",
    "        \n",
    "        # Gumbel-Max\n",
    "        (ap_g, aq_g), ate_g = get_couple(f_gumbel, p_logits, q_logits, R)\n",
    "\n",
    "        # Inverse CDF\n",
    "        (ap_c, aq_c), ate_c = get_couple(f_cdf_inv, p_logits, q_logits, R)\n",
    "\n",
    "        # Joint Predictor - gadget 1\n",
    "        (ap_m, aq_m), ate_m = sample_from_policy(model, p_logits, q_logits, R)\n",
    "\n",
    "        # Joint Predictor - gadget 2        \n",
    "        (ap_z, aq_z), ate_z = JP_Z.get_single_ate(model_z, params_z, p_logits, q_logits, r, i, p_observed=None,\n",
    "                                                  method='crn')\n",
    "\n",
    "        # Optimal Coupling\n",
    "        if not monotone:            \n",
    "            (ap_o, aq_o), ate_o = sample_from_joint(logits_pq_opt, logits_pq_opt.transpose(-2, -1), R)\n",
    "        else:\n",
    "            # The Inverse-cdf is optimal under monotone functions (see the paper). \n",
    "            (ap_o, aq_o), ate_o = (ap_c, aq_c), ate_c\n",
    "\n",
    "        ates.append((ate_i, ate_g, ate_c, ate_m, ate_z, ate_o))\n",
    "        couples.append(((ap_i, aq_i), (ap_g, aq_g), (ap_c, aq_c), (ap_m, aq_m), (ap_z, aq_z), (ap_o, aq_o)))\n",
    "\n",
    "    return ates, couples\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "n = 1000\n",
    "p_logits, q_logits = p_logits.squeeze(0), q_logits.squeeze(0)\n",
    "\n",
    "ates_monotone, couples_monotone = compare_coupling_strategies(n, p_logits, q_logits, R_monotone,\n",
    "                                                              gd1_monotone, gd2_monotone.model,\n",
    "                                                              gd2_monotone.params, monotone=True)\n",
    "\n",
    "ates_nonmonotone, couples_nonmonotone = compare_coupling_strategies(n, p_logits, q_logits, R_non_monotone,\n",
    "                                                                    gd1_non_monotone, gd2_non_monotone.model,\n",
    "                                                                    gd2_non_monotone.params, monotone=False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def joint_from_samples(couples, size=10):\n",
    "    counts = np.zeros((size,size))\n",
    "    for (ap,aq) in couples:\n",
    "        counts[ap, aq] += 1.0\n",
    "\n",
    "    joint = counts/len(couples)\n",
    "    return joint\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.ticker as plticker\n",
    "def visualize_couplings(couples, logits_1, logits_2, R, suffix=''):\n",
    "    ai, ag, ac, am, az, ao = zip(*couples)\n",
    "    k = logits_1.shape[-1]\n",
    "    RiRj = np.array([(R(i) - R(j)) ** 2 for i in range(k) for j in range(k)]).reshape(k,k)\n",
    "    independent_coupling = torch.softmax(torch.from_numpy(logits_1),-1)[:, None] * torch.softmax(torch.from_numpy(logits_2),-1)[None, :]\n",
    "    gmx = joint_from_samples(ag, k)\n",
    "    icdf = joint_from_samples(ac, k)\n",
    "    jp = joint_from_samples(am, k)\n",
    "    jpz = joint_from_samples(az, k)\n",
    "    opt = joint_from_samples(ao, k)\n",
    "\n",
    "\n",
    "    couplings = {\n",
    "        \"Independent\": independent_coupling,\n",
    "        \"Gumbel-max\": gmx,\n",
    "        \"Inverse-CDF\": icdf,\n",
    "        \"Gadget-1\": jp,\n",
    "        \"Gadget-2\": jpz,\n",
    "        \"optimal coupling\": opt,\n",
    "    }\n",
    "    \n",
    "    loc = plticker.MultipleLocator(base=2.0)\n",
    "    ncols = 2 + len(couplings)\n",
    "    _, axs = plt.subplots(nrows=1, ncols=ncols, figsize=(4 * ncols, 4))\n",
    "    axs[0].imshow(np.exp(logits_1)[:, None], vmin=0)\n",
    "    axs[1].imshow(np.exp(logits_2)[:, None], vmin=0)\n",
    "    for j, (name, coupling) in enumerate(couplings.items()):\n",
    "        axs[j + 2].xaxis.set_major_locator(loc)\n",
    "        axs[j + 2].imshow(coupling, vmin=0)\n",
    "        if name == 'Independent':\n",
    "            axs[j + 2].set_title(f\"{name}:\\n{torch.sum(coupling * torch.tensor(RiRj)).item():.2f}\",  fontsize = 18)\n",
    "\n",
    "        else:\n",
    "            axs[j + 2].set_title(f\"{name}:\\n{np.sum(coupling * RiRj):.2f}\",  fontsize = 18)\n",
    "            \n",
    "\n",
    "    plt.savefig(\"{}/toy_couples_{}_{}.pdf\".format(figpath, coupler_save_path, suffix), bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABsEAAAETCAYAAABndqhPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABLpklEQVR4nO3deZxcVZ3+8edJd3YgEEB2CCqL/BgRjQjiICCK27jOKG4ILogDKu7bKLiMuw6M46iIiggiiggujCgi4oJAQHYEAoSwhC2EACF7vr8/zm1yU6nuruo+1XX71uf9etWr0rdOfetUpZ6qW/fce64jQgAAAAAAAAAAAECdTOh2BwAAAAAAAAAAAIDcGAQDAAAAAAAAAABA7TAIBgAAAAAAAAAAgNphEAwAAAAAAAAAAAC1wyAYAAAAAAAAAAAAaodBMAAAAAAAAAAAANQOg2AAgJ5ne3/bYfuwbvelU2xfaHtet/sBjCXbJ9uODtWeZ/vCTtQGAGCs8b0GVAd5RK+ryvaLYjvRyd3uR262Dyue2/6lZbXeLsYgGABg3Ch9KX+g233B4Gw/zfZxtmd1uy/oHNtTbP+77Qts3297pe2HbF9m+4u2d+12H4FexHdla2xvYvuTxWfWQ7ZX2L7T9s9sv8q2S23nFa/pwOVR2/Ntn2v73bY3HuQxYojLR8bsyWJc6YXv1+Jz6rjBsjPE/Z5v+1vFa7GscQMekBt5HPQ+U2y/3fY5xXfkUtu32j7d9lM62F3UiO1j6jrggurp73YHAABA7TxN0rGSLpQ0r5sdQWfYfqKkX0l6iqQ/SvovSQskbaD0//8WSR+wvX1E3NWtfgJAM7b3knSOpCdI+oWk0yQ9LGkbSS+W9DNJR0n639Ld7pT00eLfUyRtLWl/SSdI+rjt10XEBU0e7kpJX22y/O+jfR6onx76ft1faV3xZEkPtXG/N0h6vaRrJd2g9JoAHUEehzRL0omS/izpu5LulvRESe+U9CrbL4yIP2TuJ+rnGKXtBSc3ue0FktxkOTrnIklTJa3sdkc6gUEwAAAAtMz2VEm/lvQkSa+KiJ83aTNF0nsldWQqQgD1YHvDiHhkjB9zS0m/VBrIem5E/LmhyWdsHyxpk4bliyPi1IZln7b9XKWBtHNs7xkRcxva3NXkfsB6+H5tycclvSMilhdHuz6ty/1BTZHHYd0vac+IuLK80PZpSjt5fFnS7C70CzURESu63YdeExFrJC3rdj86hekQAQDjlu1ZxTQox9l+aWlqlAW2v2x7vZ09bL/c9t+LdnfY/oykiYPUn2z7Y7avK9o/ZPuXtvdsaPf43Mm232X7pqL9TbbfNUjtnWz/sOjrimIaiS/bnt7Q7uSi9gzb37R9X1H7L7af1aTuJra/Y/sB20uc5tJ+xhCv4WzbPy/aL7d9o+2PN752RZ15trcuprlYZPsx2+fZ3rnU7jhJ3y/+/IPXTvt0cqnNk+owdUgPe5ukXSV9udkGAUmKiGUR8fmIuHtgmYc4P1eT90g526+xfWUxzcpc24cXbba3fabtB20/YvtU2xsOUn9z26fYXljk4ve2nz5I29fa/nNR8zHbl9j+15ZfneY1B+Zcf57T9Gu3F8/nEtt7F22eWzzukuJz4RNN6rzA9hlO080sLT6Tfuu0Eb7cbs/ic+J39jpTuvXZvqh4XruXlm9le1fb00bzPFFNbuO7snh/rbC9aZM6uxR1jm9Y3lJmBnJe5ODPth9VGoyS7Zm2/8v2LUXfFtq+3PYHm9QZbUY/qHQE2IebDIBJkiLivIj4cSvFIuKPkt6vdGRA21McFp9lu9puui6CnjLS79etbX+1+K5cVGToetsftt3XWMP2drZ/Ynux7Yed1m2f1Ozxiu+NTxTfW8tsX11k8Lgi07Ma2m/ltL46v/gsudv2ibafUGpzstJRJ5J0m9euKx433AsUEXdFxPLh2hWPM7HI1vattAcakMchRMTCxgGwYvn1Skdq7t54W5HHps8N1WZ7M9vfcNp+sqK4/kbj+qLX/uY5qHhf3u60jeFq24c0tA1JO0h6rtedLnpWcft65wTz2m0Ss5y2YTxU5Oxk2xvYnuC0/ea2IiNX2N63ocYEp+0dF9m+p3g+84usrLf+2+brZKdpQi9xmjr7UdvX2P70CF/PptkubpvnhvMGeu269kG2/+a0nnyP7RNsb9BC/x/frtVsme3DnbaNLS/+bz80SJ13Om1XWm77ZttHu8k5yIq2Y/a5wJFgAIA6eLGkf5f0LUnfk/RySR+QtEjS5wYa2X6l0hRH8yR9WtIqSYdLekljQaeNUb+R9GxJP5T0P5JmSHq7pL/Y3i8i5jTc7V2StpT0bUmPSHqdpP+2PTMiPlWq/QxJFyhNN/FtSXdJ2kPSuyXta/u5EdF4CPp5SnvcfVrSppLeJ+nXtncc2Iu+6PN5kp5Z9PlvSnvIni9pYZPn+BJJZ0maqzRV04OS9ike42mS/q3hLtOVDpH/m6SPSdpR0nuU9n7fPSJWF/W2knSE0mt/Q3HfW0p1fq+0wsv0BuPTwMbmk8bgsV4q6UilKckelPRWSd+zvULp/XWB0nvxmUpT0ixT2mjR6DfF/Y9TyujRkv5oe5+IuHagke3PKu1l/htJn5C0RtIrJf3U9tER8Y1RPp8vSOpTmj5tktKG89/aPlRpKpkTlaZle43SESa3NRxBcpikmZJOUZqabZvi+f7e9gER8SdJioi/Fz9KTlDaKP/54v7HSvpnSUeWn3dx+5slHaA0jSnqqZXvyh8ovf9ep/S9V3ZoqY2kEWVmtqRXS/pOuY6kn0rar+jb1UpTsTxFaYqmL4/i8Zp5taQVDY8/WgPrCS9ucttE25s1LFsTEQ8W/z5F0nOVvlPnZewTxp+Rfr8+VdKrJP1caX1roqQXKn3nPFHSOwYaOp3z5yJJ2ynl7Xql998flHLX6H+Uvof/IOkrkjZX+k6+rbGh02DTxUrfb98t+vJkpenRDrA9OyIWK637bqSU3fdKeqAocXWbz3s42yith/5R6bMEaAd5HAHbE5R+C97b5OYbJN2uNJUixgnbMyT9Ven98z1JV0jaU+m9dKDtvZoc1f9FpW0HA9NKHy7pdNtTIuLkYtmblKYYfUDSf5bue/8wXZqu9Bvwj0q/cwZ+B05R2ubxLElfV8reByT90vYOpT5OUtoh6mdKU2MvKWq8VdJzbD9jFEeh/VBp2t5Liuf0kNJg+r9K+qQ04tezHU8vHu87SuuYByhtZ9rd9vOLo71G4khJWyh9njwk6Y2Svmj7zoj40UAj2x9W+ry7Qmka8WlKr/dg/69j97kQEVy4cOHChcu4uCj9gA1JHyj+nlX8vUTSrFI7K+2BtqC0rE/SfKWVrM1Ky2cUX7oh6bDS8oGpLQ5u6MNGRZ0Lm/TrEUnblpZPknSp0pzK5eVXSfqHpA0bar+yST9OLpb9b0PbfyuWv6O07Ihi2aca2h5TLJ9XWjZF0j1KP7z6G9oPPPf9S8suLJZ9qKHtBxtfJ6UN9evcv+E+89IqSPffU1zavyj9uFjcZHmfpM0aLlNLt5882P978X45ufR3Ods7lJZvrjTQtUbS+xpqnKW0YXuDxscsbnNp+TOKGr8pLXt60fZzTfp3ttL5gjYsLZtX/hwY5jUbyMQVkiaVlr+sWL5S0uzS8klK55y4uKHO9Ca1t1D6XDu3yW3nFLX3VvqcWi3pzCbtBl6n/bv9/uIy+otG/125QNKlDTWt9F15dWlZu5mJ4nJQQ9sZavI916RmW483SI0NixpXD9Wuyf3mSbp2mDZXF7WbPefGyz2lNhcWy2a10ycu9bto5N+vU1X6jist/2Hxub9Vadnnivfb4Q1tjy+WX1ha9v+KZb+RNKG0/J+Kuuu8b4vvnPtUWuctls9W2vHsuNKy40b7vlfauDnU+uasxufEhUurF/I44tft34tan25y2zq/R7mMj4vSYE5I+veG5UcVyz9TWnZYsex2STNKy2cUyx5syMu8wT6jVZxfvMmykPTBhuVnKf22myNpYmn5wG+t8jYTl/tQWv7Wou1rGpav8zt1iNfpNUXbH5YzWtxWzmw7r+eg2Wz22mnteuYrGpafUCw/pMn/1f6lZftr/e1RA8vubvg/naY0sHVxadlMSUuV1omnlJZvKWlx4+OV+jyv8fl14sJ0iACAOjg7IuYN/BHp2/QPkrYsHfb9DKW97L4fEQ+U2i5W2vOu0RuVBqouLw5X36zYk3uSpN8p7SXUuIfeaRFxZ6n2CqW9m/ol/Ysk2f4npT0EfyRpckPtPyttpHxBk/78V8PfFxTXO5WWvULpR9BXG9p+U2njYNnzlTaef1/Sxg39OLdo09iPNZL+u4V+DCkiZkUER4GNXxtp/feTlI7auL/hctQoH+vsiLh94I+IuF/SjUrvxcYjPv6ktMffrCZ1vlR8LgzUuVwpxweVPiPeoLQS/oNyHopM/EJp4/k+o3w+34x19yz8U3F9SZSOLC3aXKqGXEXEkoF/O035salS5i9R2uux0eFKg92nSzpV0h1qcqRcRBwWEY6IC0fypDBuDPtdGemI3tMkPdPrTlu7v6Ttte7RUyPJzFURcX7DsqWSlkt6VrPpXkb5eI02Kq6bfYaN1kDNjRqWX6L0nVu+PH6kdUTsX+RvXgf6hPFlRN+vEbF04DvO9iSn6UU3U5odYILWPS/PK5SO0Dil4TG+2ORxX1pcnxClPccj4pqi9uOKPdtfqpTFZQ35nKc060Cz9duOiYh5Rbb2H8vHRW2QxzbZfrakrynt8Pm5xtuLPM7K/bjouFcqvc9PbFj+7WL5K5vc55vFdhZJ62xz2USjPzJ3tdKRXmV/Uhrc+lasO6POwG+tx39TRbJUenyK0Y2LbAxs12j2m6oVbyiuPxANR1s1/D2S17MdN0bE2Q3LvlB67JH6fsP/6WNKMwSVf68+X2ln629GxLJS23uUfl+sZyw/F5gOEQBQB7c2WTYw/d+mkh5Vmn5CSgNbja5vsuwpSnvyDXU4/mZKG5UH3NCkzUDtgcd/SnH9qeLSzBZNlq3zHCNiodOpfsrzRj9RaY/+hxvaLrd9q9JK54CBfnxvkD4068fd5ZWZQvl1Rm94WOtv5JXSVCzPL/69h9I0LaPVLNuLlN7njecEWVRcN3svDpbNFyhNzXmdUias5p8RA5plU9LjGxwaB8bvLwYVBjTmeFGR49ualFykhudSzJf+n5IOlrRxQ/tQ44KIB53Oofa7YtFzIuKhwZ4Daq+V70opDXS9X2n6w48Vyw5V2uhQ/gE7kszc1NggIlbYPkZpL9XbbF+vtCHi7Ij4/Ugez/ZMpZ1Wyo9zj9Zu0Gx6/sBRGmyA7YEmA39AMyP6fnU6r99HlHL6ZK0/3XR5/e+Jki5r+G5SRCyw/VDD/XYsrm9s0qcbJb2o9PcuShv431pcmmn2GbQO25OU9iQvezQiHm3WHugg8thGHovp/n+tdLTIS5r8ZsT4taOkORGxqrwwIlbZvknpSP1GrWwXGakFTd5fA78D1/lNVfqt1fib6jVK67p7av3zw2+ikdmp6FuzqUDLRvJ6tmO91770mTKa136w3xHl13a4z6muYhAMAFAHq4e4baRHHFnSNUrn3hrMcPNVD1ZXSkdr/WaQNosaFzT+OGpSb6T9+KCkKwdpc3fD3514nTH+XCtpv+J8dI//0CiOUjpfkmyvanK/9QZpirZDrY8O9p7rVOZDaSPCYPWvG+L+JyidV6us8Rw/I3k+qXPpSJ2LlObBP17p8+kRpaPiPirpwEHu+rLSv58m6S/DPRZqq6XcRMQ1tq+U9AbbH1ca3H21pN8WA0nl+7SbmceaNYqIb9k+R+kcnc9VOpfB0bbPiIiBE6m383hnFXXWeY4R8Yjt2yXtanvqwJ7Ao2V7sqSdlTZ+jOY8DuhtI/1+/ZrSeWnPUNpR4j6lqXCfrnREyVjMADTwGXKqBj/fXit5e7bSEapln1KaEgoYS+SxxTzafrrSDleLJR0QEXe13FOgfUOtzw67zcT2q5TyeanS+c3vUJpuv0/FdKN5uplF09/PhbEe0xn292rVMQgGAOgVA3uu7Nrktt2aLLtZ6fxDFzQezj6EpzRZNlB74PFvLq5Xd2DP8FslvcD2RuWjwYqNc0/UuoNrA/1Y0oF+DLWyhvHvTEn7KU2r9/E27veglI7QiIgHS8tHuzdgK56iNF1D2W5KK/MD0y3erHTi8vkR0WzvxeF8SWljQ9k9zRqO0PMkbS3pLRHx/fINtj/b7A62X6a0Ieb7Snsjf8X2RcXUOcBQfqA0De8BSie431Drb0gbbWbWERELJJ0k6STbfUrnVHid7a9GxGVtPt77NfievGcpnfvyTVp/KpqRepOkyUp7wQMjNdLv1zdJuqg0YCxJsv3kJm1vlbST7b7yDla2t9L6RxjPK6530fp7gO/S8PdcpfW/SS2uVw62rniV1h5lU+4zMNbIYwt5LAbAzlfaMeuA8jTmqI1bJe1iu7989FKxI+POav4Z/RSl89KVNW4Xkbqz3eBNSoNeBxRT+kmSGqYBH4mbJL3c9hbDHA3Wzus58Jt5pko7VtqeorR+PrdJ/fW2SZU+Uzr9fTqvuN5Fa6eXVGlZV1VpdBMAgE66XNKdkg4v5nyWJNneSNKRTdqfonQCz6ZHgtluNi3aG2xvW2ozSWlD22pJvyoW/11pz8Ijba+38d92fzGN00ico7QH0/sblr9T60/ncZ7SnokfafZ4tqfaHul0UQNTZDR9HraflGElE91zktJ0ZB+0Pdi84s2OxhqYBu2ghuWN79dO+JCLuTCkx3+wHyTp96UpXX5YXH+u2AC/jkEy/7iIuD4izm+45JwKZmDjyDqvre0XqMnc9ba3UZru9AZJRyvNU79U0o/dcD5D21vZ3tX2tIz9xfj2I0mrlKZzOlRp7+7GjRmjykyp3bTG916xMfDq4s+B75KWHy8iLm/MY6npl5SO5P6S7abnELP9AtuHNLutSdvnKh3d/Yikz7dyn4b7b1/kr3E6HvSekX6/rm5cbnu60jpoo3OUpg09tGH5h5u0/WVx/R7bj287cjq/7cHlhhGxUOmcsq+yvfd6nU42Ly1quq4YEYuafJeOaKOd7YlFtrYfyf3R88jjMHm0vafSEWCPKg0oNJveu/y4uxZTe2N8OVtp5+DG8wq/vVj+8yb3eWcxVbykx6eNP1LSQ5L+WGr3qAbZZtBBq5UG38o5sqT/GGXdgSnDv1TOaKn+gLPV+us52O/n92rwMZ1dbL+iYdnAZ8rZg9wnl98pnef3ncVAnSTJ9pZae860dYzl5wJHggEAekJErLb9Xkk/kXSp7e8obeB7i9Jcxo0/kE9Q2vPty7YPVNqT5eGi3fNU7D3UcJ+bJF1i+1tKG8NeL+mZkj4TEXcU/QjbbyrqXW37e0rTN01TOlLjVUpTm508gqf5fUlHSPqk7R0lXaw0z/W/SbpFpe/9iFhi+1ClFaEbi37MVdpDaNeiH6+UdOEI+nGZ0hRtH7e9iaQlkm6LiEuK23+vdB4mplAchyJiqe2XKA3snmX7Qkm/VTrqaSOl989rlX5glM+Zd7rSSbJPLAZBH1Q6qmMzdd4Oks6z/QulveaOVhoQ+uBAg4i4zPZxSlO8XGn7p0pTgm4l6RmSXqyGcwyNsT8rvcZftT1LaVD/aUp7M14j6Z8GGhY/vE5TmjrxtcVejo/ZfqvSUTDHS3pHqfbnlaZyPEAjyzxqJiLus/1/StMSTpH03cZB3YyZ2VnSH23/XGknkUVKe7G+U+ncDn/K+XgRcY/tlyptfPyz7bOVphp9WOloyxdKek7x+GUzbL+x+Pfkou0BSid4v0/SISPcWH+K0tSNjdOnoseM4vv1TEnvsH2G0hEZW2jt+m2jLymtn37H6Rw+1ym9h/eR9EBDf66zfaLSuuX5RUY3l3SU0k5dz9C6e/G/U+m76iLbpxRtJigd8f1ypff6cUXbgaOzv2j7NKX16msj4tqhXiPbT9XaaX73La7fZPs5xb+/HhGLi39vo7QjyB+L5wi0jDwOnUfbOyht8N5E0n9LerbtZzc0+3kxfeSAG5RmYJg1WF1U0peUtid8o9iR8O9K2xjeqnSepy81uc8DSttFBmavOFxpO8rbykdfKb333mr7M0rvjzWSftnwvsntTKVpvi8osjFR0iuUtseMWET8tMj9oUpHeP5CaZ12Z6WB6t2Lpu28nucXyz5te1Ol9eLnSNpbDZ8RJddIOrXY3nWz0rrqvyp9F54xmuc4nOK89Z9S+s3/F9unKr2uRyhtK5ut9Y/+G7vPhYjgwoULFy5cxsVF6UdBSPpA8fes4u/jmrQ9rrhtVsPyVymdA2u50g+WzygNdoWkwxra9kt6t9KgzpLicrPSxuUXNOnXYUX7m4v6N0t6zyDPZQdJ31La4LVC6YfR5Uobo7crtTs5fV03rRGSTm5YNlPSd4t6S5Q2aM8uruc1qbG70hRudxX9uFfSXyV9QtLMUrvB7t/0/0Bpg/r1Rc11+lk856bPicv4uSidJ+gopXMFPKB0voOHJM1RWnnfpcl9nqV0TqplxX1OVBp4bXyPDJXtwd6LhxX32b+07ORi2eZKR5EsVDon0QWSnjHI83qJ0pGSD2rt58T/STqyod08SRe2+Fqt17fSbevluNz3hmVPVZqrfpHSQPuFkv65sa2kTxZ139mk7jeK217d5HVar39cxt9FGb4ri9teXdwWkvYd4vFazcxg7/VNlaZevLL4DFmqtFPG8ZK2GunjtfA6zZR0rNJn1mKl76s7lTaOvKyh7bzSaxHF58jA475b0saDPEZI+tUw/bhwsP8DLr15UZvfr0obmL6stBFpmdL650eUdtpqtn67ffE+f7i4/FLSk9Tke01phoFjJc0v8na1pNdI+kpR+wkN7Tcr+nJT0ZeHlDbInSBpt4a2H1KammnlYJ9RTV6bwxqy2HiZVWo7q1h24XB1uXAZ7EIeB31d9h8mi81+h4earMNzqf5F6bfU/yqtJ60srr8habOGdgOf0QcpnT9u4L16jaTXN6n7BEk/U1qnW1N+36jJb75myxoed/8mtzXbZvJ2pW0VyyQtUPpNOnOQtk3XXwd5nSYUnxdXKK0rPlLk9NiRvJ5F252Vfv89VmT4J0o7eTT7jAil33UHSbpEaZ36Xklfl7ThcK9ZKdeHDbWsdNvJarJdp3gNbtLabWJHK03RH5L2atLn9f5PO3Fx8YAAAGCEbO+v9MPo8Ig4uaudAQAAADrI9i8lHShpoyidywjA2COPqArbhynNTnNARFzY3d70Htsh6QcRcVi3+9LI9teVBsO2ioic581uGecEAwAAAAAAwDoazyFZLHuqpBdJuoAN7sDYIY8Aqq58LrDSsq2Upom8tlsDYBLnBAMAAAAAAMD63lycQ/bXku5XOg/SEUrTh36ymx0DehB5BFB1+9v+stJ5qO9Umpr47ZI2UJoatmsYBAMAAAAAAECjKyS9UuncdzOVzm9ygaRPRcTfu9kxoAeRRwBVN1fSLUoDX5sqnXdtjqTPR8T53ewY5wQDAAAAAAAAAABA7XBOMAAAAAAAAAAAANQOg2AAAAAAAAAAAACoHQbBAAAAAAAAAAAAUDsMggEA0CG2p9m+1XbY/p8W7zPR9sds32B7ue2Ftn9me9dB2m9v+9u259peavsu27+0vV/eZwOMb+3m0fb+RduhLvuW2r/E9i9sz7P9mO1Ftq+wfYztKZ19dkC12d7F9mnFd9viIiP/sP0121u1WOM1tr9v+yrbK4sMzhqi/da2T7F9f/H9OMf2v2V7UsA4lSmPtv1623+1/YDtR2xfZ/uTtjdq0v7/2f5R8R25zPbttk+3vUf+ZwiMXyNYX51o+1u2Ly+yuNz2bbbPsL3nEPfbrcjkguI+d9r+ue0t8j4jYHxgXbX++rvdAQAAauzTkjZvtbFtSzpH0osknS3p68X9/13Sxbb3jYjrS+23lnS50vf5tyXdLGlrSW+X9AfbL4uIX+d5KsC411YeJd0g6U1Nlk+WdKKkByRdWlr+T5JWS/qupAWSpkr6Z0n/Jekltl8QETGCfgN1sK2krST9XNKdklYpZeYISYfYflpE3DdMjX+X9CxJV0m6RdIugzW0PVPSnyU9QdLXisd8vaSf2H5LRHx/dE8HGNdy5PGzkj4m6QJJn5K0UtL+xb9fbHufge+8YqDrYkmLlL4/75T0JEnvkPTKou3fsz5DYPxqd311kqTZkv4i6YeSHpG0vaTDJV1i+4URcUH5DrYPVvqteYuk/5Z0r9L35T6SNir+BnoN66o1Z36LAwCQn+2nK20g/5Ckr0r6RkQcPcx9XqG00nViRLyjtPyJkq6V9NeIOKi0/KOSPifpFRFxTmn5k5UGxM6JiFfkek7AeDWSPA5R63WSfiTpKxHxwRbaf0PFD6KIuHS49kAvKfZ2/YmkD0fEl4Zpu72kuyNiVbF3/FGSdoyIeU3afknSByW9LCJ+WSzrU9oQ/yRJO0TEo1mfDDDOtZpH2/2SFkv6h6RnRsSa0m2nSnqDpD0j4spi2beVNiI+LSKuKrU9SNLvJJ0QEcdkf0LAOJN5fXUrSfMl/S4iXlxa/gSlHb0uVfqOXDnqjgM1xrpqfTAdIgAAmRUrL9+R9BtJZ7Vx1wOK63X2+omIWyX9SdLzihWrAQPTzdzdUOceSWskLWnjsYFaGkUeB/O24vqkFtvfXlxvkuGxgbppOR8RMT8iVrVY9/WSbhnYqFDcf7XSEdYzJb14sDsCPazVPE5UOtr5nvIAWGFgnbS8DjrY+mqztkBP6sD66n2Slmn9PB+p9D34oYhYWUy/ODHD4wF1xbpqTXRkOsTNZvbFrO3yfobOu2OlHnhwtbMWBXrApL6pMXXijKw1l65crBWrl5JHYHDvlbSrpFe3eb/JxfVjTW4bWPYspb36JOm3kj4i6X9tf1Brp0P8pKRHlfYgBHrdSPO4Hts7Kg1W/zkibhykzYZKWd5I0r6SPixpoaRLRvv4wHjndH68DSRNkbSbpC8WN52b8TG2krSNpNOa3Py34vqZSnv1Aj1rpHmMiKW2L5L0QtsflvQzpWmj9lc68vnUiLi5dJfzJB0i6Ye2j1Wa8umJxeMtkPStXM8JGMdGtb5aDKJtorSddztJH1DKd2OeXyzpYUkb275S0h6S1tj+q6T3RcRlI+o9UBOsq9ZXRwbBZm03UZeet13WmnsdfEfWekCvmDpxhvbZ/tCsNS+ef0rWekCdFBvJPyXp0xExb6gToTZxXXF9oKSrSzWnKQ1+SelHjSQpIv5g+yilueP/UKpzs6S9I+KG9p8BUB+jzGMzb5FkDX0U2Pe17gaMSyQdFREPjfKxgTp4m9IergPmSXpjRPwp42NsXVzf1eS2gWXbZHw8YLwaTR7fIOlkSV8oLpIUkv5TaWessh9I2lHS+7R2456UpmObHRGNR4gBPSXT+upTJF1T+nuxpM8Xl7JdlLYF/0bSTyV9RtIsSf8h6ULbe0XEdQJ6F+uqNdWRQTAAAHrYtyTdqnRy03adqvQD5NO2l0g6X9JmSj+KNivaTGu4z/2S5hRtb5K0s9Lc0r+2/dyIYC8S9LLR5HEdxR62hyntPfvTIZp+qnjczZWOGnuqpE1H+/hATZytdB6hDSTtKellWvv9lsvA9+TyJrcta2gD9LKzNfI8Lpd0m6RTJP1fsezVSuuxy5QGwyRJERG275H0F0m/UNrA9zRJ75d0ju2DImLxKJ8LMJ7lWF+9TdLzJU2S9GRJb5Q0Q2l2gvL0bBtK6pN0WkQcNrDQ9uVKO1V+UtJrR9EPYLw7W6yr1hKDYAAAZGL7jUo/PvYbyUmGI2JRcZLwUySdWLrpj0qH4f+H0gb4gcd7u6T/VTr5+LWl5edJukJpz783juCpAOPeaPPYxMGStpX07YhoNmWpJCkirtHaPXFPt/0OSf9ne7+I+EuGfgDjVkTcqTQVmiSdbftnki6zPS0iGvdWH6mBfE5uctuUhjZAzxppHosZCv4q6YqIOKR0049t/1hpZ64zB6YNtv1ZSUdL2jUi7inanmP7EqUBtA8qreMCPSfX+mpEDOxAOVD3e0q/B38m6YWlpkuVNu6f3HD/C23PV5rWFOhZrKvW14RudwAAgDqwPVlp771zJd1j+8m2nyxph6LJjGLZxkPViYhrImJPSTtJeq6knSJif61dQfpHqflHJf2jPAA2UKNo99zRPStgfMqVxwZvLa6HmgqxmVOL6yPbvB9QexFxtaS/K51HKJeBqdWaTSMzsKzZ9DNAT2sjj/+qtJ7a7KjonyptZ3qOJNmeqHRuoj+VBsAGHu83kh4R66voUR1aX5UkRcSjks6SdLDtJ5VuGti4f8/699ICpfOKASiwrlofDIIBAJDHVKXpz16idE6ugcuFxe1vLP5+WyvFImJuRFwUEXOLRS9SOgqsfCTJNkrTWTTTL474Ru/KmkfbT5D0L5Kuiog5bfZlktI698w27wf0iqnKmI+IWKC04WDvJjcPLGs3x0CvaCWPAxvomq2D9jdcb6a0I9d6bW27WM76KnpV1vXVQepL62b60uJ62ybtt5V03wgfC6gz1lVrgJUNAADyWCLp35os31xpysLfSPqupKslyfZWSvO0zx9qarWi7bsk7S7pU8VUFwOul/RU23tHxN9K7fdROjfYr0f+dIBxLXceD5U0sbhPU7a3bNzLvfDu4vpvTW4DesJg+bB9gNL324WlZS1/Pw7hdEkfsP0vEfHLom6fpHdJekhpr3ugJ2XI4/XF9Zsl/aShzJuL68uK63slLZS0n+0dI+K2UtvXKJ3z5DIBvWnU66u2N5e0MCLWlAvY3rKo/aik60o3/VBpvfbIov5A+39RGuD+To4nBow3rKvWH4NgAABkUMzhfmbjctuzin/eEhHl2z+vtKHgAK27QnWu0omRr5cUkl4g6RVKA1r/qXUdpzTNxe9sf0tpT8GdJL1T0gpJnxrVkwLGqVx5LHmr0kmKT21y24Brbf9Z6fwLdynt/f58Sc9TOkfY8e08B6BmvllsMLhA0u1K5zt4hqRDlKZDe3+p7WDfj/tJ2q/4c3ZxfbTthyQpIj5bqvEFpY1/P7L9NaVMvk7SMyW9LSIeyfjcgPFmtHn8ldLRJC+2fZHSuqgkvUrSP0v6aURcIUkRscb2cZK+LumSYn31TklPUzq65QFJX+nEkwSqLtP66hskHWP755JuU/oNuHPRbhOl77zHN9JHxPm2T5f0uuJ356+Upl98l9J0iMdlenrAeMO6as21NAhm+4WSTlA6VP2kiPhCR3sFYFDkEai9iyW9VtJhxd83SDpK0rcjYnW5YUScY/v5SicUf4vS3kiLJJ0n6TMRceUY9RmoLdvPlrSrpB9FxKIhmp6gNGh9lNJ0GUsl3SjpY5L+u+EoTqDXnK605/mblPZwD6UNDN+W9OWImN9CjQMlHduwrLxB4vENCxGx0Pa+ShsYjpK0gdLOJYdExBkjfRJATYwqjxGx2vZBSuemfZWkLxY1bpb0YaVzHJXb/4/tBUob2d+jdPTX/ZJ+LOnYFvMPoLk/KW1sf6mkrZSm4b5X0vmSToiIvza5z6GSrlL6/fhfkhYrDcZ9PCLubtIe6AWsq9acI2LoBulQvJuU9mS9U+lQ9ddFxPWD3Wf2HlPi0vO2y9lP7XXwHZpz1TJnLQqMMyPJ44wpW8Y+2x+atR8Xzz9Fi5fdQx4BAAAAAAAAAJU1oYU2e0maGxG3RsQKpb11Xt7ZbgEYBHkEAAAAAAAAAKAFrQyCbSPpjtLfdxbL1mH7CNtzbM+5f+HqxpsB5NF2HlesXjpmnQMAAAAAAAAAoCpaGQRrSUScGBGzI2L25pv25SoLYATKeZzUN7Xb3QEAAAAAAAAAYMy1Mgh2l6TyCb62LZYBGHvkEQAAAAAAAACAFvS30OYySTvZ3lFpY/shkl7f0V4BGAx5xJiZ1Dc1pk6ckafYhGwHHismOF+tvny11vTnrJWtlGJivlqavCZbqRmT803Vut3Ex7LVuvzq5Q9ExObZCmbSN3169M+cmaXW5DuXZKkjSe7LePR/xmzLGfM4NV+IVk3J2K8pka3WzOn53hPb9OfLdlXzOGni9JgyeeMstdZMyvf9qMj3nlhd0ffqlAX5pt2P5cuz1eoFj2hRNfO48dSYtuWGWWote3RyljqS5Cn53qtr1uTL49bTF2erdd/yPK+7JK1akW99YueN7s1W69HI957YdMKqbLWq+v3YP216TJyRZ3110n0ZT+vQn3F9NeO67+opGWtNylZK2iDf59fUiSuz1VqyLF8eJ/Tl+1277JYFlcwj66vtYX21HoZaXx12M1tErLJ9tKTzJPVJ+l5EXJe5jwBaQB4xlqZOnKF9tj80S62Ylm+FdfUG+WqtmJHv18LyTfL9iHls83wrmUu3yrcyp1n5NpS/eKfrs9U6fqs52Wr1bTX39mzFMuqfOVPbvO+YLLWe9P6/ZakjSX0bZRool+SpGafvnZRv4GrJ7ltmq7Vwt3z9enSnfBsV3rBXvvfEZ59wTbZaVc3jlMkb61lPPTJLrUe3n5aljiT1Lc+3QWfRLvn2xnh01xXZaj3lC4uy1Vp90y3ZavWC8+PMSuZx2pYbar/vvCZLresvfmKWOpI0ZdeHstVa8siUbLU+vdcvstX6+i0HZKt1790bZ6t11vOPz1brb8t2yFbrDRsuzFarqt+PE2fM1Ky3vC9Lre3/+6osdSRpwmZ5BuYkac3GG2Sr9ciTN8pW6+EdMg707Zfvu3b3ze/JVuvim/J9Rm+4Sb6dKK97+WcqmUfWV9vD+mo9DLW+2tK7JSLOlXRuth4BGDHyCAAAAAAAAADA8DIezwgAAAAAAAAAAABUA4NgAAAAAAAAAAAAqB0GwQAAAAAAAAAAAFA7DIIBQA+x/ULbN9qea/sj3e4P0MvII1Ad5BGoDvIIVAd5BKqBLAKjwyAYAPQI232SviHpRZJ2k/Q627t1t1dAbyKPQHWQR6A6yCNQHeQRqAayCIxefyeK3rFymo5ZMDtzzYVZ6wE9Y8IExbTJ2WtiXNpL0tyIuFWSbP9Y0sslXd/VXgG9iTwC1UEegeogj0B1kEegGsgiMEpsyQaA3rGNpDtKf99ZLAMw9sgjUB3kEagO8ghUB3kEqoEsAqPEIBgA4HG2j7A9x/acFauXdrs7QE8r53H1kiXd7g7Q08p5XLmSPALdtM766kOsrwLdVM7jqsf4fgS6ifVVYHAMggFA77hL0nalv7ctlj0uIk6MiNkRMXtS39Qx7RzQY9rKY9/06WPaOaDHtJXHiRPJI9BB7a2vbsz6KtBBbeWxfxrfj0CHDJtFifVVYCgMggFA77hM0k62d7Q9SdIhkn7R5T4BvYo8AtVBHoHqII9AdZBHoBrIIjBK/d3uAABgbETEKttHSzpPUp+k70XEdV3uFtCTyCNQHeQRqA7yCFQHeQSqgSwCo8cgGAD0kIg4V9K53e4HAPIIVAl5BKqDPALVQR6BaiCLwOgwHSIAAAAAAAAAAABqh0EwAAAAAAAAAAAA1A6DYAAAAAAAAAAAAKgdBsEAAAAAAAAAAABQO/3d7gAAoKIi5OUru92L9fRlrDUpY63c1fLJt7/LUk3PVutc7ZatVl5zu92BpibfuURPev/fstS66Vt7ZakjSbscfUW2Wn1Tp2arpRX5PrumzV2UrdaqKZtmq9W3bGK2Wj/e4BnZak2ZkPN7o5p51JKl0t+uzlJqgzyxzi5jGrOa+/l9stV638vzvfjnPHNWtlprlizJVqsXxLw+rXj7Bllqbbh/ljLJrRtnK7XjTcuz1fran16TrdaDe+X7vJ/wSL5NUwf97phstd63z++y1fr6ojzv06Sa348T71mibT//1yy1/vHNfOurO7/z0my1JjyQ7/fQBhN3zFZLGX+n9X8733v14pftlK3W1js+kK1WT2B9tWtYX60mjgQDAAAAAAAAAABA7TAIBgAAAAAAAAAAgNphEAwAAAAAAAAAAAC1wyAYAAAAAAAAAAAAaodBMAAAAAAAAAAAANQOg2AAAAAAAAAAAACoHQbBAAAAAAAAAAAAUDsMggEAAAAAAAAAAKB2GAQDAAAAAAAAAABA7TAIBgAAAAAAAAAAgNrp70TRxcun6tybd8tc86Ks9YBeEROs1RtMzl4TAAAAAAAAAIAq40gwAAAAAAAAAAAA1A6DYAAAAAAAAAAAAKgdBsEAAAAAAAAAAABQOwyCAQAAAAAAAAAAoHYYBAMAAAAAAAAAAEDtMAgGAAAAAAAAAACA2mEQDAAAAAAAAAAAALXTP1wD29tJOkXSFpJC0okRcUKnOwZgfeQRY2pNKFasyFLKWark15ex1qSMtXJXyyffvjNLNT1brXO1W7ZavWDnIy/NVmv+mbtnq7XDobdmq5Xrs0uStOWm2UpNW7AsW60Nb1yerdY9q2dmq/W9B5+brZb0q4y1UAc7fvTibLV+/tHNs9U67+6/ZKv1op32zVZrzZIl2WpVVSxfrtU33ZKl1hYPPpSlTm53v26XbLUmPRzZap3xvG9mq3Xzii2z1frCd16brdbX9PxstSY8Muzmtzb8PmOtatr5nfnWV7e/JN/vjvnPyve56qtuylZrxv1bZKuV0yazts1W6/6H8j3Hm96c7/Mr5zYF1APrq+0Zq/XVVr6FV0l6f0RcYXtDSZfb/l1EXN/hvgFYH3kEAAAAAAAAAKAFw+7SHRELIuKK4t+PSLpB0jad7hiA9ZFHAAAAAAAAAABa09a8RrZnSdpT0iVNbjvC9hzbc1Y/XP9pF4BuazWPK1eRRwAAAAAAAABA72l5EMz2BpJ+JumYiHi48faIODEiZkfE7L6N8s23C2B97eRxYj95BAAAAAAAAAD0npYGwWxPVNrgflpEnNXZLgEYCnkEAAAAAAAAAGB4ww6C2bak70q6ISK+1vkuARgMecRo2N7O9h9sX2/7Otvv6XafgF5FHoHqII9AdZBHoDrII1ANZBEYvVaOBNtX0pskHWj7yuLy4g73C0Bz5BGjsUrS+yNiN0l7SzrK9m5d7hPQq8gjUB3kEagO8ghUB3kEqoEsAqPUP1yDiPizJI9BXwAMgzxiNCJigaQFxb8fsX2DpG0kXd/VjgE9iDwC1UEegeogj0B1kEegGsgiMHotnRMMAFAvtmdJ2lPSJV3uCtDzyCNQHeQRqA7yCFQHeQSqgSwCIzPskWAAgHqxvYGkn0k6JiIebrjtCElHSNKUCRt0oXdAb2k5j5rWhd4BvYU8AtVBHoHqII9ANQyVxeJ28ggMgiPBAKCH2J6otNJ0WkSc1Xh7RJwYEbMjYvakCVPHvoNAD2knjxM1eew7CPQQ8ghUB3kEqoM8AtUwXBYl8ggMhUEwAOgRti3pu5JuiIivdbs/QC8jj0B1kEegOsgjUB3kEagGsgiMXmemQ1w+QZo3PX9NAG2LPmvFjEnZa2Jc2lfSmyRdY/vKYtnHIuLc7nUJ6FnkEagO8ghUB3kEqoM8AtVAFoFR4pxgANAjIuLPkhjBBCqAPALVQR6B6iCPQHWQR6AayCIwehxeBQAAAAAAAAAAgNphEAwAAAAAAAAAAAC1wyAYAAAAAAAAAAAAaodBMAAAAAAAAAAAANROf7c7AACoqAhpxco8pbJUSap6Nti+jLUmZayVu1o++fbDWarp2WqhPdv967XZat14ytOz1drp0Cuy1ZqwcnW2WppazVXvqQvXZKsVN+b8NATGh4O3flq2WrtfvjRbrWufka1UT1j9wMJud6Gpbc6en63Wqm1mZqv1scPfka3WA/80JVutjxx9RrZam/Y9mq3Wey5/bbZaaM/8Zy3JVuvgax/OVuu83TfKVmvNAw9mq+Upk7PVmvJgvi0BE1bm2xKw8x/fnK2W9ImMtYDO6fX1VY4EAwAAAAAAAAAAQO0wCAYAAAAAAAAAAIDaYRAMAAAAAAAAAAAAtcMgGAAAAAAAAAAAAGqHQTAAAAAAAAAAAADUDoNgAAAAAAAAAAAAqB0GwQAAAAAAAAAAAFA7DIIBAAAAAAAAAACgdhgEAwAAAAAAAAAAQO0wCAYAAAAAAAAAAIDaYRAMAAAAAAAAAAAAtcMgGAAAAAAAAAAAAGqHQTAAAAAAAAAAAADUDoNgAAAAAAAAAAAAqB0GwQAAAAAAAAAAAFA7/Z0o6pXS1AXOXhNA+9b0W8s36cteEwAAAAAAAACAKuNIMAAAAAAAAAAAANQOg2AAAAAAAAAAAAConY5MhwgAqIEIxapVWUrlnEAzMtaq6sSeOScwnZSxVu5q+bBPTx3sdOgV2Wr1/WHrbLVWH/CPbLX6N9s0Wy1tukm2UjOueyhbren3TM1WC+hFvz9572y1tt753my1dGO+Ur1gwrRp3e5CU30PLslXbOb0bKUmPZxvDf9TP39NtlorN8t3Xo4ttn4oWy10z3m7b5St1i5zJmardePsfNn2inzv+w3P+Fu2Wsv+Za9stdZcVM3PaGC8GI/rq2w1AgAAAAAAAAAAQO0wCAYAAAAAAAAAAIDaYRAMAAAAAAAAAAAAtcMgGAAAAAAAAAAAAGqn5UEw2322/277V53sEIDhkUcAAAAAAAAAAIbWzpFg75F0Q6c6AqAt5BEjxiAqUB3kEagGsghUB3kEqoM8AtVBHoGRa2kQzPa2kl4i6aTOdgfAcMgjMmAQFagO8ghUA1kEqoM8AtVBHoHqII/ACLV6JNjxkj4kac1gDWwfYXuO7Tmrly7J0TcAzR2vNvK4ahl5xFoMogLVQR6BaiCLQHWQR6A6yCNQHeQRGJ1hB8Fsv1TSfRFx+VDtIuLEiJgdEbP7pk7P1kEAa40kj/1TyCPWcbyGGUQFMGaOF3kEquB4kUWgKo4XeQSq4niRR6Aqjhd5BEaslSPB9pX0MtvzJP1Y0oG2T+1orwAMhjxixFoZRC0fRbgilo1h74De0m4eV2r5GPYO6B2t7mBEHoHOI49AdZBHoDrIIzB6ww6CRcRHI2LbiJgl6RBJF0TEGzveMwDrIY8YpWEHUctHEU7ylG70EegVbeVxoiZ3o49AL2hpByPyCIwJ8ghUB3kEqoM8AqPU6jnBAADjHIOoQHWQR6AayCJQHeQRqA7yCFQHeQRGr7+dxhFxoaQLO9ITAG0hjwAAAAAAAAAADK6tQTAAQD0wiApUB3kEqoEsAtVBHoHqII9AdZBHYGSYDhEAAAAAAAAAAAC1wyAYAAAAAAAAAAAAaodBMAAAAAAAAAAAANQOg2AAAAAAAAAAAAConf5OFJ2wSpp2/5rsNQG0b02/9Njmece713TkkwOVtHp1ljKRpUrijLWq2q+c+jLWmpSxVu5qQNnqA+7OVmvTv2ySrdbCfRdmq5Uz29py82yl+hcuzVYL6EVbfP2v2Wrd9NW9s9XS+/KVQntiab7P1ZzrqxOW5luX2+zi+7LVWrr5FtlqrZ6c74fv/Tdulq0W6uHG2Suz1Zp+Ub51uceetzhbrf4dtstWa4PLbs9Wq/+xbbPVAnrReFxf5UgwAAAAAAAAAAAA1A6DYAAAAAAAAAAAAKgdBsEAAAAAAAAAAABQOwyCAQAAAAAAAAAAoHYYBAMAAAAAAAAAAEDtMAgGAAAAAAAAAACA2mEQDAAAAAAAAAAAALXDIBgAAAAAAAAAAABqh0EwAAAAAAAAAAAA1A6DYAAAAAAAAAAAAKgdBsEAAAAAAAAAAABQOwyCAQAAAAAAAAAAoHYYBAMAAAAAAAAAAEDtMAgGAAAAAAAAAACA2mEQDAAAAAAAAAAAALXDIBgAAAAAAAAAAABqh0EwAAAAAAAAAAAA1A6DYAAAAAAAAAAAAKid/m53AABQURHSmjXd7sV6ImMtZ6xV1X7l1Jex1qSMtXJXA8oW7rsoW601v98uW60J/7IwW634x9xstSZsMD1bLQCjM/Ve9nntljWPPZavWM5aOd3R7Q4099geG2WrtfGfpmSr9dBTcv5aANa1ZL/7s9XaZc7EbLVu2ufebLVymnx3vs8JAKMzVuurrBUDAAAAAAAAAACgdhgEAwAAAAAAAAAAQO10ZDrECatCkxetzl4TQPtiorR0q7z5iXxHxwMAAAAAAAAA0BEcCQYAAAAAAAAAAIDaYRAMAAAAAAAAAAAAtcMgGAD0ENsb2z7T9j9s32B7n273CehV5BGoDvIIVAd5BKqDPALVQBaB0enIOcEAAJV1gqTfRMS/2p4kaVq3OwT0MPIIVAd5BKqDPALVQR6BaiCLwCgwCAYAPcL2DEn7STpMkiJihaQV3ewT0KvII1Ad5BGoDvIIVAd5BKqBLAKjx3SIANA7dpR0v6Tv2/677ZNsT+92p4AeRR6B6iCPQHWQR6A6yCNQDWQRGCUGwQCgd/RLerqkb0bEnpKWSPpIuYHtI2zPsT1nhZZ3o49Ar2grjyvJI9BJ5BGoDvIIVAd5BKph2CxK5BEYSkuDYJx8D6gO8ohRuFPSnRFxSfH3mUorUo+LiBMjYnZEzJ6kyWPeQaCHtJXHieQR6CTyCFQHeQSqgzwC1TBsFiXyCAyl1SPBBk6+t6ukPSTd0LkuARgGecSIRMQ9ku6wvUux6HmSru9il4CeRR6B6iCPQHWQR6A6yCNQDWQRGL3+4Rpw8j2gOsgjMniXpNNsT5J0q6TDu9wfoJeRR6A6yCNQHeQRqA7yCFQDWQRGYdhBMK178r09JF0u6T0RsaTcyPYRko6QpMlTNs7cTQCFtvPYv/EmY95JVFdEXClpdrf7AYA8AlVCHoHqII9AdZBHoBrIIjA6rUyH2NLJ99aZd3TS9MzdBFBoO49908kjAAAAAAAAAKD3tDII1tLJ9wCMCfIIAAAAAAAAAEALhh0E4+R7QHWQRwAAAAAAAAAAWtPKOcEkTr4HVAl5BAAAAAAAAABgGC0NgnHyPaA6yCPGSkiK1Wuy1HKWKvlFxlo5n2NV+5VTX8ZakzLWAjppwvPuyFbrls/tk63WEz95WbZasWJltloARmfjm1d3uws9a8K0adlqrXnssWy1qirn67XzV5Zlq+WlD2erdf4nfpit1mbHZCsFrOfG2fnW5XzBNtlqxYF3Zau1Zu7t2WoBGJ2xWl9t5ZxgAAAAAAAAAAAAwLjCIBgAAAAAAAAAAABqh0EwAAAAAAAAAAAA1A6DYAAAAAAAAAAAAKgdBsEAAAAAAAAAAABQO/2dKOrVoUmLV2SvCWAEJq+RZi3JXxMAAAAAAAAAgArjSDAAAAAAAAAAAADUDoNgAAAAAAAAAAAAqB0GwQAAAAAAAAAAAFA7DIIBAAAAAAAAAACgdhgEAwAAAAAAAAAAQO0wCAYAAAAAAAAAAIDaYRAMAAAAAAAAAAAAtcMgGAAAAAAAAAAAAGqHQTAAAAAAAAAAAADUDoNgAAAAAAAAAAAAqB0GwQAAAAAAAAAAAFA7DIIBAAAAAAAAAACgdvq73QEAQFWFFGvyVFqdpYwkyflKZRUZa+V8jlXtV0593e4A0AU7fuzibLVu+eI+2Wo98cP5+gVgdCYtXtXtLvSsNY891u0ujCtZX6+pE7OVuvMlm2Sr9fRfHpOtlvShjLWAzokD78pWa5c5+bJ94+wV2WoBGJ2xWl/lSDAAAAAAAAAAAADUDoNgAAAAAAAAAAAAqB0GwQAAAAAAAAAAAFA7DIIBAAAAAAAAAACgdhgEAwAAAAAAAAAAQO0wCAYAPcT2e21fZ/ta26fbntLtPgG9ijwC1UEegeogj0B1kEegGsgiMDoMggFAj7C9jaR3S5odEbtL6pN0SHd7BfQm8ghUB3kEqoM8AtVBHoFqIIvA6PV3oqjXhPoeXZ69JoD2zZi8VC/e6fqsNc+YvDRrPYypfklTba+UNE3S3V3uD9DLyCNQHeQRqA7yCFQHeQSqgSwCo8CRYADQIyLiLklfkTRf0gJJiyPit93tFdCbyCNQHeQRqA7yCFQHeQSqgSwCo8cgGAD0CNubSHq5pB0lbS1puu03NrQ5wvYc23NWRt4jegGs1XYeRR6BTiGPQHWQR6A6yCNQDa1ksWhHHoFBMAgGAL3jIEm3RcT9EbFS0lmSnl1uEBEnRsTsiJg90ZO70kmgR7SXR5FHoIPII1Ad5BGoDvIIVMOwWZTIIzAUBsEAoHfMl7S37Wm2Lel5km7ocp+AXkUegeogj0B1kEegOsgjUA1kERglBsEAoEdExCWSzpR0haRrlL4DTuxqp4AeRR6B6iCPQHWQR6A6yCNQDWQRGL3+bncAADB2IuJYScd2ux8AyCNQJeQRqA7yCFQHeQSqgSwCo8ORYAAAAAAAAAAAAKidlgbBbL/X9nW2r7V9uu0pne4YgObIIwAAAAAAAAAAwxt2EMz2NpLeLWl2ROwuqU/SIZ3uGID1kUcAAAAAAAAAAFrT6nSI/ZKm2u6XNE3S3Z3rEoBhkEcAAAAAAAAAAIbRP1yDiLjL9lckzZe0VNJvI+K3je1sHyHpCEmaMnGj3P0EoJHlccMtp41tJ1EfIcWayFLKE9ZkqSNJsTpbKTlfqazyvOpJzudY1X4BGJ0nfvjibLUefMs+2Wrpu2fmqwX0oDWTOAU4es+Eq+dmq7X91dlKadGZW2erNT9bJWD8uHH2ymy1bvt8xvXVj7C+CozGWK2vtjId4iaSXi5pR0lbS5pu+42N7SLixIiYHRGzJ/VPz99TACPK49RNJo91NwEAAAAAAAAA6LpWhtoOknRbRNwfESslnSXp2Z3tFoBBkEcAAAAAAAAAAFrQyiDYfEl7255m25KeJ+mGznYLwCDIIwAAAAAAAAAALRh2ECwiLpF0pqQrJF1T3OfEDvcLQBPkEQAAAAAAAACA1vS30igijpV0bIf7AqAF5BEAAAAAAAAAgOG1Mh0iAAAAAAAAAAAAMK4wCAYAAAAAAAAAAIDaYRAMAAAAAAAAAAAAtcMgGAAAAAAAAAAAAGqHQTAAAAAAAAAAAADUTn9Hqq5ZIz+2PHtNAO3bbuJjOn6rOVlr/nXiY1nrAQAAAAAAAACQG0eCAQAAAAAAAAAAoHYYBAMAAAAAAAAAAEDtMAgGAAAAAAAAAACA2mEQDAAAAAAAAAAAALXDIBgAAAAAAAAAAABqh0EwAAAAAAAAAAAA1A6DYAAAAAAAAAAAAKgdR0T2ojOmbBn7bH9o1poXzz9Fi5fd46xFgR4we48pcel522WtudfBd2jOVcvIY83Zvl/S7S003UzSAx3uzkjQr/bQr2SHiNh8DB+vJeSxY+hXe8ijyGMH0a/2kEeRxw6iX+0hjyKPHUS/2kMeRR47iH61pzJ57B/DTgAAxpFWV+Rsz4mI2Z3uT7voV3voV7WRx86gX+2par/GGnnsDPrVnqr2a6yRx86gX+2par/GGnnsDPrVnqr2a6yRx86gX+2pUr+YDhEAAAAAAAAAAAC1wyAYAAAAAAAAAAAAaodBMADAaJ3Y7Q4Mgn61h37VQ1VfL/rVHvpVD1V9vehXe+hXPVT19aJf7aFf9VDV14t+tYd+1UNVXy/61R76NQxHRPaiM6ZsGftsf2jWmhfPP0WLl93jrEWBHjB7jylx6XnbZa2518F3aM5Vy8gjAAAAAAAAAKCyOBIMAAAAAAAAAAAAtcMgGAAAAAAAAAAAAGqHQTAAwIjYfqHtG23Ptf2RbvdHkmxvZ/sPtq+3fZ3t93S7T2W2+2z/3favut2XAbY3tn2m7X/YvsH2Pt3ukyTZfm/xf3it7dNtT+l2n6qMPLaPPLaOPLaHPLaPPLaOPLaHPLaPPLaOPLaHPLaPPLaOPLaHPLaPPLauinlkEAwA0DbbfZK+IelFknaT9Drbu3W3V5KkVZLeHxG7Sdpb0lEV6deA90i6odudaHCCpN9ExK6S9lAF+md7G0nvljQ7InaX1CfpkO72qrrI44iRxxaQx/aQxxEjjy0gj+0hjyNGHltAHttDHkeMPLaAPLaHPI4YeWxBVfPIIBgAYCT2kjQ3Im6NiBWSfizp5V3ukyJiQURcUfz7EaUVgG2626vE9raSXiLppG73ZYDtGZL2k/RdSYqIFRHxUFc7tVa/pKm2+yVNk3R3l/tTZeSxTeSxbeSxdeSxTeSxbeSxdeSxTeSxbeSxdeSxTeSxbeSxdeSxTeSxbZXLI4NgAICR2EbSHaW/71RFVk4G2J4laU9Jl3S5KwOOl/QhSWu63I+yHSXdL+n7xWH9J9me3u1ORcRdkr4iab6kBZIWR8Rvu9urSiOP7Tte5LEl5LFt5LF9x4s8toQ8to08tu94kceWkMe2kcf2HS/y2BLy2Dby2L7jRR5bUtU8MggGAKgd2xtI+pmkYyLi4Qr056WS7ouIy7vdlwb9kp4u6ZsRsaekJZK6Ph+47U2U9kTbUdLWkqbbfmN3e4WRIo8tI4/oOPLYMvKIjiOPLSOP6Djy2DLyiI4jjy0jj21gEAwAMBJ3Sdqu9Pe2xbKusz1RaYXptIg4q9v9Kewr6WW25ylNNXCg7VO72yVJaY+vOyNiYO+qM5VWorrtIEm3RcT9EbFS0lmSnt3lPlUZeWwPeWwPeWwPeWwPeWwPeWwPeWwPeWwPeWwPeWwPeWwPeWwPeWwPeWxPJfPY34miDy+/99Hzbv7yjS003UzSAy2W3WUUXQJ61uVXL3+0b6u55BG5XSZpJ9s7Kq0sHSLp9d3tkmTbSvMh3xARX+t2fwZExEclfVSSbO8v6QMR0fU9YSLiHtt32N4lIm6U9DxJ13e7X0qHze9te5qkpUr9mtPdLlUaeWwDeWwbeWwPeWwDeWwbeWwPeWwDeWwbeWwPeWwDeWwbeWwPeWwDeWxbJfPYkUEwSTdGxOzhGtme00q7gbaj7xbQk8gjsouIVbaPlnSepD5J34uI67rcLSntofMmSdfYvrJY9rGIOLd7Xaq8d0k6zfYkSbdKOrzL/VFEXGL7TElXSFol6e+STuxur6qLPNYKeRznyGOtkMdxjjzWCnkc58hjrZDHcY481gp5bJEjIn/RFjemt7vRvdW2ANYijwAAAAAAAACAXsQ5wQAAAAAAAAAAAFA7nRoEa/UQt3YOhev6YXPAOEUeAQAAAAAAAAA9pyPTIQIAAAAAAAAAAADdxHSIAAAAAAAAAAAAqB0GwQAAAAAAAAAAAFA7WQbBbM+0/TvbNxfXmzTc/kLbN9oO23fbvtL2L0q3T7Z9hu25tv9h+5bi3x9p8liH2b6/qHGl7bfleA5AXZBHAAAAAAAAAADyHQn2EUm/j4idJP2++FuSZLtP0jckvUjSEkn3S3p9RLysdP+3SlokaRdJm0i6QdJukl5ne7cmj3dGRDytuJyU6TkAdUEeAQAAAAAAAAA9L9cg2Msl/aD49w8kvaJ0216S5kbErcXfPy7aN7v/XpKukrS3pJWDtAUwNPIIAAAAAAAAAOh5uQbBtoiIBcW/75G0Rem2bSTdUfx7iqS3S3qP7Vc0abONpPmSFkvaVNKdxbJGr7Z9te0zbW+X6TkAdUEeAQAAAAAAAAA9r7/VhrbPl7Rlk5s+Xv4jIsJ2DFJmB0kHSnq+pONtXxMRt7Tah8IvJZ0eEcttv0PpiJUD26wBjGvkEQAAAAAAAACAobU8CBYRBw12m+17bW8VEQtsbyXpvtLNd0narqhxl+1tlc4xtEbSnpJuKbW5S9L2kmZIWihp22JZuR8LS3+eJOlLrT4HoC7IIwAAAAAAAAAAQ8s1HeIvJL25+PebJZ1Tuu0ySTvZ3sP2BpIOkXSRpH0lXd9w/8sk7SHpEkkTi7a/KD9QsVF/wMuUNuADWIs8AgAAAAAAAAB6niMGmymtjSL2ppJ+onTUyO2SXhMRD9qeLelISWdJ+pbSuYkelPSApNsknRQRv7A9RdIPlY5EWS1pktKRKd+LiP+0/WlJc4q2n1fa2L6qqPXOiPjHqJ8EUBPkEQAAAAAAAACATINgAAAAAAAAAAAAQJXkmg4RAAAAAAAAAAAAqAwGwQAAAAAAAAAAAFA7DIIBAAAAAAAAAACgdhgEAwAAAAAAAAAAQO0wCAYAAAAAAAAAAIDaYRAMAAAAAAAAAAAAtcMgGAAAAAAAAAAAAGrn/wMe9BeD1ohGbwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 2304x288 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "visualize_couplings(couples_monotone, p_logits.numpy(), q_logits.numpy(), R_monotone, suffix='R_monotone')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABsEAAAETCAYAAABndqhPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABOFUlEQVR4nO3deZxcVZ3+8edJdzo7YRVCAgQcFhEHkIiAjoCiKPpTwRlHXBAHRXFn3JcR1Bk3XMI4joqILCIuiKgjiiKiIgiEHRKBAGENsm8he76/P85tuKlUd1d1nUrdrvq8X696dfrWqW+dqvRT99Y9957riBAAAAAAAAAAAADQTcZ1ugMAAAAAAAAAAABAbgyCAQAAAAAAAAAAoOswCAYAAAAAAAAAAICuwyAYAAAAAAAAAAAAug6DYAAAAAAAAAAAAOg6DIIBAAAAAAAAAACg6zAIBgDoebb3sx22D+90X9rF9gW2F3W6H8D6ZPtk29Gm2otsX9CO2gAArG+s14DqII/odVXZf1HsJzq50/3IzfbhxWvbr7Ssq/eLMQgGABgzSivlD3a6Lxia7d1sH2t7dqf7gvaxPdH2O22fb/s+2yttP2z7MttftL1Tp/sI9CLWlY2xvZHtTxWfWQ/bXmH7Tts/tX2IbZfaLire08Hb47Zvt32O7ffa3nCI54hhbh9dby8WY0ovrF+Lz6ljh8rOMI97se1vFe/FstodeEBu5HHIx0y0/TbbPy/WkUtt32L7DNvPaGN30UVsv79bB1xQPf2d7gAAAOg6u0k6RtIFkhZ1siNoD9vbSfo/Sc+Q9EdJX5O0WNJUpf//f5P0QdtbR8RdneonANRje09JP5f0NEm/kHS6pEclzZR0kKSfSnqXpP8tPexOSR8r/j1R0paS9pN0vKRP2D40Is6v83RXSfpKneVXtvo60H16aP26n9K24smSHm7icW+Q9HpJ10laoPSeAG1BHoc1W9IJki6U9F1Jd0vaTtJRkg6x/dKI+EPmfqL7vF9pf8HJde57iSTXWY72+ZOkSZJWdroj7cAgGAAAABpme5KkX0l6uqRDIuJnddpMlHS0pLZMRQigO9ieFhGPrefn3ELSL5UGsvaNiAtrmnzW9oGSNqpZ/khEfL9m2Wds76s0kPZz27tHxMKaNnfVeRywDtavDfmEpLdHxPLibNfdOtwfdCnyOKL7JO0eEVeVF9o+Xekgj+MkzelAv9AlImJFp/vQayJijaRlne5HuzAdIgBgzLI9u5gG5VjbryhNjbLY9nG21znYw/arbF9ZtLvD9mcljR+i/gTbH7d9fdH+Ydu/tL17Tbsn5062/R7bNxbtb7T9niFqb2/7tKKvK4ppJI6zPaWm3clF7em2v2n73qL2X2w/t07djWx/x/b9tpc4zaW9xzDv4RzbPyvaL7d9g+1P1L53RZ1Ftrcsprl4yPYTts+1vUOp3bGSvlf8+gc/Ne3TyaU2T++GqUN62Fsl7STpuHo7BCQpIpZFxOcj4u7BZR7m+lx1/kbK2X6t7auKaVYW2n5L0WZr22faftD2Y7a/b3vaEPU3s32q7QeKXPze9rOHaPuvti8saj5h+xLb/9zwu1O/5uCc6y9ymn7ttuL1XGJ7r6LNvsXzLik+F/6jTp2X2P6R03QzS4vPpN867YQvt9u9+Jz4nb3WlG59tv9UvK5dSstn2N7J9uRWXieqyU2sK4u/rxW2N6lTZ8eiztya5Q1lZjDnRQ4utP240mCUbG9s+2u2by769oDty21/qE6dVjP6IaUzwD5SZwBMkhQR50bEDxspFhF/lPQBpTMDmp7isPgs28l23W0R9JTRrl+3tP2VYl35UJGh+bY/YruvtobtrWz/2PYjth912rZ9er3nK9Yb/1Gst5bZvqbI4LFFpmfXtJ/htL16e/FZcrftE2w/rdTmZKWzTiTpVj+1rXjsSG9QRNwVEctHalc8z/giW1s30h6oQR6HEREP1A6AFcvnK52puUvtfUUe6742VJvtTW1/w2n/yYri5zdqtxf91HeeA4q/y9uc9jFcY/t1NW1D0jaS9vXa00XPLu5f55pgfmqfxGynfRgPFzk72fZU2+Oc9t/cWmTkCtvPq6kxzml/x59s31O8ntuLrKyz/dvk+2SnaUIvcZo6+3Hb19r+zCjfz7rZLu5b5JrrBvqpbe0DbP/VaTv5HtvH257aQP+f3K9Vb5nttzjtG1te/N9+eIg6RzntV1pu+ybb73ada5AVbdfb5wJnggEAusFBkt4p6VuSTpL0KkkflPSQpM8NNrJ9sNIUR4skfUbSKklvkfTy2oJOO6N+I2kfSadJ+h9J0yW9TdJfbL8gIubVPOw9kraQ9G1Jj0k6VNJ/2944Ij5dqr2HpPOVppv4tqS7JO0q6b2Snmd734ioPQX9XKUj7j4jaRNJ/y7pV7a3HTyKvujzuZKeU/T5r0pHyJ4n6YE6r/Hlks6StFBpqqYHJe1dPMdukv6l5iFTlE6R/6ukj0vaVtL7lI5+3yUiVhf1Zkg6Uum9X1A89uZSnd8rbfAyvcHYNLiz+cT18FyvkPQOpSnJHpR0hKSTbK9Q+vs6X+lv8TlKU9IsU9ppUes3xeOPVcrouyX90fbeEXHdYCPb/6l0lPlvJP2HpDWSDpb0E9vvjohvtPh6viCpT2n6tAGlHee/tX2Y0lQyJyhNy/ZapTNMbq05g+RwSRtLOlVparaZxev9ve39I+LPkhQRVxZfSo5X2in/+eLxx0j6J0nvKL/u4v43S9pfaRpTdKdG1pWnKP39Haq03is7rNRG0qgyM0fSayR9p1xH0k8kvaDo2zVKU7E8Q2mKpuNaeL56XiNpRc3zt2pwO+GgOveNt71pzbI1EfFg8e9TJe2rtE5dlLFPGHtGu379R0mHSPqZ0vbWeEkvVVrnbCfp7YMNna758ydJWynlbb7S398flHJX63+U1sN/kPRlSZsprZNvrW3oNNh0sdL67btFX/5BaXq0/W3PiYhHlLZ9N1DK7tGS7i9KXNPk6x7JTKXt0D8qfZYAzSCPo2B7nNJ3wb/XuXuBpNuUplLEGGF7uqSLlP5+TpJ0haTdlf6WXmh7zzpn9X9Rad/B4LTSb5F0hu2JEXFysexNSlOM3i/pv0qPvW+ELk1R+g74R6XvOYPfAycq7fN4rqSvK2Xvg5J+aXubUh8HlA6I+qnS1NhLihpHSHq+7T1aOAvtNKVpey8pXtPDSoPp/yzpU9Ko389mPLt4vu8obWPur7SfaRfbLy7O9hqNd0jaXOnz5GFJb5T0Rdt3RsQPBhvZ/ojS590VStOIT1Z6v4f6f11/nwsRwY0bN27cuI2Jm9IX2JD0weL32cXvSyTNLrWz0hFoi0vL+iTdrrSRtWlp+fRipRuSDi8tH5za4sCaPmxQ1LmgTr8ekzSrtHxA0qVKcyqXl18t6W+SptXUPrhOP04ulv1vTdt/KZa/vbTsyGLZp2vavr9Yvqi0bKKke5S+ePXXtB987fuVll1QLPtwTdsP1b5PSjvq13p8zWMWpU2Qzv9NcWv+pvTl4pE6y/skbVpzm1S6/+Sh/t+Lv5eTS7+Xs71NaflmSgNdayT9e02Ns5R2bE+tfc7iPpeW71HU+E1p2bOLtp+r07+zla4XNK20bFH5c2CE92wwE1dIGigtf2WxfKWkOaXlA0rXnLi4ps6UOrU3V/pcO6fOfT8vau+l9Dm1WtKZddoNvk/7dfrvi1vrN7W+rlws6dKamlZaV15TWtZsZqK4HVDTdrrqrOfq1Gzq+YaoMa2occ1w7eo8bpGk60Zoc01Ru95rrr3dU2pzQbFsdjN94tZ9N41+/TpJpXVcaflpxef+jNKyzxV/b2+paTu3WH5Badkzi2W/kTSutPxZRd21/m6Ldc69Km3zFsvnKB14dmxp2bGt/t0r7dwcbntzdu1r4sat0Rt5HPX79s6i1mfq3LfW91FuY+OmNJgTkt5Zs/xdxfLPlpYdXiy7TdL00vLpxbIHa/KyaKjPaBXXF6+zLCR9qGb5WUrf7eZJGl9aPvhdq7zPxOU+lJYfUbR9bc3ytb6nDvM+vbZoe1o5o8V95cw2834Omc16752e2s58dc3y44vlr6vzf7Vfadl+Wnd/1OCyu2v+TycrDWxdXFq2saSlStvEE0vLt5D0SO3zlfq8qPb1tePGdIgAgG5wdkQsGvwl0tr0D5K2KJ32vYfSUXbfi4j7S20fUTryrtYblQaqLi9OV9+0OJJ7QNLvlI4Sqj1C7/SIuLNUe4XS0U39kv6fJNl+ltIRgj+QNKGm9oVKOylfUqc/X6v5/fzi5/alZa9W+hL0lZq231TaOVj2YqWd59+TtGFNP84p2tT2Y42k/26gH8OKiNkRwVlgY9cGWvfvSUpnbdxXc3tXi891dkTcNvhLRNwn6Qalv8XaMz7+rHTE3+w6db5UfC4M1rlcKccHlD4j3qC0EX5KOQ9FJn6htPN87xZfzzdj7SML/1z8vCRKZ5YWbS5VTa4iYsngv52m/NhEKfOXKB31WOstSoPdZ0j6vqQ7VOdMuYg4PCIcEReM5kVhzBhxXRnpjN7TJT3Ha09bu5+krbX22VOjyczVEXFezbKlkpZLem696V5afL5aGxQ/632GtWqw5gY1yy9RWueWb0+eaR0R+xX5W9SGPmFsGdX6NSKWDq7jbA84TS+6qdLsAOO09nV5Xq10hsapNc/xxTrP+4ri5/FROnI8Iq4taj+pOLL9FUpZXFaTz0VKsw7U275tm4hYVGRrv/X5vOga5LFJtveR9FWlAz4/V3t/kcfZuZ8XbXew0t/5CTXLv10sP7jOY75Z7GeRtNY+l43U+pm5q5XO9Cr7s9Lg1rdi7Rl1Br9rPfmdKpKl0pNTjG5YZGNwv0a971SNeEPx84NRc7ZVze+jeT+bcUNEnF2z7Aul5x6t79X8nz6hNENQ+fvqi5UOtv5mRCwrtb1H6fvFOtbn5wLTIQIAusEtdZYNTv+3iaTHlaafkNLAVq35dZY9Q+lIvuFOx99UaafyoAV12gzWHnz+ZxQ/P13c6tm8zrK1XmNEPOB0qZ/yvNHbKR3R/2hN2+W2b1Ha6Bw02I+ThuhDvX7cXd6YKZTfZ/SGR7XuTl4pTcXy4uLfuypN09Kqetl+SOnvvPaaIA8VP+v9LQ6VzZcoTc15vVImrPqfEYPqZVPSkzscagfG7ysGFQbV5vihIse31in5kGpeSzFf+n9JOlDShjXtQ7ULIh50uoba74pFz4+Ih4d6Deh6jawrpTTQ9QGl6Q8/Xiw7TGmnQ/kL7Ggyc2Ntg4hYYfv9Skep3mp7vtKOiLMj4vejeT7bGysdtFJ+nnv01A7NutcPbNFQA2z31xn4A+oZ1frV6bp+H1XK6T9o3emmy9t/20m6rGbdpIhYbPvhmsdtW/y8oU6fbpD0stLvOyrt4D+iuNVT7zNoLbYHlI4kL3s8Ih6v1x5oI/LYRB6L6f5/pXS2yMvrfGfE2LWtpHkRsaq8MCJW2b5R6Uz9Wo3sFxmtxXX+vga/B671nar0Xav2O9VrlbZ1d9e614ffSKOzfdG3elOBlo3m/WzGOu996TOllfd+qO8R5fd2pM+pjmIQDADQDVYPc99ozziypGuVrr01lJHmqx6qrpTO1vrNEG0eql1Q++WoTr3R9uNDkq4aos3dNb+3433G2HOdpBcU16N78otGcZbSeZJke1Wdx60zSFO0HW57dKi/uXZlPpR2IgxV//phHn+80nW1ymqv8TOa15M6l87U+ZPSPPhzlT6fHlM6K+5jkl44xENfWfr3bpL+MtJzoWs1lJuIuNb2VZLeYPsTSoO7r5H022IgqfyYZjPzRL1GEfEt2z9XukbnvkrXMni37R9FxOCF1Jt5vrOKOmu9xoh4zPZtknayPWnwSOBW2Z4gaQelnR+tXMcBvW2069evKl2X9kdKB0rcqzQV7rOVzihZHzMADX6GfF9DX2+vkbzto3SGatmnlaaEAtYn8thgHm0/W+mAq0ck7R8RdzXcU6B5w23PjrjPxPYhSvm8VOn65ncoTbffp2K60TzdzKLu9+fC+h7TGfH7atUxCAYA6BWDR67sVOe+nessu0np+kPn157OPoxn1Fk2WHvw+W8qfq5uw5Hht0h6ie0NymeDFTvnttPag2uD/VjShn4Mt7GGse9MSS9QmlbvE0087kEpnaEREQ+Wlrd6NGAjnqE0XUPZzkob84PTLd6kdOHy2yOi3tGLI/mS0s6GsnvqNRylF0naUtK/RcT3ynfY/s96D7D9SqUdMd9TOhr5y7b/VEydAwznFKVpePdXusD9NK27I63VzKwlIhZLOlHSibb7lK6pcKjtr0TEZU0+3wc09JG8Zyld+/JNWncqmtF6k6QJSkfBA6M12vXrmyT9qTRgLEmy/Q912t4iaXvbfeUDrGzP0LpnGC8qfu6odY8A37Hm94VK238DDW5XDrWteLWeOsum3GdgfSOPDeSxGAA7T+nArP3L05ija9wiaUfb/eWzl4oDGXdQ/c/oZyhdl66sdr+I1Jn9Bm9SGvTav5jST5JUMw34aNwo6VW2Nx/hbLBm3s/B78wbq3Rgpe2JStvnC+vUX2efVOkzpd3r00XFzx311PSSKi3rqCqNbgIA0E6XS7pT0luKOZ8lSbY3kPSOOu1PVbqAZ90zwWzXmxbtDbZnldoMKO1oWy3p/4rFVyodWfgO2+vs/LfdX0zjNBo/VzqC6QM1y4/SutN5nKt0ZOJH6z2f7Um2Rztd1OAUGXVfh+2nZ9jIROecqDQd2YdsDzWveL2zsQanQTugZnnt32s7fNjFXBjSk1/YD5D0+9KULqcVPz9X7IBfyxCZf1JEzI+I82puOaeCGdw5stZ7a/slqjN3ve2ZStOdLpD0bqV56pdK+qFrrmdoe4btnWxPzthfjG0/kLRKaTqnw5SO7q7dmdFSZkrtJtf+7RU7A68pfh1clzT8fBFxeW0eS02/pHQm95ds172GmO2X2H5dvfvqtN1X6ezuxyR9vpHH1Dx+6yJ/tdPxoPeMdv26una57SlK26C1fq40behhNcs/UqftL4uf77P95L4jp+vbHlhuGBEPKF1T9hDbe63T6WSz0qK624oR8VCddemodtrZHl9ka+vRPB49jzyOkEfbuyudAfa40oBCvem9y8+7UzG1N8aWs5UODq69rvDbiuU/q/OYo4qp4iU9OW38OyQ9LOmPpXaPa4h9Bm20WmnwrZwjS/pki3UHpwz/UjmjpfqDzlbj7+dQ35+P1tBjOjvafnXNssHPlLOHeEwuv1O6zu9RxUCdJMn2FnrqmmlrWZ+fC5wJBgDoCRGx2vbRkn4s6VLb31HawfdvSnMZ135BPl7pyLfjbL9Q6UiWR4t2L1Jx9FDNY26UdIntbyntDHu9pOdI+mxE3FH0I2y/qah3je2TlKZvmqx0psYhSlObnTyKl/k9SUdK+pTtbSVdrDTP9b9Iulml9X5ELLF9mNKG0A1FPxYqHSG0U9GPgyVdMIp+XKY0RdsnbG8kaYmkWyPikuL+3ytdh4kpFMegiFhq++VKA7tn2b5A0m+VznraQOnv51+VvmCUr5l3htJFsk8oBkEfVDqrY1O13zaSzrX9C6Wj5t6tNCD0ocEGEXGZ7WOVpni5yvZPlKYEnSFpD0kHqeYaQ+vZhUrv8Vdsz1Ya1N9N6WjGayU9a7Bh8cXrdKWpE/+1OMrxCdtHKJ0FM1fS20u1P680leP+Gl3m0WUi4l7bv1aalnCipO/WDupmzMwOkv5o+2dKB4k8pHQU61FK13b4c87ni4h7bL9CaefjhbbPVppq9FGlsy1fKun5xfOXTbf9xuLfE4q2+ytd4P1eSa8b5c76U5WmbqydPhU9poX165mS3m77R0pnZGyup7Zva31Jafv0O07X8Lle6W94b0n31/TnetsnKG1bnldkdDNJ71I6qGsPrX0U/1FK66o/2T61aDNO6YzvVyn9rR9btB08O/uLtk9X2q6+LiKuG+49sv2Pemqa3+cVP99k+/nFv78eEY8U/56pdCDIH4vXCDSMPA6fR9vbKO3w3kjSf0vax/Y+Nc1+VkwfOWiB0gwMs4eqi0r6ktL+hG8UBxJeqbSP4Qil6zx9qc5j7lfaLzI4e8VblPajvLV89pXS394Rtj+r9PexRtIva/5ucjtTaZrv84tsjJf0aqX9MaMWET8pcn+Y0hmev1Dapt1BaaB6l6JpM+/necWyz9jeRGm7+PmS9lLNZ0TJtZK+X+zvuklpW/WfldaFP2rlNY6kuG79p5W+8//F9veV3tcjlfaVzdG6Z/+tv8+FiODGjRs3btzGxE3pS0FI+mDx++zi92PrtD22uG92zfJDlK6BtVzpC8tnlQa7QtLhNW37Jb1XaVBnSXG7SWnn8kvq9Ovwov1NRf2bJL1viNeyjaRvKe3wWqH0xehypZ3RW5XanZxW13VrhKSTa5ZtLOm7Rb0lSju05xQ/F9WpsYvSFG53Ff34u6SLJP2HpI1L7YZ6fN3/A6Ud6vOLmmv1s3jNdV8Tt7FzU7pO0LuUrhVwv9L1Dh6WNE9p433HOo95rtI1qZYVjzlBaeC19m9kuGwP9bd4ePGY/UrLTi6WbaZ0FskDStckOl/SHkO8rpcrnSn5oJ76nPi1pHfUtFsk6YIG36t1+la6b50cl/tes+wfleaqf0hpoP0CSf9U21bSp4q6R9Wp+43ivtfUeZ/W6R+3sXdThnVlcd9rivtC0vOGeb5GMzPU3/omSlMvXlV8hixVOihjrqQZo32+Bt6njSUdo/SZ9YjS+upOpZ0jr6xpu6j0XkTxOTL4vO+VtOEQzxGS/m+Eflww1P8Bt968qcn1q9IOpuOUdiItU9r+/KjSQVv1tm+3Lv7OHy1uv5T0dNVZrynNMHCMpNuLvF0j6bWSvlzUflpN+02LvtxY9OVhpR1yx0vauabth5WmZlo51GdUnffm8Jos1t5ml9rOLpZdMFJdbtyGupHHId+X/UbIYr3v4aE62/Dcqn9T+i71v0rbSSuLn9+QtGlNu8HP6AOUrh83+Ld6raTX16n7NEk/VdqmW1P+u1Gd73z1ltU873517qu3z+RtSvsqlklarPSddOMh2tbdfh3ifRpXfF5cobSt+FiR02NG834WbXdQ+v73RJHhHysd5FHvMyKUvtcdIOkSpW3qv0v6uqRpI71npVwfPtyy0n0nq85+neI9uFFP7RN7t9IU/SFpzzp9Xuf/tB03F08IAABGyfZ+Sl+M3hIRJ3e0MwAAAEAb2f6lpBdK2iBK1zICsP6RR1SF7cOVZqfZPyIu6Gxveo/tkHRKRBze6b7Usv11pcGwGRGR87rZDeOaYAAAAAAAAFhL7TUki2X/KOllks5nhzuw/pBHAFVXvhZYadkMpWkir+vUAJjENcEAAAAAAACwrjcX15D9laT7lK6DdKTS9KGf6mTHgB5EHgFU3X62j1O6DvWdSlMTv03SVKWpYTuGQTAAAAAAAADUukLSwUrXvttY6fom50v6dERc2cmOAT2IPAKouoWSblYa+NpE6bpr8yR9PiLO62THuCYYAAAAAAAAAAAAug7XBAMAAAAAAAAAAEDXYRAMAAAAAAAAAAAAXYdBMAAAAAAAAAAAAHQdBsEAAMjI9jjbR9v+m+1ltu+w/RXbUxp8/FTbH7d9re3HbN9v+yLbh9t2nfYvs/172/fYXmL7Bttftr15/lcHjB0ZshhD3B6vaWfbb7T9Q9sLbT9h+3bbv7D93Pa8OmBsaTWPRY2Ni/XbwqLGfbb/YPufatp9oVhv3mt7efFc/2d7v9yvCxiLWsmj7R1tn257ge1HinXe32x/1faMIR6zd7FOvNP2Uts32/6O7e3yvzpgbGkxjzvY/oztvxbrxMdsX2X7E/UeP8y27eDtE+15lcDYkOH748ds/8T2LUWmFg3T1rbfYfvKYt34sO3f2N4r2wvCWhwRne4DAABdw/bxkt4r6WeSfi3pGZLeI+nPkg6IiDXDPHacpD9K2kfSKZL+KmmypEMl7SnpSxHxkVL7t0k6QdLlkk6XtETScyS9RdLtkp4VEUsyv0RgTGgli8Xjo2h7Qs1dKyPiR6V2EyUtlXSVpF9JulXSDEnvkLSlpMMi4vsZXhIwZmXI4zaSLpA0VdJ3Jd0oabqkf5R0bkT8sNT2Akl/k3STpIckbSHpjcVzHhYRp2V8acCY0+K26oskfUJpG/VOSaskPUtp2/NRSbtFxL2l9i9VWjfeLOkkSfdLeqakIyUtV9pWvSvzSwTGjBbz+AVJ75L0C6VMrpS0v6TXSrpG0l4RsbTU/o1DlDpW0tMl7RoR17T4koAxK9P3xwclXSFpD0mPRsTsIdp+U+n74gWSzlLa73OkpFmSDoyIC1p+QVgLg2AAAGRi+5mSrpX0s4h4TWn5eyT9t6Q3RMQPhnn83pIukjQ3Io4uLR9Q2qG3cURsWFp+g6RpkraLiGWl5f+ptIPi4Ig4O8+rA8aOVrNYtA1Jp0TE4SO065f0vIj4Y83yzSVdL2m1pBkjfWkCulWmPP5Z0mxJe0bE4lH0YarSTvgHImLnZh8PdIsceRyi7r9I+rGkj0TEl0rLz1XaKb9lRNxfWv5WSd+RdHREzB3lywHGtAzfHedIuikiHqlZPvhd8D0R8T8j9GGWpNskXRERzxn1iwHGuEzbq9tFxC3Fv6+TNLXeIJjt3SRdKek3kg6KYnDG9oZK+30elbQT3x/zYjpEAADyOVSSJc2tWf4dSU8oHYk+nA2Kn3eXF0bECqUjZ2vP6tpA0kPlAbCax3MWGHpVq1l8ku2BYgd6XRGxqnYArFj+d6UzO59W3IBe1VIebb9A0vOVzoZebHu87cnNdCAiHpf0gKSNmnkc0IWyrR9r3Fb8rM3YBpKWKZ2VWca2KtBiHiNiXu0AWGFwxoJdGujDW5T2DZ/YQFugm7W8fhwcAGvA/sXPUwYHwIrHPyzp55K2l/S8BmuhQf3tKLrpxn0xe6vxWWsuumOl7n9w9TrXQgEwvIG+STFp/PSsNZeufEQrVi8lj8C6niNpjaRLywsjYpntq4r7h3OppIclfbiYP/oSpdPi36x0Ov07atqfK+nNtr+i9MXl8eI5/kNp5/v5o38pwJjWahYH/bPSF54+2/cp7VT45BA7HOqZJWmFUq6BXtVqHg8qft5u+5eSXqaUyZskfWao6UZtb6q0Y2+GpLcpTWlz0mhfBNAlsqwfi6mAp0qaKGlnSV8s7jqnpum5kvaSdIrt45QO6tpF0lckLZD0QwG9K9f2aq1Zxc+/D9fItpUGwZZIOmOUzwV0i3blsZ4Jxc8n6tw3uGwvpWkYkUlbBsFmbzVel567Vdaaex54R9Z6QK+YNH669t76sKw1L7791Kz1gC6ypaT7I2J5nfvukrSP7YHizK51RMRDtl+pNKD149Jdj0l6TZ2pDd+nNEj2Pkn/Xlr+PUlvj4jVo3sZwJjXUhYLl0r6iaSFSkeyHyTp3ZL2tb1PcWbJkGwfpHQtv9PqnK0J9JJW87hj8fM7Stf5erOkAUkfkHSa7fER8b3yA4qzN+8rLVqqdH2/8roS6EU51o+S9FZJXy/9vkjSGyOidofd55XOhv43SW8oLT9H0qER8VgznQe6TK48Psl2n9IBkaskjTS16QslbSvp5Ih4tNHnALpU9jwO4/ri5wuVrukn6cmB6X2LX/MOrKA9g2AAAPSoyUoX+a5nWanNcBtOj0u6Tmlj6CJJGytd8PgHtl8VEb8rtV0p6XalC7f+UumooQOVdjSsVjryHehFLWcxIp5bs+hU29dI+i+lgef/GuqxtreXdJrSF6YPNNhnoFu1msdpxc/HJO0/uPPB9tmSbpH0Odun1Fw3YamkFyt9391Gaef71OJ5mH4NvSzHtqokna103ZKpknaX9EpJm9Zpt1ppXXie0vbqg0pTPL1H0g+LbduVTfQf6Ca58lg2V9Lekj4eETeM0Patxc/vNlEf6FbtyONQfi1pvqR32r5b0llF7X/XU9OYNjX1N0bGIBgAAPk8oaGv/TOx1KYu289SGvg6OiK+VVp+htLA2HdsPz0iVtsep3Qh1X5JzyvNJX2m7QckfcT2jyLivNZeEjAmtZTFYRwn6RhJL9cQg2C2t5X0e0kh6WURcV+9dkAPaTWPS4ufZ5SPvi3Onv6FpMOUzhZbULpvtdJOd0mS7RMlXSDpfNvPZqc7eliW9WNE3CnpzuLXs23/VNJltidHxOdLTU+WtI+kZ0bEYJZ/ZnuhpG8qndnJtYjQq7Jur9r+rNKsBSfU5LBe240lHSzpbxFxYaPPAXSxdn1/XEdErLL9MkmnKE0nPDil8DWSPqo0ZTBnZ2Y2rtMdAACgi9wtaVPbE+rcN1Pp9Prhjhw6WmkD6yflhRHxhKRfKR3NPrtY/HxJ/yTpp+WLqRYGH7+vgN7UahbrKnac3636R7vL9mxJf1A6Mv7FEXFts88BdKFW8zi4o/2eOvctLn5uNFwHikGx05WOrn3B8N0Fulq71o/XSLpS0jsHl9neWukszF+VBsAGsa0KZMyj7WMlfVJpWvza60jX8wal6xJxFhiQtGX9OJSIuD0i9lfax7OvpF0iYlc9ddbZ33I9FxIGwQAAyOcypXXrnuWFxcXDd5M0b4THzyx+9tW5r7/mZzNtgV7TahbrKh4/S3UuNF4MgF0gabrSANiVo3kOoAu1msfBC5TPqnPf4LJ7G+jHpOLnxg20BbpVW9aPhUlaO19sqwLDy5LHYgDsGKWzSt5a5wDJeo5QmlqfC74DSTvXj0MqBsP+FBGD1wk7SNIaSee24/l6GYNgAADk8yOlKdDeX7P8bUpzOp8+uMD2023vVNNufvHz8PJC2xtKepWkhyQtrGn7Btvja+oMPv6yZjoPdJGWsmh7kyHqflZph90va9pvo3QG2IaSXhIRl7fQd6DbtLpuPFvpemBvtD211HaGpFdLujEiFhbLNrI9UNsB21OUdvit0VODakAvanX9uEW9orb3VzrT8q+lxTcoXRPs1cW2bNnhxU+2VdHLWl0/yvanlAbATpP0bzXXx6zL9hxJu0r6ZUQ0chAJ0AtazmOrbL9Sadr90yLittz1ex1H3QAAkElEXGv7G5LebfssSedIeoak90r6o6QflJr/XunUd5eWzVW6tskXiuuD/UXpiNq3SZoh6V3FlE6KiKuL6y+8RtI8299XmqP6QEn/T2knxM/b9FKBSsuQxU/a3ktpYOt2pekND5K0v6RLJH19sKHtaUW72cXyHW3vWNOl30XEOmePAb2g1TwW1/76oKRvS/qr7ZMkDUg6qvj5ntLj95X07WL9uFBp8GxbSW9SOmvs0+xUQC/LsH78ZjEAfb6k25Sm8d5D0uuU8vaB0nM9aHtusexK29+R9KCk5ylNxXazuB4YelirebT9LkmfVtpWPU/S6+1yXPX3iPhdnac+ovhJ/oBChvWjbL+pWC5Jm0kasP3J4vfbIuK0UtvvFo+/Sun6t89XWjdeJul9WV8cJElu5CxZ2y+VdLzSaewnRsQXhms/Z9eJcem5W+XpYWHPA+/QvKuXeeSWQHdrNo/TJ24Re299WNY+XHz7qXpk2T3kEajDdp/S0UNHKu0Uv1/pqKJPRcTjpXaLJG0TEbUbTk+X9ClJL5K0udIG0VWS5kbEWTVtB4rneoOkHZTO8L5N0pmS/isilmR+ecCY0UoWbb9K6bomu0jaROlI9psk/VjSVyNiWantbEm3jtCd/SPighZfEjBmtbpuLO47RNKHJT1L6Yyui5UGtf5SavN0SR+TtI/SVGyTJT2gtEPhWxHxq/yvDhhbWlw/vlbpgK1dlXbwhdK25+8kHRcRt9c8lyW9tbg9U+kaRHcpXev22Ii4rx2vERgrWszjyZLePEz5P0bEfjXPN0npepqPSprdyJljQK/IsC/nAg19rcu18mj77cXz/IPSQV0LJZ0h6Wt1rqOJDEYcBCv+AG6U9GKlixJfJunQiJg/1GMYBAPaYzR5ZBAMAAAAAAAAANCLGrkm2J6SFkbELRGxQtIPla5LAmD9I48AAAAAAAAAADSgkUGwmZLuKP1+Z7FsLbaPtD3P9rz7Hlidq38A1tZ0Hles5ixaAAAAAAAAAEDvaWQQrCERcUJEzImIOZtt0perLIBRKOdxoG9Sp7sDAAAAAAAAAMB618gg2F2Syhf4mlUsA7D+kUcAAAAAAAAAABrQ30CbyyRtb3tbpZ3tr5P0+rb2CsBQyCPWm4G+STFp/PQ8xcZlO/FYMc75avXlq7WmP2etbKUU4/PV0oQ12UpNn5Bvqtatxj+Rrdbl1yy/PyI2y1Ywk/6JU2Jg2sZ5at2/JEud3Dw+3x9+ZKy1ZiDfDAerJmcrpeiPbLU2nprvb2Jmf75sVzWPA+OnxMQJG2apFRnXHcr3EZ21X6sHspXSwCOrstWKZcuz1eoFj+mhauZxYEpMnLhRllp+It/fRKzJd4kIT5yYrZZWZ7x0ReRbD2lcvnXtsi3zfe941rQHstW67oF88Vl+152VzGP/5Ckxfnqe7dXxf8+4vZrxT9U5v9euybjinppvBiAvz7euXT0743fkh/N9sd12i79nqzX/2pWVzOOAJ8RETel0N1AhAzvl+/yaPfB4tlo3XpPvS/Jw26sj7iGIiFW23y3pXEl9kk6KiOuz9Q5Aw8gj1qdJ46dr760Py1IrJk/IUkeSVk/NV2vF9Hx75pZvlO/L+xOb5ds4WToj47e+2fm+jB60/fxstebOmJetVt+MhbdlK5bRwLSNtdOrj85Sa5PvXpylTm79m22RrdbqLTbJVmvJ7KnZat27R75sr9g4347MQ/f+a7Zan9v8mmy1qprHiRM21J67HZWl1oqN8q2H+pbl25m2fKN8A8mPbpNv/Tjr1/dnq7V6/o3ZavWC8+LMauZx4kaaM+ddWWoNXH1rljqStOaxx7LV8g7/kK3WuIfy7bTSqoyD0tPy7ai98dhp2Wpdut/J2WrtcHKe9YYk3fyxD1Qyj+Onb6zZR/x7llpbfenSLHUkKTL+rY6bmu/vK+fnxJrdd8tWa+DWe7PVevh/8g3iP/GLfN8VfvDhL2ertcs2d1cyjxM1Rc/1izrdDVTIlqfk+/z63tZ/zlbrwC13y1ZruO3Vhr7dRMQ5ks7J1iMAo0YeAQAAAAAAAAAYWb7DUQEAAAAAAAAAAICKYBAMAAAAAAAAAAAAXYdBMAAAAAAAAAAAAHQdBsEAoIfYfqntG2wvtP3RTvcH6GXkEagO8ghUB3kEqoM8AtVAFoHWMAgGAD3Cdp+kb0h6maSdJR1qe+fO9groTeQRqA7yCFQHeQSqgzwC1UAWgdb1t6PoHSsn6/2L52Su+UDWekDPGDdOMXlC9poYk/aUtDAibpEk2z+U9CpJ8zvaK6A3kUegOsgjUB3kEagO8ghUA1kEWsSebADoHTMl3VH6/c5iGYD1jzwC1UEegeogj0B1kEegGsgi0CIGwQAAT7J9pO15tuetWL20090Belo5j6uWLel0d4CeVs7jylXkEeiktbZXV5BHoJPW2l59gjwCnbTW9qqWd7o7QKUwCAYAveMuSVuVfp9VLHtSRJwQEXMiYs5A36T12jmgxzSVx/6JU9Zr54Ae01Qex/eTR6CNmtteHSCPQBs1t706mTwCbTJiFqWa7VVlviwKMMYxCAYAveMySdvb3tb2gKTXSfpFh/sE9CryCFQHeQSqgzwC1UEegWogi0CL+jvdAQDA+hERq2y/W9K5kvoknRQR13e4W0BPIo9AdZBHoDrII1Ad5BGoBrIItI5BMADoIRFxjqRzOt0PAOQRqBLyCFQHeQSqgzwC1UAWgdYwHSIAAAAAAAAAAAC6DoNgAAAAAAAAAAAA6DoMggEAAAAAAAAAAKDrMAgGAAAAAAAAAACArtPf6Q4AACoqQl6+stO9WEdfxloDGWvlrpZPvuNdlmpKtlrnaOdstfJa2OkO1NV//xJt8t2Ls9S68btzstSRpB2OmJet1qrF92SrlXMDd+rK1RmrTc9W6cGd8r3Ki7ffNlutVz86I1utquZRS5bKF12dpdSELFXyy5mhfGsO6e9H7p2t1rS5s7PVmnDgbdlqKSJfrV7w+BPqv/CaLKU8Y4ssdSTJS5dmqxV9+bZ+bzpqq2y1nv7Dh7LV8vGPZqu1wY+elq3WQf/52my1puzvbLWqavw9SzTrcxdlqXXHx/fJUkdStj5J0rgNpmWr5c03zVZLF16VrdT8b+6ZrdbGP8v3+XXlMf+brdbLtnthtlrS9zPWqqblL3tOtloTfn1Ztlo5LfrPfNuYsz+Z53u7JN111jOz1dJe12crdaB2y1ZrfeFMMAAAAAAAAAAAAHQdBsEAAAAAAAAAAADQdRgEAwAAAAAAAAAAQNdhEAwAAAAAAAAAAABdh0EwAAAAAAAAAAAAdB0GwQAAAAAAAAAAANB1GAQDAAAAAAAAAABA12EQDAAAAAAAAAAAAF2HQTAAAAAAAAAAAAB0HQbBAAAAAAAAAAAA0HX621H0keWTdM5NO2eu+aes9YBeEeOs1VMnZK8JAAAAAAAAAECVcSYYAAAAAAAAAAAAug6DYAAAAAAAAAAAAOg6DIIBAAAAAAAAAACg6zAIBgAAAAAAAAAAgK7DIBgAAAAAAAAAAAC6DoNgAAAAAAAAAAAA6DoMggEAAAAAAAAAAKDr9I/UwPZWkk6VtLmkkHRCRBzf7o4BWBd5xHq1JhQrVmQp5SxV8uvLWGsgY63c1fLJd+zMUk3JVusc7ZytVi/Y4Yh52WrdddYzs9Wa9dobstVa88ij2Wppkw2zlZp07/JstaZMn5St1tJTZmSrtWj/VdlqAbU2PeHifMVOyFfq3LuvylbrwFl7ZKulNavz1aqqkGJVns+dVXfcmaWOJN39wX2y1Zo599Jstba9Mt9n9ENv2Ctbren7/y1brS12npqtVtyxOFutUz50arZau/9PtlKVNetzF2WrdffP8n1X2PLg+dlq5XTHf+T7zNnhqHzv/bhp07LVOvDbu2Wrde7df81Wqy/fZnRlTfj1Zdlq3fXRjOvHL+T7W539yXzbmLf9+FnZam1zyLXZalVV38475Ct2/dB3jTgIJmmVpA9ExBW2p0m63PbvIqKan/xAdyOPAAAAAAAAAAA0YMRDuiNicURcUfz7MUkLJM1sd8cArIs8AgAAAAAAAADQmKbmNbI9W9Luki6pc9+RtufZnrf60SWZugdgKI3mceUq8ggAAAAAAAAA6D0ND4LZnirpp5LeHxHrXGghIk6IiDkRMadvg3zX+QCwrmbyOL6fPAIAAAAAAAAAek9Dg2C2xyvtcD89Is5qb5cADIc8AgAAAAAAAAAwshEHwWxb0nclLYiIr7a/SwCGQh7RCttb2f6D7fm2r7f9vk73CehV5BGoDvIIVAd5BKqDPALVQBaB1jVyJtjzJL1J0gttX1XcDmpzvwDURx7RilWSPhARO0vaS9K7bO/c4T4BvYo8AtVBHoHqII9AdZBHoBrIItCi/pEaRMSFkrwe+gJgBOQRrYiIxZIWF/9+zPYCSTMlze9ox4AeRB6B6iCPQHWQR6A6yCNQDWQRaF1D1wQDAHQX27Ml7S7pkg53Beh55BGoDvIIVAd5BKqDPALVQBaB0RnxTDAAQHexPVXSTyW9PyIerbnvSElHStLEcVM70DugtzScR03uQO+A3kIegeogj0B1kEegGobLYnE/eQSGwJlgANBDbI9X2mg6PSLOqr0/Ik6IiDkRMWdg3KT130GghzSTx/GasP47CPQQ8ghUB3kEqoM8AtUwUhYl8ggMh0EwAOgRti3pu5IWRMRXO90foJeRR6A6yCNQHeQRqA7yCFQDWQRa157pEJePkxZNyV8TQNOiz1oxfSB7TYxJz5P0JknX2r6qWPbxiDinc10CehZ5BKqDPALVQR6B6iCPQDWQRaBFXBMMAHpERFwoiRFMoALII1Ad5BGoDvIIVAd5BKqBLAKt4/QqAAAAAAAAAAAAdB0GwQAAAAAAAAAAANB1GAQDAAAAAAAAAABA12EQDAAAAAAAAAAAAF2nv9MdAABUVIS0YmWeUlmqJFW9GmxfxloDGWvlrpZPvuNwlmpKtlpozsxDrs9W67afPCtbrW0OXZCtVt8DD2erNe6BbKW0obbIVmvp0yZkqzXz1/k+DW/PVglorwO33C1brYVfe062Wv9w9F+z1UJztvzyRdlq9W01K1utNffdn63Wxlc+mK1WTJuWrdbNh26SrdaKzTbMVuvNx+XbzpH+PWOt7rflwfOz1br7ZztnqzXr0Fuy1dr685dmqzVu86dlq7Vmq3y1lm8+KVutHU/aO1st6QMZa3W/mV/It36sqm1ee22nu1DXPUfvk63WFl/L9/9494s2zVZLw+ye4EwwAAAAAAAAAAAAdB0GwQAAAAAAAAAAANB1GAQDAAAAAAAAAABA12EQDAAAAAAAAAAAAF2HQTAAAAAAAAAAAAB0HQbBAAAAAAAAAAAA0HUYBAMAAAAAAAAAAEDXYRAMAAAAAAAAAAAAXYdBMAAAAAAAAAAAAHQdBsEAAAAAAAAAAADQdRgEAwAAAAAAAAAAQNdhEAwAAAAAAAAAAABdh0EwAAAAAAAAAAAAdB0GwQAAAAAAAAAAANB1GAQDAAAAAAAAAABA1+lvR1GvlCYtdvaaAJq3pt9avlFf9poAAAAAAAAAAFQZZ4IBAAAAAAAAAACg6zAIBgAAAAAAAAAAgK7TlukQAQBdIEKxalWWUjkn0IyMtao6sWfOCUwHMtbKXS0fjunpBlv/y7XZat151jOz1Zp5yPXZavVvs1W2WuPvfCBbrb5HpmSr9cCcjbLVAnrRxPvyrdP6Z2+drZZuzVeqqla89DnZag385rJsteLxJdlqLX9+vvXjxKtvz1Zr5R7bZ6v19K/+LVutrDbbOFupa7JVQrO2PHh+tlr7XZsv2+c/K9+23NLdt8lWa+Jf8uVxYPx22WrN/mS+z+iF2SplNnWS1uy+W5ZS4y68Kkud3B49dK9stTY446/ZauW0/GX5tk22+NpF2WrltPnX8/VruD0K7DUCAAAAAAAAAABA12EQDAAAAAAAAAAAAF2HQTAAAAAAAAAAAAB0HQbBAAAAAAAAAAAA0HUaHgSz3Wf7Stv/184OARgZeQQAAAAAAAAAYHjNnAn2PkkL2tURAE0hjxg1BlGB6iCPQDWQRaA6yCNQHeQRqA7yCIxeQ4NgtmdJermkE9vbHQAjIY/IgEFUoDrII1ANZBGoDvIIVAd5BKqDPAKj1OiZYHMlfVjSmqEa2D7S9jzb81YvXZKjbwDqm6sm8rhqGXnEUxhEBaqDPALVQBaB6iCPQHWQR6A6yCPQmhEHwWy/QtK9EXH5cO0i4oSImBMRc/omTcnWQQBPGU0e+yeSR6xlrkYYRAWw3swVeQSqYK7IIlAVc0UegaqYK/IIVMVckUdg1Bo5E+x5kl5pe5GkH0p6oe3vt7VXAIZCHjFqjQyils8iXBHL1mPvgN7SbB5Xavl67B3QOxo9wIg8Au1HHoHqII9AdYwmjytWMisUUDbiIFhEfCwiZkXEbEmvk3R+RLyx7T0DsA7yiBaNOIhaPotwwBM70UegVzSVx/Ga0Ik+Ar2goQOMyCOwXpBHoDrII1AdTedxYDyzQgFljV4TDAAwxjGIClQHeQSqgSwC1UEegeogj0B1kEegdf3NNI6ICyRd0JaeAGgKeQQAAAAAAAAAYGhNDYIBALoDg6hAdZBHoBrIIlAd5BGoDvIIVAd5BEaH6RABAAAAAAAAAADQdRgEAwAAAAAAAAAAQNdhEAwAAAAAAAAAAABdh0EwAAAAAAAAAAAAdJ3+dhQdt0qafN+a7DUBNG9Nv/TEZnnHu9e05ZMDlbR6dZYykaVK4oy1qtqvnPoy1hrIWCt3NaBs5iHXZ6v16Vsuz1brmO2ylVLfhtOz1cq5lbDhTRMzVgN6z6zPXZSt1o1f3DtbLX04X6mqGvjNZdlquT/fF6bVDz2Urdb48/Kt0+46ep9stbb4Wr6/e281K1utNQ8/kq3W3S/fPFst3ZCvFDrn/GdNyVbrxm/uma3WDkddmq2Wpk3LVmrcZfm27+9/a8b143fOzFcrp8eXatyFV3W6F221wRl/7XQX2m7Cr/Ntm1TVzcdlzOMHh84jZ4IBAAAAAAAAAACg6zAIBgAAAAAAAAAAgK7DIBgAAAAAAAAAAAC6DoNgAAAAAAAAAAAA6DoMggEAAAAAAAAAAKDrMAgGAAAAAAAAAACArsMgGAAAAAAAAAAAALoOg2AAAAAAAAAAAADoOgyCAQAAAAAAAAAAoOswCAYAAAAAAAAAAICuwyAYAAAAAAAAAAAAug6DYAAAAAAAAAAAAOg6DIIBAAAAAAAAAACg6zAIBgAAAAAAAAAAgK7DIBgAAAAAAAAAAAC6DoNgAAAAAAAAAAAA6DoMggEAAAAAAAAAAKDrMAgGAAAAAAAAAACArtPf6Q4AACoqQlqzptO9WEdkrOWMtarar5z6MtYayFgrdzWg7Jjt9shW68ZvPydbrWd85MZsteKJZdlqjb/zgWy1ALRm06vzbZ3ckq1SXrY1buLEPLW22zpLHUlaPT/fZ3TfZptlq6VVq7KV2vLEa7PVWn5AvnVt392PZasVi+/JVmvm92/KVuu6bJXQLXY46tJstRb9597Zam177GXZauW0yYkXd7oLyGDcrs/IVmvN1Quy1eoFj/3rXtlqPf1D+fJ46zD3cSYYAAAAAAAAAAAAug6DYAAAAAAAAAAAAOg6bZkOcdyq0ISHVmevCaB5MV5aOiNvfmJ81nIAAAAAAAAAAGTHmWAAAAAAAAAAAADoOgyCAQAAAAAAAAAAoOswCAYAPcT2hrbPtP032wts793pPgG9ijwC1UEegeogj0B1kEegGsgi0Jq2XBMMAFBZx0v6TUT8s+0BSZM73SGgh5FHoDrII1Ad5BGoDvIIVANZBFrAIBgA9Ajb0yW9QNLhkhQRKySt6GSfgF5FHoHqII9AdZBHoDrII1ANZBFoHdMhAkDv2FbSfZK+Z/tK2yfantLpTgE9ijwC1UEegeogj0B1kEegGsgi0CIGwQCgd/RLerakb0bE7pKWSPpouYHtI23Psz1vhZZ3oo9Ar2gqjyvJI9BO5BGoDrZXgepg/QhUw4hZlMgjMJyGBsG4+B5QHeQRLbhT0p0RcUnx+5lKG1JPiogTImJORMwZ0IT13kGghzSVx/HkEWgn8ghUB9urQHWwfgSqYcQsSuQRGE6jZ4INXnxvJ0m7SlrQvi4BGAF5xKhExD2S7rC9Y7HoRZLmd7BLQM8ij0B1kEegOsgjUB3kEagGsgi0rn+kBlx8D6gO8ogM3iPpdNsDkm6R9JYO9wfoZeQRqA7yCFQHeQSqgzwC1UAWgRaMOAimtS++t6ukyyW9LyKWlBvZPlLSkZI0YeKGmbsJoNB0Hvs33Gi9dxLVFRFXSZrT6X4AII9AlZBHoDrII1Ad5BGoBrIItKaR6RAbuvjeWvOODkzJ3E0Ahabz2DeFPAIAAAAAAAAAek8jg2ANXXwPwHpBHgEAAAAAAAAAaMCIg2BcfA+oDvIIAAAAAAAAAEBjGrkmmMTF94AqIY8AAAAAAAAAAIygoUEwLr4HVAd5xPoSkmL1miy1nKVKfpGxVs7XWNV+5dSXsdZAxlpAO+3w9suy1br12H2y1dr62Iuy1XJ/o8fYAWi3qOpGQEYRoTXLluUpNv/GPHWU97MwHnssX61Vq7LV6tts02y1Bi64OlutJw7cPVutyYunZasVMzbJVkv35isF1Jr9yYuz1brrZztnq7XlwfkmQcq6vboyXyk0Z83VCzrdhZ417Ud/zVbr4cP2zlZLp5w55F2NXBMMAAAAAAAAAAAAGFMYBAMAAAAAAAAAAEDXYRAMAAAAAAAAAAAAXYdBMAAAAAAAAAAAAHQdBsEAAAAAAAAAAADQdfrbUdSrQwOPrMheE8AoTFgjzV6SvyYAAAAAAAAAABXGmWAAAAAAAAAAAADoOgyCAQAAAAAAAAAAoOswCAYAAAAAAAAAAICuwyAYAAAAAAAAAAAAug6DYAAAAAAAAAAAAOg6DIIBAAAAAAAAAACg6zAIBgAAAAAAAAAAgK7DIBgAAAAAAAAAAAC6DoNgAAAAAAAAAAAA6DoMggEAAAAAAAAAAKDrMAgGAAAAAAAAAACArsMgGAAAAAAAAAAAALpOf6c7AACoqpBiTZ5Kq7OUkSQ5X6msImOtnK+xqv3Kqa/THQA6YOtjL8pW66aT98hWa8d3LshWSyvzlQJ60cSHM26AoSmxalW2Wu7Pt9smZ79WLb4nW62cr3HKNYuz1Vo9e8tstfrueSBbLWCs2PLg+dlq3fbpfbLV2vbL12arxfYq0JrJf18/IeJMMAAAAAAAAAAAAHQdBsEAAAAAAAAAAADQdRgEAwAAAAAAAAAAQNdhEAwAAAAAAAAAAABdh0EwAAAAAAAAAAAAdB0GwQCgh9g+2vb1tq+zfYbtiZ3uE9CryCNQHeQRqA7yCFQHeQSqgSwCrWEQDAB6hO2Zkt4raU5E7CKpT9LrOtsroDeRR6A6yCNQHeQRqA7yCFQDWQRa19+Ool4T6nt8efaaAJo3fcJSHbT9/Kw1fzRhadZ6WK/6JU2yvVLSZEl3d7g/QC8jj0B1kEegOsgjUB3kEagGsgi0gDPBAKBHRMRdkr4s6XZJiyU9EhG/7WyvgN5EHoHqII9AdZBHoDrII1ANZBFoHYNgANAjbG8k6VWStpW0paQptt9Y0+ZI2/Nsz1sZec/oBfCUpvMo8gi0C3kEqoM8AtVBHoFqaCSLRTvyCAyBQTAA6B0HSLo1Iu6LiJWSzpK0T7lBRJwQEXMiYs54T+hIJ4Ee0VweRR6BNiKPQHWQR6A6yCNQDSNmUSKPwHAYBAOA3nG7pL1sT7ZtSS+StKDDfQJ6FXkEqoM8AtVBHoHqII9ANZBFoEUMggFAj4iISySdKekKSdcqrQNO6GingB5FHoHqII9AdZBHoDrII1ANZBFoXX+nOwAAWH8i4hhJx3S6HwDII1Al5BGoDvIIVAd5BKqBLAKt4UwwAAAAAAAAAAAAdJ2GBsFsH237etvX2T7D9sR2dwxAfeQRAAAAAAAAAICRjTgIZnumpPdKmhMRu0jqk/S6dncMwLrIIwAAAAAAAAAAjWl0OsR+SZNs90uaLOnu9nUJwAjIIwAAAAAAAAAAI+gfqUFE3GX7y5Jul7RU0m8j4re17WwfKelISZo4foPc/QSg0eVx2haT128n0T1CijWRpZTHrclSR5JidbZScr5SWeV515Ocr7Gq/QLQmu0Pvzxbrds/tk+2WvrcqflqAT1o8q2PdLoLPatvo42y1Vr90EPZauXUv+022WqtuvW2bLXi8SX5at1xZ7Za2mpWvlpAD9rmmIuy1brnnRm3V7+RrxTQi/qXZdzJN4xGpkPcSNKrJG0raUtJU2y/sbZdRJwQEXMiYs5A/5T8PQUwqjxO2mjC+u4mAAAAAAAAAAAd18h0iAdIujUi7ouIlZLOkpRxyBxAE8gjAAAAAAAAAAANaGQQ7HZJe9mebNuSXiRpQXu7BWAI5BEAAAAAAAAAgAaMOAgWEZdIOlPSFZKuLR5zQpv7BaAO8ggAAAAAAAAAQGP6G2kUEcdIOqbNfQHQAPIIAAAAAAAAAMDIGpkOEQAAAAAAAAAAABhTGAQDAAAAAAAAAABA12EQDAAAAAAAAAAAAF2HQTAAAAAAAAAAAAB0HQbBAAAAAAAAAAAA0HX621J1zRr5ieXZawJo3lbjn9DcGfOy1rxo/BNZ6wEAAAAAAAAAkBtnggEAAAAAAAAAAKDrMAgGAAAAAAAAAACArsMgGAAAAAAAAAAAALoOg2AAAAAAAAAAAADoOgyCAQAAAAAAAAAAoOswCAYAAAAAAAAAAICuwyAYAAAAAAAAAAAAuo4jInvR6RO3iL23PixrzYtvP1WPLLvHWYsCPWDOrhPj0nO3ylpzzwPv0Lyrl5HHLmf7Pkm3NdB0U0n3t7k7o0G/mkO/km0iYrP1+HwNIY9tQ7+aQx5FHtuIfjWHPIo8thH9ag55FHlsI/rVHPIo8thG9Ks5lclj/3rsBABgDGl0Q872vIiY0+7+NIt+NYd+VRt5bA/61Zyq9mt9I4/tQb+aU9V+rW/ksT3oV3Oq2q/1jTy2B/1qTlX7tb6Rx/agX82pUr+YDhEAAAAAAAAAAABdh0EwAAAAAAAAAAAAdB0GwQAArTqh0x0YAv1qDv3qDlV9v+hXc+hXd6jq+0W/mkO/ukNV3y/61Rz61R2q+n7Rr+bQr+5Q1feLfjWHfo3AEZG96PSJW8TeWx+WtebFt5+qR5bd46xFgR4wZ9eJcem5W2WtueeBd2je1cvIIwAAAAAAAACgsjgTDAAAAAAAAAAAAF2HQTAAAAAAAAAAAAB0HQbBAACjYvultm+wvdD2RzvdH0myvZXtP9ieb/t62+/rdJ/KbPfZvtL2/3W6L4Nsb2j7TNt/s73A9t6d7pMk2T66+D+8zvYZtid2uk9VRh6bRx4bRx6bQx6bRx4bRx6bQx6bRx4bRx6bQx6bRx4bRx6bQx6bRx4bV8U8MggGAGia7T5J35D0Mkk7SzrU9s6d7ZUkaZWkD0TEzpL2kvSuivRr0PskLeh0J2ocL+k3EbGTpF1Vgf7ZninpvZLmRMQukvokva6zvaou8jhq5LEB5LE55HHUyGMDyGNzyOOokccGkMfmkMdRI48NII/NIY+jRh4bUNU8MggGABiNPSUtjIhbImKFpB9KelWH+6SIWBwRVxT/fkxpA2BmZ3uV2J4l6eWSTux0XwbZni7pBZK+K0kRsSIiHu5op57SL2mS7X5JkyXd3eH+VBl5bBJ5bBp5bBx5bBJ5bBp5bBx5bBJ5bBp5bBx5bBJ5bBp5bBx5bBJ5bFrl8sggGABgNGZKuqP0+52qyMbJINuzJe0u6ZIOd2XQXEkflrSmw/0o21bSfZK+V5zWf6LtKZ3uVETcJenLkm6XtFjSIxHx2872qtLIY/Pmijw2hDw2jTw2b67IY0PIY9PIY/Pmijw2hDw2jTw2b67IY0PIY9PIY/Pmijw2pKp5ZBAMANB1bE+V9FNJ74+IRyvQn1dIujciLu90X2r0S3q2pG9GxO6Slkjq+HzgtjdSOhJtW0lbSppi+42d7RVGizw2jDyi7chjw8gj2o48Now8ou3IY8PII9qOPDaMPDaBQTAAwGjcJWmr0u+zimUdZ3u80gbT6RFxVqf7U3iepFfaXqQ01cALbX+/s12SlI74ujMiBo+uOlNpI6rTDpB0a0TcFxErJZ0laZ8O96nKyGNzyGNzyGNzyGNzyGNzyGNzyGNzyGNzyGNzyGNzyGNzyGNzyGNzyGNzKpnH/nYUfXT53x8/96bjbmig6aaS7m+w7I4tdAnoWZdfs/zxvhkLySNyu0zS9ra3VdpYep2k13e2S5JtK82HvCAivtrp/gyKiI9J+pgk2d5P0gcjouNHwkTEPbbvsL1jRNwg6UWS5ne6X0qnze9le7KkpUr9mtfZLlUaeWwCeWwaeWwOeWwCeWwaeWwOeWwCeWwaeWwOeWwCeWwaeWwOeWwCeWxaJfPYlkEwSTdExJyRGtme10i7wbatdwvoSeQR2UXEKtvvlnSupD5JJ0XE9R3ulpSO0HmTpGttX1Us+3hEnNO5LlXeeySdbntA0i2S3tLh/igiLrF9pqQrJK2SdKWkEzrbq+oij12FPI5x5LGrkMcxjjx2FfI4xpHHrkIexzjy2FXIY4McEfmLNrgzvdmd7o22BfAU8ggAAAAAAAAA6EVcEwwAAAAAAAAAAABdp12DYI2e4tbMqXAdP20OGKPIIwAAAAAAAACg57RlOkQAAAAAAAAAAACgk5gOEQAAAAAAAAAAAF2HQTAAAAAAAAAAAAB0nSyDYLY3tv072zcVPzequf+ltm+wHbbvtn2V7V+U7p9g+0e2F9r+m+2bi39/tM5zHW77vqLGVbbfmuM1AN2CPAIAAAAAAAAAkO9MsI9K+n1EbC/p98XvkiTbfZK+IellkpZIuk/S6yPilaXHHyHpIUk7StpI0gJJO0s61PbOdZ7vRxGxW3E7MdNrALoFeQQAAAAAAAAA9Lxcg2CvknRK8e9TJL26dN+ekhZGxC3F7z8s2td7/J6Srpa0l6SVQ7QFMDzyCAAAAAAAAADoebkGwTaPiMXFv++RtHnpvpmS7ij+PVHS2yS9z/ar67SZKel2SY9I2kTSncWyWq+xfY3tM21vlek1AN2CPAIAAAAAAAAAel5/ow1tnydpizp3faL8S0SE7RiizDaSXijpxZLm2r42Im5utA+FX0o6IyKW23670hkrL2yyBjCmkUcAAAAAAAAAAIbX8CBYRBww1H22/257RkQstj1D0r2lu++StFVR4y7bs5SuMbRG0u6Sbi61uUvS1pKmS3pA0qxiWbkfD5R+PVHSlxp9DUC3II8AAAAAAAAAAAwv13SIv5D05uLfb5b089J9l0na3vautqdKep2kP0l6nqT5NY+/TNKuki6RNL5o+4vyExU79Qe9UmkHPoCnkEcAAAAAAAAAQM9zxFAzpTVRxN5E0o+Vzhq5TdJrI+JB23MkvUPSWZK+pXRtogcl3S/pVkknRsQvbE+UdJrSmSirJQ0onZlyUkT8l+3PSJpXtP280s72VUWtoyLiby2/CKBLkEcAAAAAAAAAADINggEAAAAAAAAAAABVkms6RAAAAAAAAAAAAKAyGAQDAAAAAAAAAABA12EQDAAAAAAAAAAAAF2HQTAAAAAAAAAAAAB0HQbBAAAAAAAAAAAA0HUYBAMAAAAAAAAAAEDXYRAMAAAAAAAAAAAAXef/AzR94w8KfBbGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 2304x288 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "visualize_couplings(couples_nonmonotone, p_logits.numpy(), q_logits.numpy(), R_non_monotone, suffix='R_non_monotone')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------  trial:  0 -----------\n",
      "Monotonic reward function:\n",
      " indep: 4.945008733341747, gmx: 4.317262516278003,  icdf: 0.8556931297995052, jp: 2.5089839048180056, jpz: 3.556820869445801,  opt: 0.8556931297995052\n",
      "Non-monotonic reward function:\n",
      " indep: 0.5155141654030311, gmx: 0.27721546732525076,  icdf: 0.537890074085093, jp: 0.13227421161822375, jpz: 0.13969506323337555,  opt: 0.033353026616698034\n",
      "-------  trial:  1 -----------\n",
      "Monotonic reward function:\n",
      " indep: 0.9346426465980064, gmx: 0.8546105420829483,  icdf: 0.11761825901602112, jp: 0.42962005012723403, jpz: 0.5831938982009888,  opt: 0.11761825901602112\n",
      "Non-monotonic reward function:\n",
      " indep: 0.9601617302853735, gmx: 0.6535214406187828,  icdf: 1.2309370924190854, jp: 0.3053650809559505, jpz: 0.17960374057292938,  opt: 0.08596162860123933\n",
      "-------  trial:  2 -----------\n",
      "Monotonic reward function:\n",
      " indep: 1.8068711050510027, gmx: 1.5260233996375931,  icdf: 0.41109068922261677, jp: 0.8936223986409814, jpz: 1.2781264781951904,  opt: 0.41109068922261677\n",
      "Non-monotonic reward function:\n",
      " indep: 0.7464320514487729, gmx: 0.43387548944713755,  icdf: 0.8954337832074225, jp: 0.2685935035218696, jpz: 0.198048934340477,  opt: 0.10938590499022044\n",
      "-------  trial:  3 -----------\n",
      "Monotonic reward function:\n",
      " indep: 2.722151051032917, gmx: 2.4149199385848386,  icdf: 0.84703061732004, jp: 1.6486978729859116, jpz: 2.230787992477417,  opt: 0.84703061732004\n",
      "Non-monotonic reward function:\n",
      " indep: 1.360146100445662, gmx: 0.8489671426380427,  icdf: 1.7740400690626916, jp: 0.5155387197966199, jpz: 0.4703139066696167,  opt: 0.34413860746512603\n",
      "-------  trial:  4 -----------\n",
      "Monotonic reward function:\n",
      " indep: 3.7811413928869695, gmx: 3.4387353196318826,  icdf: 0.9929884377171265, jp: 2.043034328752582, jpz: 2.888115167617798,  opt: 0.9929884377171265\n",
      "Non-monotonic reward function:\n",
      " indep: 0.23750453516659553, gmx: 0.15578766259826096,  icdf: 0.26353876827076866, jp: 0.0834774825826208, jpz: 0.07155356556177139,  opt: 0.03937486239919723\n",
      "-------  trial:  5 -----------\n",
      "Monotonic reward function:\n",
      " indep: 3.850818492947021, gmx: 3.4873245217173934,  icdf: 1.0301228395183677, jp: 2.0930390710179907, jpz: 2.8240904808044434,  opt: 1.0301228395183677\n",
      "Non-monotonic reward function:\n",
      " indep: 0.749802489440765, gmx: 0.4471540654341502,  icdf: 1.0006039688190393, jp: 0.25252284986123086, jpz: 0.20250333845615387,  opt: 0.06378512842054893\n",
      "-------  trial:  6 -----------\n",
      "Monotonic reward function:\n",
      " indep: 2.137653966891296, gmx: 1.9445716373728628,  icdf: 0.3543150213073554, jp: 0.9635957273038864, jpz: 1.4162511825561523,  opt: 0.3543150213073554\n",
      "Non-monotonic reward function:\n",
      " indep: 0.7387997766870625, gmx: 0.5805294940682612,  icdf: 0.4871678251019081, jp: 0.39027322095525474, jpz: 0.3909716308116913,  opt: 0.18446726629055923\n",
      "-------  trial:  7 -----------\n",
      "Monotonic reward function:\n",
      " indep: 3.49308644352356, gmx: 3.291803273191893,  icdf: 0.9551504915829809, jp: 1.9989202393546501, jpz: 2.7203593254089355,  opt: 0.9551504915829809\n",
      "Non-monotonic reward function:\n",
      " indep: 0.7561894095480113, gmx: 0.40496171107487094,  icdf: 0.8110921032699724, jp: 0.18887109990737178, jpz: 0.08131907135248184,  opt: 0.04767400521516415\n",
      "-------  trial:  8 -----------\n",
      "Monotonic reward function:\n",
      " indep: 2.6759320605041683, gmx: 2.5110446537153224,  icdf: 0.38617316835196597, jp: 1.2768049915989295, jpz: 1.8324733972549438,  opt: 0.38617316835196597\n",
      "Non-monotonic reward function:\n",
      " indep: 1.381363403235935, gmx: 0.8160989965439778,  icdf: 1.1825461222806628, jp: 0.27771302196032216, jpz: 0.23069524765014648,  opt: 0.11009975833023747\n",
      "-------  trial:  9 -----------\n",
      "Monotonic reward function:\n",
      " indep: 0.9875590053181937, gmx: 0.8809997633826474,  icdf: 0.13086551364733606, jp: 0.4119320823895545, jpz: 0.694908618927002,  opt: 0.13086551364733606\n",
      "Non-monotonic reward function:\n",
      " indep: 0.8994723425329494, gmx: 0.45081791815562655,  icdf: 0.9364692643855619, jp: 0.27484657854169053, jpz: 0.1907193660736084,  opt: 0.1399817781902644\n"
     ]
    }
   ],
   "source": [
    "var_mono = {'indep': [],\n",
    "            'gmx': [],\n",
    "            'icdf': [],\n",
    "            'jp': [],\n",
    "            'jpz': [],\n",
    "            'opt': []\n",
    "            }\n",
    "\n",
    "var_nonmono = {'indep': [],\n",
    "               'gmx': [],\n",
    "               'icdf': [],\n",
    "               'jp': [],\n",
    "               'jpz': [],\n",
    "               'opt': []\n",
    "               }\n",
    "for trial in range(10):\n",
    "    print('-------  trial: ', trial, '-----------')\n",
    "    torch.manual_seed(trial)\n",
    "    np.random.seed(trial)\n",
    "    RANDOM_MONOTONE_R = create_random_monotone()\n",
    "    RANDOM_NON_MONOTONE_R = create_random_non_monotone()\n",
    "\n",
    "    p_logits = (torch.arange(K, 0, step=-1).float()/4).unsqueeze(0)\n",
    "    q_logits = (torch.arange(K).float()/4 ).unsqueeze(0)\n",
    "\n",
    "    gd1_non_monotone = models.JointPredictorGivenLogits(K, 1024)\n",
    "    gd2_non_monotone = JP_Z.Gadget2Jax(s_dim=K, z_dim=20, hidden_features=[1024,1024], tmp=1.0)\n",
    "\n",
    "    train_gsm_with_gadget1(3000, gd1_non_monotone, p_logits, q_logits, R_non_monotone, tau=1.0,\n",
    "                           lr=learning_rate,verbose=False)\n",
    "    train_gsm_with_gadget2(3000, gd2_non_monotone, p_logits, q_logits, R_non_monotone, s_prime_obs=None,\n",
    "                           verbose=False, seed=trial)\n",
    "\n",
    "    gd1_monotone = models.JointPredictorGivenLogits(K, 1024) #models.FreeParamsPredictor(K)\n",
    "    gd2_monotone = JP_Z.Gadget2Jax(s_dim=K, z_dim=20, hidden_features=[1024,1024], tmp=1.0)\n",
    "\n",
    "    train_gsm_with_gadget1(3000, gd1_monotone, p_logits, q_logits, R_monotone, tau=1.0, lr=learning_rate, verbose=False)\n",
    "    train_gsm_with_gadget2(3000, gd2_monotone, p_logits, q_logits, R_monotone, s_prime_obs=None,\n",
    "                            verbose=False, seed=trial)\n",
    "\n",
    "    n = 2000\n",
    "    p_logits, q_logits = p_logits.squeeze(0), q_logits.squeeze(0)\n",
    "\n",
    "    #### plot R(a) = np.cumsum(np.random.rand(K))[a] ###\n",
    "    ates, _ = compare_coupling_strategies(n, p_logits, q_logits, R_monotone,\n",
    "                                          gd1_monotone, gd2_monotone.model,\n",
    "                                          gd2_monotone.params, monotone=True)\n",
    "    ate_i, ate_g, ate_c, ate_m, ate_z, ate_o = zip(*ates)\n",
    "\n",
    "    var_mono['indep'].append(np.var(ate_i))\n",
    "    var_mono['gmx'].append(np.var(ate_g))\n",
    "    var_mono['icdf'].append(np.var(ate_c))\n",
    "    var_mono['jp'].append(np.var(ate_m))\n",
    "    var_mono['jpz'].append(np.var(ate_z))\n",
    "    var_mono['opt'].append(np.var(ate_o))\n",
    "    print('Monotonic reward function:')\n",
    "    print(f\" indep: {var_mono['indep'][-1]}, gmx: {var_mono['gmx'][-1]},  icdf: {var_mono['icdf'][-1]},\"\n",
    "          f\" jp: {var_mono['jp'][-1]}, jpz: {var_mono['jpz'][-1]},  opt: {var_mono['opt'][-1]}\")\n",
    "    \n",
    "    \n",
    "    #### plot R(a) = sin(30*(np.random.normal(loc=np.zeros(K), scale=1)))[a] ###\n",
    "    ates, _ = compare_coupling_strategies(n, p_logits, q_logits, R_non_monotone,\n",
    "                                          gd1_non_monotone, gd2_non_monotone.model,\n",
    "                                          gd2_non_monotone.params, monotone=False)\n",
    "    ate_i, ate_g, ate_c, ate_m, ate_z, ate_o = zip(*ates)\n",
    "\n",
    "    var_nonmono['indep'].append(np.var(ate_i))\n",
    "    var_nonmono['gmx'].append(np.var(ate_g))\n",
    "    var_nonmono['icdf'].append(np.var(ate_c))\n",
    "    var_nonmono['jp'].append(np.var(ate_m))\n",
    "    var_nonmono['jpz'].append(np.var(ate_z))\n",
    "    var_nonmono['opt'].append(np.var(ate_o))\n",
    "    print('Non-monotonic reward function:')\n",
    "    print(f\" indep: {var_nonmono['indep'][-1]}, gmx: {var_nonmono['gmx'][-1]},  icdf: {var_nonmono['icdf'][-1]},\"\n",
    "          f\" jp: {var_nonmono['jp'][-1]}, jpz: {var_nonmono['jpz'][-1]},  opt: {var_nonmono['opt'][-1]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "def conv_to_np(this_list):\n",
    "    this_arr = np.array(this_list)[:, np.newaxis]\n",
    "    # Make this idempotent\n",
    "    this_arr = this_arr.squeeze()[:, np.newaxis]\n",
    "    return this_arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def plot_ate_var(variance, title='Monotonic R'):\n",
    "    gmx = conv_to_np(variance['gmx'])\n",
    "    icdf = conv_to_np(variance['icdf'])\n",
    "    #anti = conv_to_np(variance['anti'])\n",
    "    jp = conv_to_np(variance['jp'])\n",
    "    jpz = conv_to_np(variance['jpz'])\n",
    "    opt = conv_to_np(variance['opt'])\n",
    "\n",
    "    plt.rcParams.update({'font.size': 14})\n",
    "    reward = np.concatenate([gmx,\n",
    "                             icdf,\n",
    "                             #anti,\n",
    "                             jp,\n",
    "                             jpz,\n",
    "                             opt,\n",
    "                             ], axis=1)\n",
    "    reward_df = pd.DataFrame(reward, columns=['Gumbel-Max',\n",
    "                                              'Inverse-CDF',\n",
    "                                              #'Antithetic',\n",
    "                                              'Gadget-1',\n",
    "                                              'Gadget-2',\n",
    "                                              'Optimal (LP)'\n",
    "                                              ])\n",
    "\n",
    "    plt.figure(figsize=(8, 5))\n",
    "    #sns.boxplot(data=reward_df, whis=[2.5, 97.5])\n",
    "    sns.barplot(data=reward_df, ci=68)\n",
    "    plt.title(title)\n",
    "    plt.ylabel(\"ATE variance\")\n",
    "    plt.savefig(\"{}/toy_ate_var_{}_{}.pdf\".format(figpath, coupler_save_path, title))\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfwAAAFHCAYAAABEX/6lAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnsUlEQVR4nO3de5wcVZnw8d8DyHLzAgYJ4IYosIKiosYrgoiiKIIIiDfU6KuICIoIUZBV0PUWEGQVl4XdFZGFl8XLqyCIggZQEQniKogIQrgEJiRyE0iAJM/7x6nRSqdnpifp7rnU7/v59GemTp2qevp0dT1dVaeqIjORJEmT2xpjHYAkSeo9E74kSQ1gwpckqQFM+JIkNYAJX5KkBjDhS5LUACZ8SRNGRGREHDPWcUgTkQlfGiciYmaV0DIidhyizk3V+Dl9Dq+tiHhmRBwTEdPHOpbRiog5tfbOiFgcEb+LiEMjwm2jJp21xjoASStZArwduLxeGBEvAbasxo8XzwQ+DcwB5vVheesCS7s4v7uAWdX/UyjtfiKwMfDJLi5HGnP+ipXGnwuAN0fE41rK3w78Efhz/0MaHzJzSWZ2M+E/kJlnVq+vADsBtwKHRMSaXVyONOZM+NL4czawEfDawYIq+bwFOKvdBBGxXkQcFxG3RcQjEXFjRHyi9dB0dej6lIjYKyKurepeFxG7tZnncyPigoh4ICIeqg6B71gbPxM4txr8We3Q+MxanX0iYm51uPwvEXF2RExrWc7pEbEkIjaPiP8XEQ9GxMKIOL416bY7hx8RT6ze+83V+5kfEWdFxObDtHFbmbkEuAp4PPCU0U4vjWcmfGn8uYNyOP/ttbJXUxLQSgk/IgL4f8DhwE+AjwK/A74AfL3N/F8KnAz8X8rh7HWA70TEk2vz3LaK4XnAccAxwGbAxRGxU1XtMuBfq/8/D7yzel1WzWN/4NvV+COBU4A3AL+IiCktMa0B/Aj4S/U+LgU+BhzQJv76e1+/qvtR4KfAR6r3vAWw1XDTDmM6kMB9qzi9ND5lpi9fvsbBC5hJSTQvAT4APASsX407A/hV9f+1wJzadHtW0326ZX7fqMq3q5Ul8CiwVa3sOVX5wbWy71b1tq6VTQEWAXNrZftW0+7csuzHAQPAH4B1a+U7V/WPr5WdXpV9qmUev6kvqxb/MbXhY6qyN7dpzxihvecAN1bvawrwDGB2Nb/zx3p98OWr2y/38KXx6VxK0twrItYF9gL+e4i6uwPLgZNayr9cG1/3s8y8aXAgM38HPAA8Hf52+uC1wHmZeWOt3iJKcn5BRGwyQvwzgE2Af8vMxbV5zAGubhMTwGktw5cPxjSMfYHrMvPc1hGZ2cmjQLcCFlavPwJHAN8H3t3BtNKEYsKXxqHMvAe4CHgHZQ9+PeCcIapvASzIzPtaym+g/BCY3lJ+W5t53AtsWP2/cbW8G9rUu7762zrPdjENxtBuHq3TP5aZdw0T01C2pBzxWFW3A7sCrwE+WA1PARYPN5E0EXlZnjR+nUU5lP8E4OLMvLtL8102RHl0af6rYvkYLffhzLy4+v8nEXE58FtKn4RDxygmqSfcw5fGr+8DjwA7METv/MqtwCYR8cSW8n+ifMfnjXK5C4GHKee0W21T/R2c51CHzW+t/g41j9HGNJQ/A9t1aV5k5nXAN4EPtl5NIE10JnxpnMrMhymHmY8FvjdM1fMp3+UPt5QfVv394SiXu4zSY36PiNhysDwiNqKc256bmQuq4oeqv62H3ucCC4APRMQ6tXnsSDm/f/5oYhrGt4FnRcSbW0dUVy+siuMo/Sc+tjqBSeONh/SlcSwzv9VBtfMpl+MdGxFbUHq37wLsA/x7Zq7KOe6jKee1fx4RJ1Pu7vd+4EmUjnKDrqGcIjgyIp5EOfd9ZWbeEhFHUE5JXB4RZ1L6BnwYmA98aRViauc4yvs8OyJeQ+kQ+CTgdcCnKJfsjUpm3hAR5wPvi4jPVp0VpQnPPXxpgqt6o78JOAHYDfgKsD1wFPChVZzn9cDLKQn945SjDAPAqzPzslq9BZQfAhtSetmfDbyiGvctyo+DoCT4gyh3EdyhW0k0Mx+i3B3vZMp7/1fgYMq9DG4cZtKRHE/puPiR1Y1RGi+isytXJEnSROYeviRJDWDClySpAUz4kiQ1gAlfkqQGMOFLktQAk/Y6/ClTpuT06dPHOgxJkvrm6quvXpSZG7cbN2kT/vTp05k7d+5YhyFJUt9ExK1DjfOQviRJDWDClySpAUz4kiQ1gAlfkqQGMOFLktQAJnxJkhrAhC9JUgOY8CVJagATviRJDWDClySpASbtrXUnmlmzZjEwMMDUqVOZPXv2WIcjSZpkTPjjxMDAAPPnzx/rMCRJk5SH9CVJagATviRJDWDClySpAUz4kiQ1gAlfkqQGMOFLktQAJnxJkhrAhC9JUgOY8CVJagATviRJDWDClySpAUz4kiQ1gAlfkqQGMOFLktQAJnxJkhrAhC9JUgOY8CVJagATviRJDWDClySpAUz4kiQ1gAlfkqQGMOFLktQAJnxJkhrAhC9JUgOY8CVJagATviRJDdC3hB8RR0bEVRHxQEQsjIjzImK7EaaZHhHZ5rVbv+KWJGky6Oce/s7A14GXAbsAS4GLI2KjDqbdDdi09vppj2KUJGlSWqtfC8rM19aHI+KdwP3ADsB5I0z+l8wc6FVskiRNdmN5Dv/x1fLv7aDudyPi7oj4RUTs2+O4JEmadMYy4Z8E/Ba4Ypg6DwKHA/sBrwcuAc6JiP3bVY6IAyJibkTMXbhwYZfDlSRp4urbIf26iDgBeDnw8sxcNlS9zFwEfLlWNDcipgCzgDPb1D8VOBVgxowZ2dWgJUmawPq+hx8RJwJvA3bJzJtXYRZXAlt3NypJkia3vu7hR8RJwFuAV2bmH1dxNtsDd3UjnhcccUY3ZtMVj1/0V9YEblv013ET19XHvWusQ5AkdUnfEn5EnAy8E9gLuDciplajHszMB6s6XwBelJmvqobfDTwGXAMsB/YAPgR8vF9xS5I0GfRzD/+g6u8lLeXHAsdU/28KbNky/mhgC2AZ8CfgvZm50vl7SZI0tH5ehx8d1JnZMvxN4Ju9ikmSpKbwXvqSJDWACV+SpAYw4UuS1AAmfEmSGsCEL0lSA5jwJUlqABO+JEkNMCYPz5EkrbpZs2YxMDDA1KlTmT179liHownChC9JE8zAwADz588f6zA0wXhIX5KkBjDhS5LUACZ8SZIawIQvSVIDmPAlSWoAE74kSQ1gwpckqQFM+JIkNYAJX5KkBjDhS5LUACZ8SZIawIQvSVIDmPAlSWoAE74kSQ1gwpckqQFM+JIkNYAJX5KkBjDhS5LUACZ8SZIawIQvSVIDrDXWAahYvvb6K/yVJKmbTPjjxENbv2asQ5AkTWIe0pckqQFM+JIkNYAJX5KkBuhbwo+IIyPiqoh4ICIWRsR5EbFdB9M9OyIujYjFETE/Ij4VEdGPmCVJmiz6uYe/M/B14GXALsBS4OKI2GioCSLiCcBPgAXAC4GPAEcAh/U6WEmSJpO+9dLPzNfWhyPincD9wA7AeUNM9g5gPeDdmbkYuDYitgEOi4gTMjN7GbMkSZPFWJ7Df3y1/HuHqfNS4PIq2Q+6CNgMmN670CRJmlzG8jr8k4DfAlcMU2cqcEdL2YLauFu6H5YkrezSnV4x1iH8zeK11oQIFt9xx7iJ6xWXXTrWIWgEY5LwI+IE4OXAyzNzWRfnewBwAMC0adO6NVtJkia8vh/Sj4gTgbcBu2TmzSNUHwA2aSnbpDZuBZl5ambOyMwZG2+88eoHK0nSJNFxwo+IgyLiuoh4OCKeXpV9IiL2G8U8TuLvyf6PHUxyBbBjRKxTK9sVuBOY1+lyJUlquo4SfkQcChwNnArUr4GfDxzc4TxOBt4DvB24NyKmVq8NanW+EBGX1CY7C3gYOD0itouIvYFPAPbQlyRpFDrdwz8QeH9mnkS5fn7Qb4BndTiPgyg98y8B7qq9Dq/V2RTYcnAgM++n7NFvBswFTga+DJzQ4TIlSRKdd9rbAri2TfljwLqdzCAzR7w7XmbObFP2e2CnTpYhSZLa63QP/2bg+W3KXw/8oXvhSJKkXuh0D/944GsRsR7lHP5LqzvlzQLe26vgJElSd3SU8DPzGxGxFvB5yq1uv0XpKf/hzDynh/FJkqQu6PjGO5l5GnBaREwB1sjMu3sXliRJ6qaOEn5EPAtYMzN/l5mLauXPAZZmpufxJUkaxzrttHcq0O7Z9c+sxkmSpHGs04T/HODXbcqvAp7dvXAkSVIvdJrwlwFPbFO+ISveeU+SJI1DnSb8S4FPRsSagwVVr/1PApf1IjBJktQ9nfbSnwX8HLgpIn5elb0c2ADvgidJ0rjX0R5+Zt5AOY9/FrBR9fpv4LmZeX3vwpMkSd0wmuvw76IcwpckSRNMxwm/uq3u9sBTaDkykJnf7W5YkiSpmzq98c6rgbOBJ7cZncCabcolSdI40Wkv/ZOAHwJPzcw1Wl4me0mSxrlOD+lPB/bMzDt7GIskSeqRTvfwfwE8o5eBSJKk3ul0D/8U4PiI2Az4PfBYfWRm/qbbgUmSpO7pNOF/u/rb7kE5dtqTJGmc6zThP62nUUiSpJ7qKOFn5q29DkSS1JknZa7wV+rEaG68sxbwImAasHZ9XGae0eW4JElD2H/Z8rEOQRNQpzfe2QY4j3JoPyiPy12L0nnvEcCEL0nSONbpZXlfAa4Gngg8DGwLzAB+C+zTi8AkSVL3dHpI/4XAKzLzoYhYDqyVmb+JiFnAVylP0pMkSeNUpwk/KHv2AAuBzYEbgDuArXoQl6QJatasWQwMDDB16lRmz5491uFIqnSa8K8FngvcDPwa+HhELAPeD9zUo9gkTUADAwPMnz9/rMOQ1KLThP85YP3q/6MpD9L5GbAI2K8HcUmSpC7q9Dr8i2r/3wxsGxEbAfdmeiGoJEnjXcfX4bfKzHu6GYgkSeqdIRN+RPwA2D8zH6j+H1Jm7tn1yCRJUtcMt4f/F8qDcQDuqf0vSZImmCETfma+pzZ4EPBIZi7rfUiSJKnbRrzTXkSsCdwPPKP34UiSpF4YMeFXe/W30vLAHEmSNHF0ei/9zwJfjIgpq7OwiNgpIn4QEfMjIiNi5gj1p1f1Wl+7rU4ckiQ1TaeX5R1OeVLe/Ii4A3ioPjIzO72X/gaUu/adweiesLcb8L+1YS8JlCRpFDpN+N/uxsIy8wLgAoCIOH0Uk/4lMwe6EYMkSU3U6Z32ju11ICP4bkSsA9wInJiZXfkBIklSU3R6Dn+sPEg5nbAf8HrgEuCciNi/XeWIOCAi5kbE3IULF/YxTEmSxreO9vAjYm3gk8DbgGnA4+rjM3PN7ocGmbkI+HKtaG7VcXAWcGab+qcCpwLMmDHDGwVJklQZTS/9d1OS73LgCOBkyt34DupNaEO6Eti6z8uUJGlC6zTh7wccmJn/DiwDvp+ZHwY+Dezaq+CGsD1wV5+XKUnShNZpL/1NgD9U/z8IPKn6/0fAlzpdWERsAGxVDa4BTIuI7YF7MvO2iPgC8KLMfFVV/93AY8A1lCMLewAfAj7e6TIlSVLne/i3AZtV/98EvLb6/6XA4lEsbwYleV8DrAscW/3/mWr8psCWLdMcDcwFrgLeCrw3M08cxTIlSWq8Tvfwvwe8CvgVcBJwdkS8H9gcOK7ThWXmHCCGGT+zZfibwDc7nb8kSWqv0+vwj6z9/+2IuB3YAfhTZp7fq+AkdWaHr+4w1iH8zdr3rc0arMHt990+LuL6xSG/GOsQpHGh08vyplSXyAGQmVdSestLkqQJoNNz+HdGxPkR8ZbqjneSJGkC6TThvwFYRLmpzYKIOD0iXhURQ56PlyRJ40dHCT8zf1x1qNsEOADYkPIQnNsjouNOe5IkaWyM6l76mbkkM8/JzDdSboCzEDisF4FJkqTuGVXCj4j1I2L/iLiQ8nz6xwP/0pPIJElS13TaS3934B3AnpQb7fwP8JnMvKKHsUmSpC7p9MY75wLnUZ6Wd2FmLu1dSJIkqds6vpd+Zv61p5FIkqSe6bSXvslekqQJbFSd9iRJ0sRkwpckqQFM+JIkNcBqJ/yI8EeDJEnj3LDJOiLujIgn14b/OyKeUhveBHish/FJkqQuGGnvfCqwZm14T2CDljo+QEeSpHGuG4fjswvzkCRJPeT5d0mSGmCkhJ+svAfvHr0kSRPMSLfWDeDSiBi8d/66wIUR8WiH00tqmFwvWc5ycj33DaTxZKSEfWzL8Hd6FYikyeGxHbxwRxqPRkr43wDuyMzl/QhGkiT1xkjn8G8BpvQjEEmS1DsjJXyvsZckaRLwsjxJkhqgk172h0fEg8NVyMzPdCkeSZLUA50k/D2ApcOMT8CEL0nSONZJwn9FZt7d80gkSVLPdHKnPUmSNMGtdi/9iHh1l2KRJEk90smd9lbqsBcRmwPvAd4LbMGKj9CVJEnjzLB7+Jl5bGY+DBARa0bE3hHxQ2Ae8CbgFGCrnkcpSZJWy4id9iLiGcD7gHcBDwFnAa8B3pmZf+hteJIkqRuG3cOPiMuBXwEbAvtl5tMz8+hVXVhE7BQRP4iI+RGRETGzg2meHRGXRsTiarpPRYR3AJQkaRRG6rT3UuAM4MTMvLQLy9sAuBb4CLB4pMoR8QTgJ8AC4IXVdEcAh3UhFkmSGmOkhP9CymH/n0fENRHx0YiYuqoLy8wLMvOozPw20MkT+N4BrAe8OzOvrab7EnCYe/mSJHVupE5712Tmh4BNgROAPYHbq+l2j4gNexzfS4HLM7N+NOAiYDNgeo+XLUnSpNHRw3Myc0lmfiszXwlsCxwHfBQYiIgLexjfVMrh/LoFtXGSJKkDo35aXmbelJmfAP4R2A94tOtRraKIOCAi5kbE3IULF451OJIkjRur/HjczFyWmd/PzDd2M6AWA8AmLWWb1Ma1xnRqZs7IzBkbb7xxD8OSJGliWeWE3ydXADtGxDq1sl2BOyk3/5EkSR3oa8KPiA0iYvuI2L5a9rRqeFo1/gsRcUltkrOAh4HTI2K7iNgb+ARwQmb6YB9JkjrU7z38GcA11Wtdyr36rwE+U43fFNhysHJm3k/Zo98MmAucDHyZcsWAJEnq0Ii31u2mzJzDME/gy8yZbcp+D+zUu6gkSZr8xvs5fEmS1AUmfEmSGsCEL0lSA5jwJUlqABO+JEkNYMKXJKkBTPiSJDWACV+SpAYw4UuS1AAmfEmSGsCEL0lSA5jwJUlqABO+JEkNYMKXJKkBTPiSJDWACV+SpAYw4UuS1AAmfEmSGmCtsQ5A6pdZs2YxMDDA1KlTmT179liHI0l9ZcJXYwwMDDB//vyxDkPSBDAZdxBM+JIktZiMOwiew5ckqQFM+JIkNYAJX5KkBjDhS5LUACZ8SZIawIQvSVIDmPAlSWoAE74kSQ3gjXfUc7d95tljHQIAS+/ZCFiLpffcOm5imvap3491CNK48bWPnTfWIfzNfYse+tvf8RLXwV/eY7Wmdw9fkqQGMOFLktQAJnxJkhrAhC9JUgOY8CVJaoC+J/yIOCgibomIJRFxdUTsOEzdnSMi27y26WfMkiRNdH1N+BHxFuAk4PPA84BfAhdGxLQRJn0WsGntdWMv45QkabLp9x7+YcDpmXlaZl6fmYcAdwEfHGG6uzNzoPZa1vtQJUmaPPqW8CNibeAFwI9bRv0YeNkIk8+NiLsi4pKIeGVPApQkaRLr5x7+FGBNYEFL+QJg6hDTDO797wPsDdwAXDLUef+IOCAi5kbE3IULF3YnakmSJoFxfWvdzLyBkuQHXRER04EjgMvb1D8VOBVgxowZ2Y8YJUmTz/prP2GFv5NBPxP+ImAZsElL+SbAwCjmcyXw1m4FpeaYss5yYGn1V5KGtsOWe491CF3Xt4SfmY9GxNXArsC5tVG7At8Zxay2pxzql0bl8OfcN9YhSNKY6fch/ROAb0XEr4FfAAcCmwGnAETEGQCZ+a5q+FBgHnAdsDawP7AX5Zy+JEnqUF8TfmaeExFPBo6mXE9/LfD6zLy1qtJ6Pf7awHHAU4HFlMS/e2Ze0KeQJUmaFPreaS8zvw58fYhxO7cMzwZm9yEsSZImNe+lL0lSA5jwJUlqABO+JEkNYMKXJKkBTPiSJDWACV+SpAYw4UuS1AAmfEmSGsCEL0lSA5jwJUlqABO+JEkNYMKXJKkBTPiSJDWACV+SpAYw4UuS1AAmfEmSGsCEL0lSA5jwJUlqABO+JEkNYMKXJKkBTPiSJDWACV+SpAYw4UuS1AAmfEmSGsCEL0lSA5jwJUlqABO+JEkNYMKXJKkBTPiSJDWACV+SpAYw4UuS1AAmfEmSGsCEL0lSA5jwJUlqgL4n/Ig4KCJuiYglEXF1ROw4Qv1XVPWWRMTNEXFgv2KVJGmy6GvCj4i3ACcBnweeB/wSuDAipg1R/2nABVW95wFfAL4aEfv0J2JJkiaHfu/hHwacnpmnZeb1mXkIcBfwwSHqHwjcmZmHVPVPA74JHN6neCVJmhT6lvAjYm3gBcCPW0b9GHjZEJO9tE39i4AZEfG47kYoSdLk1c89/CnAmsCClvIFwNQhppk6RP21qvlJkqQOrDXWAXRTRBwAHFANPhgRN4xlPKtgCrBorIMYFMe/e6xD6IVx1cZ8OsY6gl4ZN+0cH7aN+yImZTuPqzY+5ISOqm0x1Ih+JvxFwDJgk5byTYCBIaYZGKL+Utp8CJl5KnDq6oU5diJibmbOGOs4JjPbuD9s596zjXtvsrVx3w7pZ+ajwNXAri2jdqX0wm/niiHqz83Mx7oboSRJk1e/e+mfAMyMiPdFxLYRcRKwGXAKQEScERFn1OqfAmweEV+p6r8PmAkc3+e4JUma0Pp6Dj8zz4mIJwNHA5sC1wKvz8xbqyrTWurfEhGvB06kXLp3J/DhzPxOH8Pupwl7OmICsY37w3buPdu49yZVG0dmjnUMkiSpx7yXviRJDWDC77GImBkRD3ZhPnMi4mvdiElaXRGxb0R4eLCHbOPVExE7R0RGRM/v2VItZ98R6jwuIm6IiJ26uNxfj+ZW841I+BGxSUScGBE3Vg/huTsifhkRh0TEBmMdX7dExLxqxXtnm3FXVuN6dlviiDg9Is7v1fzHg4jYOyJ+GhH3RcRDEfH7iPhcRDylGj+zaueMiGVVvbn1OrV5HVOrW3/t1cV4J9S638mGs6r3rIj4dvVArYyIY/oQ3lCxTNY2fn9EXB4R91br8c8i4uVdjGPziDg1Iu6IiEcjYn5EnBYRT12Fec1rs237JaWv2F+6EvDqO4Byq/jLBguG+yxqP1gGXwsj4sKIeG6t2meBL0ZER7l80if8iJgO/AbYDfhn4PnAiykP8HkVsOeYBdcbtwPvrRdExHbAdoyfFb9nqls492renwPOBX4LvAF4JvAR4Gms+DyIhykbmqdS1rWvUNazayNi25bZ3lDVrb8u7FK805m86/56wDxKB+BbxiqISd7GOwPnALtQ3tMNwEURsfXqzjjKg9HmUrZL7wa2AvYHngVcVbXrasnMRzNzIMdBR7WICODDwH+uwuTPomwXdgc2BH4UEU+sxl0APB54XUdzysxJ/aJsPG8H1h9i/GDHxQT2bRk3Dzi8NpyUDfv3KRv1PwGvpGzYLwIeoiSD59emmQk8COxR1V8C/Ax4esuy9qDcp2AJZQP2OWDt2vg5wNdGeK/zgC9W83h6rfxEykOHWt/P/sBVwF+BuynJbPPa+H+m3PzoKbWysykbuLXbLP904Pz6/5SEOB+4F/gGsF41/gDKbZLXbJnHWcAPRtEu84BjgP8C7gPOrco/BdwKPFK9hzPqnzkwC/gzsBj4PbD/CG37ourzP2yI8U+qf95txm9Qff6X1sqOAa4dB+v+YcDvqvV3PvAfg++nVvddVXs+XH2uHwKypc6R1Wf6IHAG8GlgXkud9wB/qD7PPwEfBdaofZZZe83r8H1eCxzTq3a0jVf43gwAh3Sh3S6o2mG9lvL1qvIf1srmUC7RPomyHbkXOK72nua0vKesyneuhqdUwzOrdnsd8MeqnX8APBHYF7gRuB/4FrBubfm7AZdXy72Hsq3ftiXulfJHy/gZwPI2n/mQ07XGX5XtUJW9tlb2X8CZHbX7WHxJ+vUCnlw18ic6qNtpwp8PvA3YmpKcFlQrwBuBf6pW5N/VppkJPEb5NbsD5TG/l1F+GAxuDF4LPFB9Ubek/Ii4ATi+ZaXvJOEfTkncn63K1gYWAju1eT/vBV4PPJ2S0H4GXFYbv2a1og8m8XdVX5Jthlj+6ayY8O8HTgO2BV5DSchHVuM3pGyQdqtNvwFlg7jfKNplXlVnFmUvYWtgn6psd8qlnjOAg2vTfK6az26UvfO3V8vdfZi2PYmysXjcCJ/BTNok/GrcYdU6tHE1fAw9SviMbt0/lLIXNx14BSUxfas2/sXVvD5JWcc/QDlalLU6b60+z/dVdY6sPv95tTrvpzwdc9+q3fegJJCDq/EbV+3zPspzNDbu8L2OScJvUhtX0/4D5Ts87I/jDuazUfVejxpi/Cer8RtWw3MoOyVfBbYB9qve92G1+d0OHFu9p6lV+c6snPAfAy6mPMjtpZRLvS8GzgOeQ9nG3At8rBbPPtVr66rO/wA3seKOx0gJ/6PAjW3KR5vwn1+VvaFWdiBwe0dt3+8vST9f1ZcogTe1lN9B2Xg/CJwyVMPTPuF/oTa8HS17fUOsZAnsUKuzBeU2w6+uhi8D/rll2XtV8Q3+KJhD5wn/dcBtlFM2+w6uaK3vp83021SxPrUl1vuA2ZQk+sFhpj+dFRP+7dT24CnJ/+La8HdZcaO3P+WLvM4o2mUecF5LncMoCX2l5AysT9mr37Gl/CvABcO8twuA/+1gnZvJ0Al/t6p9X1QNH1OtBw/WXtf1e90fIs5H+Pse1FnAT1rq/AcrJqMrWudHedLlvNrwbcA7W+ocCvyh5Ts25IZziHjHKuE3po2r6Y6r3tsTetFutfFvavmezKEcqYhanaOBO2rD82jZtjH0tvgZtTrHV9/BelI9nWo7NkR861fTvLzTNqVsXy5tU95xwqf8wPw+ZTtcP+q6J+UH0lojtf2kP4c/hB2B7YFfA+uMctrf1f4ffJLf79uU1TtoLa+WBUCWGw3dSTkHDOXX5icj4sHBF2UDsD5tniQYEUfV60bEtJYqF1EOv+0K/B/KIZ+VRMTzI+L7EXFrRPyVchQCajdAqmL9CHAEZe//39rNawh/yMxlteE7WbFdzgT2ioj1quF3AN/JzCXVcKftMpcVnUv5XG+JiP+MiDdHxD9U455ZjftRy3w/SDmKQNUxZnDcddV03XgyyOA8slb2Z8q6OPh6fReWM5yV1v2I2CUiflJ1nvor5YfY2vy9jbelJJu61uFtqK3jlSsH/4mIjYF/BP69pd2/SNXu7UTEtJZ1/ahRvNexMunaOCI+QjnqsHdmPjDMe++VX2WV3SpXUO7C+oRRzueRzKw/VG0BMJCZi1rK/radiogtI+KsiPhzRDxQjV+DlhvFjWBdytGZVTGv+hwXUdaTN2fm3bXxiynblhFz2aR6Wl4bN1E2rtvUCzPzFoCIeLhezMob9ce1mWf9Hv45TFnrj6lkaGtQDked22bcwjZlp1AOKw26c4UFZS6PiG8CRwEvoST9FUTE+pQfBhcD76Scw59COYTf2vFtJ8ov2n+MiH/IzEeGeS91rc87SFZslx9SHoT0xoi4BHg15TD+oE7b5aEVFpJ5e0Q8g9Jp6tXAl4FPR8SLa8vfg7I31C7e91G+oPWyPwE7RsTaWZ4LsSqeSXXetFb2aGbetIrzG05H635EbEH5HE6j9Hv4C+Ww4dmsvB6sjsF2P5Chn53Rzp2U5Dnonm4F1AWNaOOIOJTSG/x1mdn6g2NVDLbbM4HvtRk/+D3pxfdiactwMvJ26nzKkY0PUE7pLqX0kRjNZ7eIcjp3VbyS8pksHOLH1kbAkswc8fLvSb2Hn5l/oRzyOriDy2MWUnpCAuVSm/rwalqDco58cN7TKM8QuL4q+g3lvPhNbV6tKyiZec9IdSh79TtSDhPe2Wb8NpQEf1RmXpaZf2TFve/BWPem7HnvQunc8oXRvPHhVD8czq3m/xbKucY5tSqjapeWeS/JzB9m5keBF1J6uu5A+aI+AmzRZp63VtPOby3j70cWDm63vIh40nDxVOvfgZTDeu1+xHXVKNb9GZQN10cz84rM/BNl3ay7nvLDsa51+I+Udq772zqfmQsoiWXLdp9nbZrHKH1HBqdb2lJ33CT8JrRxRBxGSfa7Z+bPh3mPHava7SLgoNrRvcHlrUfprHhhy2f94qqn+6CXUC5xG0yAj9bfU7dEuRX8NsDnM/PizLye0it+tDvL1wDP6PTyuRa3ZOafhzmysh1lWzmiyb6HD3AQ8Avg6ijX6v4v5RfaC4DnUr6wAD8FPhQRv6TszX6eVT8E02op8JXqsNhiSq/56yh71wCfAc6PiFspe+5LKR/iizJz1qosMDNvjnLDicVDVLmNkvgOjoiTKYeKPluvEBGbU/ZKjsrMy6Jc3//TiLggMy9eaY6r5kzgEkoHo7Mzc3lt3Cq1S0TMpKzbV1LOo76FspG7MTP/GhHHA8dXG5DLKJ0FXwIsz/KI5ZVk5pURMRs4Lsp1wt+h/Op/GuUIyk2UoxFVCDF4qPaJlI30x6v/+3mZVifr/o2UH6SHRsR3Ke1waMt8/hX4ZUQcCXybcm7xTS11TgK+ERFXUY4SvYlyrvbeWp1PA1+NiPsofSIeR9nT3TwzB39IzgNeFRGXUg6/1qf/m+ryy8FTYusAUyNie0r/iV7sGQ5lMrfxEZQOrvsDf6qt04sz8/4R2mUkB1OOQlwcEUdT2mjLannByj+sN6NsQ78OPJtyivFfauPnUY7AnVm9p249w/5eyt75+yPidmBzSl+GYXc42vgZZT19DqXDdt30at2tu3kU894R+FFHNUc6yT8ZXpTzZCdRNsqPUJLAVZRero+v6mxGucTmQcp51X1o32lv39rwlKps51rZYMe37fLvHUUepPTiv7Fa/qXAVi0xvobyJX6Y0iljLiv2LJ9Dh532Oh1PSYR/pvyw+TXlcHpSNjYB/ITyo6TeWeazlMNaT24z/9NpuSyvZfwxtPRKr5Yzr1ruc9rMc6R2Wek9Uzr2XUHpbPhQ9Vm/oWWZh/D3vf2F1XvdtYN16c3VZ3F/Ne9rKRupwZ73M/n75UHLq3q/ofyAfMpI7TFG6/6Hq890MeXH135V/NNr83kP5UfiYsr35GBWvmTsKMqpocFLxr4IXN9S521VeyyhbEx/Dry1Nn4PyvfkMYa5ZIzS2z3bvOb0sj0b1sbzhmjj07vUbv9I2aGYX8VyJ6Wj4lNb6s2hnMb8GuU7fS/lNF29Q/BLKD+2lgy2GUNcltcy78Nb26Bq07m14V0o3/Ml1d/XVu0/s1Zn2E57VZ2zgeNaytq1b1Lu87FC/EPMc3PK0Y2nDrfswZcPz5HUExHxPUrP4T3GOpbJqgltHBFzKD+M255Omygi4lmUPf2tsksdHyPiOOCJmXlAJ/WbcEhfUo9V514/SDm0uJRyhOyN1V91gW08sWXmdVFu//s0ytGIbribcmlhR9zDl7TaImJdys1Lnke5wuFG4EuZedaYBjaJNLWNJ8se/nhgwpckqQEm9WV5kiSpMOFLktQAJnxJkhrAhC9JUgOY8CVJagATviRJDfD/AUh8buxFUfeXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfsAAAFHCAYAAACmg+XcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApH0lEQVR4nO3deZgdVbnv8e+PYBgdkMSEwRDFgcmROCIQhyiCIAKCA2L0KAciKOZAjiBKwPEERDgqF4lKBIXLEfUqCIKoYRZIOA4gMigJEOiQyCCBhCTkvX+s1VCp7N1dnd5Dd+X3eZ56uveqVVXvrl273l2rVlUpIjAzM7P6Wq/bAZiZmVl7OdmbmZnVnJO9mZlZzTnZm5mZ1ZyTvZmZWc052ZuZmdWck72ZrbMkzZM0q9txmLWbk73ZWpA0WVJIelLSuAbjL5Y0rwuhDWuStpQ0XdKrux3LQEmalbeJKGwbd0g6SdKG3Y7P1m3rdzsAs2FuJHAccFi3A6mJLYETgHnAHzuwvJcDq1o4vxXAx/P/zwXeC3wB2Bb4cAuXYzYgPrI3G5w/Ah9rdHRvQ19EPBkRK1o4y1UR8aM8fAd4F/AH4IOSxrRwOWYD4mRvNjhfy38/319FSSMkfV7SXbmJ9x5JMyRtVKo3T9KvJb1F0o2Slkn6h6RDqgQkaWJuRv6QpBMkLZC0RNLPJG0maaSkUyUtzOU/bBBDy2OVNF7SBZL+KWmppJsk7VuMG7gpvzy70Bw+vVBnd0lXSnpc0qP5dMlOpeVMz9O9PDetP5Lrni1p4wbxzyqVjZR0vKS/5fe+UNIvJO1YZf0XRbof+TWAgBcPdHqzVnGyNxuce4AfUO3o/rvAl4E/AZ8FrgCOAX4qSaW6LwIuBH4D/AfwMDBrgAlnGvBu4OvA94F9gbOAmcAOwInAT4FDgGPbGaukFwDXAXsCZ+blCfiZpA/marcBX8z/nwV8JA8/y/N4a45jC2A6cAowAbhW0ssavP//Czw7L+t/gMmkUwRNSVoPuAj4EvBnYCrwDdLpmp37mrYP4/Pfh9dyerPBiwgPHjwMcCAljgDeCLwQeBL4bmH8xcC8wutX5vpnl+YzPZe/p1A2L5ftVigbDSwDTqkQ28Q8/V+BkYXy80jnpy8DVCi/DrivnbECp+Z6EwtlG+UYHwCelcsm5HqTG7yvm4HFwOaFspcCy4ELG8T5g9L0PwMWl8rmAbMafK7HNFi+ymWl8bPy+x6Vh21JP35WAX/pb3oPHto5+MjebJAi4l6eObrfpkm1vfLfU0vl3wSeKozvdUdEXFVYxiLgdgbWFHxuRCwvvL6BdDR9dkREqXxLSRu0Mda9gJsjYnah3lLgDGAs8Nq+3oikLYDXAD+MiH8W5nEn8EtgD0kjSpPNLL2+Gthc0nP6WNQBpCPw08sjSuusmQ2ARXm4i9T6cDXw3orTm7WFk71Za3yVdETY7Nz9Nnn8HcXCiHiUdGQ7vlT/ngbzeBjYrPeFpLGlYaNS/fI8Hs1/721QLuB57Yo1z/P2BvVuy3/L8yzr/RHVbB6bkI6m+4qrtxl9M5rblvTjZXkfdfqyApiUh48BtwIvAB5fy/mZtYSTvVkL5KP77wOT+zi6H4inmpQXz5c/UBoOqjiPKvMeiFbPr1W6EdeqiLgiD7OAt5NaLr7bxmWa9cvJ3qx1eo/uj28wbj4pyazWkSw3KW9BOnc8UJNKw2VrMY9G2hHrfNI17WXb5b+982zW1D0//202j8dJ5/MH6+/AyySNbMG8iIiFpA5+75X0xlbM02xtONmbtUhE3Ad8D/gozzQ79/pV/ntUqfwzwAhSh76BLu+K0vDAQOfRRMtjzdO8VtKuvQX5rnKHAz3A3Fzc29y9WlN7fm83A4dIen5hHtsC+wCXRkSzI/mBuDAv+zPlEQ2uQqjqO6T3Vb7iwaxjfAc9s9b6GvBvwE48czRKRPxZ0veBf5P0XOB3pE5pHwd+DVzShVgbalOs/wV8EPiVpP8mdWA7mHQJ4IcjYmWu93fSufXDJS0BHgNuiYhbgKOBy4HrJc0ENgQ+ReoB3+99Dio6N8c1Q9IE4Kq8nLcCF+TxAxIRD0v6AXCEpB0j4tYWxWpWmY/szVqocHTfyL+TmvhfA5xGurvaKcB+Q7CndktjjYgHgV1IPxamkJK/gP0j4vxCvRWka+uXkY6Izyf1kCcifk86XbEQOIl0H4GbgbdExGqdCddWbh14D+k6+51JVyAcA6zkmdaHtfFN0iV4nxtsjGZrQ0NvH2NmZmat5CN7MzOzmnOyNzMzqzknezMzs5pzsjczM6s5J3szM7Oaq+119qNGjYrx48d3OwwzM7OOmTt37uKIGF0ur22yHz9+PHPmzOl2GGZmZh0jaX6jcjfjm5mZ1ZyTvZmZWc052ZuZmdWck72ZmVnNOdmbmZnVnJO9mZlZzTnZm5mZ1ZyTvZmZWc052ZuZmdWck72ZmVnN1fZ2uWZl06ZNo6enh7FjxzJjxoxuh2Nm1jFO9rbO6OnpYcGCBd0Ow8ys49yMb2ZmVnNO9mZmZjXnZG9mZlZzTvZmZmY152RvZmZWc072ZmZmNdfRZC9pN0m/lLRAUkiaXGGaV0i6UtLSPN0XJakD4ZqZmdVCp4/sNwVuAT4DLO2vsqTnAL8BFgKvy9MdA0xtY4xmZma10tGb6kTEJcAlAJJmVZjkw8DGwEcjYilwi6TtgKmSTo2IaFuwZmZmNTHUz9m/Cbg6J/pelwFbAuO7EpGZmdkwM9ST/VhSE37RwsI4MzMz68dQT/YDIulQSXMkzVm0aFG3wzEzMxsShnqy7wHGlMrGFMatJiLOiogJETFh9OjRbQ/OzMxsOBjqyf56YFdJGxbKJgH3A/O6EpGZmdkw0+nr7DeV9GpJr87LHpdfj8vjvybpt4VJzgOeAGZJ2knSfsDnAPfENzMzq6jTR/YTgP/Nw0bAifn/k/L4LYBteytHxKOkI/ktgTnAd4BvAKd2LmQzM7PhrdPX2c8Gmt79LiImNyj7C7Bb+6IyMzOrt6F+zt7MzMwGycnezMys5pzszczMas7J3szMrOac7M3MzGrOyd7MzKzmnOzNzMxqzsnezMys5pzszczMas7J3szMrOac7M3MzGrOyd7MzKzmnOzNzMxqzsnezMys5pzszczMaq6jz7O3ddM9J72i2yEAsPKh5wPrs/Kh+UMmpnFf/Eu3QzCzdYCP7M3MzGrOyd7MzKzmnOzNzMxqzsnezMys5pzszczMas7J3szMrOac7M3MzGrOyd7MzKzmnOzNzMxqzsnezMys5pzszczMas7J3szMrOac7M3MzGrOyd7MzKzmnOzNzMxqzsnezMys5pzszczMas7J3szMrOac7M3MzGrOyd7MzKzmnOzNzMxqzsnezMys5pzszczMaq7jyV7SFEl3S1omaa6kXfup/yFJf5T0hKQeST+SNLZT8ZqZmQ13HU32kg4CTge+CrwGuA64VNK4JvV3Ac4FfgjsCOwL7AD8uBPxmpmZ1UGnj+ynArMiYmZE3BYRRwIPAIc3qf8m4L6I+GZE3B0RfwC+BbyhQ/GamZkNex1L9pJGAjsDl5dGXQ68uclk1wJbSNpbySjgA8Al7YvUzMysXjp5ZD8KGAEsLJUvBBqeg4+I60nJ/cfAcmARIOCjjepLOlTSHElzFi1a1Kq4zczMhrUh3Rtf0g6kZvsvkVoF9iD9MPhuo/oRcVZETIiICaNHj+5coGZmZkPY+h1c1mLgKWBMqXwM0NNkmmOBGyPi5Pz6z5IeB66WdFxE3NeeUM3MzOqjY0f2EbEcmAtMKo2aROqV38jGpB8IRb2vh3SrhA09ozZcxZiNVjJqw1XdDsXMrKM6eWQPcCpwrqQbSZ3vDgO2BM4EkHQOQEQckutfBMyUdDhwGbAFcBpwc0Tc09nQbbg7+pWPdDsEM7Ou6Giyj4gLJG0OHE9K3LcAe0bE/FxlXKn+LEnPBo4AvgE8CvwO+M/ORW1mZja8dfrInog4AzijybiJDcq+ReqkZ2ZmZmvB573NzMxqzsnezMys5pzszczMas7J3szMrOac7M3MzGrOyd7MzKzmnOzNzMxqzsnezMys5pzszczMas7J3szMrOac7M3MzGqucrKXNEXSrZKekPTiXPY5SQe2LzwzMzMbrErJXtJRpCfVnQWoMGoB6Yl0ZmZmNkRVPbI/DPhkRJwOrCyU3wzs2PKozMzMrGWqJvttSM+eL1sBbNS6cMzMzKzVqib7fwCvbVC+J/DX1oVjZmZmrbZ+xXqnAN+WtDHpnP2bJH0EmAZ8vF3BmZmZ2eBVSvYRcbak9YGvAhsD5wL3A5+OiAvaGJ+ZmZkNUtUjeyJiJjBT0ihgvYh4sH1hmZmZWatUSvaSdgRGRMSfI2JxofyVwMqI8Hl7MzOzIapqB72zgJ0alO+Qx5mZmdkQVTXZvxK4sUH5TcArWheOmZmZtVrVZP8U8NwG5Zux+h31zMzMbIipmuyvBD4vaURvQe6d/3ngqnYEZmZmZq1RtTf+NOAa4C5J1+SytwCbAru1IzAzMzNrjUpH9hFxO+m8/XnA8/PwY+BVEXFb+8IzMzOzwRrIdfYPkJrtzczMbBipnOzzrXJfDbyAUotARPystWGZmZlZq1S9qc47gPOBzRuMDmBEg3IzMzMbAqr2xj8d+BWwdUSsVxqc6M3MzIawqs3444F9IuL+NsZiZmZmbVD1yP5a4OXtDMTMzMzao+qR/ZnAKZK2BP4CrCiOjIibWx2YmZmZtUbVZH9h/tvooTfuoGdmZjaEVU32L2prFGZmZtY2lZJ9RMxvdyBmZmbWHgO5qc76wOuBccDI4riIOKfFcZmZmVmLVL2pznbARaTmfJEeebs+qaPek4CTvZmZ2RBV9dK704C5pGfaPwFsD0wA/gjsP5AFSpoi6W5JyyTNlbRrP/VHSjopT/OkpHskfXogyzQzM1uXVW3Gfx2we0Q8LmkVsH5E3CxpGvAt0hPx+iXpINLd+KaQHpk7BbhU0g4RcU+Tyf4vsDVwKHAnMAbYqGLcZmZm67yqyV6kI3qARcBWwO3AfcBLBrC8qcCsiJiZXx8paQ/gcODYNRYqvRN4O7BtRCzOxfMGsDwzM7N1XtVm/FuAV+X/bwT+U9LuwInAXVVmIGkksDNweWnU5cCbm0y2L3ATMFXSfZLulPTfkjatGLeZmdk6r+qR/VeATfL/x5MeivN7YDFwYMV5jCLdfGdhqXwh8I4m07wYeAupE+D+wPNIpw22BA6ouFwzM7N1WtXr7C8r/P8PYHtJzwcejohoV3CklocAPhQRjwJIOgK4TNKYiFjth4OkQ0nn9hk3blwbwzIzMxs+qjbjryEiHhpgol9MumRvTKl8DNDTZJoHgAW9iT67Lf9dI5tHxFkRMSEiJowePXoAoZmZmdVX0yN7Sb8EDo6If+X/m4qIffpbUEQslzQXmAT8pDBqEvDTJpNdC7xf0qYRsSSXvSz/9V39zMzMKuirGf+fpCZ0gIcK/w/GqcC5km4kJfLDSOffzwSQdA5ARByS658HfAE4W9J00jn704ELI+LBFsRjZmZWe02TfUR8rPByCvBkRDw1mIVFxAWSNid18tuC1Mt/z8K998eV6i+R9A5Sp7ybgIeB/wd8bjBxDEXTpk2jp6eHsWPHMmPGjG6HY2ZmNdJvBz1JI4BHSZfe/XWwC4yIM4Azmoyb2KDsduCdg13uUNfT08OCBQu6HYaZmdVQvx308tH8fEoPvzEzM7PhoWpv/C8BX5c0qp3BmJmZWetVvanO0aQn3i2QdB/weHFkRFS6N76ZmZl1XtVkf2FbozAzM7O2qXoHvRPbHYiZmZm1x1rfQc/MzMyGh0rJXtJISSdKukPSMklPFYd2B2lmZmZrbyC98T8KfANYBRwDfId0l70p7QnNzMzMWqFqsj8QOCwivkt6mM0vIuLTwAmke9ubmZnZEFU12Y/hmbvnLSHdox7g16wDd7czMzMbzqom+3tID6wBuAt4V/7/TcDSVgdlZmZmrVM12f8ceHv+/3TgREl3A7OA77UhLjMzM2uRqtfZH1v4/0JJ9wK7AHdExMXtCs7MzMwGr1KylzQqIhb3vo6IG4Ab2haVmZmZtUzVZvz7JV0s6SBJG7Y1IjMzM2upqsn+PcBi4CxgoaRZkt4uSe0LzczMzFqhUrKPiMsjYjLpErxDgc2AS4B7JZ3cvvDMzMxssKo+9Q6AiFgGXABcIGl74DxgKumOesPOzsec0+0QnvbsxY8xArhn8WNDJq65Jx/S7RDMzKwFBvQgHEmbSDpY0qXAn4BnA19uS2RmZmbWElV74+8FfBjYh3QTnf8BToqI69sYm5mZmbVA1Wb8nwAXAR8ELo2Ile0LyczMzFqparIfExGPtTUSMzMza4uqvfGd6M3MzIapAXXQMzMzs+HHyd7MzKzmnOzNzMxqbtDJXpJ/MJiZmQ1hfSZqSfdL2rzw+seSXlB4PQZY0cb4zMzMbJD6OyofC4wovN4H2LRUxw/DMTMzG8Ja0QQfLZiHmZmZtYnPt5uZmdVcf8k+WPPI3UfyZmZmw0h/t8sVcKWk3nvhbwRcKml5xenNzMysy/pL1ieWXv+0XYGYmZlZe/SX7M8G7ouIVZ0IxszMzFqvv3P2dwOjOhHIum7VyE14aoPnsGrkJt0OxczMaqbKOXvrgMdf+s5uh2BmZjXlS+/MzMxqrkpv+qMlLemrQkSc1KJ4zMzMrMWqJPu9gZV9jA+gcrKXNAU4BtgCuBU4KiKurjDdW4DZwN8iYqeqyzMzM1vXVUn2u0fEg61YmKSDgNOBKcA1+e+lknaIiHv6mG4z4Bzgt8BWrYjFzMxsXVHlDnqtNBWYFREzI+K2iDgSeAA4vJ/pvg/8ELi+xfGYmZnVXn/Jvt/e+JLeUWVBkkYCOwOXl0ZdDry5j+mmAGOAL1dZjpmZma2uv2R/IrBG5zxJW0k6XtI/gMsqLmsU6XG5C0vlC0mP0l2DpFcAJwAHR8RT/S1A0qGS5kias2jRoophmZmZ1VufyT4iToyIJwAkjZC0n6RfAfOA9wFnAi9pR2CSNgAuAI6OiLurTBMRZ0XEhIiYMHr06HaEZWZmNuz020FP0suBTwCHAI8D5wHvBD4SEX8dwLIWA0+RmuSLxgA9DepvAWwPnC3p7Fy2XgpJK4E9I6J8SsDMzMxK+jyyl3Q18AdgM+DAiHhxRBy/NguKiOXAXGBSadQk4LoGkywAXgG8ujCcCdyV/280jZmZmZX0d2T/JuA7wFkRcWsLlncqcK6kG4FrgcOALUlJHEnnAETEIRGxArilOLGkB4EnI2K1cjMzM2uuv2T/OlIT/jWS5pGudT9/bRcWERdI2hw4ntRMfwupOX5+rjJubedtZmZmjfXXQe9/I+JTpMR8KrAPcG+ebq98s5sBiYgzImJ8RGwQETtHxFWFcRMjYmIf00733fPMzMwGptKDcCJiWUScGxFvJXWaOxn4LNAj6dJ2BmhmZmaDM+Cn3kXEXRHxOeCFwIHA8pZHZWZmZi1T5d74DeWb3PwiD2ZmZjZE+Xn2ZmZmNedkb2ZmVnNO9mZmZjXnZG9mZlZzTvZmZmY152RvZmZWc072ZmZmNedkb2ZmVnNO9mZmZjXnZG9mZlZzTvZmZmY152RvZmZWc2v9IBwzs0amTZtGT08PY8eOZcaMGd0Ox8xwsjezFuvp6WHBggXdDsPMCtyMb2ZmVnNO9mZmZjXnZG9mZlZzTvZmZmY152RvZmZWc072ZmZmNedkb2ZmVnNO9mZmZjXnm+qY1cAu39ql2yE8beQjI1mP9bj3kXuHRFzXHnltt0Mw6zof2ZuZmdWck72ZmVnNuRnfzGyY8cOGbKCc7M3Mhhk/bMgGys34ZmZmNedkb2ZmVnNO9mZmZjXnZG9mZlZzTvZmZmY152RvZmZWc770zsxaKjYOVrGK2Di6HYqZZU72ZtZSK3ZZ0e0QzKyk4834kqZIulvSMklzJe3aR939JF0uaZGkxyTdIGmfTsZrZmY23HX0yF7SQcDpwBTgmvz3Ukk7RMQ9DSbZHfgdcDzwEPBh4OeSJkbE1R0K28yMK3fbvdshPG3p+iNAYul99w2ZuHa/6spuh2B96HQz/lRgVkTMzK+PlLQHcDhwbLlyRHymVHSipL2AfQEnezMzswo61owvaSSwM3B5adTlwJsHMKtnAw+3Ki4zM7O66+Q5+1HACGBhqXwhMLbKDCR9CtgaOLfJ+EMlzZE0Z9GiRYOJ1czMrDaGzXX2kvYHTgY+FBHzG9WJiLMiYkJETBg9enRnAzQzMxuiOpnsFwNPAWNK5WOAnr4mlHQA6Wj+kIi4qD3hmZmZ1VPHkn1ELAfmApNKoyYB1zWbTtKBpEQ/OSIubF+EZmZm9dTp3vinAudKuhG4FjgM2BI4E0DSOQARcUh+/QFSoj8auEpS77n95RHxUIdjNzMzG5Y6muwj4gJJm5Oum98CuAXYs3AOflxpksNIMZ6Wh15XAhPbGauZmVlddPx2uRFxBnBGk3ET+3ptZmZmA+d745uZDTPPi1jtr1l/nOzNzIaZg59a1e0QbJgZNtfZm5mZ2dpxsjczM6s5J3szM7Oac7I3MzOrOSd7MzOzmnNvfDMzs5Jp06bR09PD2LFjmTFjRrfDGTQnezMzs5Kenh4WLFjQ7TBaxs34ZmZmNedkb2ZmVnNO9mZmZjXnZG9mZlZzTvZmZmY15974ZmY2JHz7Py7qdghPe2Tx40//HSpxHfGNvdd6Wh/Zm5mZ1ZyTvZmZWc052ZuZmdWck72ZmVnNOdmbmZnVnJO9mZlZzTnZm5mZ1ZyvszczMyvZZORzVvs73DnZm5mZleyy7X7dDqGl3IxvZmZWc072ZmZmNedkb2ZmVnNO9mZmZjXnZG9mZlZzTvZmZmY152RvZmZWc072ZmZmNedkb2ZmVnNO9mZmZjXnZG9mZlZzTvZmZmY152RvZmZWc072ZmZmNdfxZC9piqS7JS2TNFfSrv3U3z3XWybpH5IO61SsZmZmddDRZC/pIOB04KvAa4DrgEsljWtS/0XAJbnea4CvAd+StH9nIjYzMxv+On1kPxWYFREzI+K2iDgSeAA4vEn9w4D7I+LIXH8m8EPg6A7Fa2ZmNux1LNlLGgnsDFxeGnU58OYmk72pQf3LgAmSntXaCM3MzOqpk0f2o4ARwMJS+UJgbJNpxjapv36en5mZmfVj/W4H0EqSDgUOzS+XSLq9m/GshVHA4m4H0UunfLTbIbTDkFrHnKBuR9AuQ2Y969Nexx2hWq7nIbWOjzy1UrVtGhV2MtkvBp4CxpTKxwA9TabpaVJ/JQ0+gIg4CzhrcGF2j6Q5ETGh23HUmddxZ3g9t5/XcfvVaR13rBk/IpYDc4FJpVGTSL3tG7m+Sf05EbGitRGamZnVU6d7458KTJb0CUnbSzod2BI4E0DSOZLOKdQ/E9hK0mm5/ieAycApHY7bzMxs2OroOfuIuEDS5sDxwBbALcCeETE/VxlXqn+3pD2Bb5Iuz7sf+HRE/LSDYXfSsD0FMYx4HXeG13P7eR23X23WsSKi2zGYmZlZG/ne+GZmZjXnZN9mkiZLWtKC+cyW9O1WxGQ2WJIOkORmwTbyOl57kiZKCkltvx9LXs4B/dR5lqTbJe3WwuXeOJBbx68TyV7SGEnflHRnfqDOg5Kuk3SkpE27HV+rSJqXN7yPNBh3Qx7XtlsNS5ol6eJ2zX8okLSfpN9JekTS45L+Iukrkl6Qx0/O6zkkPZXrzSnWKcxreqFucdi3hfEOq22/yo4z19tR0oX54VghaXoHwmsWS13X8SclXS3p4bwd/17SW1oUw1aSzpJ0n6TlkhZImilp67WY17wG+7XrSP3C/tmKeFvgUNKt36/qLejrcyj8WOkdFkm6VNKrCtW+BHxdUqU8XvtkL2k8cDOwB/AF4LXAG0gP43k7sE/XgmuPe4GPFwsk7QTsxNDZ8Nsm35a5XfP+CvAT4I/Ae4AdgM8AL2L15zs8QdrRbE3a1k4jbWe3SNq+NNvbc93icGmL4h1Pfbf9jYF5pM6+d3criJqv44nABcDbSO/pduAySS8dzEyVHnA2h7RP+ijwEuBgYEfgprxOByUilkdETwyBTmmSBHwa+P5aTL4jaZ+wF7AZ8GtJz83jLgGeDby70pwiotYDacd5L7BJk/G9nRQDOKA0bh5wdOF1kHbqvyDt0O8A3kraqV8GPE5KBK8tTDMZWALsnesvA34PvLi0rL1J9yFYRtp5fQUYWRg/G/h2P+91HvD1PI8XF8q/SXqAUPn9HAzcBDwGPEhKZFsVxn+BdGOjFxTKzift3EY2WP4s4OLi/6RkuAB4GDgb2DiPP5R06+MRpXmcB/xyAOtlHjAd+AHwCPCTXP5FYD7wZH4P5xQ/c2Aa8HdgKfAX4OB+1u3r8+c/tcn45xU/7wbjN82f/5WFsunALUNg258K/DlvvwuA7/W+n0LdQ/L6fCJ/rp8ColTn2PyZLgHOAU4A5pXqfAz4a/487wA+C6xX+CyjMMyr+D5vAaa3az16Ha/2vekBjhzkOrskr4ONS+Ub5/JfFcpmky7BPp20D3kYOLnwfmaX3k/k8on59aj8enJeZ+8G/pbX8S+B5wIHAHcCjwLnAhsVlr8HcHVe7kOk/fz2pbjXyB2l8ROAVQ0+76bTlePPZbvksncVyn4A/KjSeu/GF6RTA7B5Xsmfq1C3arJfAHwQeCkpMS3MG8B7gZflDfnPhWkmAytIv2R3IT2q9yrSj4LeHcG7gH/lL+m2pB8QtwOnlDb6Ksn+aFLS/lIuGwksAnZr8H4+DuwJvJiUzH4PXFUYPyJv6L0J/JD8JdmuyfJnsXqyfxSYCWwPvJOUjI/N4zcj7Yz2KEy/KWlneOAA1su8XGca6QjhpcD+uWwv0uWcE4AjCtN8Jc9nD9JR+YfycvfqY92eTtpZPKufz2AyDZJ9Hjc1b0Oj8+vptCnZM7Bt/yjS0dt4YHdSUjq3MP4NeV6fJ23j/05qJYpCnQ/kz/MTuc6x+fOfV6jzSdJTLg/I631vUvI4Io8fndfPJ0jPxRhd8b12JdmvS+s4T7sB6Tvc5w/jfubx/Pw+j2sy/vN5/Gb59WzSwci3gO2AA/N7nlqY373Aifn9jM3lE1kz2a8AriA9kO1NpEu5rwAuAl5J2r88DPxHIZ798/DSXOd/gLtY/YCjv2T/WeDOBuUDTfavzWXvKZQdBtxbad13+gvSySF/gQJ4X6n8PtKOewlwZrMVT+Nk/7XC650oHe012cgC2KVQZxvSrYPfkV9fBXyhtOx9c3y9PwhmUz3Zvxu4h3Sa5oDeDa38fhpMv12OdetSrI8AM0gJ9PA+pp/F6sn+XgpH7qTEf0Xh9c9YfYd3MOmLvOEA1ss84KJSnamkZL5GYgY2IR3N71oqPw24pI/3dgnwpwrb3GSaJ/s98vp9fX49PW8HSwrDrZ3e9pvE+STPHD2dB/ymVOd7rJ6Iri/Pj/TEynmF1/cAHynVOQr4a+k71nTH2STebiX7dWYd5+lOzu/tOa1eZ4Xx7yt9R2aTWidUqHM8cF/h9TxK+zWa74dfXqhzSv7+FRPqLPI+rEl8m+Rp3lJ1fZL2LVc2KK+c7Ek/LH9B2gcXW1r3If04Wr+/dV/7c/ZN7Aq8GrgR2HCA0/658H/vE/n+0qCs2BlrVV4WAJFuInQ/6ZwvpF+an5e0pHcgffk3ocETASUdV6wraVypymWkJrdJwL+RmnrWIOm1kn4hab6kx0itD1C4uVGO9TPAMaSj/v/TaF5N/DUiniq8vp/V18uPgH0lbZxffxj4aUQsy6+rrpc5rO4npM/1bknfl/R+SRvkcTvkcb8uzfdwUusBuSNM77hb83SteMpH7zyiUPZ30rbYO+zZguX0ZY1tX9LbJP0md5Z6jPQjbCTPrOPtSYmmqPx6OwrbeHZD7z+SRgMvBL5bWu9fJ6/3RiSNK23rxw3gvXZL7daxpM+QWhv2i4h/9fHe2+EPkTNbdj3pzqrPGeB8noyI4sPRFgI9EbG4VPb0PkrStpLOk/R3Sf/K49ejdAO4fmxEapFZG/PyZ7iYtI28PyIeLIxfStqv9JvHavXUuwbuIu1YtysWRsTdAJKeKBaz5g79WQ3mWbwnf/RRVv4hFTS3HqkZ6icNxi1qUHYmqTmp1/2rLShilaQfAscBbyQl/NVI2oT0o+AK4COkc/ajSM325U5uu5F+zb5Q0gYR8WQf76Wo/PyCYPX18ivSQ43eK+m3wDtITfe9qq6Xx1dbSMS9kl5O6iD1DuAbwAmS3lBY/t6ko6BG8X6C9AUtlt0B7CppZKTnPKyNHcjnSQtlyyPirrWcX18qbfuStiF9DjNJ/Rz+SWouPJ81t4PB6F3vh9H8WRiN3E9KnL0ealVALbBOrGNJR5F6fr87Iso/Ngaqd53tAPy8wfje70g7vhMrS6+D/vdRF5NaM/6ddAp3Jak/xEA+t8Wk07dr462kz2NRkx9ZzweWRUS/l3fX+sg+Iv5JauY6osIlMItIvR6BdDlN8fUgrUc6J94773GkZwLclotuJp0Hv6vBUN5AiYiH+qtDOprfldQ0eH+D8duRkvtxEXFVRPyN1Y+6e2Pdj3TE/TZSZ5avDeSN9yX/aPhJnv9BpHOLswtVBrReSvNeFhG/iojPAq8j9WrdhfRFfRLYpsE85+dpF5TLeKZF4YhGy5P0vL7iydvfYaTmvEY/4FpqANv+BNKO67MRcX1E3EHaNotuI/1oLCq//htpPRc9vc1HxEJSUtm20edZmGYFqa9I73QrS3WHTLJfF9axpKmkRL9XRFzTx3usJK+zy4AphRa93mVtTOqUeGnpc35D7tHe642ky9h6k9/y4vtpFaVbu28HfDUiroiI20i93wd6kPy/wMurXiJXcndE/L2P1pSdSPvJftX9yB5gCnAtMFfpWtw/kX6d7Qy8ivRlBfgd8ClJ15GOYr/K2je9lK0ETstNYUtJveNvJR1VA5wEXCxpPumIfSXpQ3x9RExbmwVGxD+UbiixtEmVe0hJ7whJ3yE1EX2pWEHSVqSjkeMi4iql6/d/J+mSiLhijTmunR8BvyV1Jjo/IlYVxq3VepE0mbRt30A6b3oQaQd3Z0Q8JukU4JS8A7mK1DHwjcCqSI9JXkNE3CBpBnCy0rXAPyX94n8RqeXkLlIrRA5Bvc2zzyXtoP8z/9/JS7GqbPt3kn6MHiXpZ6T1cFRpPv8NXCfpWOBC0vnE95XqnA6cLekmUuvQ+0jnZx8u1DkB+JakR0h9IJ5FOsLdKiJ6f0TOA94u6UpSs2tx+qflSyx7T4NtCIyV9GpSf4l2HBU2U+d1fAypM+vBwB2FbXppRDzaz3rpyxGklocrJB1PWj/b5mWJNX9Qb0naf54BvIJ0SvHLhfHzSK1uP8rvp1XPn3+YdFT+SUn3AluR+i30eaDRwO9J2+grSR2zi8bn7bboHwOY967AryvV7O+kfh0G0nmx00k75CdJCeAmUm/WZ+c6W5Iuo1lCOo+6P4076B1QeD0ql00slPV2ctspnukYsoTUW//OvPwrgZeUYnwn6Qv8BKkTxhxW70E+m4od9KqOJyXBv5N+1NxIakIP0o5GwG9IP0iKnWO+RGrO2rzB/GdRuvSuNH46pd7neTnz8nJf2WCe/a2XNd4zqRPf9aSOhY/nz/o9pWUeyTNH+Yvye51UYVt6f/4sHs3zvoW0k+rtYT+ZZy4DWpXr3Uz68fiC/tZHl7b9T+fPdCnph9eBOf7xhfl8jPQDcSnpe3IEa14WdhzpdFDvZWFfB24r1flgXh/LSDvTa4APFMbvTfqerKCPy8JIvdqjwTC7netzHVvH85qs41ktWGcvJB1ILMhx3E/qkLh1qd5s0mnLb5O+zw+TTssVO/6+kfQja1nv+qLJpXeleR9dfv95fc4pvH4b6Tu+LP99V173kwt1+uygl+ucD5xcKmu0boN0D4/V4m8yz61IrRpb97Xs3sEPwjGztpD0c1Iv4b27HUtd1X0dS5pN+kHc8PTZcCFpR9IR/kuiRR0cJZ0MPDciDq1Sf11oxjezNsvnWw8nNSmuJLWMvTf/tRbwOh6+IuJWpVv6vojUCtEKD5IuH6zER/ZmNmiSNiLdnOQ1pCsZ7gT+KyLO62pgNbIuruO6HNkPBU72ZmZmNVfrS+/MzMzMyd7MzKz2nOzNzMxqzsnezMys5pzszczMas7J3szMrOb+Pwh1BGWaJyI1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_ate_var(var_mono, title='Monotonic R')\n",
    "plot_ate_var(var_nonmono, title='Non-monotonic R')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "monotone:\n",
      "indep: 2.733486489809488 0.3913439631449939\n",
      "gmx: 2.4667295565595384 0.35167200462405224\n",
      "icdf: 0.6081048167483316 0.10889778343050012\n",
      "opt: 0.6081048167483316 0.10889778343050012\n",
      "gadget-1: 1.4268250666989728 0.22174222145468936\n",
      "gadget-2: 2.0025127 0.3010683485305203\n"
     ]
    }
   ],
   "source": [
    "print(\"monotone:\")\n",
    "sqrt_length = np.sqrt(len(var_mono['gmx']))\n",
    "print(\"indep:\", np.mean(var_mono['indep']), np.std(var_mono['indep'])/sqrt_length)\n",
    "print(\"gmx:\", np.mean(var_mono['gmx']), np.std(var_mono['gmx'])/sqrt_length)\n",
    "print(\"icdf:\", np.mean(var_mono['icdf']), np.std(var_mono['icdf'])/sqrt_length)\n",
    "print(\"opt:\", np.mean(var_mono['opt']), np.std(var_mono['opt'])/sqrt_length)\n",
    "print(\"gadget-1:\", np.mean(var_mono['jp']), np.std(var_mono['jp'])/sqrt_length)\n",
    "print(\"gadget-2:\", np.mean(var_mono['jpz']), np.std(var_mono['jpz'])/sqrt_length)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "non-monotone:\n",
      "indep: 0.834538600419416 0.10410340659394125\n",
      "gmx: 0.5068929387904362 0.06618207663718108\n",
      "icdf: 0.9119719070902205 0.12953292939956165\n",
      "opt: 0.11582219665192553 0.02799775553045251\n",
      "gadget-1: 0.2689475769701155 0.03697646188685072\n",
      "gadget-2: 0.21554239 0.037788736229982856\n"
     ]
    }
   ],
   "source": [
    "print(\"non-monotone:\")\n",
    "sqrt_length = np.sqrt(len(var_nonmono['gmx']))\n",
    "print(\"indep:\", np.mean(var_nonmono['indep']), np.std(var_nonmono['indep'])/sqrt_length)\n",
    "print(\"gmx:\", np.mean(var_nonmono['gmx']), np.std(var_nonmono['gmx'])/sqrt_length)\n",
    "print(\"icdf:\", np.mean(var_nonmono['icdf']), np.std(var_nonmono['icdf'])/sqrt_length)\n",
    "print(\"opt:\", np.mean(var_nonmono['opt']), np.std(var_nonmono['opt'])/sqrt_length)\n",
    "print(\"gadget-1:\", np.mean(var_nonmono['jp']), np.std(var_nonmono['jp'])/sqrt_length)\n",
    "print(\"gadget-2:\", np.mean(var_nonmono['jpz']), np.std(var_nonmono['jpz'])/sqrt_length)\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Joint Predictor.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}