izmailovpavel/mode_connectivity_example.ipynb

## mode_connectivity_example.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "env: CUDA_VISIBLE_DEVICES=0\n"
     ]
    }
   ],
   "source": [
    "%env CUDA_VISIBLE_DEVICES=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import torchvision\n",
    "import torch.nn.functional as F\n",
    "import tqdm\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LeNet(torch.nn.Module):\n",
    "    \n",
    "    def __init__(self):\n",
    "        super(LeNet, self).__init__()\n",
    "        self.conv1 = torch.nn.Conv2d(1, 32, 5, padding=2)\n",
    "        self.conv2 = torch.nn.Conv2d(32, 64, 5)\n",
    "        self.fc1   = torch.nn.Linear(64*5*5, 100)\n",
    "        self.fc2   = torch.nn.Linear(100, 10)\n",
    "        self.num_flat_features = 64*5*5\n",
    "    \n",
    "    def forward(self, x):\n",
    "        x = F.max_pool2d(F.relu(self.conv1(x)), (2, 2))\n",
    "        x = F.max_pool2d(F.relu(self.conv2(x)), (2, 2))\n",
    "        x = x.view(-1, self.num_flat_features)\n",
    "        x = F.relu(self.fc1(x))\n",
    "        x = self.fc2(x)\n",
    "        return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform = torchvision.transforms.ToTensor()\n",
    "\n",
    "trainset = torchvision.datasets.MNIST(root=\"/datasets/\", train=True, transform=transform)\n",
    "testset = torchvision.datasets.MNIST(root=\"/datasets/\", train=False, transform=transform)\n",
    "\n",
    "trainloader = torch.utils.data.DataLoader(trainset, batch_size=32, shuffle=True)\n",
    "testloader = torch.utils.data.DataLoader(testset, batch_size=32, shuffle=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = LeNet().cuda()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "criterion = torch.nn.CrossEntropyLoss()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_epochs = 10\n",
    "optimizer = torch.optim.SGD(model.parameters(), lr=0.01, momentum=0.9)\n",
    "scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=n_epochs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate(model, loader):\n",
    "    criterion = torch.nn.CrossEntropyLoss(reduction=\"sum\")\n",
    "    loss = 0.\n",
    "    with torch.no_grad():\n",
    "        correct, total = 0, 0\n",
    "        for (x, y) in loader:\n",
    "            x, y = x.cuda(), y.cuda()\n",
    "            logits = model(x)\n",
    "            preds = torch.argmax(logits, axis=1)\n",
    "            correct += (preds == y).sum().item()\n",
    "            loss += criterion(logits, y).item()\n",
    "            total += len(x)\n",
    "        return correct / total, loss / total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.1135, 11805508.464)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evaluate(model, testloader)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:07<00:00, 246.29it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 347.23it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9844\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:09<00:00, 202.46it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 343.35it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9867\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:07<00:00, 256.00it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 346.95it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9897\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:12<00:00, 149.80it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 342.27it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9917\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:07<00:00, 256.36it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 344.25it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9928\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:07<00:00, 255.57it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 347.69it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9927\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:08<00:00, 217.34it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 343.95it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9931\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:07<00:00, 256.49it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 347.81it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9936\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:07<00:00, 256.48it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 347.29it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9933\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:07<00:00, 256.88it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 347.28it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9933\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "for epoch in range(n_epochs):\n",
    "    total_loss = 0.\n",
    "    pbar = tqdm.tqdm(trainloader)\n",
    "    for (x, y) in pbar:\n",
    "        x, y = x.cuda(), y.cuda()\n",
    "        optimizer.zero_grad()\n",
    "        logits = model(x)\n",
    "        loss = criterion(logits, y)\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        total_loss += loss.item()\n",
    "    scheduler.step()\n",
    "    acc, _ = evaluate(model, testloader)\n",
    "    print(f\"Test accuracy: {acc}\")\n",
    "    pbar.set_description(f\"Epoch {epoch}: Loss {loss:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "# torch.save(model.state_dict(), \"mnist_model1.pt\")\n",
    "# torch.save(model.state_dict(), \"mnist_model2.pt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mode Connectivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = LeNet().cuda()\n",
    "model.load_state_dict(torch.load(\"mnist_model1.pt\"))\n",
    "weights_1 = torch.nn.utils.parameters_to_vector(model.parameters())\n",
    "\n",
    "model.load_state_dict(torch.load(\"mnist_model2.pt\"))\n",
    "weights_2 = torch.nn.utils.parameters_to_vector(model.parameters())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_bezier_curve(weights_1, weights_2):\n",
    "    def curve(theta, t):\n",
    "        return (1 - t)**2 * weights_1 + 2 * t * (1 - t) * theta + weights_2 * t**2\n",
    "    theta = torch.from_numpy((weights_1 + weights_2).cpu().detach().numpy() / 2)\n",
    "    return curve, theta\n",
    "\n",
    "curve, theta = get_bezier_curve(weights_1, weights_2)\n",
    "theta = theta.cuda()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 313/313 [00:01<00:00, 299.77it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9933\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 313/313 [00:00<00:00, 344.51it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9933\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 313/313 [00:00<00:00, 346.92it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9553\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "w = curve(theta, 0)\n",
    "torch.nn.utils.vector_to_parameters(w, model.parameters())\n",
    "print(evaluate(model, testloader))\n",
    "\n",
    "w = curve(theta, 1)\n",
    "torch.nn.utils.vector_to_parameters(w, model.parameters())\n",
    "print(evaluate(model, testloader))\n",
    "\n",
    "w = curve(theta, 0.5)\n",
    "torch.nn.utils.vector_to_parameters(w, model.parameters())\n",
    "print(evaluate(model, testloader))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_epochs = 5\n",
    "optimizer = torch.optim.Adam([theta], lr=1.e-4)\n",
    "scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=n_epochs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([213206])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grads.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:12<00:00, 151.85it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 343.60it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9891\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:08<00:00, 232.31it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 349.25it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9899\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:14<00:00, 133.84it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 344.55it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9879\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:13<00:00, 135.04it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 342.73it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.989\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1875/1875 [00:13<00:00, 138.75it/s]\n",
      "100%|██████████| 313/313 [00:00<00:00, 319.25it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9878\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "for epoch in range(n_epochs):\n",
    "    total_loss = 0.\n",
    "    pbar = tqdm.tqdm(trainloader)\n",
    "    for (x, y) in pbar:\n",
    "        x, y = x.cuda(), y.cuda()\n",
    "        optimizer.zero_grad()\n",
    "        \n",
    "        t = torch.rand(1).cuda()\n",
    "        curve_t_fn = lambda theta: curve(theta, t)\n",
    "        w = curve_t_fn(theta)\n",
    "        \n",
    "        torch.nn.utils.vector_to_parameters(w, model.parameters())\n",
    "        logits = model(x)\n",
    "        loss = criterion(logits, y)\n",
    "        loss.backward()\n",
    "        \n",
    "        # chain rule\n",
    "        grads = torch.nn.utils.parameters_to_vector([p.grad for p in model.parameters()])\n",
    "        theta.grad = torch.autograd.functional.vjp(curve_t_fn, theta, grads)[1]\n",
    "        \n",
    "        optimizer.step()\n",
    "        total_loss += loss.item()\n",
    "    scheduler.step()\n",
    "    acc, _ = evaluate(model, testloader)\n",
    "    print(f\"Test accuracy: {acc}\")\n",
    "    pbar.set_description(f\"Epoch {epoch}: Loss {loss:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_surface(weights_2, weights_1, theta, N, model, loader):\n",
    "    u = weights_2 - weights_1\n",
    "    v = theta - weights_1\n",
    "    v = (v - u * (u @ v) / (u @ u))\n",
    "    us = torch.linspace(-.5, 1.5, N)\n",
    "    vs = torch.linspace(-.2, 1.2, N)\n",
    "    grid_us, grid_vs = torch.meshgrid(us, vs, indexing='ij')\n",
    "    accs = torch.zeros_like(grid_us)\n",
    "    losses = torch.zeros_like(grid_us)\n",
    "    with torch.no_grad():\n",
    "        for i in tqdm.tqdm(range(N)):\n",
    "            for j in range(N):\n",
    "                x, y = grid_us[i][j], grid_vs[i][j]\n",
    "                w = x * u + y * v + weights_1\n",
    "                torch.nn.utils.vector_to_parameters(w, model.parameters())\n",
    "                acc, loss = evaluate(model, loader)\n",
    "                accs[i, j] = acc\n",
    "                losses[i, j] = loss\n",
    "\n",
    "    return grid_us, grid_vs, accs, losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import islice\n",
    "import cmocean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 30/30 [00:23<00:00,  1.27it/s]\n"
     ]
    }
   ],
   "source": [
    "grid_us, grid_vs, accs, losses = loss_surface(\n",
    "        weights_2, weights_1, theta, N=30, model=model, loader=list(islice(testloader, 40)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f5d243948b0>]"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD8CAYAAABq6S8VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABkWklEQVR4nO29d5gkV3nv/3krdJqeHDbnXYXdFSAQwQKMMEmArzAmmGjAJJtg/8AJLr6GC7YljA22r7BBxphggkG2dcUFDCIIMCCQMApIWkm72tXm3cmxU1W9vz9O9UzP7ITenZ6Z7pnzeZ56uvI5XV39rVPvec/7iqpisVgslsbFWekKWCwWi2VxWCG3WCyWBscKucVisTQ4VsgtFoulwbFCbrFYLA2OFXKLxWJpcGoi5CLySRE5KyK/mGP7K0XkbhG5R0R+JCKPrkW5FovFUo9UoYkiIn8nIgdjbXxsxbbXiMhD8fSaasqrVYv8U8DV82w/DDxNVS8DPgDcUKNyLRaLpR75FPNr4nOBPfH0JuAfAESkA3gv8ETgCcB7RaR9ocJqIuSq+n1gYJ7tP1LVwXjxNmBzLcq1WCyWemQhTQReAHxGDbcBbSKyAXgOcIuqDsSaeQvzPxAA8GpR6fPk9cDXZ9sgIm/CPJ1IuO7j1rU0L2e9Vi+qSCkHgJdKrXBlLGuRoFAAjcyC46OOC+KASE3Of2xwqE9Vuxdzjqt+eZMODOar2veeXwzcC1TufIOqno+lYRNwrGL5eLxurvXzsqxCLiJPxwj5U2bbHl+IGwC2drTrHz/jquWr3GqllEN67wOayfZsIN3WidToz2OxVEsUhhTHRhg5cwo0RAAFI+bJVnA8cFzU8eL5yskHx533/G+78aZHFlvHgcE8X/v351e175aLPptX1SsWW2atWDYhF5FHAZ8Anquq/ctV7pomN4AMPgxA25adJDJNK1why1rFcV1Sre2kWtuJopDS+BiFsRHyoyOQHyKW9VltvZOCn+lGk82QyBqBb2xOAFsqljfH604AV81Yf+tCJ1uWqyEiW4F/B16tqg8uR5lrGlVk5DgyfgYVh84dF+H6/krXymIBwHFcks2tJJtbadkwtV5V0ShCw5AoDNDIfIaFAuOD/TB+Bmf8jNkXgUzXlLA3HjcDbxORL2I6NodV9ZSIfAP4i4oOzmcD717oZDURchH5AuYp0iUixzG9rj6Aqn4M+FOgE/j7+LU+qKfXklVFWELO3oNohDoe3bsvRsQOF7DUPyKCuC64Li6JaduautejUUQpn6OUG2e8vxcmenEmeleotvNThSZ+DXgecBCYAF4XbxsQkQ8At8ener+qztdpCtRIyFX15QtsfwPwhlqUZZmH4hjS9wCgNK/fTKp1Qa8li6VhEMchkWkikWmiqbMHVSUo5ClNjGP6BOuHKjRRgbfOse2TwCfPp7yGNzRZYsZ7keFHAKF92278VHqla2SxLCkigp9K23sdK+SNT1hCzv4C0RAVh65dF+O4q/Nnvf/MAJ+9/QAP943gew4vvGwnL7hsJ65jTUeWtc3q/MevBVQh148MHQEg09lDprNn1boW3nT3Ib743w8RqpJwHUphxI13HaJ/PM8br9y/0tWzWFYU25RpRAqjyKn/xhk6AuLQsX0PTV3rVq2I3396gC/feRAEnrRtHX/2/F/i0nXthFHEdx86QSEIV7qKFsuKYlvkjURQQHrvM2YUhOYNW0g2t65aAS/z2dsPUAwjNrRkeMfTL8cR4fl7t/Nw3wiRKgPjeTa0Wh95y9rFCnkjEIXI2CkYOw3EZpSObmQN2IaPD43xyOAoKc/lpZfvwYkfWr5rRvqJCIXQtsgtaxsr5PVMhR1cABWXzh171tTgnh8+fJIwUkoa8fk7HuBff/4QAIVSSBhFiAjuKn8jsVgWwgp5vVIcR/oOICgqDm1bduCnMytdq6VDFQhNYCUNzUTIXSfOEKoZvt07fm5AI1dgXWoCAgUnBeLXLBCTxdIoWCGvRyb6Ym+UVWQHj0oQTUCYgygHpUHK8TWmPs/lzMg4AFtbXK7cnJxc/9WDOUaLSk+TS6J0DEqVRwl4rUbYnaSZJBFH23OA8qfFsjqwQl5PRAFy5p7G9glXBS1COGGm4hmMUFeKtQCC39SCOC7iOIjEn87UJ+IQyQAQ8rhd63jd03YBwmgh4MYDPybhOjzv8m20bt1MWCoQlQqEpSKF4QEIhoGhKios0+f9NpDk1APASdpWvqXuaTCVWMUURpD+hwBtTJ/wqASlfiicYKZoJ7JtuIk0bjKNm0id18Mp4blQCMmVFNc3LfJ//a8j5VPzq4/ZguP5OF4C0iZ+faZzI2CCMEVBMZ4C0MgEZdJo2rxGZrk0MTbjTaGScis/CW4zeC22VW+pG6yQrzQaISMnTKRChPatuxrHFq6RafnmjgBlzxGHdOdGvGQGN5FatGfNlo4mBsaL/PBgH1cfH+Lek0Pc9HMTd/8Fl2+mOTV3x6+I4PrJyQdAtahGRKUSUWBa+FFQICoVKU0MY0TeRODD7wS/wwi7FXXLCmKFfCUpTSC995sOTceje9fFjeFSGE5AqQ+KZ+MVQrK1m2RzB26ithmIXvlLOzhwapjRfIk/vvG/KQYRDsL+za381lN317SsMiIObiKJm0gy8zGhcct9/MxR8wZSikPr+10Vot5Ab1KWVYEV8pVAFcZOI6MnAKF183YSTXWe1k4jI1r5o0yZHlyy67fipZuXzAz02G0dvP+Fj+bD37if3tE8Xdkkr3nyTp61b8OKxFgRcUg0tZDYud+EVc2NMn7mmHmwlfrMTn43JLrBbZA3K0vDY4V8uQny8ejMCBWXrl0X1XeHpgZQ7I1t3wBCunMjiWz7stX7cds7+dybZ80OuKKI45BoaiWxs9WI+sQI42ePQ6nXTMmNkFhvzS6WJaeOFWSVoQoTvcjwUQCa128m2dJWvx2aUdF4nBRjezAO2Q078FJN9VvnFUQch0S2jUS2jSgMGH7kABROQuE0NF0Krg21utYQkauBvwVc4BOqet2M7dswcce7gQHgVap6PN72QaCcQPQDqvqv85VlhXw5CItxqNmo/lOvRUUYuw8I4hUuzZt24SWtEFWL43q079xPcWyI8bOPwPi9kNgAyfUg8ycRtqwORMQFPgo8C5P14nYRuVlV76vY7a+Az6jqp0XkV4BrgVeLyPOBxwKPAZLArSLydVUdmas8K+RLTW4QGTwEQHbdJlKt7fXZoo0CKJ6abIEnmjtJtXWft8eHZYpEtg0v3cTwIw/E1/Y0pLeD12E7RFc/TwAOqurDAHFuzhcAlUK+F3hnPP9d4KaK9d9X1QAIRORu4GrgS3MVZo13S0UUIqd+jjN4yISa3XER6baO+hPxqAj5YzB2ZyziLi1bLqGpe7MV8RrguD7tO/fTvHGXWZE7DKM/h3B8ZStmWWo2Accqlo/H6yq5C/j1eP6FQLOIdMbrrxaRjIh0AU8HtsxXmG2RLwUVuTMzHd1k6jFWeJiD8QNM+X+7tGzehZuwJpSlwEtladtxGcXRASb6jsP4/YALmV3WZbFO0GJE8fhYtbt3icgdFcs3qOoN51nkHwDXi8hrge8DJ4BQVb8pIo8HfgT0Aj9m6o86K1bIa4kqjJ1CRk8CQtuWnSQydRYnOxiDiYco3xfJlk6SrdaEshyICMmWThLZNnKDZygM98LEg4BAosf4oTsZK+qNQZ+qXjHP9hNMb0VvjtdNoqoniVvkIpIFXqSqQ/G2Pwf+PN72eeDB+SpjhbxWVA7uETeOk1InHVuq8QjMh4EIgFTbOpKtXfXt+rhKEccl07mRdPv6KZfFSQ8hgeQGY0d3azu4yrKs3A7sEZEdGAF/GfCKyh1is8mAqkbAuzEeLOWO0jZV7ReRRwGPAr45X2E1+ReLyCeBXwXOquo5CRTF2BX+FngeMAG8VlX/uxZlrzgaIaOnYOwUAC0bt5LIttSPKSUcg3Fj5in7gCebO0xQqgYijJTbD/dx8Mwou9c18/gdXbhOnVzjC2Smy2JpfJiJvpOx2+JJwIHUJjNq1Hq7NBSqGojI24BvYNwPP6mq94rI+4E7VPVm4CrgWhFRjGnlrfHhPvCDWENGMG6JwcwyKqlVc+xTwPXAZ+bY/lxgTzw9EfiH+LOxiW3hU63wOhrco5Fp4RXM6NFM91YS2Tr2W5+HMFLe9eX/5sCpEfKlkJTvcsmGFq57yWMbXszLOK5nzFwtnURBkeLYELmB06YjOn8cUttMbJcG/P3WKqr6NeBrM9b9acX8jcCNsxyXx3iuVE1NVEdVvy8i2+fZ5QUYf0kFbhORNhHZoKqnalH+slPPQ+w1NDFQJkdiurRuu7R+zDwXwO2H+zhwaoRcydj1c6WQ+0+NcPvhPp60q3uFa1d7HC9Bqq2HVFsPpdwYY6cOQ/4I5B+B9E7w2qygW6axXM3HuVxxpgm5iLwJeBNAe6ZOvSfCEnL2nooh9nViC49KRsCL5UtaxyMxVTGZIOKMQMyYNJyalyQHzwySL03vtC+UQg6eHV2VQl6Jn87StmM/pfFhxs8ehdwhwIHMHvDqpPFgWXHqxA5giN13bgDY2tE+d9qYlaIwivQ/CCjZdRtJtdaBX3hUNANNJiMRujRv2omXXOGATapAAFoA8uZTC6AjzJcRaDZ2dxZJ+ZCryAKU9B1296wNIRMxMd39ptbYffEETDwAONB0Mbh15hllWXaWS8gXdMWpa4I80nt/3FIUOrbtxkut8BuDBlA4VRELxaVl8+6ah5Gtri6KEetxk85NB5grOYOXaka8JI6XMNmAxGQCEsekX5M4HZsJ5ytEQYGnto5z6d0Pct+JHIUSJH24dL3D4zcfh2ACJM7ViV/xWQdvSTVmyn2xncJIH7mBU1P+6Nl94CRWuoqWFWK5hPxm4G3xMNUnAsMNYR8v5ZC+A7GAQ6a9k6audSvr8aEah5M9Eq9wadlyEa6/jH9iLYJOxMLdR9mlcQoHP9uJEwu24yURN3FBby+un8L1U1z/lifx4wP9PHBihN3dPo/f7hOO9ZqHxpwNfJNSDmmL07e1gtSpye48EMch1dZDsqWT/FAv+aEzMHY3ZC4ymYssa45auR9+AeNK0yUix4H3YlxoUNWPYXpunwccxLgfvq4W5S4ZpQnjjaIhCmQ6usi0d+F4KxzoKpyIR2NGgBMHs1omE4pGEJ2C6OyMDYKf7cJNZHATmQsW7IVwHeEpe7t4yt6uqZVt61FV0IgoLKHxFIUlNArQsESQGwEdNA/A6BQg4KwDp9205BsYcVzSHetJZNsYOf6gGVyU3BSHzq2zfhHLklIrr5WXL7BdmfKRrF+KY0j/g6YjE8h09pBu71x5l0INYPQXlCMSZro2k2heRvt8NAGh6RsAh2TretxEBieRNqaQFUREQFxcxwV/bmGOwhJBbpjC0CmITpsJAWc9OG0NLepuIkXb9n0MHbnfeCsVTkPzZSB11QVmWULsLw2xCeX+KQHvWke6rXPlvVFUIRgwgZaAZEsXqfb1y1evaALCg5jh/EK6aydeqjE7GB3XJ5HtIpHtmiHqpypa6l2x6SUV290bx84ujkvbjn3Gdt5/EkbvgvRu8FtXumqWZWBtC3kUmLgo48Zc0NS1nnR7nYx6jIowdi9GRB2aN+1enpjgGhlTRHiMsvHZb+ok2boecVbH7XKuqA9RGDoNUe9se4PTGbfYywJfn9dBREi1duMlM4yePAS5hyDnQnYvODaWzmqmPu/IpWZGth4cj66de1behALxiMzT8RBtSHdsJNnatfRmFC2ajsuo7AUjJNs24mfq5MG2RBhR7yaR7UZV0bBIVMoTlgpEQZ5gYmgWgRdwN4O012Wr3Us10bZjP4Xh2LNl7B6b2GKVUwfKtczEvuBmWL1Dx9adK+9KCBVmlCOYlrBLy5Y9Sx+VUCcgKJtPABzSXdtxk9mV95FfZkQkdo1M4pVviY6tscCXYoHPURw5E7+xHANccLeD1FcoWhHj2ZLItjN8tJzY4hT4PZDosomhVxlrR8jDInL23tgTRWjesIVkc2t9iFUwFoczjQAhu2EXfjq7dOWpgo5BOBUN0c92k8h24XjWF3kmRuATOF4CL91CormHqJSjND5AabwfwkMYG3tPbIapHzOG45nEFkF+nNGTh6F01kw4kNpsQufWqanIUj2r/xfUKI6LYkwVmc4eMh3d8YCTFWaaHRwy3VtIZJcwFZwq6DCERyjbvxMt60lku1a1+aTWiMiku2WybSNBboT8wDFjlorOAA64W2P/9TpoKGDMLe079xOFAcWxQXL9pyB/1Ey4kNkNbrZu6ms5P1avkEehsYOPnJiMTti5Y/fyDpyZCw2Ni1gcF6UcIGnJxFQjM3AmPE45nG2ybTN+U/uKuw82OiIOfqYNP9NGFJYojQ9SHDkdPywF3B0gLXUjkI7rkWrtJtnSRVjMURgZoDjaHw/5F+OD7mXNsH/bUm8YVtcvpQpBDskNmFY4oOKY6ISZJTRVVIsqlAbiUZllO/gSjsrUEKJ+iMrREIRUxza8dJ2YlFYZjuuTbOkh0dxNkBsyrfTwYaYEvRnq5MEpInjJDF53hkznBorjw0z0noht6ZN7GdOLmzURF50VHhBnmZPGF/KwCIURZPgYaGjEG0Bc2rZsx0/XSadOMBqnWItHZW7chZdaomBHWjCeFpPeFmu3A3MlEBH8TDteuo3SxACFwROxoAPSEQ9AqiNRd1ySzR0kmzvQKCQoTBDkJ8gPnjXhIEr9wCNM2tW9divqdUbjCLlGEORNAKtSzsznh5DY1qtAqqWNRCaL35TFXenh9GWiPIzdT3lQzZLZwVVN7JPwEFOxT1wyPTtxE3XyMFtjiAiJpk78TDthfoxSbohgYgDCgXgPF7yL6mpUqTgufroZP91Mun0dqkpUylMcGyY/dDa2qR9lMs+o1xLb1m0fy0xE5GpMZjQX+ISqXjdj+zZMerduYACTCeh4vO0vgecDDnAL8HvxCPlZqV8hj0Jk5AQEOcgPA0pZ+sy3ERCHpq5uEpksbjJVX63NGdEJU+3rSbUuQSerajyA5yiTHZjN3fhN1gOlXhBx8NIteOkWtH1zhagPQnA/pnN0Nzj1F47WdOymSXekSbWvIyzmKU2MmNb6ZJ5RAMf4qbvNsX29Pt42Voo47+ZHgWdh8i/cLiI3q+p9Fbv9FSbhzqdF5FeAa4FXi8iVwJMxuToB/gt4GnDrXOXVrZBLWIjzYArJbDNeMoWbSOIlkriJZH14ncyGRlDqi1suAC6tWy+ufcAtDeMBPKeY6sDchJ9ptx4odUylqEetGymN9VIcPWti2YROhS29jholMcaunsZLpk1rPYoI8uOUcqMUhvsmB7EZHOOvLj5IwphixI9NMm5dfr8a8wTgoKo+DBBHfn0BUCnke4F3xvPfBW6K5xVIAQlMCE8fOMM81K2Qe8kU3RftaxyvClUIBuO4KCa41JJEJ9RSbP+eagmlO7fjpprr643EsiCO65Fs3UCiuYfSeD+F4VNTPunudpDWuhY8cRz8TDN+pplM50aiMCDIjRHkxyiMDFQkO5n1aDP5naYF72bASa/o941KEblTE9Xu3iUid1Qs3xAnxikzW1a0mXmK7wJ+HWN+eSHQLCKdqvpjEfkuJoOaANer6v3zVaZuhRxoDBE/Z0Sm0LRuO36mpbbCqnkIHqRyBGamZ5e1f68CxHFJNPfgZ7sojQ9SGDoB4WGMoG8xHaR1LOhlHNcjkW0jkW0j07UZAI1CojAgCuIQw0Fpcrk0PgKlXjNNnQUS3bG4N5nWfH1+9z5VvWKR5/gD4HoReS3wfUyynVBEdgOXYhLwANwiIk9V1R/MdaK6FvK6ZjYB79mG31Rj175oPI5AGI/AbOok0dyN49XP6EFLbRBxSGQ78Zs6CHLD5AeOxn0fxxpK0CsRx4QYnivUhOlMLRAUcoSFibglP9OK4EJqE3itjRT8a8GsaKp6EtMiR0SywItUdUhE3gjcpqpj8bavA78EWCGvGRrGNvCpwTVLIuBagKCcRILJFpvj1ok3jmXJMO6LbXjpVsL8KLn+I1OC7l28KrIclTGdqSmTorC5nUzXJlQjwmKesJAjKExQHB2s6HMSkzzD76x3F8jbgT0isgMj4C8DXlG5g4h0AQOqGgHvxniwgHELeqOIXIsxrTwN+Jv5CrNCXi1RIe6lL9v9HJp6ti6BgEfxUO/TACRaNpDIdtoOzDWIiOClW8huuiwW9MPm4e7uBqcx48JXg4hjBislMyTpJNO1mahUoDQxSm7gNBSOmwkX0ttNS73OUNVARN4GfAPjfvhJVb1XRN4P3KGqN2Oyql0rIooxrZST79wI/ApwD6a1+J+q+pX5yrNCvhDhBIw/wJRt2qV54068VK07MWfGQXFp2nCxbYFbJgW9af2ljJ8+EJvatpoAXWuAylZ7qq2bsJinMDpgPGVyh1a6enOiql/DpLmsXPenFfM3YkR75nEh8ObzKcsK+VxMG4kJydZuUq3dS5O3U/MQPEA5+mG6exdesg5CCljqCsdLkN24j7GT9xlTixZNqroGs5svFjeRItO5kXTHBoLcKGOnjy580CrHCnklqhAMQ24qvGuqfT3JliXK26lhnD/SmGuSbRvxm5YhiYSlYRHHJbtpH2Mn7ovvnV7w9q/JATimL6GF9p37gdMrXZ0VxQo5zOIDLqQ7N5JsXqLsOJOjMR+JV6wWM4oCJUyS6BRmdLGl1og4ZDftozh61kRaDO4Bb68ZcGNZk6xdIVc1yYVL/RUdmEsYCwVMC1wHp4WTzXTvwk3W39DshSmLdh4kDzIIRJNv+SYqhAPaDdrE1CA1Sy0QEZIt63C8hHFTDH4BzuY4gbS9zmuNmgh5FcFhtgKfBtrifd4VdwQsP1EBin0mLyblGDQuTeu21H4QD8St79G4E7PcYVqOB95Rp2YUjacwniLzKSFQnEO0BT/dgeOnQTyi0hiliQHEOTO1j7YCTbGwL1UbovyAKYAU4/qOxJvWg7awmh4ofqYdx08zceYhiI6bkMV1FgPdsvQs+t9UZXCYPwG+pKr/ICJ7MT252xdbdtVoFNu+DzMVGdAh07UJv6kNx62x+aQs3tEQaP/kaj+OhOckMnUm4DlwykG34pBkc1SvLNpeuh3Hz+B4aRwvde4o3FQrieZNRGGRqDhGYeQkyDAiw/F5JBbVFGgKSGKe8dVSFuwiSAEj3COATqt7ub4A4pxE9RREW4BGfAuaHddPkd20v2IQURwy11kHTkddRVe0LA21aBZVExxGgZZ4vhU4yVKiClHOuA6G4xVDgIVU+zqSzR21jwxYzoOpgyaZwyQuqc4teKnmOgw5UATnMCIRquClWkA8xHERcSH+nFr2EMc7r+/huAmcdAdeusMkMQ5yhMUximNnpwk7xOJO1gj7pLhHpp5SFm0TCXN2wXbwM504XhLx0jhucrKPQ1UJC8MUho8i7lFUHYi2x2U0PlODiFpmST0n4GyKRd2OR1iN1ELIqwkO8z7gmyLydkxT6JmznUhE3gS8CaCrqcpWhEYVoh3bvCdb3WVcsuu34qVrHFhK1WSh18GKJA6mvPoVb4ASOIcpm3q8TA9+pnvJBx2JCOJncPwMflOPEfYoQIMcUZAjCvIE+RHEGZ31+Mm3gVQL4iZx3IT59BLIAqP8RAQv1YabbCGY6KM4dhqch0HbQbtYLd1FM1PPBRNDJhhXdNxM7jaQdmt2WWUs1937cuBTqvrXIvJLwGdFZH88NHWSOHrYDQC7ulrPDaKuGg9dH4dwbA7Rdki2duMmTLhNx0/WXrzJGbNJdJYpO7tDqmOLEZl6DbFLCNJvJsDPdOJnepAV8pYREVO26+MmzQtbsjUOtBTkiYKceRtwEzhuEsRd9G8p4uA39eClO5jofcDY+2UQtMeI+irytHFc38Smz3YRFSeY6D0Ue0odixNarJ6h/mudWgj5gsFhgNcDVwPEIRpTQBcwX5xL4+URjpupcIq5RNtLZnCTadMyW7IM9KU4/vcZpol3+2aTA7Ouh9BHsWCdRQRUHdKde+o28JY4Lm6iCTexdHZscTya1u0jCvLk+g8izllUe0E3rLoOURHBTTaR3XQZpfEBCkPHzVB/pyceUFTP966lGmoh5AsGh8EEgXkG8CkRuRTjZNzLfIQ5GP15xQohkW3HSzXhJjO4iWXKCKTF2NbYF69wSLZtxEu3Lc0goZoQdwTKBDAR26KNgKfadxnPEgsAjpeiad1+wsIo+aEjUx2iq1TQE9lOvHQr46cOmDfKqLchYp9b5mfRSlRlcJjfB/5RRN6BUZnXzpd/DsDxfVJt6/BSGdxkZvlF85z43y5N6+u1FVsW7nGMcI+c48+dbNuOm7DD/ufCTTaT6dlPWByhMHR0VQu643o0b95PUBgn13sojn3ugHcJSD3e35aFqIk6VhEc5j5MDrqqcRyXdMf6WlSverQI0TBEJ5mM/53tIpHtrtP8l9NNJlDhvZFdj5vIIm6N+whWMSaVWSvurILeGfvB1+N9cGF4ZXPLWB+F4ZMQ3AdON0gbSJNtoTcQ9WobWB40AB0307SOS6nz+N/n2rz95g24fpMV7hpwrqAfA+lDnL7Y/319LOqNf51FhERzN16mjfFTD8TeV2Wrpwvu5nhw0dqWinpnbf06WoyFeyz29a607jgkWtfjp1pw/HodQBGBDIGcmRTwZNuOJe0UXMuUBd1b10oUlggLQxRHTyPOqbiV3gPaxvkNZKpPHNenefN+NAoJ8qME+RGCicp4QE48wKgFWNncmpZzWd1CroGJ8R2Ngg5xjnC39OAmsriJTB27DAJoLOCnYwGXWMCtzXu5cFwfJ9ONl+4iKo6RH3ok9nQ5C9phJurx7e38EMed9EPX9i1EpRxBfoTiyFmITpkJMd4uTpdtqc9DFaFLtmGyAnUDA8CrVPV4vC3EJJYAOKqq18xX1ur7FWaxcxsckq0bcJNNOH66QcwPakYyymlE1Ah463acRLZB6r/6MK58zTSt209UmiA38DDIgJlwINqGccpqfExChwxuIkOyZb15K8mPkh88MSXqTldsV18d37lWVBm65K+Az6jqp0XkV4BrgVfH23Kq+phqy1sdQq75WLxPMdPO7aVbG0i4y6iJGyKnZgh4jUemWhaF42eMoIdFgok+ShN9iHs4Hv6/EciyGuzoZRzXx2nqwG/qICzlmDhzKB5b0Qc44O4EyVqzi6Ga0CV7gXfG898FbrrQwhpTyFVjW/fIuZ2ULevx0q24dWvnno+ZNnAh0bINN7kEURktNcNxEySaN+I39RDkBiiOnUbc43HEx464YzTJahJ110/TvHk/UViiNN5PceRMnIJOwN0ShwGoZ3PluUQlZfxUqdrdu0TkjorlG+KR6WWqCV1yF/DrGPPLC4FmEelU1X4gFZ8/AK5T1Zvmq0zjCPlkLO8TzGYy8dKtdeoiWA1B7IXSNyXgrdtwbQu8oRDHM8P/M90V3i4DiDMQe7t0r5rO0TKO65NsWU+iuYdgYtCYXcKjwFFwNpq8oqvTjt6nqlcs8hx/AFwvIq/FJF8+wdTAlW2qekJEdgLfEZF7VHXOBKX1f4U1V/H6BiD4TR24yWa8VBZx6v8rzE0xtq8OWi+UVUSlt4tGAUF+mOLoyYrO0fa4c7RRGx7nYmLYdOJlOggLY+T6jph+qugkxuyyJR49unoeYguwYOgSVT2JaZEjIlngRao6FG87EX8+LCK3ApcDjSjkIZTuYqr17ZLp3lGHsbzPlzjolnOMyu+W6tiF4zWiOcgyH+J4cXCyTqIgR67/0FSgLhyIdrC6BF3wUs00b76MsJgjyA1SHO2bltYQd2vsm95YppfzZMHQJSLSBQzEwQPfjfFgQUTagQlVLcT7PBn4y/kKq18h1yKQIdGyHr9piZIfLwuxcMtEHHUwilvf4Ge68DJdOO7q+SNb5sbx0nHnaIkg10dpvBecQ6sulG4ZN5HGTaRJtGwgKk5QmhikNN4fhwQApAOcOA2gJOLP1SHuVYYuuQq4VkQUY1p5a3z4pcDHRSTChOO8boa3yznU7Z3j+Cmymy5rwNZ3xJRwDzA9JZqYxAd+k7F/13XERMtS4bg+iewG/HQXE30PVoTS7Y5NLqtDzMqUoy+6ySaSbZsIC2OUJgYJJgYgHJjtiNgMUxb3VBwyoLGuSxWhS24EbpzluB8Bl51PWXUr5IjUqYiXU4yVprLWUAIZpTJN2mQey0wXTqLJDJ9vaHu+pdaI688IpdtbEUp3dYQAmEnZ9OKlmtH2LWhYIgqLaFAkCktoUKQ0MWQG8p0TVy8eXSrNIBnr5liBVZYFUcoRBZEhYPr9M3WvmTyW4vi4fgbHb7ItbktVTIbSLY6RHzwchwA4DdFmjC/66kREEC9hvM0qgi6mOkwfYTmDVFTMERRGKY31x2NFTsV7OuBsAKd52eteb1ghnxXFJPMdNu5jky1sB7+pK04xFk+OX6dvDpZGw01kTSjdwhCF4eOIeyweXLSLtfhXLWeQctI+XrqFVNumOObNmDHPjA9CdOLcfDNrkLV3d8yLAuPgHI9HVIKJybI5HpTTWDY6S+NhTA/tuMlWglw/xdFT4ByMvVtsrHAT86YdP9NOqn0LUVAkKIwC/Qseu5qxQj5JDpyjcUZ5wW/ehJdstXZty4pgkih34/gZ8gOHTKLoaDNgzQiVOF6ChNe50tVYcaxKUQLn4VjAwc9uxEt32Na3pS5w/SbSXZeQ63vADPuPumJXRWvOs0yxhtVKjV+3cxCI8DI9ZLr34We6rIhb6grHTZDp2Yeqgzh94DyINQxbKlmjijUGzgOIcxZwSHddTCK73nqZWOoWEYdMzz787AYgisW8uNLVstQJa8y0UmlGMaFh3WTLSlfKYqkKETF2czdFfuiwGRUabQVsbJ61zhppkUcgfZNmFL9pHZmefVbELQ2Jm2wm3XkxJlPP0XgE8czBM5a1xCoWcsVEFxwC50HE6cWYUS7Bb1pn7eCWhsbxkmS69wIO4pwB5wFgBCvoa5OaqJmIXC0iD4jIQRF51xz7vFRE7hORe0Xk87Uod4pYtBkBOWtuaucA4h5CHDMKLNm6g6Z1+22AKsuqQRzX2M2bN5ll90RsOy+sbMUsy86ibeTV5KYTkT2YMI1PVtVBEelZbLmmtT08S2AqMDFOOnC8DI6XRryUHX1pWZWICH66Ey/VQZAfoDhywvica2fspmjfPNcCtejsrCY33RuBj6rqIICqnr2woiITnGoylyWYGCcdOH4ad1K018LNq+DkgQiiDNaveG0zKejJViZ6DyBOP6r9dhDRGqEWQl5NbrqLAETkh5jYvO9T1f+s7vQK5ONclkOTqdD8pvW4qXYc119s/RuECJwCODnwhgGdFh6X/DZs68sijjc9AJd7PI7XsroSWFims1z/fA/Ygwmk/nLgH0WkbeZOIvImEblDRO4YGo/ToDkPIO6ROPKgQ7JtJ5me/fhNPatcxCNwJsDrh9RhSB1BkqfAGwLAS3XiN23FS60HFFKPxGF1LZapAFxTfueHrHfLMrNQ36GIbBORb4vI3SJyq4hsrtj2GhF5KJ5es1BZtWiRL5ibDtNK/4mqloDDIvIgRthvr9wpzkJ9A8ClWxIqzhnT+s5uwku1rZ0BO1KE5PFpcc29VBeO14TjZZDKvIc+iJuiOHYEkseguD42tVjWOmW/czfZSq7vQcz/qReiPdi3t6Wlmr5D4K+Az6jqp0XkV4BrgVeLSAfwXuAKzJP3Z/Gxg3OVV4tfczI3nYgkMLnpbp6xz02Y1ng5T91FwMPzndTxkqTa99C07jL8TOfaEXFnApLHAfAzW0m17SXdvg8/vR7Xb54u4jGunyXZchEgkDgN7hC25WUpUx7ib9xuy6NC7dvbEjPZd6iqRaDcd1jJXuA78fx3K7Y/B7hFVQdi8b4FuHq+whbdIq8yN903gGeLyH1ACPyhqs4fd1IcHD+92Oo1Fs4YJM4CQqp1D+JUb9N03ASptkvIDx1AEgOoDgACYTNEPqgL6pkJF9s5urYQESPkXorC0CNmcJz2gLZh7gdLVILRM1W3bbtE5I6K5Rtii0KZavoO7wJ+Hfhb4IVAs4h0znHspvkqU5Mh+lXkplPgnfFkmZVwSsTbLr0gzxsRl1TbXsLiAFEwTlgcBXcEmfErV2Y1IkpPCbzWStwdCDMsTiBC07GrSdDV3BeyvHjJVpzOi+LUcmdRPRu7KrYD9jqfB32qesUiz/EHwPUi8lpM8uUTmIbuebPGYq3UMSnzAE4271yU+6SI4CU7IdkJTSZdFhqiUQnVEhoFaFQCLREUho0ph9qnPyzb9il1xqJeza0WgDsBfj9lr5xycg/ym6s8h2UhyqnlotIEuYGHTRRQ6QdciLZhE1jUhAX7DlX1JKZFjohkgRep6pCInCA2RVcce+t8hdl/Rj0gxh/cS3bieLU1J4kIiBcnyJh+bj8z2UmOajhLstsLQ6MSYWmEIN+LJPoqRL0DoqbYvFMmAHccfONRURZvL9WF6zcTFocJCgOQOgpBKwRtWFNAbXD8jBH0oECQ66U0MYC4Dxt3Vu0BbcFKxAUz2XeIEfCXAa+o3CHuLxxQ1QgzYPKT8aZvAH8hIu3x8rPj7XNif6UVRyFpwgh46RoMeL1ARNyamc3F8XC8NF6qB40KRKURSrleJNEP9BuhCFrAG0HEPDxUBS/VjZtoRZzk5Ehcx2vCTXZRGDlo/Oe9YQjajahbz4ua4HhJEs2b8ZvWEeSHKI6ejj1czgAO6HrQZuz1rp4q+w6vAq4V8yf4PvDW+NgBEfkAU15971fT6TUnVshXGteImZ/ZMqtHSiNjkuemcNwUXqqHKCyL+lnEH47FuwfHb8Vx536dd9wE6fa9RGGewsghxB9EvcHYbNOC7bitDeL4+Jlu/Ew3UZAnyA9SGu9FnJNTJq5oEyZsrr3mC1FF3+GNwI1zHPtJplroC2KFfEUJwDctVMdf/SF1HTeJ43bjpbpRjc67L8BxU6Tb9xEFExRGDyOJeBh60AFR0ky21VgTHC9FIrsBv2k9UWmcID9EkBtA3GNTpjLAdMSUPaHiSV2MtKSwgr88WCFfSVLGXzzZsnvNBfVaTIeu42VIte0lCsYojh1F/Km3TlUgzBpRD5uwt/jiEBHcRBY3kSXRvJGwMEpUGkWjkCA/CoxDbB6beQtPtuKtaWbJsXf5SuFMIBIZ08I8ZgXL7IgIrt9Mun0fGgVEYQ4Nc5RyveCOId6Yaa0jUFxn3Cxt63BRiDh4qVZItQKQbJ3aZryjIlQDNArRsERUHKWUG5gyzWhX7OZoZafW2Cu6IkSQOI0quMmula5MwyOOh+s0g9+MlzIdxlFYICwOEuT7kKS51gQdEDRjvV5qj/GOchFcc3l9INWK37yJqDRGfvARxOlDtQ/j5rgVY3qx1AIr5CuBZ6I4+k3b10jI3eXHcZM46fWmk7U0QnH8BOIPoN4A4MStdGvDXWqMaaY5dnPME+T6KU30I+7hOCrjRiCL/R0WhxXyZUXBG4g9NhxcP7vSFapbVJXDx8bobE/S2pzg5NkJvvrdmbHY4JefsI6Ld8zdUSzi4CbaSCfaiMI8YWGQoNCPJE+ZVnqp09jUbSt9yXG8FInmTcbNMTdAcex0HGaX2G+9FStJF4a9astGCKmjcUIMh1TbJStdoboiipSHjozwk7v6uP3ufn56dx99gwWu+8PLeclzt3P05DjXffwX5xy3d08rF+9o4We/6OffvnGUJz66iyc8uosN3ecOrHLcFE5mA156HWFxmNLEySnPl3IrXX1QB9NCtK3EpUAcD7+pBy/TTVgYoTB8bCpcgDW7XBBWyJeDOCwtgJfeiJfsWOEK1Q+FYkgy4TIwVOB5bzCB4Db0pHnK43p4/KM6ecoV6wB4/GVd3P3V/zHt2ChS0inTkn7k5Dhfu/UE//rVIwBs2ZDhOU/dyNtefQnNTdNjiIg4eMl2vGQ7UZAjLA4QFAZNvPeY6YNcY0GPUsa1LkrFHjG2Fb8YRAQv1YqXaiUKcgQT/aZz1D1sBo1FGzHZjewDdSGskC81zoQJLQskm3fgeE0rXKH64M77B7j+swcA+MRfXElXR4rr3/sELru4jU3rMue4Y7qu0JSe+3b99Wdv5QXP2MKBh4f56d193PbzXr566wne+Vt7ASgFEb53bn+E46VxvE14aeMvrRrE3hfh1CchYXEsTq2nsUdMHybsQFcs6ravYzE4XppEy2b87HqC/KAZXeqeiM0unaBpTIgJK1mzYa9KzVHTAndz4A2a5QsIS7tauf3uPq7/7AP818/O0tbi81sv3oOqIiI892nzRupcENcV9u1pY9+eNl73ot3k8gHJhEspiHjub32bJ13exW+//CI2rz/3YSri4ibmGZQVH6KqaJgnLA0bj5hEr0nWUFxn4shYFoU4Hn6mGy/dRVgcpTB0FKSfsk+AiQOTBdKgKYwJxr4ZWSGvJRWZfaAcP6QLL9VjvVOAL3/9CO/60M/pbE/yR2/axyuv2UE2s3ShU9Mpc3vn8yG/dHk3N/7nI3z5a4/wa8/awltfdQlbN56/8IoI4qXjWDLr0HCCwugRSJyBwiYTdteyaEwUzxa8dftRjYhKOaLSOMWxsyBjiIxO7qs1C7/cuFghrxVSqsjsswnHyyKOje9cKIac7c+zZUMTz37KRsYmAl72/O2TIrscNGd9PvCOx/CWV13MP/7rQ3zx/x3ma7ee4N8+ehUXzePxshBG1JtItV5EfvgBSJ6Ik2DbFmItMZ5HTbiJJvwmM05Ao8CIezBBcax3hWu48thmYk0ITL5MzHB7N9FuRRw4fGyUl7zte7z2j39IsRTR2pzgdS/avawiXsmG7jR/+rZH8e3PPpvXvGgXu7c1A/GoxEUgjk+yeZdZSB3FptlbesTxcJPN+E3raFq3f6Wrs+JYIV800bSkEI5r3aYA/uObR7nmzd/lxNkJ3v3my0j49XOrbehO8wev34fjCCfPTvDit32Pex4YICyOYEJDnz+OlybRtMWE5U09wnQxVyAyb22SNx3g7ig441jRt9QCa1pZFNHknzaR3Ybj2ez1uXzAn3zkTm665RhPeFQnH37P42f16a4X+ocKnO7L85K3f593vrqTVz2/nXTbxRf0RuUmWuNBR71o6ki8dvaAUmVM6IBWE47XprSzXCD100xqOKZaXommzbh+80pXqC7wfYeTZyb4vddcwr/89VPrWsQBLruonZs/diVPvTzDB/+5j7dee4KTR+8lCnIXdD4v1WOGngOun8VLtsUd3uvw0pvwm7aSyO4k0XwRftM2QMAbRlLHIHXYttJXESJytYg8ICIHReRds2zfKiLfFZGfi8jdIvK8eP12EcmJyJ3x9LGFyrIt8gvCiLiI4qU34CbaVrpCK85td/Zy6a5WWpsTfPavn4LnNkYbQVXJcIy//eMNfOErLh/63IP83ef7ed/vHCLRtAU30brwSSoQEdLte6vb2U3E0RtLhMVBk3AjeSZupbfFWZBsx2kjIiZLzEeBZwHHgdtF5GZVva9itz8BvqSq/yAiezFJKLbH2w6p6mOqLc8K+YXgjk2GoPWSnStdmxXnWz88xdvf/1N+7VlbuPYPHtswIm6IMLlChVc8v4krHr2Vres9RKA4fpz0eQr5hSCOj5fqwU12EwWjFMeOgTdkpqA9zlNqXewajCcAB1X1YQAR+SLwAqBSyBUou021AicvtLCa/OMWeoWo2O9FIqIickUtyl0ZQvB7URXcZPdKV2bF+c/vn+Ct7/sJl+5q5d2/3XjeAyIuyZY9gAPuKBdtS5Jt7SanG/m9vx7h8LHRBc9Ru7oIrt9Cun0fyZbdgCD+IKSOAOGy1WOtEgTCYJ9f1QR0icgdFdObZpxuE3CsYvl4vK6S9wGvEpHjmNb42yu27YhNLt8TkacuVPdFC3nFK8Rzgb3Ay+PXhJn7NQO/B/xksWWuKGUPlZaday6rz0y+8p3j/O77b+dRl7Tz6Q89mZZsY45cddwk6fa9pNr2kW7fh5/eQP+Ix8/vG+Dl7/gBBx8ZWYE6mbR2XnoTk/0xUlz2eljmpE9Vr6iYbriAc7wc+JSqbgaeB3xWzMjBU8BWVb0ceCfweRGZd8BDLVrkk68QqloEyq8QM/kA8EEgX4MyVwYnH5tUOnHc+u7EW2oKxZAPf/I+Hre/g0/95ZPPCUzViFQ+mHdva+HzHzYNoZe/4wc88PDwitTJS7aTbN5pFpLH485QSwNwAthSsbw5XlfJ64EvAajqjzHxBrpUtaAmJCeq+jPgEHDRfIXVQsgXfIUQkccCW1T1qzUob4VQSJgY1uUsNGsVVSWZcPnch5/CP1175bzBrBqZPdtb+MJHnkrCc3jlO/+L+w4OrUg9HC9DqvViQJDkGfAGsJ4tdc/twB4R2SEiCeBlwM0z9jkKPANARC7FCHmviHTHlg5EZCewB3h4vsKWvFcqflX4MPD7Vez7prLNaXC0zl4jvUFENB70sXY9CT7/lcP8z7/+OVGkbOzJkFmlIl5mx5ZmPv+Rp7JnezOtK2g6Escn1XYpqg7iD8V28wsbvGRZelQ1AN4GfAO4H+Odcq+IvF9Erol3+33gjSJyF/AF4LVqhhn/MnC3iNwJ3Aj8tqoOnFNIBbX4Fy70CtEM7AdujV9d1wM3i8g1qnpH5YliO9MNAPt2ttZPk8MdQ/whVB0c/8JjczQ6//n9E/yvj9zJ05+0jlIQkUysjQfatk1ZPv+RpyIiqCrDoyXaWpZf1EUcUm2XEhYHKE2cMnbzwmY7kKhOUdWvYToxK9f9acX8fcCTZznu34B/O5+yatEin/cVQlWHVbVLVber6nbgNuAcEa9PQvMa659FVUi1XbJmOzjvPzTMH1z3My7f285H3/fENSPiZcq/+/v+7i5e/o4fMDZRWrF6eMlOEtntgMYxfmzLfK2zaCGv8hWiwQjB6zeDfvwhTDzxi9dsKNr+oQJv/pPbaM36/P3/ftKaE/FKnv3UjRx6ZJTfv/ZnRNHKvTROxfSxKeksNRoQtNArxIz1V9WizKUhAH8Y3LKHgkPCBsLi4COj5AoB/3TtlfR0ru1r8eTH9vA/33IZH7j+bv7mU/dPZiBaTlTVhM0FKGzECrlldfdUVU0A/hC4ZX9hh2TLLhzXJgkAeOKju/j+55+zYuFn643XvHAnDzw8zEf/5QEu3tnC86/avKzlh8Uhk8S71GETWVgAK+QmD2OiPDLWIdmyG8dtzIEttebzNx8mjJRXvWCHFfEKRIT3/e6j6R8s0Na8vPdKFBYpTZwAJI7FYrGs9eiH7sikiCdbdpNu32tFPObHP+/lfX93F9/76WkWmXdhVZJMuNzw57/Ekx9nxhQE4dJ3OKoqhZGHzEJhM9akYimzRoVcIXUESZhM6KnWS9e8HbySY6fGefv//ik7tmT5yHsej+NYwZiPT3zpIV79+/9FsbS0Yh4W+kziilL3ebochqbfxx3FerisTtagkIdGxCXCTXaRatuLOGvXC2MmuXzAm//XbUSqfPzPnrQqht4vNeu70/z07n7e93d3LVkZUZCjlDtjYp2H2eoOkrwZOJR6BEn0I4les+z1YwV9dbG2DJ9SMAlyMQmS3UT7Cleo/vjBHWd56MgIn7z2SrZvqlIw1ji/+vTN3H9wmI994UGueuI6nv2UjTUvIwrGEYnzizo5iObJRiUlSJ5AJDIhJZKduMl2NAoojh0Fbxi8Uchvr3k9LSvD2hFyZwwSZ4E4t6ZNyzYrz37KRr79mWezdWPTSlelofi9117KD24/w3s+/HMet7+TzrbaepO4yU7E8SiOH0eSp1EV43o4zWtluveVm+zGS3VNhZRwId2+l1LujElH5+QgWtvB31YLDWpaiW1+qcNm8s+aEZjuiIkOJwUgwAQWUvAGkORZygN7rIify+hYiTt+0Q9gRfwCSPgOH3r3FRQKEXfc3Vfz84sIbqKNVNtevPR6QJHUCWMqkYK5/1NHwR3BS3aQar0YP71u1rhAXqrbdGAnTmGDb60OGqhFrkakE72YjC6YVgmAOwacm+C27G1h9nXiIfYN+uxaYj7w0bu5+dvHuPVzz2F9nefZrFcu3tHC97/wnCWNwyLi4CW7cBPtBPk+gnyvEXTMPV7N+AcRh0R2G6XxR9ByFiJLQ1PnQq7gFExvuzsaCzJ4qS7cRNs0TxNVBQ3QqIRGAaqlyXnHa8JNtK3ZOCkL8a0fnuLfvnGUt77qYivii6Qs4t/58Sku3d22ZMmnRVz89Dq8ZAdhcQjHazqvN03Xb6aoDniDpvPUBt5qaOpXyCWE1OFJ8QYHv2kLjpedVZBFBMRHHHtDng/9QwX+54d/zt7drbzt1ZesdHVWBf1DBX73A7dzxf5O/vmDVy5pA8Lk+7ywlIOp1j1mqH/yeNzxaRs6jUod2xlCQPAzm0i1XUq6fS+u32xb1TVEVfnTv7mT0bESf/Wux5Hw6/h2aCA625K86837+cEdZ/nCV46sdHXmRBwfP7Pe+KbbzEMNTd3+cx3P5Cx0E+1rOpHDUvPYfR384Rv2cfFOO9y7lrzymh085XE9XPuxe3jkxNhKV2dO3ESn6WtKnI2dBCy1YqGk9CLyERG5M54eFJGhim1hxbaZmYXOoW6F3LL0iAivf8kefuslu1e6KqsOEeG6P3osruvwhx/8GWFYn94hIkKyZZdZSJ6IR3/WZ10biWqS0qvqO1T1Mar6GOD/AP9esTlX3qaqC4YDt0K+BlFVfv8v7uA/vz8zF6yllmzoTvNn73wML33edpw6/qc5bmoqJ2h59KczgRX0RVFtUvoyL8eke7sg6rez07JkfOU7x7npW8d4zN6Ola7KqudXn768IW4vFJMTdC9hcYjSxImKQUcbQNdGHKJS5HBmtOqBXF0iUpnl7IY4VWWZ2ZLSP3G2E4nINmAH8J2K1an4/AFwnareNF9lrJCvMcZzAdd97BdcdnEbr7xmx0pXZ83w+a8c5u4Dg1z3h49d6arMiUkj146baCUsDlKaOIWkTsaCvgnURgatoE9Vr6jRuV4G3KiqYcW6bap6QkR2At8RkXtU9dBcJ6jjFz7LUvAPn3uAM/153vv2R9uohsvI0EiRL3/9EX7432dXuioLYgYddZJquxQv1YPJDXo8HkW6MrlKG5CFktJX8jJmmFVU9UT8+TBwK3D5fIVZIV9DnOrN8U9fPsgLn72Fy61ZZVl5/Ut2s3VDEx+4/m5KQWNEHhRx8VI9pFovwUt1AZFJ9uz3gRRXunr1zrxJ6cuIyCVAO/DjinXtIpKM57uAJwP3zVeYFfI1xPquFH/7vx7PH71x/0pXZc2RTLj8z7fs56Ejo3zu5sMrXZ3zQhwPP70+7hB1wB1BUsdNnCMnt9LVq0vOIyn9y4Avqk5L33IpcIeI3AV8F2Mjn1fIrY18jRCGiuvKkoRYtVTHM6/cwFOv6OFvP3U/L3rO1oaL9S6OT7p9LxoFhMUBgvxZNHEa8ravZTaqSUqvqu+b5bgfAZedT1lWyNcAxVLEi956Ky//1R28wnZwrhgiwp++7VGcHcg3nIhXIk6FbBTXrVxFLJNYIV8DfPrfD3HfwWE29NiAWCvNzq3N7NzaDEy9JTUaUViglDsLOPMnuLAsGzWxkVcxFPWdInKfiNwtIt+O/SYty0DvQJ7rP3uApz9pHU9/0vqVro4l5v989gCv/oP/QheR2VpVicICQWGAsDhEWBojCvNx9M+lGcxjEkAfNAv5xvCRXwssukVeMRT1WRin99tF5OYZxvmfA1eo6oSI/A7wl8BvLLZsy8J86B/vpVAMec9bHrXSVbFUsL4rxU/u6uPmbx/nBc/csvABMRoFRMEYUTBGUBg6Jwb/5H6TOi64fhZxPBwvg+O3XnBMflUlLPQiomixE/tCXz/U4peYHIoKICLloaiTQq6q363Y/zbgVTUo17IAJ05PcNMtx3j9S3ezY7PNv1lPvOg52/jczYe57uO/4BlXriebmd1mrhoRBRNGuPP9TCVVARAjqFE88lJCMxGCBGbeHScsjZnjioOonsBLdeImOhZMQDGzHvmhA3EeUAfClkVeATVhADRhY6HXgFoIedVDUWNeD3x9tg0i8ibgTQAb11l77mLZtD7Dv//9VWzfbFO31RuOI7z37Y/mxW/7Hv/wuQf5wzfum7Zdo4Cg0EuQ758m3ATtaJiOc3XOTIk1S0GlqY3q5CFxhiDfT1joR1VING3G8ZvnbaVrVCQ//JA5R6kDgtZzy66aOHl04oxp2StQXG9t7YtkWf3IReRVwBXAh2bbrqo3qOoVqnpFR2ttk9euNfIFM9p3/0Vtc7b2LCvL5Xs7+LVnbeEz/3GI0bGpEZOqIfnhA4SFfkDQwjqT+CG/w6Rl0xTnL6RiEi3nt0N+K1pqB5TSxDHyQ/eh0ewjNlWV/PCDJmZ50LZIEQe8QSR5GlC89EZzrsRp64++SGoh5FUNRRWRZwLvAa5RVRv4eAkJwohr3vwdPvzJeccQWOqAP37Tfr76iWfQnK182DpM/TUVI5y1bHN55oGQ32EeEhC3uOfYO07WLP6QGQTkDmMSv1wAaowAiaZNeMmOOLEFoHZs4mKoxdVbcCiqiFwOfBwj4vUfbKLBuemWYxw6Osa+PW0rXRXLAvR0pti60Zi+ykP3RYR0+14S2V1Mtli9IWofVlYgaoKgDZEotqXP2EMEP72OVNtevPQmTKjbfkg9EsdeyZ9fvcJmVIXi+Ak0CilNnDJBudS+gS+GRQt5lUNRPwRkgS9Xm/HCcmGUgojrP3uAfXvaePZTNqx0dSxVoKq89X0/4T1//fNp6x0vTartEkAQf8CIJ0sQpyVoQxWKY0dQnf38JpBWO+n2fSSyu/CS7UCEpE4aQXdHqqybQMGMLs4PHzCt8VJXjb7I2qUm/kMLDUVV1WfWohzLwtz4n49w7NQEn/iLR9v8pg2CiLCxJ82n/v0Qv/OKi9ixpblim2vihBd6zZD41COQ30JtXf8cKK5HkqcJC31xxMN59vbSON4mvPR6wuKwCXeb6EO1D8JmCFrmb2FrEsIWxBsxnZ2h7YxfLNYwtcJ89F8OsO+5N/P8N357ct13bzvFZc+7mf3PvZmv3Wq6G6JIuebN32Hfc2/mr//p3lnPpar8840Hecyl7Vz1RDt0upF488svJplw+bvPHADg/3zmAJc8+yae+Zu3mDjhqR5+fF8Tj3vZQzzmNd/m379nAm+pKq/409t4zG9+k/d+4hcXXoEog6pDKXeWKKyuC0vExUt2kGrbSyK7ExNQaxRJnYht6aPM2UovdRgRD9qwMrR4rEf/CjM2HpAvhBx6xNgnVZU///tfMJE3nUkHj44Am7j1J6c5fHwMAR63v3PWc4kIn//wUxkcKdrWeIPR1Z7kN39tJzf860P8zisuJl8MKQXKsVMmu72q8sF/PEShZOzRh8+cBbeT2+8u8dDxMXxPeNK+2e+LqslvhtRRCiMHSbXtrfoeEhHEy1QE1BqilDuNJHpR7TU+50HLjMQUThxsy96ntcA+CleY9tYEniuUgogoUn70815O9eYmY3D0D5rW0fWffYCJXEhHW4KnPeHc1nYYKqpKV0eKPdsXO1jDshK88Tf20JT2+PgXH6SjNYHvCWGkBGHE7Xf3c+z0OF58X/QNhkiij4/ffB+5Qkg66fGsWe6L88ODUiciSlQauaAziOPhpbpIte0jkd3BuWFvxyv3XmR9LWVsi3yFaWtJ4PsOEDGeC/jIJ++b9AEH6BsscODQMA8cHiadcvnd37wEEeHAoWH++74B7n1wCM8TNvZk+NaPTvHPH7zS+o03KO2tSf7h/U9i/0WtfPMHp/A8BxFlfCLgw/98H/nClJlicDzF8f4m7nwgRzIhvPGFG8Ad4as/7udbtw9wsrfA5p4kL3/Oeq64tNoHu4I/CEAUjOMmWi/4u5hWelNFK32QUu4MkjyDlloh6MAKee2wQr7CtGZ9XEdQz+H2u/u4/+AwqaTLlg0ZHjoySv9Qgb///AMUihEtWZ9rnrmV0fESz3/jVJ5Wz4VsU4K9u1utiDc4Vz62G4DmrIcbp+L72S/6uefAIKmky47NWe4/NMzQaIlP/scgpVBJJYQXP8fjljse4Y+vPz15rnsfHudbtw/wZ29dxzVXVSfmquClN+AmapdByrTSu3GTnWaYvz+MeiOQ3wq4NSun3hCRq4G/xXzJT6jqdTO2vxbj0Vced3O9qn4i3vYa4E/i9X+mqp+erywr5CtMa0sCEXAd4S//8V7yxYhdW5t5xpUbeOjIKMdOjXP3gUF8z+ENL91DwneIImX3tmbWdab48Z29IMLQSJF3vO7Slf46lhpw6Ogo1338XqJIcR3hQ/94H4VixJaNTTz3aRs58PAwp3tzHHxkBNcRXvVru2jv2Utrx2ngNH/yO/u4eGcL/3LzEb71w1P83RdHeMk180XNqECcJetfEXFIt+8lKAxQmjhp3CmL682I01XWOq8ymCDAv6rq22Yc2wG8FzMKXoGfxccOzlWeFfIVprU5YUbNOXDkxDippMM7XncJx0/ncB3hdG8e3xNc1+HVL9gJQCrp8o1/fibFUsTeq/8vQaA86TFdPHaxnV2WumBdV4qhkSKFYkQm7XL4+CjJhMP/99pLGBkL8Fyhb7BAwndwHOH1L96DOC7PuHITh77zwsnzFAO49SdnGB4pIk79tHy9pAnYVRg9jCRPxy6IzcYNcfWI+oLBBOfhOcAtqjoQH3sLcDUzEjRXYjs7V5jWrE8YKWGghEFEc5PPs5+yieYmb7LDM1J46XO3zRjGPZ03vGT3clXZssRkMz6vvGYHYaQUShFhpGTSHs9/+ub4vjB/2yCIeN7TNtHVkTrnHIViyHUf/wWFYsSTHze/X/hK4HhNpNr24me2MOm2mDxtOkRTR+JO0eVNUl2KlOP56iagS0TuqJjeNON0swUT3DRLsS+K8zTcKCLlUCfVHjuJFfIVprUlQakUkSuEJHyHt7/6ElxXaM76JHzz83iu8MaX7ZnzHAnf4eKdF94xZak/XvfiXQAUixG+5/CWV12M5zo0N/n4nnnA+77DW1918TnHloKIN77nxzx0ZJSWrM///r1HL2vdq0XEwU20km7fa0S9aStesg2IkOQZI+ipI+COccGxXZaOvnKAv3i64QLO8RVgu6o+CrgFmNcOPh/WtLLCZFIuIkI24+I4wouuNsmT2luSFIoh6aTL0564jo09s4f5FKEh04VZ5qejNYnrGPdDR4TfeN52wHg55Qsh6ZTL4/Z3ThsFCkbE3/wnt/HDn/WSTrl88W9/mXVd9R8SWsTB9Vtw/Ra89CaTPKM0QlAYRBImPJOqGH/0KB3HYK/rduiCwQRVtb9i8ROYhDvlY6+aceyt8xVmhXyFERG+8c/PYHQ8YH1XilTS2DKf+Ogu/uMfnk4YKnu2N59z3NhEidf+0Q8nM8HkCyHFUjTZirc0NiLCTR+7im/+1ylefPVWMmnzV718bwc3Xn/VZId3JUEY8Zb3/oTv/fQMyYTD5z/8VHZtzTZcblARwfWbcf1mvPRGNMwRBWOUcr3gDSMyPBWfPWyq18QUk8EEMcL8MuAVlTuIyAZVPRUvXoOJVQUmbtVfiEh7vPxs4N3zFWaFvA7Ytunc7D2uK1y6a25zyS+//BsMj5YQgSBUnvXab5FKutz79WvmPMbSWOzd3cbe3W3T1jmOsP+itln3v+mbx/jRf59FBArFiBe+5VbAvPXd87XGvC/Ko0YdL4OX6kE1rMiYNADu2Jzp7lYSVQ1EpBxM0AU+WQ4mCNyhqjcDvxsHFgyAAeC18bEDIvIBzMMA4P3ljs+5sELegESRksuHiBj7uIjgJl3amuuyZWJZJJ+96WGGRoq8/TcvmXe/VMolDHVy9GeZdKp+PFYWi4g72Vr30ya6p2oIzB1PfaWoIpjgu5mjpa2qnwQ+WW1ZVsgbkO/8+DTFUsRH3nMF1zyj+sS9lsbk3oeGuPnbx3jFNTvobJs7quCvPn0zv/r0tZfZ3rhsr22sQXUFUIlQt4S6RdQroF4e9XNoYsJMyXE0NYZ6s0eh+8SXHmLTugzPu2pejyTLKuENL91NoRjx+ZsPr3RVLHWKbZEvEyohJHKQGQbR6sY8KOhINxJMtcJUldf8+i4zlNq1z+G1wO5tLfzKk9bzmZsO8cbf2DPZIW6xlLFCvoSoE8TiPQLE4q3g+p2IJEEEwQHE+BHGuRklVvli/hC09KJD65HI/FQiwnOfZlvia403/MYeXvGOH/Dv3zzKK/7HjpWujqXOsEJeYxSF1JgRb9HySly/C8dtxnGqz03op3ZSyj8MbafRgY0cPTHBzd85zmteuJOWbGLhEywzqgoaxp1PkYnAhAmva14v4k+ieN05Z5jrxFNbZj2uAmHyQbg4dwaZ9mlOJedum7a+cv84YXL8gJ65PLWuOp7wqE5e/YKd7Np6riuqxWKFvNZIFJtP4mUFL7EOx2lCnIXF1whchGqJKByd2uCG/PONh/jiVw/z0uduo+Vcj8UlJwoLRME4aIBGIaoBYbEct9oIbD26gtUrsz+Tpi6g6zeZtzZxQVz+5M0bQISwOAziTq43n0sX7MpS/1ghrzGiLjq4AfwC+HlIThCUzpiNCo7XhuM2ITioluIpMIH8yy34iocAKjDSzdBAxJf/8xFe8Mwtyz5SLwqLFEYOAtGkUE+JkJiRduqCukShEpXyRGFAWBg3rWlV85U0fmNR0IjJT1VFA0U1fpBFFY338jFlqknYXm6QT7aYK9aVl52p7Y7jok5ofJYlbk07FfvP1Mf4JOXG+OS1iNe7iQzieIjr4bg+4nqx0Gr8G1dMcy07ecJSOQnD1EPy7EDALbeN8crntU2r0vSHQlwPPzun0J/7RqTTjpeKN4f53yK0omydtr7yLUyn/aCV3/nccs4tm8kbRrV8c0STjZ4F39LWAFbIlwBRF4oZKGbQ8XZwAkgUID1MFAwRhUNTO1eKQDENkTs1BYlJ2/i/3HyAfCHk9S+ZO+ZKrdEoIMj3EhTikcRBKxq2GNGe4fAUBXly/VNiPznyTpsBjygUwkKR4vAgURhNfu/yfzoKISyZzygUNIQwnBL8sjhFFXGUKhO+C1NVcspVk2hSbMUx68UBcdQsu4LjKuKWzPoK3ahEHAFH8NLNOJ4XTw7ixuJMBIQgITBOVBgDdNq5ytfDS7UgbhLHTSBeGsdLIVJdp7WifPP7R7j2M71cumUTl1+UNW+AhOZTIlMXicDJEZbGyhdvzjel2TRwqRr21ertfOUv9BazVrFCvsQIApEPeR/yWdO69IpmYyzYssCNWCiGfOY/DnHVE9dx0Y6lT+OmGhLk+wjyvfEax+RznOV2qRRwc3AnGnWgkRAWxgjyo5TGz1AWFFWhVEhTGMkTBqAhlFtfiqBOB8WiMj40RljSBa9N1d8JxfWETEsTXtLF9UAIkGCYyhaxEXhwPMHPpHA9B8dTRCYoBsOz1MbU3c924vgZXK8Dx08ijosSAgWQIlAEGSLIm47v6W82gp/pMNnp/QziJucwkwgvfNoWPvpvD/Opr5zk8ndcfl7ff7bzzb9v5Wf5baFy/bl9A9OLqWhZz1Pe/OVXljXb61GZB6s692qlJkJeRSaMJPAZ4HFAP/AbqnqkFmU3GoJAUF2HZxgq3/vpaX5yVx+b12V4/RKHqlWNCAsDJnFuuYlb2Dgjaa4hCgrk+h9iuoC3E+YnyPXfP7UeUBzCoIuxswNoqEAB/G7wm9FEkmI+pDgxRm5oEMH0C6g4aPN61E2Wm9FzTxDbaOImfBR/VqyTsEg4fpbRgfFpliskRbazGy+Vwk86oAXCcIIwd5ZSrsBsQpLuaMdLefHDoAA6RGmsl+kIiZZ1+Ok2HL8pLnDdZMlK0VwHyYMMUpronyHuDn5T/HDws5PxxDMpj994xhY+cfNhrv3M/fzS/k6e+pjuyWxCc3M+D8Q5hHfZLBjnJ/yWGgh5lZkwXg8MqupuEXkZ8EHgNxZb9momDJXX/tEPufPAALm8iXb3D597kCc+urvmAZBU1WQ+nzgx6Z2h+Y2g5z5wzhXwDtBOwmKJibMH4vUCTg9KluJEwETvCczzG0huRL0O8mMTjJ45DRqWvTJBXKKWzZBqBXfujmFVJQpDoiAkCooggut5uF4CcWY3UyhAy2Zjsw9yUBxDimOQG2S87/TUPuKQ7ewh0XQpbjKFEEKYg3ACohyU+skNVIa9cEi1b8BPN+EmvFjY8xCdpjhiJhASreuNqHsJc31ImklbQHuYFHfJAXmQIUrjvRVmKodky0bwm/nZA4Mo8LlvHOU/bj3BZbtbueFdV1Qh5pbVSi1a5NVkwngB8L54/kbgehERnd0HLcZ0bFRrP1xtfO+np7nzwAATOROHeSIXcuf9A3zvp6f5lV/aULNyomDCZGoR0/GkhXLqreloFDLRe4DJuNCxgEehMn6qYr2zmUjbKIz0kx88Eh8tkNpOSDO54UEmBh4qO+lBppso1QbJ5qnWNRAUi+RGhgmLJcIwIDcae/DEt8xckmXE2GzNNDfjeL4R+USCVDaL47rgZ8DPoE090A4alqA0jhRGYfws432nGe8DRUi3tpFoasbPdJlj09tBAwjGIByBYi/5wdPkJ5NwOaQ71uFnLsb1HYiGIDpJcfgUxeFTgEOydT1epg3HrYyNE4t7+eGp64EIJQ8yCjJAcfQ4P/hFnvsPD08eNVEIuefgMD+4s5erHlt/CSQsy0MthHy2bBYzEwRO7hNHBRsGOoG+yp3iLBtvAtjQ5ZEfurf617kavPZd8CkWOvACGkp33m1a4pXkCiF33fMgT7qob46j5ilqjg3iCBoqhWElHI+IgunBh2TaW64STPgURiDInQE5Pbm9lE8y0R+gegrlJAKEoUNRuyjllfyp00AcsVNcos49ELvXlQmKBXIjowyfPYPEgl0pzGSy4HrguKjrgeuaZRSCAMKKaXyUidFRKHvMVJyrpauLZFOWRDptWvCuD24bmmqD1i1oWIT8MDJynNzwIPnhQSYfHyKkmltwXBdxUzjOdsQF1yngkMMJh8gNnCI3UI5OKoCL4yqp1gReMkdh+CT5oZNTHb6VP5CCm0hR7qEVp+zJkcLx4f4jE+SL02+4XCHknnsf4vEbT09bb5051g511dkZZ9m4AeCSzQktDhZnF6BK97I5ttUMqZxdvlfXizd4pBPCRGHq35hKCBet9wjHg2n76nRFOJdp3ns6fTlQisNF48VV3nfSQ0wnvUpUoTQeEeRKU2bo2LMkPyqExRJRJIShEIYOhQmHMHSAYbTsRta0Hs10gWdanapKkM+TGxlhpK93UrwRQTt6INsK/uIHPmkUQTEP46Mw3M9Iby/S2ztpSmnt6SbZlMVPpUwno5uApm60qRs0QovjSGEEgjwUhsmPmA7S2e+GBI6jJFIRrhcZTxkxHailQhFxHFxfSWbA8So8SkQnvfyiIFfpjTdtzNHWTETSh3xpqsRUAnZ2upTGw+n367QKyqyzVSO1/2tZakcthHzBTBgV+xwXEQ9oZdJoOjuFUTj4n7nZN1ZzR8m0jwtmrv/CnPufV5/TXGIAGyJlV5vHQ/0lCgEkPdjd5rHhVMjDZ8amH13hmjvZz1/R8T9tuaLE6evNfFTuI4yESMWIduzVpvioGvfAKIIoktgdMFarRDM4Hrg+2pwCLw1eynRcRhFhqUSYKxKWxikVC4z1DyDlx5AI2rkOmlpqIt7TcBxIZczUuQ7CEM2Pw8Q4jA4ycuYMcKai9S80tbbgej6O5+H6Pq7XjpP2cNt3TnqUaPmJFgXTpigKCDWEoNzpGkIhHjg12UyOcCqEWuL4O+X7x7hClh9sU2K+UX12ZAs8PFKiGCoJF7ZnfXp6kxzqm1Lb2R4CZea9RWc4hsjM1kG1gj7Z8VED5jjXrO7ra5RaCPmCmTCAm4HXAD8GXgx8Z377OISSpHd4jhCt5/XOuAy/8pz1mbF+vkEYquese9Ml3dw7MM7xsTybsyn2dTZz9qzM8rSobL7N9hnvM99XmO04ccARcCV+aEj8oHBRhEgdM5gnitAoIip/lkLCUoHx4bOUnzIzSy+3hrVrPWSawVs4lnoURYRBQFgKzGdQIgjM24njuDiug+O6OK6LW7nsONM7QV3XPDCaWqB7AxqUIDcO+QljosmNMT40NOdVm3omTl3bdDaL4zqIUy7PQxxBYvOIZLphctmZGnyEgFRcI8G8mcQiLhpR7vEUAVT5nSsj7j3Tz4mRMTa1NLGvu51xmaF2WvHmtdD/Zd7tVfx/zuf/uFBLp5pzzXqOs9XXYZmowpvvncAbMIkleoHfUtVH4m0hcE+861FVnTczyKKFvMpMGP8EfFZEDmIyYbxsofNGEZwdqLARL/LNsDrO48zzNNWl4g8+GaNDzl2WScEsz89cB1vbO9jabkwcI6HGujjV/FYUjdSYANTM5yfGpte1BsbS873mk63bTBY8H/V8I9aeb+zanj/tD6mqhEFAqVCkVCwQFIuUikXGR0ZrZ+wVIZPN4iUSeAkf3zefXiKB29wGzW0zvoSiZZt7EJgRS0EQt8KjqSk/Tm5sej2X6h4tl9AOtLsejBc4OX56vkOWiWq+cbW/40Lnqv8meJXefD8HrlDVCRH5HUzOzrI3X05VH1NteTWxkVeRCSMPvOR8zhkGAYXcRC2qVx+cx3vgvK3AuY6YeVAq7kicfH+v0h5l3vMrls26STs3YkwV5UnOnVcxLfXJFnqkUy32QkAUFYmikKBYolQsMDE6NrtYx/VwslnEdcFzzafrIp43OW8ukBqTSSyu5U/CeD4M0aDExNgcZVWUl4lb15OtfMdFysteyrS4yy3tik/HEdPaJn7YloeTawRR5XLFumm2r7JtbOZ6mPbGdk5Le1onyLnrZlmcZ+X5s9wt8mrOs/Is6M2nqt+t2P824FUXWlhddXZWErk+fan2Bfer2l63EOdzT59zvlla5GI+yy1siQVy+roZ62UqlsWUHVYBmWyFaxy7pNwRqRqhUaVoxp9heT4kNz5OzalhKxkR3NYWxPcnJzzv/IJAiXnAVPVTq+kM0FKABqXJT4KAaCI3v9jXghqIUCqTpvK+k3NmphW4uMJmOXw5O/4blGq8+Sp5PfD1iuWUiNyBMbtcp6o3zVdY3Qp5qVhkpH/e/lDLQlQIhqRSVR4zy590Zo+ZVLbczzUJle3r4kzZ2aVsb4/dOIzbn3t+Yl0jjGdK3Lpn/lG25YfmtFa+KkQ6vdUdTT1gp1reVHgJVba2y5/l+RllzvAsOmcnhfzEHI4AlgumGCknclVf165YaMvcEHvdnTci8irgCuBpFau3qeoJEdkJfEdE7lHVQ3Odo26FXJJJkrt3rXQ1Lhid9odl6pV6cr78h9XJ5XLHoFYuT7bgmSac5ZasINMFsv5fORsKKUfSqrK1b1khjhxc7hL7VPWKebZX482HiDwTeA/wNFWdzO2oqifiz4dF5FbgcqDxhLzRmRTUeYTVCoPFsmpZ0JtPRC4HPg5crapnK9a3AxOqWhCRLuDJmI7QObFCbrFYLDWmSm++DwFZ4Mtxw6/sZngp8HERiTDBma+b4e1yDlbILRaLZQmowpvvmXMc9yPgsvMpa21GpLJYLJZVhBVyi8ViaXCskFssFkuDY4XcYrFYGhwr5BaLxdLgWCG3WCyWBscKucVisTQ4VsgtFoulwbFCbrFYLA2OFXKLxWJpcKyQWywWS4NjhdxisVgaHCvkFovF0uBYIbdYLJYGxwq5xWKxNDhWyC0Wi6XBWZSQi0iHiNwiIg/Fn+ekvReRx4jIj0XkXhG5W0R+YzFlWiwWSyMgIleLyAMiclBE3jXL9qSI/Gu8/Scisj1e3yki3xWRMRG5vpqyFtsifxfwbVXdA3w7Xp7JBPCbqroPuBr4GxFpW2S5FovFUreIiAt8FHgusBd4uYjsnbHb64FBVd0NfAT4YLw+D/wv4A+qLW+xQv4C4NPx/KeBX5u5g6o+qKoPxfMngbNA9yLLtVgslnrmCcBBVX1YVYvAFzF6WUmlft4IPENERFXHVfW/MIJeFYvN2blOVU/F86eBdfPtLCJPABLAoTm2vwl4U7xY+Oy3v/mLRdavFnQBfbYOQH3Uox7qAPVRj3qoA9RHPS5e7AkGRke+8dlvf7Oryt1TInJHxfINqnpDxfIm4FjF8nHgiTPOMblPnKx5GOjkAq7lgkIuIt8C1s+y6T2VC6qqIqLznGcD8FngNaoazbZPfCFuiPe/Q1WvWKh+S0091KMe6lAv9aiHOtRLPeqhDvVSjxmiekGo6tW1qMtKsKCQz5XpGUBEzojIBlU9FQv12Tn2awG+CrxHVW+74NpaLBZLY3AC2FKxvDleN9s+x0XEA1qB/gspbLE28puB18TzrwH+78wdRCQB/AfwGVW9cZHlWSwWSyNwO7BHRHbEGvgyjF5WUqmfLwa+o6pzWjXmY7FCfh3wLBF5CHhmvIyIXCEin4j3eSnwy8BrReTOeHpMFee+YeFdloV6qEc91AHqox71UAeoj3rUQx2gPupRD3WYRFUD4G3AN4D7gS+p6r0i8n4RuSbe7Z+AThE5CLyTCq8/ETkCfBijm8dn8XiZhlzgA8BisVgsdYId2WmxWCwNjhVyi8ViaXBWVMirGeIf7xdW2Ndvrli/Ix7aejAe6ppYqnrMF2pARD4lIofPsw+gfOwFDeONt707Xv+AiDzn/L951XV4p4jcF3/vb4vItopts/42S1SP14pIb0V5b6jY9pr493tIRF4z89ga1uEjFeU/KCJDFdtqci1E5JMiclZEZh1HIYa/i+t4t4g8tmJbTa5DlfV4ZVz+PSLyIxF5dMW2I/H6O2URroFV1OEqERmuuO5/WrFt3t9yVaGqKzYBfwm8K55/F/DBOfYbm2P9l4CXxfMfA35nqeoBXATsiec3AqeAtnj5U8CLL6BcFzM4aidmoNRdwN4Z+7wF+Fg8/zLgX+P5vfH+SWBHfB53ierwdCATz/9OuQ7z/TZLVI/XAtfPcmwH8HD82R7Pty9FHWbs/3bgk0twLX4ZeCzwizm2Pw/4OiDAk4Cf1PI6nEc9riyfHzMU/ScV244AXctwLa4C/t9if8tGn1batLLgEP+5EBEBfgUztPW8jz/feujShBq44GG88fovqmpBVQ8DB+Pz1bwOqvpdVZ2IF2/D+MTWmmquxVw8B7hFVQdUdRC4BRPXZ6nr8HLgCxdQzryo6veBgXl2eQHGnVfVjMtoEzOOo1bXoap6qOqP4nJgie6LKq7FXCzmfmo4VlrIqx3inxKRO0TkNhH5tXhdJzCkxs0HzBDYTUtcD2DOUAN/Hr9mfkREklWWO9sw3pnfYdowXqA8jLeaY2tVh0pej2kNlpntt7kQqq3Hi+LrfKOIlAdcLPu1iM1LO4DvVKyu1bVYiLnqWavrcCHMvC8U+KaI/ExM6I2l5JdE5C4R+bqI7IvXreS1WHYWG2tlQaQ2Q/y3qeoJEdkJfEdE7sEI2nLXY65QA+/GPAASGH/WPwbefz71awRE5FXAFcDTKlaf89uo6qyxdGrAV4AvqGpBRN6MeVP5lSUqayFeBtyoqmHFuuW8FnWDiDwdI+RPqVj9lPha9AC3iMiBuHVda/4bc93HROR5wE3AniUop65Z8ha5qj5TVffPMv1f4EwsjGWBnHWIv6qeiD8fBm4FLscMZW0TM7QVZh8CW9N6yByhBlT1VPyaWwD+mepNHOczjBeZPoy3mmNrVQdE5JmYh9418fcE5vxtLoQF66Gq/RVlfwJ43Pl8h1rUoYKXMcOsUsNrsRBz1bNW16FqRORRmN/iBao6Oby84lqcxYzsvhCz34Ko6oiqjsXzXwN8EeliBa7FirKSBnrgQ0zvZPzLWfZpB5LxfBfwEHGnBfBlpnd2vmUJ65HAxFz//2bZtiH+FOBvgOuqLNfDdEjtYKpDZt+Mfd7K9M7OL8Xz+5je2fkwF9bZWU0dLseYkfZU+9ssUT02VMy/ELgtnu8ADsf1aY/nO5aiDvF+l2A682QprkV8ju3M3cH3fKZ3dv60ltfhPOqxFdM3c+WM9U1Ac8X8j4Crl6gO68u/A+ZhcTS+LlX9lqtlWtnCja332/FN/63yTYd5ff9EPH8lcE/8Q9wDvL7i+J3AT+Ob6cvlP9IS1eNVQAm4s2J6TLztO3HdfgH8C5A9j7KfBzyIEcr3xOvej2n5AqTi73Yw/q47K459T3zcA8BzF/E7LFSHbwFnKr73zQv9NktUj2uBe+PyvgtcUnHsb8XX6CDwuqWqQ7z8PmY8rGt5LTAt/VPx/XYcY7b4beC34+2CSVpwKC7rilpfhyrr8QlgsOK+uKPif3lXPN1bvo5LVIe3VdwTt1HxUJntt1ytkx2ib7FYLA3OSnutWCwWi2WRWCG3WCyWBscKucVisTQ4VsgtFoulwbFCbrFYLA2OFXKLxWJpcKyQWywWS4Pz/wObBbENv+Z1BwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "levels = [0.01, 0.25, 0.5, 0.75, 0.85, 0.90, 0.95, 0.98, 0.99,  1.,]\n",
    "\n",
    "xs = np.linspace(1, 250, len(levels))\n",
    "cmap_colors = [cmocean.cm.solar(int(x)) for x in xs]\n",
    "\n",
    "plt.contour(grid_us, grid_vs, accs, zorder=1, levels=levels,\n",
    "            colors=[cmap_colors[i] for i in range(len(levels))])\n",
    "plt.contourf(grid_us, grid_vs, accs,  zorder=0, alpha=0.8, levels=levels,\n",
    "             colors=[cmap_colors[i] for i in range(len(levels))])\n",
    "plt.colorbar()\n",
    "\n",
    "w1 = np.array([0., 0.])\n",
    "w2 = np.array([1., 0.])\n",
    "th = np.array([0.5, 1.])\n",
    "curve_xy = np.stack([w1 * t**2 + w2 * (1 - t)**2 + 2 * th * t * (1 - t) for t in np.linspace(0, 1, 20)])\n",
    "\n",
    "plt.plot([w1[0], w2[0], th[0]], [w1[1], w2[1], th[1]], 'o', ms=5, mec='black', mfc= 'black', mfcalt='black')\n",
    "plt.plot([w1[0]], [w1[1] + 0.12], marker='$w_1$', mec='black', mfc='black', mfcalt='black',  ms=20)\n",
    "plt.plot([w2[0]], [w2[1] + 0.12], marker='$w_2$', mec='black', mfc='black', mfcalt='black',  ms=20)\n",
    "plt.plot([th[0]], [th[1] + 0.12], marker=r'$\\theta$', mec='black', mfc='black', mfcalt='black',  ms=12)\n",
    "\n",
    "plt.plot(curve_xy[:, 0], curve_xy[:, 1], \"--k\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 2)"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "curve_xy.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py38",
   "language": "python",
   "name": "py38"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}