Semi-supervise learning with (250,125) adversarial examples with spatial feature.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%reload_ext autoreload\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# always import gbm_algos first !\n",
    "import xgboost, lightgbm, catboost\n",
    "\n",
    "# Core\n",
    "import numpy as np\n",
    "import scipy.io as io\n",
    "from torch.autograd import grad\n",
    "from torch.utils.data import DataLoader, Dataset\n",
    "\n",
    "# Sklearn\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from mlens.ensemble import SuperLearner\n",
    "\n",
    "# Let's do facy optimizers\n",
    "from optimizers import Lookahead, AdamGC, SGDGC\n",
    "from madgrad import MADGRAD\n",
    "from lbfgsnew import LBFGSNew\n",
    "# Modify at /usr/local/lib/python3.9/site-packages/torch_lr_finder/lr_finder.py\n",
    "from torch_lr_finder import LRFinder\n",
    "from onecyclelr import OneCycleLR\n",
    "import pcgrad\n",
    "from pytorch_stats_loss import torch_wasserstein_loss, torch_energy_loss\n",
    "from geomloss import SamplesLoss\n",
    "from utils import *\n",
    "\n",
    "from tqdm import trange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training with 250 samples\n",
      "Training with 125 unsup samples\n"
     ]
    }
   ],
   "source": [
    "DATA_PATH = \"./experimental_data/burgers_shock.mat\"\n",
    "data = io.loadmat(DATA_PATH)\n",
    "\n",
    "t = data['t'].flatten()[:,None]\n",
    "x = data['x'].flatten()[:,None]\n",
    "Exact = np.real(data['usol']).T\n",
    "\n",
    "X, T = np.meshgrid(x,t)\n",
    "\n",
    "X_star = np.hstack((X.flatten()[:,None], T.flatten()[:,None]))\n",
    "u_star = Exact.flatten()[:,None]              \n",
    "\n",
    "# Doman bounds\n",
    "lb = X_star.min(0)\n",
    "ub = X_star.max(0)\n",
    "\n",
    "N = 250\n",
    "print(f\"Training with {N} samples\")\n",
    "idx = np.random.choice(X_star.shape[0], N, replace=False)\n",
    "X_u_train = X_star[idx, :]\n",
    "u_train = u_star[idx,:]\n",
    "\n",
    "# Unsup data\n",
    "N_res = N//2\n",
    "idx_res = np.array(range(X_star.shape[0]-1))[~idx]\n",
    "idx_res = np.random.choice(idx_res.shape[0], N_res, replace=True)\n",
    "X_res = X_star[idx_res, :]\n",
    "print(f\"Training with {N_res} unsup samples\")\n",
    "X_u_train = np.vstack([X_u_train, X_res])\n",
    "u_train = np.vstack([u_train, torch.rand(X_res.shape[0], 1) - 1000])\n",
    "# del X_res\n",
    "\n",
    "# Convert to torch.tensor\n",
    "X_u_train = torch.tensor(X_u_train).float().requires_grad_(True)\n",
    "u_train = torch.tensor(u_train).float().requires_grad_(True)\n",
    "X_star = torch.tensor(X_star).float().requires_grad_(True)\n",
    "u_star = torch.tensor(u_star).float().requires_grad_(True)\n",
    "# lb and ub are used in adversarial training\n",
    "scaling_factor = 1.0\n",
    "lb = scaling_factor*to_tensor(lb, False)\n",
    "ub = scaling_factor*to_tensor(ub, False)\n",
    "\n",
    "feature_names=('uf', 'u_x', 'u_xx', 'u_xxx', 'x')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Network(nn.Module):\n",
    "    def __init__(self, model, index2features=('uf', 'u_x',  'u_xx', 'u_tt', 'u_xt', 'u_tx')):\n",
    "        super(Network, self).__init__()\n",
    "        # pls init the self.model before\n",
    "        self.model = model\n",
    "        # For tracking, the default tup is for the burgers' equation.\n",
    "        self.index2features = index2features\n",
    "        print(\"Considering\", self.index2features)\n",
    "        self.diff_flag = diff_flag(self.index2features)\n",
    "        self.uf = None\n",
    "        \n",
    "    def xavier_init(self, m):\n",
    "        if type(m) == nn.Linear:\n",
    "            torch.nn.init.xavier_uniform_(m.weight)\n",
    "            m.bias.data.fill_(0.01)\n",
    "        \n",
    "    def forward(self, x, t):\n",
    "        self.uf = self.model(torch.cat([x, t], dim=1))\n",
    "        return self.uf\n",
    "    \n",
    "    def get_selector_data(self, x, t):\n",
    "        uf = self.forward(x, t)\n",
    "        u_t = self.gradients(uf, t)[0]\n",
    "        \n",
    "        ### PDE Loss calculation ###\n",
    "        # Without calling grad\n",
    "        derivatives = []\n",
    "        for t in self.diff_flag[0]:\n",
    "            if t=='uf': derivatives.append(uf)\n",
    "            elif t=='x': derivatives.append(x)\n",
    "        # With calling grad\n",
    "        for t in self.diff_flag[1]:\n",
    "            out = uf\n",
    "            for c in t:\n",
    "                if c=='x': out = self.gradients(out, x)[0]\n",
    "                elif c=='t': out = self.gradients(out, t)[0]\n",
    "            derivatives.append(out)\n",
    "        \n",
    "        return torch.cat(derivatives, dim=1), u_t\n",
    "    \n",
    "    def gradients(self, func, x):\n",
    "        return grad(func, x, create_graph=True, retain_graph=True, grad_outputs=torch.ones(func.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SeclectorNetwork(nn.Module):\n",
    "    def __init__(self, X_train_dim, bn=None):\n",
    "        super().__init__()\n",
    "        # Nonlinear model, Training with PDE reg.\n",
    "        self.nonlinear_model = TorchMLP(dimensions=[X_train_dim, 50, 50, 1], activation_function=nn.Tanh, bn=bn, dropout=nn.Dropout(p=0.1), inp_drop=False)\n",
    "        \n",
    "    def xavier_init(self, m):\n",
    "        if type(m) == nn.Linear:\n",
    "            torch.nn.init.xavier_uniform_(m.weight)\n",
    "            m.bias.data.fill_(0.01)\n",
    "        \n",
    "    def forward(self, inn):\n",
    "        ut_approx = self.nonlinear_model(inn)\n",
    "        return ut_approx\n",
    "    \n",
    "    def loss(self, X_input, y_input):\n",
    "        ut_approx = self.forward(X_input)\n",
    "        mse_loss = F.mse_loss(ut_approx, y_input, reduction='mean')\n",
    "        return mse_loss\n",
    "\n",
    "class AttentionSelectorNetwork(nn.Module):\n",
    "    def __init__(self, layers, bn=None, reg_intensity=0.3):\n",
    "        super(AttentionSelectorNetwork, self).__init__()\n",
    "        # Nonlinear model, Training with PDE reg.\n",
    "        assert len(layers) > 1\n",
    "        self.linear1 = nn.Linear(layers[0], layers[0])\n",
    "        self.nonlinear_model = TorchMLP(dimensions=layers, activation_function=nn.Tanh, bn=bn, dropout=nn.Dropout(p=0.1))\n",
    "        self.latest_weighted_features = None\n",
    "        self.th = 0.5\n",
    "        self.reg_intensity = reg_intensity\n",
    "        \n",
    "    def xavier_init(self, m):\n",
    "        if type(m) == nn.Linear:\n",
    "            torch.nn.init.xavier_uniform_(m.weight)\n",
    "            m.bias.data.fill_(0.01)\n",
    "        \n",
    "    def forward(self, inn):\n",
    "        return self.nonlinear_model(inn*self.weighted_features(inn))\n",
    "    \n",
    "    def weighted_features(self, inn):\n",
    "        self.latest_weighted_features =  torch.sigmoid(self.linear1(inn)).mean(axis=0)\n",
    "        return self.latest_weighted_features\n",
    "    \n",
    "    def loss(self, X_input, y_input):\n",
    "        ut_approx = self.forward(X_input)\n",
    "        mse_loss = F.mse_loss(ut_approx, y_input, reduction='mean')\n",
    "        return mse_loss+self.reg_intensity*torch.norm(F.relu(self.latest_weighted_features-self.th), p=0)\n",
    "\n",
    "class SemiSupModel(nn.Module):\n",
    "    def __init__(self, network, selector, normalize_derivative_features=False, mini=None, maxi=None):\n",
    "        super(SemiSupModel, self).__init__()\n",
    "        self.network = network\n",
    "        self.selector = selector\n",
    "        self.normalize_derivative_features = normalize_derivative_features\n",
    "        self.mini = mini\n",
    "        self.maxi = maxi\n",
    "    def forward(self, X_u_train):\n",
    "        X_selector, y_selector = self.network.get_selector_data(*dimension_slicing(X_u_train))\n",
    "        if self.normalize_derivative_features:\n",
    "            X_selector = (X_selector-self.mini)/(self.maxi-self.mini)\n",
    "        unsup_loss = self.selector.loss(X_selector, y_selector)\n",
    "        return self.network.uf, unsup_loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using old implementation of TorchMLP. See models.py for more new model-related source code.\n",
      "Considering ('uf', 'u_x', 'u_xx', 'u_xxx', 'x')\n",
      "Using old implementation of TorchMLP. See models.py for more new model-related source code.\n"
     ]
    }
   ],
   "source": [
    "# network = Network(model=TorchMLP(dimensions=[6, 50, 50, 50 ,50, 50, 1], bn=nn.BatchNorm1d))\n",
    "# selector = SeclectorNetwork(X_train_dim=6, bn=nn.LayerNorm)\n",
    "\n",
    "### Version without normalized derivatives ###\n",
    "# semisup_model = SemiSupModel(network=Network(model=TorchMLP(dimensions=[2, 50, 50, 50 ,50, 50, 1], activation_function=nn.Tanh, bn=nn.LayerNorm, dropout=None)),\n",
    "#                              selector=SeclectorNetwork(X_train_dim=6, bn=nn.LayerNorm),\n",
    "#                              normalize_derivative_features=False, \n",
    "#                              mini=None, \n",
    "#                              maxi=None)\n",
    "\n",
    "### Version with normalized derivatives ###\n",
    "pretrained_state_dict = torch.load('./saved_path_inverse_burger/lbfgsnew_results/semisup_model_with_LayerNormDropout_without_physical_reg_trained250labeledsamples_trained0unlabeledsamples_2.2e-03.pth')\n",
    "network_state_dict = None\n",
    "use_pretrained_weights = False\n",
    "lets_pretrain = True\n",
    "\n",
    "# semisup_model = SemiSupModel(network=Network(\n",
    "#                                     model=TorchMLP(dimensions=[2, 50, 50, 50 ,50, 50, 1],\n",
    "#                                                    activation_function=nn.Tanh,\n",
    "#                                                    bn=nn.LayerNorm, dropout=None),\n",
    "#                                     index2features=feature_names),\n",
    "#                             selector=SeclectorNetwork(X_train_dim=len(feature_names), bn=nn.LayerNorm),\n",
    "#                             normalize_derivative_features=True,\n",
    "#                             mini=None,\n",
    "#                             maxi=None)\n",
    "\n",
    "semisup_model = SemiSupModel(network=Network(\n",
    "                                    model=TorchMLP(dimensions=[2, 50, 50, 50 ,50, 50, 1],\n",
    "                                                   activation_function=nn.Tanh,\n",
    "                                                   bn=nn.LayerNorm, dropout=None),\n",
    "                                    index2features=feature_names),\n",
    "                            selector=AttentionSelectorNetwork([len(feature_names), 50, 50, 1], bn=nn.LayerNorm),\n",
    "                            normalize_derivative_features=True,\n",
    "                            mini=None,\n",
    "                            maxi=None)\n",
    "\n",
    "if use_pretrained_weights:\n",
    "    print(\"Use pretrained weights\")\n",
    "    semisup_model.load_state_dict(pretrained_state_dict)\n",
    "    network_state_dict = semisup_model.network.state_dict()\n",
    "    semisup_model.eval()\n",
    "    referenced_derivatives, _ = semisup_model.network.get_selector_data(*dimension_slicing(X_star))\n",
    "    semisup_model.mini = torch.min(referenced_derivatives, axis=0)[0].detach().requires_grad_(False)\n",
    "    semisup_model.maxi = torch.max(referenced_derivatives, axis=0)[0].detach().requires_grad_(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### You can pretrain the semisup_model if pretrained weights are not provided"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pretraining\n",
      "Epoch 0:  4.308073130232515e-06\n",
      "Test MSE: 9.5e-04\n",
      "Epoch 100:  3.050152486139268e-07\n",
      "Test MSE: 7.0e-04\n",
      "Computing derivatives features\n"
     ]
    }
   ],
   "source": [
    "if lets_pretrain:\n",
    "    print(\"Pretraining\")\n",
    "    pretraining_optimizer = LBFGSNew(semisup_model.network.parameters(), \n",
    "                                     lr=1e-1, max_iter=300, \n",
    "                                     max_eval=int(300*1.25), history_size=150, \n",
    "                                     line_search_fn=True, batch_mode=False)\n",
    "\n",
    "    semisup_model.train()\n",
    "    for i in range(200):\n",
    "        def pretraining_closure():\n",
    "            global N, X_u_train, u_train\n",
    "            if torch.is_grad_enabled():\n",
    "                pretraining_optimizer.zero_grad()\n",
    "            # Only focusing on first [:N, :] elements\n",
    "            mse_loss = F.mse_loss(semisup_model.network(*dimension_slicing(X_u_train))[:N, :], u_train[:N, :])\n",
    "            if mse_loss.requires_grad:\n",
    "                mse_loss.backward(retain_graph=False)\n",
    "            return mse_loss\n",
    "\n",
    "        pretraining_optimizer.step(pretraining_closure)\n",
    "\n",
    "        l = pretraining_closure()\n",
    "        if (i % 100) == 0:\n",
    "            curr_loss = l.item()\n",
    "            print(\"Epoch {}: \".format(i), curr_loss)\n",
    "\n",
    "            # Sneak on the test performance...\n",
    "            semisup_model.network.eval()\n",
    "            test_performance = F.mse_loss(semisup_model.network(*dimension_slicing(X_star)).detach(), u_star).item()\n",
    "            string_test_performance = scientific2string(test_performance)\n",
    "            print('Test MSE:', string_test_performance)\n",
    "    \n",
    "    print(\"Computing derivatives features\")\n",
    "    semisup_model.eval()\n",
    "    referenced_derivatives, _ = semisup_model.network.get_selector_data(*dimension_slicing(X_star))\n",
    "    semisup_model.mini = torch.min(referenced_derivatives, axis=0)[0].detach().requires_grad_(False)\n",
    "    semisup_model.maxi = torch.max(referenced_derivatives, axis=0)[0].detach().requires_grad_(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pcgrad_closure():\n",
    "    global N, X_u_train, u_train\n",
    "    uf, unsup_loss = semisup_model(X_u_train)\n",
    "    losses = [F.mse_loss(uf[:N, :], u_train), unsup_loss]\n",
    "    updated_grads = []\n",
    "    \n",
    "    for i in range(2):\n",
    "        optimizer1.zero_grad()\n",
    "        losses[i].backward(retain_graph=True)\n",
    "\n",
    "        g_task = []\n",
    "        for param in semisup_model.parameters():\n",
    "            if param.grad is not None:\n",
    "                g_task.append(Variable(param.grad.clone(), requires_grad=False))\n",
    "            else:\n",
    "                g_task.append(Variable(torch.zeros(param.shape), requires_grad=False))\n",
    "        # appending the gradients from each task\n",
    "        updated_grads.append(g_task)\n",
    "\n",
    "    updated_grads = list(pcgrad.pc_grad_update(updated_grads))[0]\n",
    "    for idx, param in enumerate(semisup_model.parameters()):\n",
    "        param.grad = (updated_grads[0][idx]+updated_grads[1][idx])\n",
    "        \n",
    "    return sum(losses)\n",
    "\n",
    "def closure():\n",
    "    global N, X_u_train, u_train\n",
    "    if torch.is_grad_enabled():\n",
    "        optimizer2.zero_grad()\n",
    "    mse_loss = F.mse_loss(semisup_model.network(*dimension_slicing(X_u_train))[:N, :], u_train)\n",
    "    if mse_loss.requires_grad:\n",
    "        mse_loss.backward(retain_graph=True)\n",
    "    return mse_loss\n",
    "\n",
    "def selector_closure():\n",
    "    global ders, dys\n",
    "    if torch.is_grad_enabled():\n",
    "        optimizer3.zero_grad()\n",
    "    unsup_loss = semisup_model.selector.loss(ders, dys)\n",
    "    if unsup_loss.requires_grad:\n",
    "        unsup_loss.backward(retain_graph=True)\n",
    "    return unsup_loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Learning rate finding\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "13bcc90eb3854f24af00cf431c1c3e08",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=300.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stopping early, the loss has diverged\n",
      "\n",
      "Learning rate search finished. See the graph with {finder_name}.plot()\n",
      "LR suggestion: steepest gradient\n",
      "Suggested LR: 4.29E-07\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAt1UlEQVR4nO3deXhU5d3/8fc3ISSBrEBYA4SAIPsWdlm0gopV27pXbV1xqVofrU9r26dqW39P+9jaurRu1eK+oVZEFFuVXYEAguxrgLBIgJCF7Mn9+2MGBExCtslJZj6v6zpXZs6cM+d7MySfOec+5z7mnENEREJXmNcFiIiItxQEIiIhTkEgIhLiFAQiIiFOQSAiEuIUBCIiIa6F1wXUVrt27VxKSorXZYiINCvLly8/4JxLquy1ZhcEKSkppKene12GiEizYmY7qnpNh4ZEREKcgkBEJMQpCEREQlyz6yMQkdorLS0lMzOToqIir0uRAIuKiiI5OZmIiIgar6MgEAkBmZmZxMbGkpKSgpl5XY4EiHOOgwcPkpmZSY8ePWq8ng4NiYSAoqIi2rZtqxAIcmZG27Zta73nFzJBUFxWztvLM9Gw2xKqFAKhoS6fc8gEwb9W7uaet1bxzordXpci0vQ5B198Ae++6/sZoC9Qf/3rXykoKAjIe9fU4cOH+fvf/95o20tJSeHAgQMAjB07ts7vM336dPbs2dMgNYVMEFwyvCsjU9rwwMy17D5c6HU5Ik3X7NnQrRtMngzXXuv72a2bb34DC5YgKCsrq9N6ixcvrvM2FQR1EB5m/PmywVQ4x22vrKCgpG4fnEhQmz0bLrkEMjMhPx9yc30/MzN98+sYBkeOHOH8889n8ODBDBgwgDfeeIPHHnuMPXv2cOaZZ3LmmWcC8PHHHzNmzBiGDRvGpZdeSn5+PgDLly9n4sSJDB8+nHPOOYe9e/cCMGnSJH76058yZMgQBgwYwNKlS49t7/rrr2fkyJEMHTqU9957D4C1a9cycuRIhgwZwqBBg9i8eTO/+MUv2Lp1K0OGDOHee+/9Vu2/+93v6NOnD2eccQZXXnklf/rTn45t+6677iItLY1HH32U999/n1GjRjF06FDOPvtsvv76awAOHjzIlClT6N+/PzfeeOMJh6djYmKOPX744YcZMWIEgwYN4v777wcgIyODvn37ctNNN9G/f3+mTJlCYWEhM2bMID09nauuuoohQ4ZQWFjPL7fOuWY1DR8+3NXHnDV7XY9fzHI/em6JKykrr9d7iTQX69atO/VCFRXOdeninO9AUOVTcrJvuVqaMWOGu/HGG489P3z4sHPOue7du7usrCznnHNZWVlu/PjxLj8/3znn3B/+8Af34IMPupKSEjdmzBi3f/9+55xzr7/+urvuuuucc85NnDjx2PvOmzfP9e/f3znn3H333edeeukl55xz2dnZ7rTTTnP5+fnu9ttvdy+//LJzzrni4mJXUFDgtm/ffmy9ky1dutQNHjzYFRYWutzcXNerVy/38MMPH9v2rbfeemzZQ4cOuQr/v82zzz7r7r77buecc3fccYd78MEHnXPOzZo1ywHH2ty6dWvnnHNz5sxxN910k6uoqHDl5eXu/PPPd/PmzXPbt2934eHhbuXKlc455y699NJj7Zo4caJbtmxZpXVX9nkD6a6Kv6shd/rolP4d+f33BvLLd7/iHwu2c+uknl6XJNI0LFkCOTnVL3P4MCxdCqNG1eqtBw4cyD333MPPf/5zvvvd7zJ+/PhvLfPFF1+wbt06xo0bB0BJSQljxoxh48aNrFmzhsmTJwNQXl5Op06djq135ZVXAjBhwgRyc3M5fPgwH3/8MTNnzjz27b2oqIidO3cyZswYHnroITIzM/nBD37AaaedVm3dixYt4qKLLiIqKoqoqCguuOCCE16//PLLjz3OzMzk8ssvZ+/evZSUlBw7fXP+/Pm88847AJx//vkkJiZ+azsff/wxH3/8MUOHDgUgPz+fzZs3061bN3r06MGQIUMAGD58OBkZGdXWXBchFwQAPxzVjXmb9vPX/2zi3AEd6dGutdcliXhv714IO8XR4rAwqMNx6d69e7NixQpmz57Nr3/9a77zne/wm9/85oRlnHNMnjyZ11577YT5X331Ff379+fzzz+v9L1PPkvGzHDO8fbbb9OnT58TXuvbty+jRo3igw8+YOrUqTz99NOkpqbWuj1HtW79zd+OO+64g7vvvpsLL7yQuXPn8sADD9T4fZxz3Hfffdx8880nzM/IyCAyMvLY8/Dw8PofBqpEyPQRnOy3Fw2gZYswfvv+Wq9LEWkaOnWCiorql6mogM6da/3We/bsoVWrVlx99dXce++9rFixAoDY2Fjy8vIAGD16NIsWLWLLli2A7zj/pk2b6NOnD1lZWceCoLS0lLVrv/m9feONNwBYuHAh8fHxxMfHc8455/D4448fOx6/cuVKALZt20Zqaip33nknF110EatXrz6hhpONGzeO999/n6KiIvLz85k1a1aVbczJyaFLly4AvPDCC8fmT5gwgVdffRWADz/8kOzs7G+te8455/D8888f6xPZvXs3+/fvr/bftLq6aysk9wgAOsRFccvEnjw8ZyNrducwoEu81yWJeGvUKIiP93UOVyUhAUaOrPVbf/XVV9x7772EhYURERHBk08+CcC0adM499xz6dy5M5999hnTp0/nyiuvpLi4GIDf//739O7dmxkzZnDnnXeSk5NDWVkZd911F/379wd8QyoMHTqU0tJSnn/+eQD+53/+h7vuuotBgwZRUVFBjx49mDVrFm+++SYvvfQSERERdOzYkV/+8pe0adOGcePGMWDAAM477zwefvjhY3WPGDGCCy+8kEGDBtGhQwcGDhxIfHzlfyseeOABLr30UhITEznrrLPYvn07APfffz9XXnkl/fv3Z+zYsXTr1u1b606ZMoX169czZswYwNeJ/PLLLxMeHl7lv+m1117LLbfcQnR0NJ9//jnR0dG1/Vi+UVXnQVOd6ttZfLycwhI34DcfudteXt5g7ynSFNWos9g55z74wLno6Mo7iqOjfa83IdV1mDaUvLw855xzR44cccOHD3fLlzf9vxe17SwO2UNDAHFREVwzpjuz1+xlW1Y134JEQsXUqTBjBiQnQ0wMxMX5fiYn++ZPnep1hY1u2rRpDBkyhGHDhnHxxRczbNgwr0tqcAE7NGRmUcB8INK/nRnOuftPWiYSeBEYDhwELnfOZQSqpspcf0YPnlu4nafmbeX/LhncmJsWaZqmToWdO31nB+3Z4+sTGDkSmuAQFXPnzg34No4e3w9mgdwjKAbOcs4NBoYA55rZ6JOWuQHIds71Av4C/DGA9VSqXUwkV4zoyrsrd7NHVxyL+Jj5+gy+/33fzyYYAtJwAhYE/sNSR4+3RPinkwcsuQg42r0+A/iOeTAy1k0TUnEOnl2wrbE3LdJonAZcDAl1+ZwD2kdgZuFm9iWwH/i3c27JSYt0AXYBOOfKgBygbSBrqkxyYisuGtKF15fu4mB+cWNvXiTgoqKiOHjwoMIgyDn//QiioqJqtV5ATx91zpUDQ8wsAXjXzAY459bU9n3MbBowDaj01KuGcOukVN5Zmcn0xRncM6XPqVcQaUaSk5PJzMwkKyvL61IkwI7eoaw2GuU6AufcYTP7DDgXOD4IdgNdgUwzawHE4+s0Pnn9Z4BnANLS0gLylaZX+1jO6deRFxZnMG1CKrFRNb/Nm0hTFxERUas7VkloCdihITNL8u8JYGbRwGRgw0mLzQR+7H98CfCp83Df9bYze5JbVMYrS3Z6VYKISKMLZB9BJ+AzM1sNLMPXRzDLzH5rZhf6l3kOaGtmW4C7gV8EsJ5TGpScwPjT2vGPBdspKi33shQRkUYTsENDzrnVwNBK5v/muMdFwKWBqqEubp3Ukx8+u4S3lmdyzejuXpcjIhJwIX1lcWXGpLZlaLcEnp63lbLyUwzAJSISBBQEJzEzbpvUi8zsQt5f3TC3gRMRacoUBJX4zunt6dMhlifnbqWiQuddi0hwUxBUIizMuHVSTzZ9nc/cTdWPCS4i0twpCKpw/qBOtI+N5KXPd3hdiohIQCkIqhARHsaVI7sxd1MWOw8WeF2OiEjAKAiqceXIboSZ8cpS7RWISPBSEFSjY3wUU/p14M1lu3SBmYgELQXBKVwzujvZBaXM/mqv16WIiASEguAUxvRsS2pSa176QoeHRCQ4KQhOwcy4ZnR3Vu7IZuvM/8C778IXX/hu5y0iEgQaZRjq5u6y/V9x7pM3kPhYAbRsARUVkJAATz8dkjfzFpHgoj2CU5k9m9ZXXUGnvANEFRVAbi7k50NmJlxyCcye7XWFIiL1oiCojnMwbRoUVnFT+8JCuPlmHSYSkWZNQVCdJUsgJ6f6ZQ4fhqVLG6UcEZFAUBBUZ+9eCDvFP1FYGOzRKKUi0nwpCKrTqZOvY7g6FRXQuXPj1CMiEgAKguqMGgXx8dUvk5AAI0c2SjkiIoGgIKiOGTzzDERHV/56dLTvFFKzxq1LRKQBKQhOZepUmDEDkpMhJgbi4iiOasXe2Hbkv/yariMQkWZPF5TVxNSpsHOn7+ygPXvY0zKeM+cX8Mt2fZnmdW0iIvWkIKgpM1+fAdADGJn5OS9+voMbzkglPEyHhkSk+dKhoTq6dmwKmdmFfLpBt7IUkeZNQVBHU/p1oFN8FC8szvC6FBGRelEQ1FGL8DCuHt2dhVsOsGV/ntfliIjUmYKgHq4Y0ZWWLcJ4YbHuVSAizZeCoB7axkRywaDOvL0ik9yiUq/LERGpEwVBPV07NoWCknJmpGd6XYqISJ0oCOppYHI8w7ol8OLnGVRUaDhqEWl+FAQN4MdjU8g4WMC8zVlelyIiUmsKggZw3oBOJMVG6lRSEWmWFAQNoGWLMK4a1Y25G7PYfuCI1+WIiNRKwILAzLqa2Wdmts7M1prZTytZZpKZ5ZjZl/7pN4GqJ9B+OKobLcPDmL5ou9eliIjUSiD3CMqAe5xz/YDRwE/MrF8lyy1wzg3xT78NYD0B1T42iouGdObN9Eyyj5R4XY6ISI0FLAicc3udcyv8j/OA9UCXQG2vKbhpQiqFpeW8skQXmIlI89EofQRmlgIMBZZU8vIYM1tlZh+aWf8q1p9mZulmlp6V1XTPzOndIZZJfZKYvngHRaXlXpcjIlIjAQ8CM4sB3gbucs7lnvTyCqC7c24w8Djwr8rewzn3jHMuzTmXlpSUFNB662va+FQO5Bfz3pe7vS5FRKRGAhoEZhaBLwRecc69c/Lrzrlc51y+//FsIMLM2gWypkAb07Mt/TrF8eyC7brATESahUCeNWTAc8B659wjVSzT0b8cZjbSX8/BQNXUGMyMaRNS2bI/n3mbmu5hLBGRowK5RzAOuAY467jTQ6ea2S1mdot/mUuANWa2CngMuMI51+y/Rp8/qBOd4qN4Zv42r0sRETmlgN2q0jm3EKj2Ho7OuSeAJwJVg1ciwsO4flwPHpq9nq8ycxiYHO91SSIiVdKVxQFyxciuxEa24NkF2isQkaZNQRAgsVERXDGyKx98tZe9OYVelyMiUiUFQQD9aEwKzjle+WKn16WIiFRJQRBAXdu04jt9O/Dq0p26wExEmiwFQYBdNzaFQ0dKmLV6r9eliIhUSkEQYGN6tqV3hxj+uWg7QXBmrIgEIQVBgJkZPx6bwto9uSzfke11OSIi36IgaATfH9qFuKgWTNcdzESkCVIQNIJWLVtwxchufLhmH7sOFXhdjojICRQEjeSGM3oQbsaT87Z6XYqIyAkUBI2kQ1wUl6QlMyM9k305RV6XIyJyjIKgEd06sSflzvHcQg07ISJNh4KgEXVt04qpAzvx+tJd5BWVel2OiAigIGh0N43vQV5xGW8s2+V1KSIigIKg0Q1KTmBkjzb8c1EGZeUVXpcjIqIg8MJN41PZfbiQ2Wv2eV2KiIiCwAvfOb09qe1a848F2zTshIh4TkHggbAw4/ozerA6M4el2w95XY6IhDgFgUcuHpZMYqsInl2w3etSRCTEKQg8Et0ynGtGd+eTDV+zLSvf63JEJIQpCDx0zZgUIsLDeG6h9gpExDsKAg8lxUby/SFdmLE8kwP5xV6XIyIhSkHgsWkTUykpr+Cfi7RXICLeUBB4rGdSDFMHdOLFxTvIKdSwEyLS+BQETcCtk3qSV1zGy1/s8LoUEQlBCoImYECXeCb1SeK5hdspLCn3uhwRCTEKgibi9jN7cehICa8t3el1KSISYhQETURaShtG9mjDM/O3UVymvQIRaTwKgibkjrN6sS+3iDc1RLWINCIFQRNyRq92jEhJ5InPtlBUqr0CEWkcCoImxMy4e3Ifvs4t5tUl6isQkcZRoyAws9ZmFuZ/3NvMLjSziFOs09XMPjOzdWa21sx+WskyZmaPmdkWM1ttZsPq1ozgMaZnW8aktuXvc7fqDCIRaRQ13SOYD0SZWRfgY+AaYPop1ikD7nHO9QNGAz8xs34nLXMecJp/mgY8WcN6gtrdU3pzIL+Yl77I8LoUEQkBNQ0Cc84VAD8A/u6cuxToX90Kzrm9zrkV/sd5wHqgy0mLXQS86Hy+ABLMrFOtWhCERqS0Yfxp7Xhq3jaOFJd5XY6IBLkaB4GZjQGuAj7wzwuv6UbMLAUYCiw56aUuwPGnyGTy7bDAzKaZWbqZpWdlZdV0s83a3ZN7c+hICdMXZ3hdiogEuZoGwV3AfcC7zrm1ZpYKfFaTFc0sBngbuMs5l1uXIp1zzzjn0pxzaUlJSXV5i2ZnaLdEzjq9Pc/M30ZekcYgEpHAqVEQOOfmOecudM790d9pfMA5d+ep1vN3KL8NvOKce6eSRXYDXY97nuyfJ8B/nd2bnMJS/rkow+tSRCSI1fSsoVfNLM7MWgNrgHVmdu8p1jHgOWC9c+6RKhabCfzIf/bQaCDHObe3FvUHtYHJ8Uzu14FnF2wjp0B7BSISGDU9NNTPf1jne8CHQA98Zw5VZ5x/mbPM7Ev/NNXMbjGzW/zLzAa2AVuAZ4HbatuAYPdfZ/cmr6iM5xZu87oUEQlSLWq4XIT/MM/3gCecc6Vm5qpbwTm3ELBTLOOAn9SwhpDUr3McUwd25PlFGfx4bAptYyK9LklEgkxN9wieBjKA1sB8M+sO1KnjV2rv7sm9KSwt5/FPt3hdiogEoZp2Fj/mnOvinJvqP+d/B3BmgGsTv17tY7l8RFde/mIH27LyvS5HRIJMTTuL483skaPn8pvZn/HtHUgjuevs04hsEcYfP9rgdSkiEmRqemjoeSAPuMw/5QL/DFRR8m3tY6O4ZWJP5qz9mmUZh7wuR0SCSE2DoKdz7n7n3Db/9CCQGsjC5NtuHJ9Kh7hIfv/Benz97CIi9VfTICg0szOOPjGzcUBhYEqSqkS3DOeeKX1Ytesws1brcgsRaRg1DYJbgL+ZWYaZZQBPADcHrCqp0sXDkjm9Yyz/N2eDbmkpIg2ipmcNrXLODQYGAYOcc0OBswJamVQqPMz45dS+7DpUyIuLd3hdjogEgVrdocw5l3vcwHF3B6AeqYEJvZOY0DuJxz/dTPaREq/LEZFmrj63qqz2qmEJrF9N7Ut+cRmPfrLZ61JEpJmrTxDotBUP9ekYy+UjuvHyFzvYqovMRKQeqg0CM8szs9xKpjygcyPVKFW4e3JvIluE8b+zdZGZiNRdtUHgnIt1zsVVMsU652o6YJ0ESFJsJLed2Yv/rP+axVsPeF2OiDRT9Tk0JE3ADWf0oEtCNL+ftZ7yCh2tE5HaUxA0c1ER4fz3uX1YtzeXd1Zkel2OiDRDCoIgcOHgzgzumsDDczZSUFLmdTki0swoCIKAmfE/5/dlf16x7m8sIrWmIAgSaSltOLtve56at5XDBbrITERqTkEQRO4953Tyi8v422e6k5mI1JyCIIj06RjLZcO78vyiDNbszvG6HBFpJhQEQeaXU/vSpnVL7p2xmtLyCq/LEZFmQEEQZOJbRfDQ9wawfm8uT83d6nU5ItIMKAiC0JT+HblgcGce+3Qzm77O87ocEWniFARB6oEL+hEbFcG9b62iTIeIRKQaCoIg1TYmkgcu7M+qzByeX7Td63JEpAlTEASxCwZ1YnK/Dvz5401sP3DE63JEpIlSEAQxM+P33xtAy/AwfvPeGpzToHQi8m0KgiDXIS6Kn53ThwWbDzBz1R6vyxGRJkhBEAKuHt2dQcnxPPj+OvbnFXldjog0MQqCEBAeZjxy2WCOFJfxs7dWU6H7FojIcRQEIaJX+1h+/d1+zN+UxfTFGV6XIyJNSMCCwMyeN7P9ZramitcnmVmOmX3pn34TqFrE5+pR3Ti7b3v+8OEG1u/N9bocEWkiArlHMB049xTLLHDODfFPvw1gLYLvLKI/XjyIuOgIfvr6SopKy70uSUSagIAFgXNuPnAoUO8vddM2JpI/XTqITV/n84cPN3hdjog0AV73EYwxs1Vm9qGZ9fe4lpAxqU97rhuXwvTFGXy2Yb/X5YiIx7wMghVAd+fcYOBx4F9VLWhm08ws3czSs7KyGqu+oPbzc0/n9I6x3DtjFVl5xV6XIyIe8iwInHO5zrl8/+PZQISZtati2Wecc2nOubSkpKRGrTNYRUWE8+gVQ8ktKuO/Z6zSVcciIcyzIDCzjmZm/scj/bUc9KqeUNSnYyy/mtqXzzbqlFKRUNYiUG9sZq8Bk4B2ZpYJ3A9EADjnngIuAW41szKgELjC6Wtpo/vRmO7M35TF/87ewOjUtvTtFOd1SSLSyKy5/e1NS0tz6enpXpcRVA7mF3PeowuIj45g5u1nEN0y3OuSRKSBmdly51xaZa95fdaQNAFtYyJ55LIhbN6fz+8+WOd1OSLSyBQEAsAZp7Xjlok9eXXJTt5enul1OSLSiBQEcszPpvRmdGobfvnuV6zdk+N1OSLSSBQEckyL8DCe+OEwElu15JaXl3O4oMTrkkSkESgI5ATtYiL5+9XD2JdTxF1vfKkhq0VCgIJAvmVYt0Tuv6A/czdm8egnm70uR0QCTEEglbpqVDcuHpbMo59s5tMNX3tdjogEkIJAKmVmPPT9AfTvHMddr3/JjoNHvC5JRAJEQSBViooI56mrh2Nm3PzScvKLy7wuSUQCQEEg1eraphVP/HAom/fnc/urKygrr/C6JBFpYAoCOaXxpyXxu4sGMHdjFr+dtU4jlYoEmYANOifB5YejupFx8AjPzN9GStvWXH9GD69LEpEGoiCQGvvFuaez82ABv/tgHV3btGJyvw5elyQiDUCHhqTGwsKMv1w+hEFd4rnztZV8lalhKESCgYJAaiW6ZTjP/jiNNq1bct30pTqtVCQIKAik1trHRvHiDSMpr3Bc89xS3fNYpJlTEEid9EyK4flrR5CVV8y1/1xKXlGp1yWJSB0pCKTOhnZL5O9XD2PjvjxueXk5xWXlXpckInWgIJB6ObNPe/548SAWbTnI3W+s0gVnIs2QTh+Vert4eDKHjpTw0Oz1lJRX8PiVQ4mK0H2PRZoL7RFIg7hpQioPXNCPf6/7muunL9O4RCLNiIJAGsy143rwl8sHs2T7Ia569guyj+gOZyLNgYJAGtT3hybz1NXDWb8vj8ue/px9OUVelyQip6AgkAY3uV8Hpl83gj2HC7nkqcW66EykiVMQSECM7dmO16aN5khxGRc/+Tlf7jrsdUkiUgUFgQTMoOQE3rplDNEtw7j86c+ZuWqP1yWJSCUUBBJQvdrH8t5PzmBwcgJ3vraSR/69iYoK3c9ApClREEjAtWndkpdvHMWlw5N57JPN3PHaSgpLdBWySFOhIJBG0bJFGP93ySB+NbUvs9fs5Xt/W8Tmr/O8LktEUBBIIzIzbpqQyvTrRnIgv5gLnljIa0t36taXIh5TEEijm9g7iQ9/Op607m24752vuP3VleQUavRSEa8ELAjM7Hkz229ma6p43czsMTPbYmarzWxYoGqRpqd9XBQvXj+S/z63Dx+t3cfURxewcPMBr8sSCUmB3COYDpxbzevnAaf5p2nAkwGsRZqgsDDjtkm9eOuWMbRsEcbVzy3hZ2+t4nCBhqYQaUwBCwLn3HzgUDWLXAS86Hy+ABLMrFOg6pGma1i3RD786Xhum9STd1fu5uxH5vH+qj3qOxBpJF72EXQBdh33PNM/T0JQVEQ4/33u6bx/+xl0io/mjtdWcv30ZWQc0PAUIoHWLDqLzWyamaWbWXpWVpbX5UgA9escx7u3jeXX5/dlWUY2U/4ynz/N2ajrDkQCyMsg2A10Pe55sn/etzjnnnHOpTnn0pKSkhqlOPFOi/Awbhyfyqf3TGTqwI488dkWzvrzXN5enqmrkiVk5RQE7sw6L4NgJvAj/9lDo4Ec59xeD+uRJqZ9XBR/vWIob948hqTYSO55axUX/m0h8zdlqf9AQkpJWQWT/zKPh+dsCMj7B+xWlWb2GjAJaGdmmcD9QASAc+4pYDYwFdgCFADXBaoWad5G9mjDv24bx8xVe/i/jzbwo+eXcnrHWG4cn8qFgzvTskWzOMIpUmezVu9hf14xI1LaBOT9rbl9s0pLS3Pp6elelyEeKS4r570v9/Dcgu1s/DqPDnGR3DQ+lStHdqN1pG7BLcHHOcd3H19IcVkF//6vCZhZnd7HzJY759Iqe01fpaRZiWwRzmVpXfnorvG8cP1IerRrze8/WM+4P37Kw3M26CwjCTpLth9i7Z5crh/Xo84hcCr6CiXNkpkxsXcSE3snsWJnNk/O3cqTc7fyt8+2MiIlkUvTunL+wE7aS5Bm76M1+4iOCOcHwwJ3dr0ODUnQ2JdTxDsrM5mRnsm2A0eIiWzBBYM7c8WIrgxKjg/YtymRQLrg8YW0jgzn9Wlj6vU+1R0a0tclCRod46O4bVIvbp3Yk+U7snl92S7eXZnJa0t3cnrHWC4Y3JnzBnQkNSnG61JFaqSwpJz1e3O5eWJqQLejIJCgY2akpbQhLaUNv7mgH++v2sNb6Zk8PGcjD8/ZSO8OMZzbvyPf6duB/p3jaBGurjJpmlZnHqaswjGsW2JAt6MgkKAWFxXBVaO6c9Wo7uw5XMictfv4aM0+nvhsC499uoVWLcMZ2i2BtO5tSEtJZGi3RGLUryDHcc6xfEc2M5ZnsmJnNpnZhZSVOyrc0Qk6x0cxpX9HhnRNYECXeHomtW6QQ5HLd2YDMDTAQaA+AglJB/KLWbz1IOkZh0jPyGb9vlycgzDzDXNxNBhGpLShQ1yU1+WKR95dmcnjn2xh24EjtG4ZzpiebenWpjWREWGEGYSZYcC6vXnM35xFSVkFAEmxkYzt2ZYxqW3967SqUzDc+EI627Ly+fRnk+rdFvURiJykXUwkFw7uzIWDOwOQW1TKyp2HWZ5xiGUZ2by+bCfTF2cA0LVN9AnB0CsphrAwdTwHuzlr93H3m6sY0Dmehy8ZxNRTnIVWUlbB9gNHWLkzm8VbD7Joy0He+3IPALGRLejVIYbe7WM5rUMMQ7omMKxbYrX/j8orHCt3ZjOpT/sGb9vJtEcgUonS8grW7cllmX+PIX1HNgfyiwGIj45gePdE0lISSevehkHJ8URFhHtcsTSUt9J3MWft1yzacoDeHWJ44+Yxdfp8nXNszcrni22H2Lgvj83789j8dT4Hj/jut9ElIZrvDe3M94d2oVf72G+tP2v1Hm5/dSVPXjWM8wbWf4T+6vYIFAQiNeCcY8fBAtJ3ZJOecYhlGYfYmuW7eK1leBgDk+NJ655IWkobhndPpE3rlh5XLHWxZncOFzyxkOTEaPp1iuN3Fw2gfQMfGjyQX8zCzQd4d+VuFmzOosLBwC7xnDugI2f37UDvDr6z2qY+tpDisnL+/V8TCW+APVAFgUgAHDpSwnJ/MKTvyGZ15mFKy32/T10SounZPoaeSa3p1T6Gnkm+qV1MS13P0EQ557js6c/ZmnWEz342ifjoiIBvc39eEe+v2st7X+5mdWYO4DsU2TWxFYu3HuRPlw7mkuHJDbItBYFIIygqLWd1Zg7pO3yHArZm5bN1/xEKS7+5l0JcVAt6to+ha2IrkhOj6ZQQTWR4GEdKyjiQX0xWXjFlFY6E6Jb07hBDqj88kmIjiYlsoRAJkF2HCvh/s9fz4Zp9/OEHA7liZLdGr2FfThGfbtjP3I372XmogMRWLXnxhpFENNDpzQoCEY9UVDj25RaxZX++Lxiy8tmWdYTM7EL2HC6k7Lj7K7QIM9rFRBIeZhw6UnJCgABEtgijfVwkvZJi6Nspjr6d4ujdIZaUdq2IbNG8+igyswtYut23J7Vmdw5Hissor3CUljvKKxytWoaTmhRD+7hIuia2YmzPtrRp3ZK4qAjiWzXsN/WD+cWc89cFHCku49ZJPbn9zF5BeTKAzhoS8UhYmNE5IZrOCdFM6H3iTZXKKxwH84spKa8gKiKcNq1aHvsDVFHh2JVdwK5DhWTlF5GVV8yB/BL25hSxaV8eCzYfOBYi4WFGtzat6JkUQ6/230yd4qMoLa9gf14xuYWlRLYIJzkxmi4J0Z78oSuvcPxn/de8sDiDxVsPAr6zaQZ3TaBrYitahBvhYUZEWBi5RaVszcpn5c7sY52rR/XuEMNp7WPp3rYVg5LjGZScQKf4qDrtLTnn+PW/1pBbWMq/fjKOfp3jGqStzY2CQMQj4WFWZUdkWJjRvW1rurdtXenrxWXlbNmf/61p3qb9x/opqtIyPIyk2EjaxUaSFBNJUmwkXRKi6NqmFcmJrejaJpqkmMgGPQy161ABV/1jCTsPFdA5Pop7z+nDWae3p3eH2FN2hGblFbN0+yEKSsrYn1fMsoxDrNuby5y1+46FYVJsJIOT4+nfOZ7+nePo1zmOLgnRp2zDW8sz+XDNPn5+7ukhGwKgIBBpliJbhPv/6MWfML+0vIKdhwrYsj+frLxiIsJ9h5sSWrWkuLScnYcK2H7wCFl5vv6IzOwCVuzM5tBJ37qjIsJ8oZAYTdc2rejqDwhfULSqdUfqAzPXcjC/mCevGsbkfh1qNaxHUmwk5w/69umTRaW+cXhWZ+awatdhVmUe5pMN+zl6tDs+OoL+neMY3DWBwckJDOmaQMf4b4J3ze4cfv2vNYzt2ZZpEwI7lk9TpyAQCSIR4WHHzlCqzNgq1isoKSMzu5BM/+GoXYcKjh2aSt+RTV5R2QnLx0W1OCEgTg6L48+7/2jNXj7ZsJ9fTe3bIOfDHxUVEc7QboknDL9QUFLGhn15rNuTy7q9uazZncM/Fmw7tpfUIS6SHu1aYxgrdmbTtnVLHr9yaIOcntmcqbNYRE4pp6DUHwzfBMTR55nZhRT7h1Y4Kik2kq6J0bQID2Pp9kP06RDLrDvPaLAzYGrj6J6Db68hh8zsAopKK0hLSeSa0d1DZjRadRaLSL3Et4ogvlU8A7rEf+u1igrHgfzibwLiuLDILijh7sm9+dGY7p6EAFS+5yAnUhCISL2E+Tu928dFMby719VIXWggdhGREKcgEBEJcQoCEZEQpyAQEQlxCgIRkRCnIBARCXEKAhGREKcgEBEJcc1uiAkzywM21nH1eCCnHsuc/Fptnh99fPy8dsCBU1Zd+zpPtUxl82tSe1WP1Y7GbcfJz/V/q3Y11mSZYGxHd+dcUiXL+Mbjbk4TkF6PdZ+pzzInv1ab50cfnzTPk7ZUNr8mtVfTJrWjEduh/1tqR33aUdkUaoeG3q/nMie/Vpvn71exTF3Vpy2Vza9J7dU9riu1o/btOPm5/m/VvJaaLhPM7fiW5nhoKN1VMYJecxMsbVE7mp5gaYva0Tia4x7BM14X0ICCpS1qR9MTLG1ROxpBs9sjEBGRhtUc9whERKQBKQhEREKcgkBEJMQFVRCY2Xgze8rM/mFmi72up67MLMzMHjKzx83sx17XUx9mNsnMFvg/l0le11MfZtbazNLN7Lte11JXZtbX/1nMMLNbva6nPszse2b2rJm9YWZTvK6nrsws1cyeM7MZXtXQZILAzJ43s/1mtuak+eea2UYz22Jmv6juPZxzC5xztwCzgBcCWW9VGqIdwEVAMlAKZAaq1lNpoLY4IB+IwqO2NFA7AH4OvBmYKk+tgX5H1vt/Ry4DxgWy3uo0UFv+5Zy7CbgFuDyQ9ValgdqxzTl3Q2ArPYW6Xu3W0BMwARgGrDluXjiwFUgFWgKrgH7AQHx/7I+f2h+33ptAbHNtB/AL4Gb/ujOa82cChPnX6wC80ozbMRm4ArgW+G5zbYd/nQuBD4EfNuf/W8et92dgWBC0w7Pf9SZz83rn3HwzSzlp9khgi3NuG4CZvQ5c5Jz7X6DS3XMz6wbkOOfyAllvVRqiHWaWCZT4n5YHsNxqNdRn4pcNRAak0FNooM9kEtAa3y90oZnNds5VBLLukzXU5+GcmwnMNLMPgFcDWHKVGugzMeAPwIfOuRUBLrlSDfw74pkmEwRV6ALsOu55JjDqFOvcAPwzYBXVTW3b8Q7wuJmNB+YHsrA6qFVbzOwHwDlAAvBEQCurnVq1wzn3KwAzuxY40NghUI3afh6TgB/gC+XZgSysDmr7e3IHcDYQb2a9nHNPBbK4WqjtZ9IWeAgYamb3+QOjUTX1IKg159z9XtdQX865AnyB1uw5597BF2xBwTk33esa6sM5NxeY63EZDcI59xjwmNd11Jdz7iC+fg7PNJnO4irsBroe9zzZP6+5CZZ2QPC0Re1oeoKlLc2uHU09CJYBp5lZDzNria+zbqbHNdVFsLQDgqctakfTEyxtaX7t8KqXupLe99eAvXxzyuQN/vlTgU34euF/5XWdodKOYGqL2tH0pmBpS7C0Q4POiYiEuKZ+aEhERAJMQSAiEuIUBCIiIU5BICIS4hQEIiIhTkEgIhLiFAQSNMwsv5G316j3vDCzBDO7rTG3KaFBQSBSBTOrdiwu59zYRt5mAqAgkAanIJCgZmY9zewjM1vuv1Pa6f75F5jZEjNbaWb/MbMO/vkPmNlLZrYIeMn//Hkzm2tm28zszuPeO9//c5L/9RlmtsHMXvEPkYyZTfXPW25mj5nZrEpqvNbMZprZp8AnZhZjZp+Y2Qoz+8rMLvIv+gegp5l9aWYP+9e918yWmdlqM3swkP+WEsS8vrRZk6aGmoD8SuZ9ApzmfzwK+NT/OBGOXVl/I/Bn/+MHgOVA9HHPF+MbtrkdcBCIOH57wCQgB9/gYmHA58AZ+O7Ktgvo4V/uNWBWJTVei294gjb+5y2AOP/jdsAWwIAUTrwByhTgGf9rYfhudDLB689BU/Obgm4YapGjzCwGGAu85f+CDt/cHCcZeMPMOuG7i9T241ad6ZwrPO75B865YqDYzPbju9vaybfdXOqcy/Rv90t8f7TzgW3OuaPv/RowrYpy/+2cO3S0dOD/mdkEoALf+PYdKllnin9a6X8eA5xG07uHhTRxCgIJZmHAYefckEpeexx4xDk303+zlgeOe+3IScsWH/e4nMp/b2qyTHWO3+ZVQBIw3DlXamYZ+PYuTmbA/zrnnq7ltkROoD4CCVrOuVxgu5ldCr5bG5rZYP/L8XwzRvyPA1TCRiD1uFsZ1vQG6/HAfn8InAl098/PA2KPW24OcL1/zwcz62Jm7etftoQa7RFIMGnlv9/zUY/g+3b9pJn9GogAXsd3M/EH8B0yygY+BXo0dDHOuUL/6Z4fmdkRfOPU18QrwPtm9hWQDmzwv99BM1tkZmvw3af3XjPrC3zuP/SVD1wN7G/otkhw0zDUIgFkZjHOuXz/WUR/AzY75/7idV0ix9OhIZHAusnfebwW3yEfHc+XJkd7BCIiIU57BCIiIU5BICIS4hQEIiIhTkEgIhLiFAQiIiFOQSAiEuL+P0JdY53qljhrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Deleted the fake labels used in Learning rate finder\n"
     ]
    }
   ],
   "source": [
    "params = semisup_model.parameters()\n",
    "\n",
    "### For SGD and Adam ###\n",
    "learning_rate1, learning_rate2 = 1e-7, 1e-1\n",
    "\n",
    "### For LBFGS (a good choice already!!!) ###\n",
    "# print(\"Using LBFGS's learning rate set\")\n",
    "# learning_rate1, learning_rate2 = 8e-2, 5e-2 # (1e-1, 5e-2) is also OK!\n",
    "\n",
    "choice = 'MAD'; auto_lr = True\n",
    "if choice == 'LBFGS':\n",
    "    optimizer1 = torch.optim.LBFGS(params, lr=learning_rate1, \n",
    "                                   max_iter=100, max_eval=125, \n",
    "                                  history_size=120, line_search_fn='strong_wolfe')\n",
    "if choice == 'Adam':\n",
    "    optimizer1 = AdamGC(params, lr=learning_rate1, use_gc=True, gc_conv_only=False, gc_loc=False)\n",
    "if choice == 'SGD':\n",
    "    optimizer1 = SGDGC(params, lr=learning_rate1, use_gc=True, nesterov=True, momentum=0.95)\n",
    "if choice == 'MAD':\n",
    "    optimizer1 = MADGRAD(params, lr=learning_rate1, momentum=0.9)\n",
    "\n",
    "if choice != 'LBFGS' and auto_lr:\n",
    "    print('Learning rate finding')\n",
    "    bs = 4000; bs = X_u_train.shape[0] if bs>X_u_train.shape[0] else bs\n",
    "    criterion = LadderLoss(return_list=True)\n",
    "    trainloader = get_dataloader(X_u_train, u_train, bs=bs)\n",
    "    \n",
    "    lr_finder = LRFinder(semisup_model, optimizer=optimizer1, \n",
    "                         closure=pcgrad_update, criterion=criterion, device=\"cpu\")\n",
    "    lr_finder.range_test(trainloader, val_loader=None, end_lr=100, num_iter=300)\n",
    "    \n",
    "    # to inspect the loss-learning rate graph\n",
    "    suggested_lr, _ = lr_finder.plot()\n",
    "    # To prevent divergence during the second stage training.\n",
    "    # suggested_lr = min(suggested_lr, 5e-3)\n",
    "    lr_finder.reset(); plt.show()\n",
    "\n",
    "else:\n",
    "    lr_finder = None\n",
    "    suggested_lr = None\n",
    "    \n",
    "print(\"Deleted the fake labels used in Learning rate finder\")\n",
    "u_train = u_train[:N, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the learing_rate to the suggested one.\n",
    "# suggested_lr = 1e-5\n",
    "\n",
    "if lr_finder and suggested_lr:\n",
    "    optimizer1 = lr_finder.optimizer\n",
    "\n",
    "for g in optimizer1.param_groups:\n",
    "    g['lr'] = suggested_lr\n",
    "\n",
    "epochs1 = 1000; epochs2 = 10000;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using old implementation of TorchMLP. See models.py for more new model-related source code.\n"
     ]
    }
   ],
   "source": [
    "# Setting up the generator\n",
    "generator = TorchMLP([2, 50, 50, 2])\n",
    "# generator_training_epochs indicates how string the generator is\n",
    "adv_f = 100; generator_training_epochs = 300; generator_training_limit = epochs1-100\n",
    "# I can use the Learning rate finder to find a good lr for the generator optim  as well\n",
    "generator_optimizer = torch.optim.SGD(generator.parameters(), lr=3e-4, momentum=0.95)\n",
    "sinkhorn_loss = SamplesLoss(\"sinkhorn\", p=2, blur=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiments with the generator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 0/300 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:24<00:00, 12.49it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -3.7287700176239014\n",
      "Semi-supervised solver loss @Epoch 0:  3.4609358310699463\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:28, 10.57it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:24<00:00, 12.12it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: 0.10761624574661255\n",
      "Semi-supervised solver loss @Epoch 100:  0.8972504138946533\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 1/300 [00:00<00:32,  9.24it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:26<00:00, 11.36it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: 0.22298726439476013\n",
      "Semi-supervised solver loss @Epoch 200:  0.7774209976196289\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 1/300 [00:00<00:31,  9.41it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:24<00:00, 12.46it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.1296619176864624\n",
      "Semi-supervised solver loss @Epoch 300:  0.7753944396972656\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:23, 12.87it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:24<00:00, 12.19it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.18763139843940735\n",
      "Semi-supervised solver loss @Epoch 400:  0.758278489112854\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:24, 12.41it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:25<00:00, 11.55it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.4919699430465698\n",
      "Semi-supervised solver loss @Epoch 500:  0.7504470944404602\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:20, 14.43it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:25<00:00, 11.83it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.5540295839309692\n",
      "Semi-supervised solver loss @Epoch 600:  0.7530421614646912\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:24, 12.39it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:25<00:00, 11.64it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.5718031525611877\n",
      "Semi-supervised solver loss @Epoch 700:  0.716842770576477\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:23, 12.77it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:22<00:00, 13.15it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.6112722754478455\n",
      "Semi-supervised solver loss @Epoch 800:  0.7253022789955139\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:22, 13.07it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:24<00:00, 12.48it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.6140658259391785\n",
      "Semi-supervised solver loss @Epoch 900:  0.7153809070587158\n"
     ]
    }
   ],
   "source": [
    "curr_loss = 1000; F_print = 10 if choice == 'LBFGS' else 100\n",
    "\n",
    "# Stage I\n",
    "# How long should I pretrain selector part of the model?\n",
    "for i in range(epochs1):\n",
    "    if i%adv_f==0 and i<=generator_training_limit:\n",
    "        best_generator_loss = 1000; best_generator_state_dict = None\n",
    "        o_tensor = X_u_train[:N, :] # or X_u_train ?\n",
    "        print(\"Training the generator for \")\n",
    "        for _ in trange(generator_training_epochs):\n",
    "            semisup_model.eval()\n",
    "            generator.train()\n",
    "            generator_optimizer.zero_grad()\n",
    "            X_gen = scale_to_range(generator(o_tensor), lb, ub)\n",
    "            unsup_loss = semisup_model(X_gen)[1]\n",
    "#                 d_loss = F.mse_loss(X_gen, o_tensor)\n",
    "            d_loss = F.mse_loss(X_gen, o_tensor)\n",
    "            generator_loss = d_loss-unsup_loss\n",
    "            generator_loss.backward(retain_graph=True)\n",
    "            generator_optimizer.step()\n",
    "            # Saving the best_generator_state_dict\n",
    "            if generator_loss.item() < best_generator_loss:\n",
    "                best_generator_loss = generator_loss.item()\n",
    "                best_generator_state_dict = generator.state_dict()\n",
    "\n",
    "        print(\"The best generator loss:\", best_generator_loss)\n",
    "        generator.load_state_dict(best_generator_state_dict)\n",
    "        generator.eval()\n",
    "        X_gen = scale_to_range(generator(o_tensor), lb, ub)\n",
    "        if N_res<X_gen.shape[0]: X_gen = sampling_from_rows(X_gen, N_res)\n",
    "        X_u_train = torch.cat([X_u_train[:N, :], X_gen], dim=0).detach().requires_grad_(True)\n",
    "\n",
    "    semisup_model.train()\n",
    "    generator_optimizer.zero_grad()\n",
    "    optimizer1.step(pcgrad_closure)\n",
    "    l = pcgrad_closure()\n",
    "    if (i % F_print) == 0:\n",
    "        if l.item() != curr_loss:\n",
    "            curr_loss = l.item()\n",
    "        else:\n",
    "            print(\"Semi-supervised solver loss @Epoch {}: \".format(i), curr_loss)\n",
    "            print(\"Finishing the first stage\")\n",
    "            break\n",
    "        print(\"Semi-supervised solver loss @Epoch {}: \".format(i), curr_loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 0:  1.6481116063005175e-06\n",
      "Epoch 100:  2.6807390440808376e-07\n",
      "Epoch 200:  1.9208675894333282e-07\n",
      "Epoch 300:  1.5326980928875855e-07\n",
      "Epoch 400:  1.3145856314622506e-07\n",
      "Epoch 500:  1.0589058518917227e-07\n",
      "Epoch 600:  9.419993318715569e-08\n",
      "Epoch 700:  8.683706909096145e-08\n",
      "Epoch 800:  7.822136893764764e-08\n",
      "Epoch 900:  6.842613942126263e-08\n",
      "Epoch 1000:  6.083710246684859e-08\n",
      "Epoch 1100:  5.4225040457822615e-08\n",
      "Epoch 1200:  5.054048202168815e-08\n",
      "Epoch 1300:  4.666467745551017e-08\n",
      "Epoch 1400:  4.061818970058084e-08\n",
      "Epoch 1500:  3.744245447023786e-08\n",
      "Epoch 1600:  3.399074088861198e-08\n",
      "Epoch 1700:  3.1948403034220973e-08\n",
      "Epoch 1800:  3.0352030222502435e-08\n",
      "Epoch 1900:  2.77018248340255e-08\n",
      "Epoch 2000:  2.5318604102153586e-08\n",
      "Epoch 2100:  2.4422529776302326e-08\n",
      "Epoch 2200:  2.279122846005066e-08\n",
      "Epoch 2300:  2.1638648206590005e-08\n",
      "Epoch 2400:  1.9212235358168073e-08\n",
      "Epoch 2500:  1.8598921514012545e-08\n",
      "Epoch 2600:  1.7606833324634863e-08\n",
      "Epoch 2700:  1.659081405591678e-08\n",
      "Epoch 2800:  1.6141104453026855e-08\n",
      "Epoch 2900:  2.3311349295340733e-08\n",
      "Epoch 3000:  1.4466602138440976e-08\n",
      "Epoch 3100:  1.4304887940852495e-08\n",
      "Epoch 3200:  1.3806793930370986e-08\n",
      "Epoch 3300:  1.3120138753208721e-08\n",
      "Epoch 3400:  1.301831531463904e-08\n",
      "Epoch 3500:  1.2827553241834266e-08\n",
      "Epoch 3600:  1.2114295344645143e-08\n",
      "Epoch 3700:  1.184445341806395e-08\n",
      "Epoch 3800:  1.1455215442879307e-08\n",
      "Epoch 3900:  1.0453792498310577e-08\n",
      "Epoch 4000:  1.0018262663891164e-08\n",
      "Epoch 4100:  9.911864218281607e-09\n",
      "Epoch 4200:  9.203590778383841e-09\n",
      "Epoch 4300:  8.415875107914417e-09\n",
      "Epoch 4400:  7.962842829556394e-09\n",
      "Epoch 4500:  7.685143188496113e-09\n",
      "Epoch 4600:  7.511879118737852e-09\n",
      "Epoch 4700:  7.42566719225124e-09\n",
      "Epoch 4800:  7.161819137735392e-09\n",
      "Epoch 4900:  6.8120118434933374e-09\n"
     ]
    }
   ],
   "source": [
    "optimizer2 = LBFGSNew(semisup_model.network.parameters(), \n",
    "                              lr=learning_rate2, max_iter=300, max_eval=int(300*1.25),\n",
    "                              history_size=150, line_search_fn=True, batch_mode=False)\n",
    "\n",
    "curr_loss = 1000\n",
    "# Stage II: Train semisup_model.network\n",
    "semisup_model.network.train()\n",
    "semisup_model.selector.eval()\n",
    "for i in range(5000):\n",
    "    optimizer2.step(closure)\n",
    "    l = closure()\n",
    "    if (i % 100) == 0:\n",
    "        if l.item() != curr_loss:\n",
    "            curr_loss = l.item()\n",
    "        print(\"Epoch {}: \".format(i), curr_loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Testing\n",
      "Test MSE: 4.3e-04\n"
     ]
    }
   ],
   "source": [
    "print(\"Testing\")\n",
    "semisup_model.network.eval()\n",
    "# Compare btw the two semi-supervise learning?\n",
    "test_performance = F.mse_loss(semisup_model.network(*dimension_slicing(X_star)).detach(), u_star).item()\n",
    "string_test_performance = scientific2string(test_performance)\n",
    "print('Test MSE:', string_test_performance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.42464846, 0.3930947 , 0.43438876, 0.5195853 , 0.5051994 ],\n",
       "      dtype=float32)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(semisup_model.selector.latest_weighted_features).detach().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 0:  1.4170188903808594\n",
      "[4 3 1 2 0]\n",
      "Epoch 5:  1.0771796703338623\n",
      "[4 3 1 2 0]\n",
      "Epoch 10:  1.0343698263168335\n",
      "[4 3 1 2 0]\n",
      "Epoch 15:  1.0276494026184082\n",
      "[4 3 1 2 0]\n",
      "Epoch 20:  1.0300183296203613\n",
      "[4 3 1 2 0]\n",
      "Epoch 25:  1.022573709487915\n",
      "[4 3 1 2 0]\n",
      "Epoch 30:  1.0247873067855835\n",
      "[4 3 1 2 0]\n",
      "Epoch 35:  1.042427897453308\n",
      "[4 3 1 2 0]\n",
      "Epoch 40:  1.01158607006073\n",
      "[4 3 1 2 0]\n",
      "Epoch 45:  1.0314855575561523\n",
      "[4 3 1 2 0]\n",
      "Epoch 50:  1.0298184156417847\n",
      "[4 3 1 2 0]\n",
      "Epoch 55:  1.0229662656784058\n",
      "[4 3 1 2 0]\n",
      "Epoch 60:  1.0415446758270264\n",
      "[4 3 1 2 0]\n",
      "Epoch 65:  1.0460270643234253\n",
      "[4 3 1 2 0]\n",
      "Epoch 70:  1.0263696908950806\n",
      "[4 3 1 2 0]\n",
      "Epoch 75:  1.0206340551376343\n",
      "[4 3 1 2 0]\n",
      "Epoch 80:  1.0002954006195068\n",
      "[4 3 1 2 0]\n",
      "Epoch 85:  1.013137698173523\n",
      "[4 3 1 2 0]\n",
      "Epoch 90:  1.0339939594268799\n",
      "[4 3 1 2 0]\n",
      "Epoch 95:  1.039125680923462\n",
      "[4 3 1 2 0]\n"
     ]
    }
   ],
   "source": [
    "ders, dys = semisup_model.network.get_selector_data(*dimension_slicing(X_u_train))\n",
    "\n",
    "optimizer3 = LBFGSNew(semisup_model.selector.parameters(), \n",
    "                      lr=learning_rate2, max_iter=300, max_eval=int(300*1.25),\n",
    "                      history_size=150, line_search_fn=True, batch_mode=False)\n",
    "\n",
    "# Stage II: Train semisup_model.selector\n",
    "semisup_model.network.eval()\n",
    "semisup_model.selector.train()\n",
    "for i in range(100):\n",
    "    optimizer3.step(selector_closure)\n",
    "    l = selector_closure()\n",
    "    if (i % 5) == 0:\n",
    "        print(\"Epoch {}: \".format(i), l.item())\n",
    "        print(np.argsort((semisup_model.selector.latest_weighted_features).detach().numpy()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.8624258 , 0.6867128 , 0.8042834 , 0.36450225, 0.10327309],\n",
       "      dtype=float32)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "feature_importance = (semisup_model.selector.latest_weighted_features).detach().numpy()\n",
    "feature_importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEQCAYAAAC9VHPBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAi5ElEQVR4nO3de7xf053/8dc7QYO6JulUJYQSqv2loaFFq4rWtUxJRas6VCfaqopWZ8pQQWt0XEYoHaFEQ92rTVxKqdAhrSREEHWZoIlekCpa1/D5/bHWkW9OzmXl5Oz9zTnn/Xw8vo+z99q3z97f7/mu71p77bUUEZiZWd/Vr9kBmJlZczkjMDPr45wRmJn1cc4IzMz6OGcEZmZ9nDMCM7M+zhmBmVkf12lGIGl1Sf3y9HBJe0taufrQzMysDursgTJJs4CPAesAdwEzgNcj4sDqwzMzs6qVVA0pIl4G9gXOi4jPAu+vNiwzM6tLUUYgaVvgQOCGnNa/upDMzKxOJRnBOOAY4LqIeEjSxsDtlUZlZma16fQewdsrSqvlKiIzM+tFSloNbStpLvD7PP9BSedVHpmZmdWipGroLGBXYCFARNwP7FBhTGZmVqOVSlaKiPmSGpPerCaczg0aNCiGDRvWrMObmfVIs2bNei4iBre1rCQjmC9pOyDyg2RHAg93Z4DLYtiwYcycObNZhzcz65EkPdXespKqoa8AhwPrA08DI/O8mZn1Ap2WCCLiOdIzBGZm1guVtBq6RNLaDfPrSLqo0qjMzKw2JVVDIyLiby0zEfE8sGVlEZmZWa1KMoJ+ktZpmZG0LoWtjczMbMVX8oV+BjBd0tWAgNHA9yuNyszMalNys/gnuSvqT+SkfSNibrVhmZlZXUqreH4PPN+yvqQNIuIPlUVlZma16TQjkHQEcALwF9ITxQICGFFtaGZmVoeSEsGRwGYRsbDqYCq3ZDcZ3auwF1czsxVNSauh+cALVQdiZmbNUVIimAdMk3QD8FpLYkScWVlUZmZWm5KM4A/5tUp+mZlZL1LSfPTEOgIxM7PmKGk1NBj4N+D9wICW9IjYqcK4zMysJiU3iy8jPUewEXAi8CQwo8KYzMysRiUZwcCI+DHwRkTcERFfAlwaMDPrJUpuFr+R//5J0p7AH4F1qwvJzMzqVJIRfE/SWsC3gHOANYFxVQZlZmb1KckIno+IF0gPlX0CQNL2lUZlZma1KblHcE5h2lIk7SbpEUmPS/pOG8s3kHS7pPskzZG0R8l+zcys+7RbIpC0LbAdMFjSNxsWrQn072zHkvoD5wKfBBYAMyRNadWF9XHAVRHxI0lbADcCw5b5LMzMrMs6KhGsAryTlFms0fB6kTQ4TWe2AR6PiHkR8TpwBbBPq3WClLEArEW6EW1mZjVqt0QQEXdI+l/SmMVdebp4fVKHdS0WAB9utc544Jbc1fXqwC5dOI7ZsquqJ1r3Qms9UIf3CCLiTeA9FR7/c8CkiBgC7AFMlrRUTJLGSpopaeazzz5bYThmZn1PSauh2ZKmAFcD/2hJjIifdbLd08DQhvkhOa3RocBueX/TJQ0ABgHPNK4UEROBiQCjRo3yTy4zs25UkhEMABay5NPEAXSWEcwANpW0ESkDOAD4fKt1/gDsDEyS9L58LP/kNzOrUUnvo4d0ZccRsUjS14GbSa2MLoqIhySdBMyMiCmkh9QukHQUKXM5OMKVrGZmdSrpfXQI6bmBlofIfgMcGRELOts2Im4kNQltTPtuw/Tchv2amVkTlDxQdjEwhXTT+D3A1JxmZma9QElGMDgiLo6IRfk1CRhccVxmZlaTkoxgoaQvSOqfX18g3Tw2M7NeoCQj+BKwP/Dn/BoNdOkGspmZrXhKWg09BexdQyxmZtYEnZYIJG0saaqkZyU9I+kXkjauIzgzM6teSdXQT4GrgPVIrYauBi6vMigzM6tPSUawWkRMbmg1dCnpCWAzM+sFSrqYuCkPKnMF6enfMcCNktYFiIi/VhifmZlVrCQj2D//PaxV+gGkjMH3C8zMerCSVkMb1RGImZk1R0lfQ/2BPUlDSL69fkScWV1YZmZWl5KqoanAq8ADwFvVhmNmZnUryQiGRMSIyiMxM7OmKGk+epOkT1UeiZmZNUVJieC3wHV5LOE3AAEREWtWGpmZmdWiJCM4E9gWeMCjh5mZ9T4lVUPzgQedCZiZ9U4lJYJ5wDRJNwGvtSS6+aiZWe9QkhE8kV+r5JeZmfUiJU8Wn1hHIGZm1hztZgSSzoqIcZKmkvoUWkJEeLAaM7NeoKMSweT89/Q6ArGaSNXt2+0JzHqkdjOCiJiV/95RXzhmZla3kuajZmbWizkjMDPr44ozAkmrVRmImZk1R6cZgaTtJM0Ffp/nPyjpvMojMzOzWpSUCP4b2BVYCBAR9wM7VBmUmZnVp6hqKCLmt0p6s4JYzMysCUq6mJgvaTsgJK0MHAk8XG1YZmZWl5ISwVeAw4H1gaeBkXnezMx6gQ5LBHng+gkRcWBN8ZiZWc06LBFExJvAhpLc66iZWS9VOh7BXZKmAP9oSfR4BGZmvUNJRvB/+dUPWKPacMzMrG4ej8DMrI/rNCOQdDttj0ewUyURmZlZrUqqho5umB4A7AcsqiYcMzOrW0nV0KxWSXdJuqdk55J2AyYA/YELI+LUNtbZHxhPKnXcHxGfL9m3mZl1j5KqoXUbZvsBHwLWKtiuP3Au8ElgATBD0pSImNuwzqbAMcD2EfG8pHctY/xmZracSqqGZpF+rYtUJfQEcGjBdtsAj0fEPABJVwD7AHMb1vlX4NyIeB4gIp4pD93MzLpDSUbwvoh4tTFB0jsKtlsfaOysbgHw4VbrDM/7u4tUfTQ+In5ZsG8zM+smJX0N3d1G2vRuOv5KwKbAjsDngAskrd16JUljJc2UNPPZZ5/tpkObmRl0UCKQ9G7Sr/pVJW1JqhoCWBMoGa3saWBow/yQnNZoAfC7iHgDeELSo6SMYUbjShExEZgIMGrUqKWaspqZWdd1VDW0K3Aw6Qu8sTuJl4BjC/Y9A9hU0kakDOAAoHWLoJ+TSgIXSxpEqiqaVxK4mZl1j3Yzgoi4BLhE0n4Rce2y7jgiFkn6OnAzqf7/ooh4SNJJwMyImJKXfSoPhfkm8O2IWNilMzEzsy5RROc1LZL2BN5PeqAMgIg4qcK42jVq1KiYOXNm1zaWOl+nqwqu4wrB1yCp6jr0pGtgfYqkWRExqq1lJYPX/w8wBjiCdJ/gs8CG3RqhmZk1TUmroe0i4ovA87kDum3JzT7NzKznK8kIXsl/X5b0HuANYL3qQjIzszqVPFB2fW7bfxpwL+kp4wurDMrMzOpT0uncyXnyWknXAwMi4oVqwzIzs7qU3CxeTdLxki6IiNeAd0naq4bYzMysBiX3CC4GXiPdJIb0cNj3KovIzMxqVZIRvDci/ot0k5iIeJnF3U2YmVkPV5IRvC5pVfJwlZLeSyohmJlZL1DSaugE4JfAUEmXAduT+iAyM7NeoKPeR7ePiLuAO4F9gY+QqoSOjIjnaorPzMwq1lGJ4GzSsJTTI2Ir4IZ6QjIzszp1lBG8IWkiMETS2a0XRsQ3qgvLzMzq0lFGsBewC2lcgln1hGNmZnXraDyC54ArJD0cEffXGJOZmdWo0+ajzgTMzHq3kucIzMysF3NGYGbWx3X6QJmkI0n9Db1E6n56S+A7EXFLxbE13Y5tpO0PfA14GdhjiZXT2gcffDAHH3wwzz33HKNHj15q+69+9auMGTOG+fPnc9BBBy21/Fvf+haf/vSneeSRRzjssMOWWn7cccexyy67MHv2bMaNG7fU8lNOOYXtttuOu+++m2OPPXap5WcBI4FbabvDqPOBzYCpwBltLJ8MDAWuBH7UeuGOO3LNNdcwaNAgJk2axKRJk5ba/sYbb2S11VbjvPPO46qrrlpq+bRp0wA4/fTTuf7665dYtuqqq3LTTTcBcPLJJ3PbbbctsXzgwIFce20aXvuYY45h+vTpSywfMmQIl156KQDjgNmtjj0cmJinxwKPtlo+knT9AL4ALGi1fFvgP/P0fvvtx8KFSw6/vfPOO3P88ccDsPvuu/PKK68ssXyvvfbi6KOPBmDH/HlqtP/++/O1r32Nl19+mT322GOp5Sv8Z++ssxg5ciS33nor3/ve0p++888/n80224ypU6dyxhlLf/omT57M0KFDufLKK/nRj5b69PWcz964ccyePXuJ5cOHD2fixPTpGzt2LI8+uuSnb+TIkZx11llLxdxdSp4s/lJETJC0K7AOcBDp+6DHZQQav4wbXLx00h3vh8O3AV4HLluc/vGuh1WrLQ8jDSv0f6RHBVvZ/NPAIOAR4O6ll2+wL7AW8CAwY8llPeUaAEz4MPDnJdPuGAgX7J1npgALWy1/N0zYPc9cC7zYavnQxRmBWU/S6eD1kuZExAhJE4BpEXGdpPsiYst6QlzS8gxerxOr6ysvTugZg5b7GiRVXYeedA2sb1muweuBWZJuIdWE3CxpDeCt7gzQzMyap6Rq6FBS9ei8iHhZ0kDgkEqjMjOz2pQMVfmWpCeA4ZIG1BCTmZnVqKTV0JeBI4EhpIYWHwGmAztVGpmZmdWi5B7BkcDWwFMR8QlS89G/VRmUmZnVpyQjeDUiXgWQ9I6I+D2pqbmZmfUCJTeLF0haG/g58CtJzwNPVRmUmZnVp+Rm8Wfy5HhJt5MeJ/plpVGZmVltSkoEb4uIO6oKxMzMmsOdzpmZ9XHOCMzM+rhOMwJJq0vql6eHS9pb0srVh2ZmZnUoKRHcCQyQtD6px9GDgElVBmVmZvUpyQgUES8D+wLnRcRngfdXG5aZmdWlKCOQtC1wIHBDTutfXUhmZlan0i4mjgGui4iHJG0M3F5tWGZmVpeS5wj+KSJaxm0iIuZJ+k2FMZmZWY1KSgTHFKaZmVkP1G6JQNLupFHJ1pd0dsOiNYFFJTuXtBswgXRP4cKIOLWd9fYDrgG2joiujUNpZmZd0lHV0B+BmcDewKyG9JeAozrbsaT+wLnAJ4EFwAxJUyJibqv11iDdh/jdsoVuZmbdod2MICLuB+6XdFlEFJUAWtkGeDwi5gFIugLYB5jbar2TgR8A3+7CMczMbDmV3CN4TNK81q+C7dYH5jfML8hpb5O0FTA0Im7AzMyaoqTV0KiG6QHAZ4F1l/fAuduKM4GDC9YdC4wF2GCDDZb30GZm1qDTEkFELGx4PR0RZwF7Fuz7aWBow/yQnNZiDeADwDRJT5LGQp4iqTHjaYlhYkSMiohRgwcPLji0mZmVKhm8fquG2X6kEkJJSWIGsKmkjUgZwAHA51sWRsQLwKCG40wDjnarITOzepV8oZ/RML0IeBLYv7ONImKRpK8DN5Oaj16Un0w+CZgZEVO6EK+ZmXWzkqEqP9HVnUfEjcCNrdK+2866O3b1OGZm1nUl4xEMlHS2pHslzZI0QdLAOoIzM7PqlTQfvQJ4FtgPGJ2nr6wyKDMzq0/JPYL1IuLkhvnvSRpTVUBmZlavkozgFkkHAFfl+dGkG8Bm1oPpRFWy3zghKtmvVaekauhfgZ8Cr+XXFcBhkl6S9GKVwZmZWfVKWg2tUUcgZmbWHCWthm4rSTMzs56po/EIBgCrAYMkrQO0VCiuSavO48zMrOfqqGroMGAc8B7g3ob0F4EfVhiTmZnVqKPxCCYAEyQdERHn1BiTmZnVqKT56AuSvtg6MSJ+UkE8ZmZWs5KMYOuG6QHAzqSqImcEZma9QEnz0SMa5yWtTXqWwMzMeoGSB8pa+wewUXcHYmZmzVEyMM1UoOWZ8f7A+1jc3YSZmfVwJfcITm+YXgQ8FRELKorHzMxqVjJm8R3A70ljDK8DvF51UGZmVp+SLib2B+4BPksaovJ3kkZXHZiZmdWjpGroP4CtI+IZAEmDgVuBa6oMzMzM6lHSaqhfSyaQLSzczszMeoCSEsEvJd0MXJ7nx9BqQHozM+u5Sh4o+7akfYGP5qSJEXFdtWGZmVldSkoERMTPgJ9VHIuZmTWB6/rNzPo4ZwRmZn3cMmUEktaRNKKqYMzMrH4lD5RNk7SmpHVJ3U9fIOnM6kMzM7M6lJQI1oqIF4F9gZ9ExIeBXaoNy8zM6lKSEawkaT1S9xLXVxyPmZnVrCQjOAm4GXg8ImZI2hh4rNqwzMysLiUPlF0NXN0wPw/Yr8qgzMysPu1mBJJOIA1I8/eI8M1hM7NeqqMSwZP578s1xGFmZk3S7j2CiLgEuBT4cH3hmJlZ3Tq8WRwRbwLb1xSLmZk1QUmnc7MlTSHdMP5HS2LuiM7MzHq4koxgAGkwmp0a0gL3Rmpm1iuUNB89pI5AzMysOUr6Ghou6TZJD+b5EZKOqz40MzOrQ8mTxRcAxwBvAETEHOCAkp1L2k3SI5Iel/SdNpZ/U9JcSXNyZrPhsgRvZmbLryQjWC0i7mmVtqizjST1B84Fdge2AD4naYtWq90HjIqIEcA1wH8VxGNmZt2oJCN4TtJ7STeIkTQa+FPBdtuQ+ieaFxGvA1cA+zSuEBG3R0TLA2u/BYYUR25mZt2ipNXQ4cBEYHNJTwNPAAcWbLc+ML9hfgEdP5x2KHBTwX7NzKwblWQET0XELpJWB/pFxEvdHYSkLwCjgI+3s3wsMBZggw026O7Dm5n1aSVVQ49JOg3YYBkzgaeBoQ3zQ3LaEiTtAvwHsHdEvNbWjiJiYkSMiohRgwcPXoYQzMysMyUZwQeBR4EfS/qtpLGS1izYbgawqaSNJK1Camk0pXEFSVsC55MygWeWMXYzM+sGnWYEEfFSRFwQEdsB/w6cAPxJ0iWSNulgu0XA10mD2jwMXBURD0k6SdLeebXTgHcCV0tq6crCzMxq1Ok9gtwMdE/gEGAYcAZwGfAx4EZgeHvbRsSNeZ3GtO82THvsYzOzJiu5WfwYcDtwWkTc3ZB+jaQdqgnLzMzqUpIRjIiIv7e1ICK+0c3xmJlZzUoygkWSDgfeT+qJFICI+FJlUZmZWW1KWg1NBt4N7ArcQWoG2u3PEpiZWXOUZASbRMTxwD/y8JV74uErzcx6jZKM4I3892+SPgCsBbyrupDMzKxOJfcIJkpaBzie9EDYO4HvdryJmZn1FCUjlF2YJ+8ANq42HDMzq1u7GYGkb3a0YUSc2f3hmJlZ3ToqEaxRWxRmZtY07WYEEXFinYGYmVlzlLQaMjOzXswZgZlZH+eMwMysj3OrITOzPq6k1dBmwNYsHl3s08A9VQZlZmb16bTVkKQ7ga1axiuWNB64oZbozMysciX3CP4JeL1h/vWcZmZmvUBJX0M/Ae6RdF2e/2dgUlUBmZnVRSeqkv3GCVHJfqtS0tfQ9yXdRBqjGOCQiLiv2rDMzKwuHbUaWjMiXpS0LvBkfrUsWzci/lp9eGZmVrWOSgQ/BfYCZgEBNJahAvdEambWK3TUamiv/Hej+sIxM7O6ldwsRtLewA55dlpEXF9dSGZmVqdOm49KOhU4EpibX0dKOqXqwMzMrB4lJYI9gJER8RaApEuA+4BjqwzMzMzqUdrp3NoN02tVEIeZmTVJSYngP4H7JN1Oajm0A/CdSqMyM7PalDxQdrmkaaSO5wD+PSL+XGlUZmZWm9KqocH570rAdpL2rSgeMzOrWaclAkkXASOAh4C3cnIAP6swLjMzq0nJPYKPRMQWlUdiZmZNUVI1NF2SMwIzs16qtBvq6ZL+DLxGajkUETGi0sjMzKwWJRnBj4GDgAdYfI/AzMx6iZKM4NmImNL5amZm1hOVZAT3SfopMJVUNQRARLjVkJlZL1CSEaxKygA+1ZDm5qNmZr1EyZPFh9QRiJmZNUfJA2Vnt5H8AjAzIn7Ryba7AROA/sCFEXFqq+XvILVK+hCwEBgTEU+WhW5mZt2h5DmCAcBI4LH8GgEMAQ6VdFZ7G0nqD5wL7A5sAXyujecRDgWej4hNgP8GfrCM8ZuZ2XIquUcwAtg+It4EkPQj4DfAR0lNStuzDfB4RMzL210B7EMa3KbFPsD4PH0N8ENJiohYlpMwM7OuKykRrAO8s2F+dWDdnDG81vYmAKwPzG+YX5DT2lwnIhaRqpwGFsRkZmbdpKRE8F/A7NwVdct4BKdIWh24tcLY3iZpLDA2z/5d0iN1HBcYBDxXsqLGq+JQmsbXwNcAfA2g51+DDdtbUNJq6MeSbiRV9QAcGxF/zNPf7mDTp4GhDfNDclpb6yyQtBJp9LOFbcQwEZjYWazdTdLMiBhV93FXJL4GvgbgawC9+xq0WzUkafP8dytgPVIVznzg3TmtMzOATSVtJGkV4ACg9RPKU4B/ydOjgV/7/oCZWb06KhF8k1Qdc0YbywLYqaMdR8QiSV8HbiY1H70oIh6SdBKp6ekUUj9GkyU9DvyVlFmYmVmN5B/g7ZM0NldL9Vm+Br4G4GsAvfsatJsRSNoamN8yPrGkLwL7AU8B4yPir7VFaWZmlemo+ej5wOsAknYATiU9BfwCTbhxa2Zm1egoI+jf8Kt/DDAxIq6NiOOBTaoPbcUhaXNJsyXdJ+m9zY7HzKw7dZgR5CadADsDv25YVvL8QW/yz8A1EbFlRPxfs4MxM+tOHWUElwN3SPoF8AqpWwkkbUKqHup1JA2T9GDD/NGS7gHGAV+VdHvTgqtIO+c8vo31VpI0Q9KOef4/JX2/tkC72fKet6S1JD0iabOcfrmkf60p/C6r6rwlbSjpMUmDJPWT9BtJn2q9355G0taS5kgaIGl1SQ9J+kCz4+pu7f6yj4jvS7qN9AzBLQ3t+/sBR9QR3Arixvz37xFxelMjaaLcHPhg4BpJRwC7AR9ublTVa++8I+L13Dx6kqQJwDoRcUEzY+1OXTlvST8AfgTcA8yNiFuaFH63iYgZkqYA3yONzXJpRDzYyWY9TodVPBHx2zbSHq0uHFuR5edAJgPXA9tGxOvNjqkO7Z13RPxK0mdJvex+sJkxVmFZzzsiLszpXyH1WNxbnER6QPZV4BtNjqUSJZ3O9SWLWPKaDGhWIDVa1nP+f8DfgHdVFVBNlvu8JfUD3ge8TOqcsSeo7LwlrUbqSgaW7KiypxtIOp816KXfCc4IlvQX4F2SBuZBc/ZqdkA1KD5nSfsC65I6HjxH0tr1hFiJ7jjvo4CHgc8DF0taudqQu0WV5/0D4DLgu0CvqSYjNaU/nnRuvXLMlL7W+qdDEfFG7gLjHlKHeL9vckiVKz1nSYNIz5LsHBHzJf2QNPrcv7S1/opuec9b0inAl4FtIuIlSXcCxwEn1HMGXVPVeUv6NbA1eewSSftJOiQiLq7lxCqSH6R9IyJ+qjTY1t2SdoqIX3e2bU/iLibMzPo4Vw2ZmfVxrhqypUg6F9i+VfKEnl7M74zPewm9/rxtMVcNmZn1ca4aMjPr45wRmJn1cc4ImkDSm7k30wclXZ0fxCnddqSkPRrm95b0nU62OTg3/+ts30/mZoPFJF0oaYtl2aZh22Nbzd/dlf20sd9vSHpY0mVd2HaYpM93RxzLeNxxuali1ce5UdLarfscarXONElNHZu3MT5JoySd3c56nX5mW3/OljGO0yV1OBpjb+CMoDleiYiREfEB0pgPXynZSKk32JHA2xlBREyJiFMribLzePpHxJcjYm4Xd7HEP2hEbNcNYQF8DfhkRBzYhW2HkR6UWia5jXmX5Pf1S8BPu7qPUhGxR0T8rerjdKeImBkRy9O1Q5czAuAcoMMfWr2BM4Lm+w2wiaRPS/qd0pgHt0r6JwBJ4yVNlnQXMJnU78mYXKIY0/hrv719tCc/XXpL7lHxQkANy74g6Z58nPNbvugk/V3SGZLuB7Zt+fUo6SuSTmvYvjGun0ualY8zNqedCqya939Zy77z3ysk7dmwr0mSRkvqL+k0pV4x50g6rI1z+h9gY+AmSUcp9Rh5UT6X+yTtk9cbptRD5r351ZIJnQp8LMd1VOvSlKTrtbhHztbXYqlrll+TcunvAUlHtfFW7ATcGxGL8n6nSZqgxaXGbXL6uvlazpH0W0kjcvrH87otY2asIWk9SXc27ONjed3GX9ArSbpMqfR0jdoomUr6lKTp+RpdLWmpriMkbZI/b/fn9d4r6Z2SbsvzD7S67g9LuiB/Hm6RtGpe9qG8j/uBwxv2v6Ok6/N0R5/Z0s9Z8fsUEU8BAyW9u433rfeICL9qfpF6MoXUfPcXwFdJfba0tOL6MnBGnh4PzAJWzfMHAz9s2Nfb8x3sY4ltGrY9G/hunt4TCGAQqR+ZqcDKedl5wBfzdAD7N+xjGjAKGAw83pB+E/DRPL1u/rsq8CAwsPE6tHFdPgNckqdXAebnbccCx+X0dwAzgY3aOK8ngUF5+hTgC3l6beBRYHVgNWBATt8UmJmndwSub+v65vnrgR1bX4v2rhnwIeBXDduv3Ua8JwJHtLqmF+TpHYAH8/Q5wAl5eidgdp6eSnqiF1KfOCsB3wL+I6f1B9ZovDakkk80bHcRcHSr93QQcCewek7/d/LnpVX8vwM+k6cH5Gu7ErBmThsEPE760h5G6u9oZF52VcP7MwfYIU+f1nDeb78ntPOZLf2cdeV9InWXsV+zvzeqfPk5guZYVdLsPP0b4MfAZsCVktYjffk90bD+lIh4pWC/QzrYR1t2APYFiIgbJD2f03cm/WPMkATpH+uZvOxN4NrWO4qIZyXNk/QR4DFgc+CuvPgbkj6Tp4eSvngXdhDXTaTuDN5B6v74zoh4Ral/+xGSRuf11sr76ug8PwXsLenoPD8A2AD4I/BDSSPzOQ3vYB/tabwW7V2zqcDGks4BbgDa6pp5PVLfPY0uB4iIOyWtqdTPz0dJ44YTEb/Ov47XJF3nM/Mv3p9FxAJJM4CLlPoB+nlEzG7juPMjouU9upTUs2ZjV+sfAbYA7srntAowvXEHktYA1o+I63Jcr+b0lYFTlIa5fQtYH2gpoT7REM8sYFg+v7Uj4s6cPhnYvY2Y2/vMQtnnrCvv0zPAe9qIpddwRtAcr0TEyMaE/AE8MyKm5KqH8Q2L/1G43472sSxE+kV+TBvLXo2IN9vZ7gpgf1L/NddFROQ4diF1Y/yypGl00oNjRLya19uVNEzqFQ1xHRERNy/juewXEY8skZgGY/kLqRvlfqQuhtvSUW+djdei3Wsm6YP5XL5Cuj5farXKKyx9TVo/4NPuAz8RcaqkG0j3ju6StGvOQHYg/WqeJOnMiPjJMh5DpF/Jn2vv2B04kFRK/FCk/o2eZPE5vtaw3pukL+Plsgyfs668TwNI71Gv5XsEK461SJ2AQccdub1E6g53efbR4k7yjVFJu7O4S+HbgNGS3pWXrStpw4L9XQfsA3yOxV/eawHP53/OzUm/Mlu8ofZ77LwSOAT4GPDLnHYzaaS4lXNcwyWt3klMNwNHKP/8k7RlQ1x/ioi3gINI1Sew9PV9EhipNOrWUGCbdo7T5jXL9fH9IuJaUqd0W7Wx7cMsPQ74mLyfjwIvRMQLpNLjgTl9R+C5iHhR0nsj4oGI+AGp3/zN8/v1l0iDxlzYznE3kLRtnv488L+tlv8W2F5pVEKU7rcsUXKKiJeABZL+Oa/zjnyvYS3gmZwJfALo8PMT6Qb23/L50nKebWjvM1v6OevK+zScVNXUazkjWHGMB66WNAt4roP1bge2yDe6xnRxHy1OBHaQ9BCpuP0HgEitgI4DbpE0B/gVqfqiQxHxPOlLbcOIuCcn/5J0U/Jh0o3YxsGOJgJz1HYzz1uAjwO3xuIBcC4E5gL3KjUtPJ/OS7UnAyvn4zyU5yHVDf9LvjG5OYtLXXOAN/NNy6NI1S5P5OOeDdzbzrm3d83WB6blqsBLgbZKWTeRqjwavSrpPuB/gENz2njgQ3n/p7I4sx+Xb3LOAd7I+9sRuD/vYwypp9jWHgEOz+/NOqTRxRrP6VnSPZLL876n52vV2kGkapk5wN3Au0ldNo+S9ACpDr6kJ99DgHPztVI767T5maXwc7as71POQDYh3Y/qtdzFhNkKQNJ1wL9FxGO5WuPoiOjVXz49Qb7nsFVEHN/sWKrkEoHZiuE7FJS6rHYrAWc0O4iquURgZtbHuURgZtbHOSMwM+vjnBGYmfVxzgjMzPo4ZwRmZn2cMwIzsz7u/wMkP47nPSdfvwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plt.bar(feature_names, feature_importance)\n",
    "# plt.xlabel(\"Partial derivative features (possible candidates)\")\n",
    "# plt.ylabel(\"Feature importance from the sigmoid layer\")\n",
    "# plt.show()\n",
    "\n",
    "def plot_feature_importance_neural_net(feature_importances, feature_names, threshold=0.5, save_path=None):\n",
    "    # split it up\n",
    "    above_threshold = np.maximum(feature_importance - threshold, 0)\n",
    "    below_threshold = np.minimum(feature_importance, threshold)\n",
    "\n",
    "    # and plot it\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.bar(feature_names, below_threshold, 0.35, color=\"g\")\n",
    "    ax.bar(feature_names, above_threshold, 0.35, color=\"r\",\n",
    "            bottom=below_threshold)\n",
    "    # horizontal line indicating the threshold\n",
    "    ax.plot([0., 4.5], [threshold, threshold], \"k--\")\n",
    "    plt.xlabel(\"Partial derivative features (possible candidates)\")\n",
    "    plt.ylabel(\"Sigmoid layer's outputs as feature importances\")\n",
    "    \n",
    "    if save_path is not None: fig.savefig(save_path, dpi=200)\n",
    "\n",
    "plot_feature_importance_neural_net(feature_importance, feature_names, save_path='./visualization/feature_importances_selector.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving\n"
     ]
    }
   ],
   "source": [
    "# BEST-2000: 1e-06 (LBFGS)\n",
    "# For saving the plain neural networks.\n",
    "print(\"Saving\")\n",
    "MODEL_PATH = './saved_path_inverse_burger/export_publication/semisup_model_with_LayerNormDropout_without_physical_reg_trained250labeledsamples_trained125unlabeledsamples_'+string_test_performance+'.pth'\n",
    "torch.save(semisup_model.state_dict(), MODEL_PATH)\n",
    "\n",
    "# torch.save(semisup_model.state_dict(), \"./saved_path_inverse_burger/semisup_model_with_LayerNormDropout_without_physical_reg_trained2000labeledsamples_trained1000unlabeledsamples.pth\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading the best model and testing\n",
    "# semisup_model.load_state_dict(torch.load(MODEL_PATH), strict=False)\n",
    "# semisup_model.eval()\n",
    "# F.mse_loss(semisup_model.network(*dimension_slicing(X_star)).detach(), u_star)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# derivatives_test, dynamics_test = semisup_model.network.get_selector_data(*dimension_slicing(X_star))\n",
    "# derivatives_train, dynamics_train = semisup_model.network.get_selector_data(*dimension_slicing(X_u_train))\n",
    "\n",
    "# derivatives_test, dynamics_test = to_numpy(derivatives_test), to_numpy(dynamics_test)\n",
    "# derivatives_train, dynamics_train = to_numpy(derivatives_train), to_numpy(dynamics_train)\n",
    "\n",
    "# np.save(\"./saved_path_inverse_burger/lbfgsnew_results/data/derivatives-2000-V1.npy\", derivatives_train)\n",
    "# np.save(\"./saved_path_inverse_burger/lbfgsnew_results/data/dynamics-2000-V1.npy\", dynamics_train)\n",
    "# np.save(\"./saved_path_inverse_burger/lbfgsnew_results/data/derivatives-25600-V1.npy\", derivatives_test)\n",
    "# np.save(\"./saved_path_inverse_burger/lbfgsnew_results/data/dynamics-25600-V1.npy\", dynamics_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABMRUlEQVR4nO29e3wV5bno/30SEkjUEiEokIDS1l2vFCRoW2hrRettC3jZePn1FPfW6t5ttbr3gWLrRmpbjbhPtfYqRXf1tFtAsRiP9lAr0h7YtQUEQWytVm2bIApIsJpAQvL8/lhrxcnKzFpzeeeysub7+fAha9bMmnfemXmf932uoqqkpKSkpJQvFXE3ICUlJSUlXlJBkJKSklLmpIIgJSUlpcxJBUFKSkpKmZMKgpSUlJQyZ0jcDfBDfX29Hn300XE3IyUlJaWk2LRp025VHZW/vSQFwdFHH83GjRvjbkZKSkpKSSEif7bbnqqGUlJSUsqcVBCkpKSklDmpIEhJSUkpc1JBkJKSklLmpIIgJSUlpcwx4jUkIvcBfw+8qaon2nwvwLeBc4EO4ApVfTb73Vzgpuyu31DV+020KSU4qza3ccfqF9nR3snYuhrmnfUhZk9uiLtZfSS9fSnlQ6k/i6bcR38MfBd4wOH7c4Bjsv9OBX4AnCoiI4CbgSZAgU0i0qKqew21K8Unqza3ceMj2+js7gGgrb2TGx/ZBpCIBzzp7UspH1ZtbmPeQ8/R3ZvJ5NzW3sm8h54DSudZNKIaUtVfA28V2GUW8IBmeAaoE5ExwFnAk6r6VnbwfxI420SbBiOrNrcxrXkNExY8zrTmNaza3Bbaue5Y/WLfIJujs7uHO1a/GNo5vZD09qWUD4tatvcJgRzdvcqilu0xtcg7UdkIGoC/Wj63Zrc5bR+AiFwtIhtFZOOuXbtCa2hSyc2A29o7Ud6bAYclDHa0d3raHjVO7Whr7wxdSKaER5STHVO0d3Z72p5ESiayWFWXAEsAmpqayq6aTqEZsHX56UdXaXfM2Loa2mwG27F1NWYuKCBO7YPkqIlKXW8cNam6Lz6iWhG0AeMsnxuz25y2p+ThZobuZ9XgdMynjh1FTVVlv31rqiqZd9aHTFxOYOad9aEB7bMSt5po1eY25j38XL9+nffwcyUxw42LUlX3HV5b5Wl7EolqRdACfFFElpExFu9T1ddFZDVwq4gcnt3v08CNEbXJF3HN8tzM0N2uGqw4HfP0H3Zx24UnJWJGW6jP71j9ouPKIE411tce2053T57euEf52mPbB93s1tQ74VYdGedKy+7cN59/AvMefq7f/a6qFG4+/4RI2mQCU+6jDwKnAfUi0krGE6gKQFV/CDxBxnX0ZTLuo/+Y/e4tEfk6sCH7U7eoaiGjc6z4XbqaeHDnnfWhfueGgTN0P3r9QsfMntwQ+6BVrM9nT25gWvOaxKmx9nbY64edtpcqJtU5biY7caqPnM5924UnccfFH07EpMkvRgSBql5W5HsFvuDw3X3AfSbaETZ+ZtymHlzrDNjpYfOj14/aFuBVKLrpczdCMiUc/LwTTri5jybP55VC516/4HTf50+CLalkjMVJwM+M2+SDW2yG7mdAjHIQ9SMU3fS5GyEZNXU1VbZeI3U1paM3doNJ7zI39zFObzbj5966go6fL2Rmx06adCSLK+bQ0j49FgN5Kgg84Gf2HOWD62dAjHIQ9SMU3fZ5EtRYVhbNPKFfkBFAVYWwaGbp6I3dYHpFWew++jrf1hXw1C2wrxWGN8KMhTBxjue2Gb3WrSvgseuo7e4EgUbZTXPVUuiGlu7pkaxwrKSCwAN+Zs9Rq178DIhRDaJ+hGKpqn2SuEoJg6jvj+fzZQdcurPP2L6/Zj6DZ2Fg9FqfuuW9NmWplS7mD1lBS9f0yB0dUkHgAT8vd6kOZGHgRyiW8oAaRMAmQW/shqjvj+fz2Qy4dHdmtnsUBEavdV+r7eaxsifzf8SODqkg8IjXl7sUBrKoBh2/QjFpap+wMeFgUCqCxA+engeHAdd2uwsVkrFncXhjZnWSxw4dGctEMRUEEZDkgSxKd7xSEIom8TsYB3UwiPKeJj4a2GHAZXhjv48bWu7hxGf/nRoOZDYEUCG5YsbC/ioroEOrWVr9GW4776TUayglWqJ2x0uyUDRJkAEyqINBKPfUYbYc5vNjYlWz4QPXcuKmm6iRrvfap9U8/4FrmWo5z9RNi6mRA/0P9qlCckXuNy19WjtjIYvCOJcLUkFQ5iQ9uVypEmSADOpgEIabo5PBdUf7IcHPZSNkVvVMM7LSuP6FY5jSfRXzh6xgrOxhh45k8cE5bHrhGNbPzOxzx+oX+X/stv8BJ9WSCSbOCUfI+CCtUFbmOA0uSUkuV6oEGYzt8ih50Rsbv6cFDK6Bz5UTMvv+CmifkNny+BIjeYd2tHfS0jud6V138/4DP2V619209Pb3ytnR3skOrbf/gTwV0mAlFQQeKcU0uYUIOuik2BNkgJw9uYHbLjyJhroaBGioq+G2C93rjY3f0wIG18DnchAyV3X9xHZ3r6saN/dhbF0Niw/OoUOr++3TydCMCqwMSFVDHki8YcwH5WbAjYqgbsNBbCnG72kBg2vgczm5UVbssd/ucVXj5j5k9umCbvpUSK8zkh1T5jM1IaqbsJFMGqDSoqmpSTdu3Oj7eKsRqq62ClXY19ld9CF2Sm7WUFfD+gWn+25PyuCklFw4C7Y130YAUFUD598dXMd954m2QqajZgxT3rlrwADuZWWUw819KKV7FQQR2aSqTQO2l5sgyK8vmk+hh23CgsexO0qAV5vP89WelBQvhDFg5a90weY9MJSmYQAFhMyqnmn+rzWs9pY4ToKgbFRDuRfIKXd9jkKeHUmv2pUyuAlLNenKwyksDxcbN8rcoD0bn9dlMK1EPoN15VAWgsBuxlMIJ4NUmi4iOgbrCxeEsHz2Y3chNi1kDKaVsDIYbYQ5ykIQ2L1AhXCa4aeG1WgYzC9cEArVaAZ8q0PiXukaF/pe0kp4IM5aCGFPjExVKDsb+DZQCSxV1ea87+8EPpX9WAscoap12e96gG3Z7/6iqjNNtMmKl5lNsRl+uUTGxkmcL1ySqRShx8amVykSSB0S50o3FKHvMq2EV+JaOUUxMQocRyAilcD3gHOA44HLROR46z6qeoOqTlLVScB3gEcsX3fmvgtDCID7mY1Xf+1yIsr4idhVFQnFTgj0bS+kDinE1hXMXnsWL1ReyjPDvsSsinWRvgehFKyfsTBjcLZSVRM4JiCu4MtQ+igPEwFlpwAvq+orqtoFLANmFdj/MuBBA+d1jV3Qi2T/b6ir4a5LJvFa83mBys0NZnIzkrb2TpT3ZiRhCYOyjnbeuiLjUrmoLvP/1hV9XzU4XH9DXY0/dYglqldQRrOLbx/yn6w/d3dk74EroV+gT2yZOCfj2jp8HCCZ/w24usYVfBnFxMiEaqgBsK7DWoFT7XYUkaOACcAay+ZhIrIROAg0q+oqh2OvBq4GGD9+vKcGprr9YEStqilbo3wR9U7BflnrQx0SklHVC1b7xMyKddmArt28KaNg67uZnfyovELwcoprHInChhO1sfhS4GFVtY4qR6lqm4i8H1gjIttU9U/5B6rqEmAJZOIIvJ441e37p9CMJAwjVtkKboeBeecjX+GZnmmF+6VyYFrjouqQkIyqXsgJtzN7fkVz1VJqs1lCR7Mrcz1DamIXVlbiGEeimBiZEARtwDjL58bsNjsuBb5g3aCqbdn/XxGRtcBkYIAgSDGLlwHcaUZSV1sVmhGrLAW3wwB8hO7u16+2/VLAH9+RkIyqXshdy0ce/SK1dPX/srtzoBDI4UFYlborchQTo8CRxSIyBPgjMIOMANgAXK6q2/P2Oxb4v8AEzZ5URA4HOlT1gIjUA78BZqnqC4XOGTTFRLnjKpLUxf5Dh1TQ3tk9YP805YZPHNIttPbWM73rbvP9GmbqCK8sqgPbuH0Hho+DG54fuD3PhXbDB67lsxuOCpaqYhBFKTtFFgc2FqvqQeCLwGrg98AKVd0uIreIiNUL6FJgmfaXPMcBG0XkOeBpMjaCgkIgyZRKZlKvXghO2TD32QgBCNm7x6vhMC78tNPG26VDq1l8MDPo7GjvNHv9IRlVc3h6H5xWITUj3HsA2aS0PvHZf+fMnl/1282Tx41DmuzEPnc+KbtcQ2HhdZYdJ6ZyJkWehM9uBltZDdWHQufe5MzWgsy0t65g5yNf4Qjd3VdEpaV3OgBXHPo7Fsk9eeoSgaZ/gr//lvnrCIDn96FQn4G7GXmRFZUV18+6w2/2rUhKbLUQ2oogJUMUvr6mMOWeGbk7nZ0xtacLOt8iUbM1vz79ABPn8MysX3FC77K+IiqQ6df5VcttdOYKG++L/5rzcHofrl++xX51kF2ddNSMoRehtbeeRXoNq3qmwcQ5rDptNdOGPcKEN25n2hP19qsLp5TWMjCltetnvZBBfRCtFlJBYIhSCoIyNYAHLaDiGTcGQrcDbpgE9MZx6tfazp0OR2ho15xT73zpKzeyc9EHUZcqqULPvVMcyobX9rK3oxuyWoq3Orq48ZFt3LRqm7s4Fgf10uuM7PfZ07PupLIa3hhM4CeMssg1FAVx52vxgkkvhEi9e5y8XPKJ0P3RFgPeOLb96hQrAKFcc069c2bPr7itaul7Xj0ufPmd3occA+JQtq7gxGf/va+AfKPsprlqKXTDg7/9+ICoats4lhn2LrQ7TppPwws1/p51h99kxkJ45Gr7Y+J+/nyQrggMUWolH2dPbmD9gtN5tZQiqu1SB9gRd53ZkFIcZI4X++9CuOacemf+kBV9/v19FJn52r0P+fRbNTx1CzUc6Pd9rXQxf8gKx9QaA1YdDsbvqTOv8f+sFzKoF1otlBjpiiAgVh/lYVUVVAj0aiYR2EVTytAXPkzyfeVrDocDf4Nei/eSiQE3KH58+p3IN0ZO+AS8+mv6uVqGdM25gXas7LbfocDM17rqdFoZ9FstF9DvOyXbs11tm44oLmQMLrRaKDFSQRCAfM+Izu7evu96VFm5qY2mo0akwsAk+S96Ur02TAxIdiknOnZnvIRe+kXo15xT7+zQehrthEGRmW9OveXkQdRvteygTnudkVx26jhWbmqLPuWITf8ffPRavtGynfvfOYWxdfXcddLXmPqn7yTv+fNI6j4aACf3SStpcFWKb4q5LoaM1UZgTf8AOLrDOkXxFo3utXEf7WQoz5/8dabOvCae6GAX7qhJdRF3ouxLVYaBG4+gJHoNDRasg8PwmipEoL2juyTTCNgScy6g99Q71dz4NtxY/RBHshtxmPkWy5tf8H7YqNNqZixkanZ7LClHXLijDpY6GakgCEAxz4jcPinmyR90rKkuBk1Fs4TkAsr04enAbQX3DZylNqy6yH5x6P8d2t8ddTBM9lKvoQAU84xIstdQqVOs/GhSg/k8EZb3UUhEEksTZYqRIik/cgyGyV66IghAvj/+oFRPJJSw1HKJylRZxPsoUW3FfyyN6+soUq/BeH/k9X9HzWgWvnsRLb0f69slypKeic4+GgdJMRanxEcYhvpSyheVxLb6aZOnYwoYz1edtjqS/ohD+Jq8107G4lQQpJQkdi+HFT8vSuRJ9AJQsK3n7o7NpdbrQOmpzx1TVQvThj1SMvfOKyafy9RrKEYim0Uk1ac+BMJQy5VSviinNjW9/SQ89p/eSzsawqt3j6c+L2A83/FG6dw7r5RKzeKUAhRzqTNGEf3pYMRp0MkJ3huWb/EkFJx03EpmVha3Dt6KU1tvrH7IV2nHuOwNnuwKBSJ5xz5ROrm+vBJFHrPUayggxYpvBE1P7bq4xyDKhBiEnOAtmqnShkJeYF5+Jwqcclsdifd0EEH6LCifOnaU++0F8v6UWq4vL0V7org2I4JARM4WkRdF5GURWWDz/RUisktEtmT/XWX5bq6IvJT9N9dEe6LCzQsUZFnn6QUNEHy0oeUedi76IL03D2fnog+yoeWeosfERZiC15r+2Y4kuaQ6paoWH4nQ4qyl8fQfdnnazsQ5majqRe2Z/y0BZ5GmRA+AV8EbxbUFVg2JSCXwPeBMoBXYICItNiUnl6vqF/OOHQHcDDSRWYFvyh67N2i7osBNAE2QZZ2nAB2fwUcbWu7hxE03USNdIDCaXQzfdBMbgKkzrynaxihVCnZqthuWb+H65VtoyJ47qD41p25yquKWJJ2zrWqs0kZ9UlEFXe9mjK02tqM4bSMmzx1W9LHpZ9xP4F3YkdUmVgSnAC+r6iuq2gUsA2a5PPYs4ElVfSs7+D8JnG2gTZHg5iEOsqzz9JL4DD4a9+wdGSFgbZ90Me7ZO4q2L2qVgt0LlBusc+ceXlNle6xXfaqpKm5RsmpzG9OeqOdL7/4jOxmFIpmavyIFq7jFea1J7+cwnnFP73VEAXQmBEEDYJ2Ktma35XORiGwVkYdFZJzHYxGRq0Vko4hs3LXLYdkYMW4e4iDLOk8vic9C5EeofV8eoQ66ZgtRqxSKzRI7u3sQwYg+tRR1zrkB69He6Xxk/7c5vmcZHQzNlPO0kmc7ivNak97PYTzjrt/rCEthRuU19BjwoKoeEJFrgPvJJC9xjaouAZZAJo7AfBO9M++sDxVPr4v/ZZ3b3+/DR66WN2UUoxkoDN6Uep4psiSOWqXgJrdTe0c3d14yKfBS3mQVtyhwGrCGOZW3tNiOPF+rQTflpPdzGM+46/e6kAOIYU9AE4KgDRhn+dyY3daHqlqrRy8FFluOPS3v2LUG2hQJYT/EUbwkfz15HsNzNoIsnVrNr4/6F24u4vYadXlOuxfI7tym9KmxZLz0ieOA1TuSxoritQRcX6uNm3LnI1/k+df2urIp2ZHkfg7jGXf9XkeYfTZwZLGIDAH+CMwgM7BvAC5X1e2Wfcao6uvZvy8AvqyqH8kaizcBJ2d3fRaYoqpvFTpnGllslg0t9zDu2Ts4QnfzptTz15Pncf0LxxSNZowjzUHOcNfW3onQP8407hQLceIUfXrZsGdYqD8cIOifn/INfwO3Q5qHNq1nw+xfD7q+L/SMQ8grmRDqUYQWWayqB0Xki8BqoBK4T1W3i8gtwEZVbQGuE5GZwEHgLeCK7LFvicjXyQgPgFuKCYHYGYTRu1NnXgPZQWF09t+O/37cdl/rzDN/ZjP30N8xv2o5tY/uzBRaD6FvrLPHpCVdixMndcPP5eO823WQ+UNWMFb2sENHsvjgHDa9cAzrZ/o4kcNsdAx7BkVe/nycZu9A+IGiEZbCTHMNecGmipJTpaZSwm5Adao165jfJCl9E5WgTuCEwO4+3rB8i0N2Hni1+TzvJylQtevjXXf7+00HkizoI8tLZfg5S3MNmSBC401UOKXAuGhKg7c6sUnom6jSbCQ0nYedrt1JoPvWcc9YSOcjX6SGA32bcjn6TdqGIkvN4pPIHCUiKtZTlikmvIR39yPm0oFh4ORt8vQfdnlze42obwreu6jSbJRQOg/j7pkT5/D8yV+nTevpVaG1t54F3VfxZOUnjbp8Ruqa7MNXP+nxD14puxWBm+hUx8EuAaUDTVNoZuPJmyOCvik6S4xKUJfQhCAMz7OpM69h1bi/7/ebtxlW20Q243ZY3W14bS/Xv3CMY595du1OOGUnCNxEp4LD8jNC401UGHOPi6BvnGaJ1y/fwh2rX+TJmtHUdr4+8EDTgrrEJgRhuGeG7fIZmWuyw+pu7KbFtB24G7AfF5Ie/+CVslMNuYlOdVx++ozeDYJvNZZLjKkOJs5hw0lfYyej6FVhJ6PYcNLXjPZNoXvX1t7Jwncv4mDlsP5fhCGoS6yWsFfCfubcEFnEcQEvKCt248LsyQ2sX3A6rzafx/oFp5esEIAyXBG4iU4tKCxcGm9MeDxEYTAzNbNZtbmNGzccRWf3t/u21Wyo5LZxbf1+K0i/FLt3D3d9jEOrh7Bo+MpwvXmK1BIuhqs+iMkrKSlG2shm3A6rux06cuC2BCUcNE3ZuY8WK3EIwV3ATAVaxVE60e9AXaitVnfUIEFgbu6dAK9e/m7iXDtzuHo2InLFNeI2XOrY9HUnQ/ly15W09E7vt+tg6AMn99GyUw3l55yXvO9NLD9NeTxEncsnSKZFpzblfiM3uORPO7z0S7F6AQBzD/1dZIm6/ODq2YjAK8npXjutuBI/G/abpdNG3fv8yV/nycpPMrNiHeuqr+OVoZezfuh13HX8S2FeQayUnWoIwo9ONTWAR53Lx22edLs+c2prpUjBGTx465fcvXOaWc+vWg6dyY31cPVsROCV5HSvK0XosdESJNotMmhcR566dyrwAPdw4rP39sVLNLCbhm03w9GHD8oAybJbEeQThsHHlI9x1Cl63QxSTjPJTx07yratdoNKPn4GGaf03rUusm3Giatnw9H7SI3lpHe61z2q0aeFDppzP4QV1NQ/fadf0JyJ3/RMhGmoy14QhIGpATzq8ntuBimvAWiF1DgQbJCxFeI+SjVGiatnw84rKYehwcDpXlvvWyQlH00MdmGsoJIQKxJh4GJZqobCxqTHQ5Qpet0EyfgJQMv/zZzBuGgAnx8SHuvh6tno55VkE69gQNVV6F5HmhbaRGoSA3Ed+erOyGJSChGhMCpPQRCB3i3JOdadcDNIebVbRB54E9C1MwrsbFQ3LN/Sv29yeutFdQw0sRN4MEhMQJSJwS6g8LdzmV1YfRHNVUsZ0rPf128awUnASUVmDDP4TJed+2hismSWKHHUIBisuOrLEHLSJwpT1xdgcufk+nzFob9j0SEhx6QUwm6syuFzzErdR3OUUMKwJOLabhFR0e1SxpUr6SCPYjZ2fRPnsOq01Uwb9ggT3ridaU/Uu46IdlJ33v/OKRlhtKg983/UE8Wca6tUDvzO8JhlRDUkImcD3yZTmGapqjbnff+vwFVkCtPsAv5JVf+c/a4H2Jbd9S+q6qdchnvC1rslME+9aYqqvZKSpjnh98KVK6lhVdeqzW1seXwJV3X9hLEVe9hfM5rac2J0rTV0fUEioqN207bF6VmdOAceudr+GIO2gsCCQEQqge8BZwKtwAYRaVHVFyy7bQaaVLVDRP6FTM3iS7LfdarqpKDtcI2D3q21dySXNK8JpidNygAYN+VUmyAArgcgQznpV21uY93Pvs8tsoTaikzpytrO1zn46LWZgSBOYRDw3G5jYOyIPZNosWc1giSHJlRDpwAvq+orqtoFLANmWXdQ1adVtSP78RkyRerjwWYpmius4SWS1pbBpHYKotopM9c7v0QdJ3LH6he5nmXUWuoXAxmDaFL6xedzFySIM2o37QEUe1YjUA+aUA01AFZx1QqcWmD/K4GfWz4PE5GNZNRGzaq6yu4gEbkauBpg/Pjx/ltrWYr27mtlR2+mhmsur4jbWYQtSRgATVBshlJM5ZKENM0lcC+i9tzZ0d7J2KG77b9MQr8EWMUFVe/E6uVX7FmNwBMuUvdREfkM0AR80rL5KFVtE5H3A2tEZJuq/in/WFVdAiyBjNdQoIZkl6IfWPC4bT1X33lVkjAAmqDYDKXYy5oEX/4SuRdRDkBj62rY0VFPo9gIgyT0SwCVYuzqnSC4eVZDLllpQjXUBoyzfG7MbuuHiJwBfBWYqap9sduq2pb9/xVgLTDZQJtcYbTc3NYV0PXuwO2l6OFRaIbiRuUSQ92GASTM2yYpOf7v4lI6tLrf9oOVw5LxjAZYxcWu3glCAp5VEyuCDcAxIjKBjAC4FLjcuoOITAbuAc5W1Tct2w8HOlT1gIjUA9PIGJIjwdgswsnft2YEnHN7YoyTrik0Q3H7skZUdNuRBAWWJSvH/+dZ/PiQ5HgNWQm4iivFIE4gEc+qkYAyETkXuIuM++h9qvpNEbkF2KiqLSLyS+AkIBez/RdVnSkiHyMjIHrJrE7uUtV7i50vUEBZHtbQ8uE1VYhAe0e3N33tYAv6KRR055T2oFSvNQLiqCtRkpRAsGcY2YqjxCmgzIiNQFWfAJ7I27bQ8vcZDsf9NxkBETsK7Ovsdl+/2EoJGCY9UWyGErf+PyJMvfRR15UoWRIwMy5EUlZ2YVCeuYay5N9Yp6IpRW9yiRgmPeGk2kn4y2oKYy/91hX8ZthXOEJ3sUPr+3moJTrHf1zErVIsQJBYhaRT1oLA7sbm42rWlgQvmShJ4Mtqesle7KV3XXf4sesYTScINMpumquWQjc8WflJ4x4tJau2SHgEeI7BvLIra0Hg5ga6mrWVySw5qYSxZC/00rs+n42HVa108ZXqhzh91heNDtIlq7YogQjwHIlIRRES5Zd0zkKxG+jJg2jiHGMJqpLgajgAtxGf2f10UR07F32QL33lxtCvwVSNaCuFXItdn8/BRjSa3cYH5zD6IBJKIAI8R9SR4FFS1oLA7sbmitnH5YccpIB8aLitImXZT1BGs4vbqpYy5e0nQ72GMJbshV561+eLsFpayaotSsjRwlesQolk4S1r1VBiinNYSKRBym3Ep4MqZP6QFbR0TQ/tGsJYshd6Nu5Y/aK780VoO0qi2sKVzaLEHC0GxCpsXQF3OqiES0jtVdaCAJIXhJLImZ3bWZvDfmNlD4DtQGWCsNILOD0brs8Xoe0oaSkWitksckKi6e3zaa6+t3+heMPC0kiskB3FBvokZOF1SdkLgqSRxJmd61mbw347dCQAlSIDvjNB1Cs7T+eLyMMqsj5w6eFTzGaRExJtTEe74MtVKxgrexDDwjJfILV3dvd9F9Sg3vHzhdQWGuhLSO2VCoKEkbSZHeBexWGzXy7FN0BPiGVRo17ZJW0lCRG0yYOqo9DKNl9ItPROp+XA9Eyk9Q1mI62LuYj7Vbuu2tzGzI6d7xkVreQGeh9qr7hcgFNBkDCSaLdwreLIft75yFc4QnezQ/un+G4YBG52gwnPg46DqqPj5ws584n6fr9TaGXrR/3pd4B0o1L1o3a9Y/WLNOnIwplcPdqI4nQBTgVBAknibNO1imPiHJ7pmZa8VU1KP3wNOg4qjWEdO2k70Nnvdy6a0sDKTW22z4BrY3uQtlp+s5htyo/adUd7J4sr5tBctbRfoZ8OraY2N9B7tBFZVy8zK9Yxf8gKxspu3nx0FFTemvg01Ckp/SjplMBlgq+4AweVRs4GZP2dp/+wy/EZ8OqPHyRGwu5cbs9biLF1NbT0TmdB91W09tbTq0Jrbz2Lqz7ff8D2EF+UW5nMrFhHc9VSGit2UyEwml327toGSVcECadU0wYkclWT0ocv77QiNqD833F6BryqP4OWobSey5TXUM6W19I9nZaujOqzpqqS287zn0Mzt3qZP2TFgHKiYXsbpYIgwZRs2oCUZOM3EZ6NqmPxuxfRcuCUAbsWU7d4mSjEXobSxlNq9uRMX5icpN11/EuM3bSYsURfTjQVBAkmkcFlKaVN0ER4ebaiSZvbqAnZHhSrJ10BT6nZk+eYew+3rmDqtptBCqxyQgyySwVBgonSuyKlTDCcCC8KL7ewzuHqXYkqKMzuPFZCzmZsRBCIyNnAt8lUKFuqqs153w8FHgCmAHuAS1T1tex3NwJXAj3Adaq62kSbBgNel8SpKinFkT71ho1fO8ES4UVhDzJ9DtfvSlRBYYV+b/i40LMZB/YaEpFK4HvAOcDxwGUicnzeblcCe1X1g8CdwO3ZY48nU+P4BOBs4PvZ30vBe7bDks1AmRIu/ZIGOhBybp+kZdR1/a6YShxYLPmc43nGBc5m7AYT7qOnAC+r6iuq2gUsA2bl7TMLuD/798PADBGR7PZlqnpAVV8FXs7+XjwkLFOgVzfMROYpSomfmNUOScyo6/pdmbEw0z9WvPaXm+y9Js4TABOqoQbAOtVoBU512kdVD4rIPmBkdvszecfajnIicjVwNcD48eMNNDuPhGYKjNK7ImWQErPawWn2veXxJcxeuzKWYk6u3xUTiQPd2BliLm5VMsZiVV0CLAFoamoyn7SmhDIFOhGHd0USjdNJbFOsOOa8yaodQsZu9j2zYh3zu5fCvqy/fMQTL0/vStDEgU4quXwBHWMJWBOqoTZgnOVzY3ab7T4iMgQYTsZo7ObYaCihTIFORB3Rm8QlfxLbFDsxqx3sVqQFg6YiILJ3ZesK7DPTkaiaCyZWBBuAY0RkAplB/FLg8rx9WoC5wG+Ai4E1qqoi0gL8l4h8CxgLHAP8zkCbvFNiBTKciDKiN4lxDsWMgGW5UohZ7WA3+87VqBhAhBMv1++Ky9Tbtjx1C2CnwJDIBLEbAguCrM7/i8BqMu6j96nqdhG5Bdioqi3AvcD/FpGXgbfICAuy+60AXgAOAl9QVeecsWESYTWpwUJUxmkvqh6nc+dWBmXrWmundggywFkodn/s4gD2y2hqO18f+GNJm3i5tB069oGjYNNEqZyN2AhU9QngibxtCy1/7wf+weHYbwLfNNGOQBiYNZWbbjoK47TX2AinNlWKJG71EiuGnCPc3p+BJR5vKY2JlwvbYcE+KGSfSRBp9lErHjIF5lOOummvcQ5+8Bob4dQmp6I4ZetaW2iA84Dv2JWJc+D8u7MDomT+P//u/vV+k+DK7cJ2WLAPYrbPuKVkvIaSThL15WHjNvQ/yErJq/rJqU1ec+APegw5RxS6P0Xvu5OXTJJcuV3YDgs+ozHbZ9ySCgJDxBnMFadKqpjBLWjaC9fqJ4u+e/bwRmafO/BlS4vlWHAY4Fp7R3JJ8xrXz5DT/RleU+X/vifJlduF7bDoMxqjW6hbUtWQlQDLUaeZpdEZp037kq6SCpr2wpX6yUXkZlosJ48ZCzlYOazfplxtAS/PkNP9EcH/fS+yWok0XUUxFRbRqEjDJl0R5Ai4HA09mMuhfVv0Gjq7+2flSJJKKuhKyZX6yeUMMi2W8x6reqaxrvsqrmcZY2XPgPrSbp8hp/tzw/Ittvu7uu8F1DGxJFYs4nE1e3gjDVOv5foXjgm8Ko9rdS/qYERLMk1NTbpx40Yjv5Xr+OUdn6Oxwq4Qtfvoy5tWbePB3/6VHlUqRbjs1HF8Y7b/ikX9uPNEh6V8PdO77h6wXYBXm88zc24f5Pq1UL3YBlMP+qI6HH21F7UH++0oCejO6WUQmda8pmgt3yDPkNPvN9TVsH7B6YUPzp/0QEYdc/7dTHui3v/vmqJA+4KogPKFHGSrnhlcuYrIJlVtyt9e1qohq1plrASrCrRqcxsrN7X1eaf0qLJyU5u5ZatDO8ZW2AfmxGkEtfZrIYypsUxliIwTN4nJCrBqcxvzHn6un4pw3sPPOfatm5l5XW2VhwvoTyB1SQF1TCISKxryuMonzuzBZS0IrB2/Q+vtd3I5mIR+Ex3asb9mdOL0k3Z94YSRPgrTRS8qN8aAg8vXHttOd0//VVF3j/K1x7bb7u9mohBEWRDYJuPgym1t98yKdayrvo5Xhl7Ob4Z9KToX05DS0cQp5MpaEFg7ePHBOXRodf8dKquh611Xg0DoN9FhsKs955bEGUG9XnPgPnJh0PNFwFm6J3wMLlaj6d6Obtt9nLbbzdgHnLrT/li3zJ7cwPoFp/Nq83msX3C6kWcy1+6ZFetorlpKY8VuKgRGsyu8e5NPSCvQSBxOHChrQWDt4Jbe6SzovorW3np6EagZkZkSdb6Fm0Eg9JtoN9h9+HJ46hZmP3oC64dex6uXv2vshQuCk0qhwiH3lpE+ChAM6EhIKgBbPA4u+d5iXrHO2J1IYoxFrt1fqX4ovqR1Ia1A7YSzkFHzhe0dVdaCIL/jW3qnc6Z+j5ZZ26H6EOjNmxEVeNAicSGzDnYzFsJz/xXNbNUjTiqFoUMq3PVRCUWVGsPj4OJW/VZX46znz83Y77pkUuLUi4WYPbmB0QSz6QXCwArUzgU2XzgL77lAhO0WXtbuowVdEx/1NghEUcS7H0kKusnDSaWwv7uXOy+ZVLiPSiyq1BgeI1DdqNOqKoRFM08oul/kz64J4s4WHCBIrJgL7OzJDbZeV2G6hZelIHDlZufjQYvUTz3B9RMKRVoW7aOIBJyrZyDqjLQeBpdCyfV6VT0P5iUXY1HC2YLdpKOJ2nBcdqoh15G4SU8WlWCXyUBqsggEnOtnIIgKIGT1llMf/685H3ZvnE2KCs4PYTkIRICbQT5qw3HZrQhcJ4dLYLIo6yx27qEXcVPlDxnSs/+9HRIiqFyrGuwCqCJY8ntKEOhHBRCBeiuwOidJKji/lEAOHzvc5M+Kuuxs2UUWT1jwuFMMaqyRuMWwizq8uPq/ueWQldR27nQlqBJVL8EpOvPDl2eM4IajNq2E/gw4RIFHVSPYFaXQxkGK2wjiMN5Xp8jiQCsCERkBLAeOBl4D5qjq3rx9JgE/AN4H9ADfVNXl2e9+DHwS2Jfd/QpV3RKkTcWIophKGNjNYh/u+hi/qZ3B+kXFQ+tjydFSCCdbwEu/yAz6Ia7EQn8GEmy/6aMU2jhIcbuai9JuE1Q1tAB4SlWbRWRB9vOX8/bpAD6rqi+JyFhgk4isVtX27PfzVPXhgO1wTdRLLlMENR4lrl5CoYEo5CV/wWfARPnGuD1a3FAKbRzEJM04H9RYPAu4P/v3/cDs/B1U9Y+q+lL27x3Am8CogOf1TammI/ZtPMoaBP9f5wWsq76OmRXr+n0dW4WuGI3djs9A5XozkcQ2jgYHK4ex6N2LokmdjItUzXE5Q5SygXoQE3RFcKSq5ipQ7wSOLLSziJwCVAN/smz+pogsBJ4CFqjqAYdjrwauBhg/fnygRidNGrvB10rGooevEGiU3TRXLYVu+tINx6YSS6L7nynX1TxHg46a0Sx89yIe7sqkCzemlnNYvbhSA8bhDDEYDNSDlKLGYhH5JTDa5quvAverap1l372qerjD74wB1gJzVfUZy7adZITDEuBPqlo0RtxkGupSwrPxqEjqatMpbj1jQg3jAydj3QuVlyIhpLMOlJLZiaSnarYjNVDHjm9jsaqeUeBH3xCRMar6enZQf9Nhv/cBjwNfzQmB7G/nVhMHROQ/gf9ZrD3ljOeVjFPqatlDpQgXTYl5ZRST+5+TveSNyvpM8rJ8AqqrQgkOKrB62dF+u/nzmSA1UCeWoDaCFmBu9u+5wKP5O4hINfAz4IF8o3BWeCAiQsa+kE4LTOIwgO3QkebrJZQQTgPibV3/MEBv3slQvrTr/EB6/VCCgwoMqnFmsSxIwoIg/ZS8jLRMZoQEFQTNwJki8hJwRvYzItIkIkuz+8wBPgFcISJbsv8mZb/7qYhsA7YB9cA3ArYndErqQbAxCObq0kJ0RS+ShtOAuPF9Z/ZFqypCm9bz5a4rebR3eqCkX6EkJCwwqCa2hm6CovX91PpOen3wIJRdQFkQfJWSi0kPnn/+3vbWAXVpIfmBdCaxls+0ZnaEgffRqtefWbGO+UNWMFZ286aMYvSFt7KqZ5one43x4KAi5RKTEDxo24bK9YmI1vdjtwnF1hMxoQSUlRueffFD9JJw/aJn9fAfd3iIY1cXRES+EFfeS/NrVzt5h0UINFct7ct9P5pdHHz0WtZ1X0Vb18cAd15Axj3Vinj9xO0Z5+i5dOE0ZhczDEcwefJjt3F7TBKEsFfKLulcEDw/PCEVNvGzRE2KuiAu1ZqdEM8JAbsEbTkBOX/IigEFUIb07Od6lvXbFoeabVXPNKYduJsJ+3/KtAN3s6pnWqTnL4Tv0q0RVYXzY0dxc0ypqo9SQeCA3YBV7EHIP0ZD8pLw85IVDaSLINDH1UsSUju8CvGc4Bwr9gVQxsoe178VBkkfcHx7SkVUFc7PxMjNMXEWoA9CqhqywWlZe9GUBlZuarMN6rI7ZsfQkTTYDSQxuSM6qgsiCvQpqloLsR1e8wvl+unNR0fZupTu0JGuf8stXlQKiUsZkofvfE4mJk8uVEt+sre6OSbOAvRBSAWBDU4v2dN/2MVtF55k+yBMa14z4Jjbu+fwH9X3UI1le0WVby+J3EDhZN73PRBFVAym6EsSYjv8RGbPntwAlbcOMMoerBzGXb2X9ts3qJrNa1LApA84vnN6Bc2B5GEyUciO4iSUi9leSjWpZaoasqHQS5ar85pf/MPxBdS8iu3iUMG9CFZVgB0XV/83T8rn/alUIgr0KapjDbEdvnNM2RRAGTLrO0y/4PNG81V5VSkkNlYgi+/+DupiakC1FETtlhRbnFfSFYENfqS63THzh6ygWg7237Gny9cMt1Cx8isO/R036VKGdGaL1HhVqdQcDp1vDdxuONCn6Cwx5IyYvj1pbCKgZ2M2fbcfG0bSs+j66u+gOZAMTCaCqN1Ksv4zqSCwxc9LZneMnUER8DXDdRoQBFh0yErYt7//F25VKltXQNc7A7cHUGE5UfQlSWIiuojwa8MotQHHFUFSjxiYTDitup225xO3664fUkFggylD0n4ZTW3n6wN39jHDLThQBJkFPXVLZpWSz9DDQgn0KfiSJKA8aFw+4H5tGKU24ISOgclEpQg9NoG2lT7VuqVAKggc8POSDThm6y3GZrgFB4q1AWZBDsKit2MvH29eE/0sM8Y6tHFWcRvUM/woMTCZsBMChbYPBlJBECYGZ7gFB4rKALMgh6X0Dh3peSAsqYhKi4thR81oFndfwo/fOWXAblG6ZJbrDN/4cxNwMtHgsPpuSIghPgxSQRA2Bme4jgNFEIFjs5S2S0xX7MUsOpuOO+eSlTwXw9rO15mv3+etiq5+eZhyJMUlM1ZCun+Jq6VNvIb4uCZTqSAYLPgVOBYh4pSYzs1AWNDTIlcCMimVqWxcDGuli/lDVtDSNVAQJMUlMzZCDPRLYmBcXGq6OIViKghSjCSmK+j+GFHAmmsKFOzJJ2kumbEQwv2zZoK1I+5VWBxqOieheP3yLdyx+sVQhVEqCBJIIj1XiqgGQvNqCoMCdhErdllJS5HAz5Ph+2eXzj2fclyFFRJ+Ya8OAkUWi8gIEXlSRF7K/u9Ur7jHUpSmxbJ9goj8VkReFpHl2WpmZU2cycQco0Fzqp0CGSELRlS6rEwVWWbSIgV7aqoqueuSSbZZSUsNI8+T4cpihYIjoXxXYcWEX5jJ64KuCBYAT6lqs4gsyH7+ss1+nao6yWb77cCdqrpMRH4IXAn8IGCbEod1RlZXW4Uq7Ovstp2dxa0ztV0S31lcNRDUqylS/WiecT3nNfTYgVMSuQoIMqM38jwZDvQrNPNNYv9HxbyzPsS8h5+ju8fZTTUslVlQQTALOC379/3AWuwFwQCydYpPBy63HL+IQSYI8ge4vR3dfd/ZDXaJTCbmUjUQxKspcgFoMa7XknnwFpk/S2CCCkgjz5PhQD8nNWIpVfoKjSKhCmGpzIIKgiNVNRc6uxM40mG/YSKyETgINKvqKmAk0K6quWQ8rYDjky0iVwNXA4wfPz5gs6Oj2DLYutwLJbOoCUzkACri1ZRIAZgAggpIY9kwDbpBl0KepDAotrK7Y/WLdPc6S4Iw+6ioIBCRXwKjbb76qvWDqqqIOF3FUaraJiLvB9ZkC9bv89JQVV0CLIFMzWIvx8aJm4EsN8tzEhheHoBQDM0R5AAq1fS9YRNUQCZx0C21KGoT75SblV2cKrOigkBVz3D6TkTeEJExqvq6iIwB3nT4jbbs/6+IyFpgMrASqBORIdlVQSOQjPJKBnEa4KxUijgKgbqaKkTgBhcuZKHp2V2qBoK8MGENWCUV7WxDUAEZxaDrp49LJYra1DvlZmUXp8osqGqoBZgLNGf/fzR/h6wnUYeqHhCRemAasDi7gngauBhY5nR8qWM3wFmpqaosqDo6cLDX9UMYqp69iGog6AsTxoCVxKhVr5gQkGEOuib6OMnC2tQ75WZlF+fqLaggaAZWiMiVwJ+BOQAi0gT8s6peBRwH3CMivWTcVZtV9YXs8V8GlonIN4DNwL0B2+OKKB+8/AHOzmvIKbDGbqVQ6CH0okYw3QcmXphAA5ZNnMMdq+sTF7XqlaSrUYLedz+CJMr315Ttys3KLs57HUgQqOoeYIbN9o3AVdm//xs4yeH4V4CBmb5CJI5ZopsBzm4m4LRScHoI3aoRwuiDWI29DikQmt79R9oo/dxBSVGj2A3AQe+7V0ES9ftrynb1qWNH8ZNn/mK73Upc97rsSlV6LQnYx9YVmRKQfkpBFsEpkMsp26HTQ+i2TJ7vPiiA79KJJvrVIQXCjdUP+WtTkgjxufOCU2BaXW2V7f5u+9irIAnj2S2EqdKTT/9hl6ftUVN2KSZ8zWBCTLqVw2kmkJv9zKxYx/whKxgre9gvozO1DmyKcUPxpWUYs3df+k1T/eoQ53AkuwesrOL2mPFECM+dX7WK0wA8dEhFoD72OuOOeuVpSl3j1L629k6mxVH3I4+yEwS+lnoxJU3LPRhbHl/C/O6l1Eqmklht5+uOA4KtQMnTn8899CLb3PtBZsq+XhhT/eoQ5yDDG7nttJMSq18viuHnLohaxWkg29fZzZ2XTIrMWywON2MT6ppC3oNJcGIoO0Hga+YaY9K02ZMbmL12JezLKyfppSZx3qzypsof8k71QR7u+ljfbiZmyrMnN2TyEuWEztrGTHoJpzaa6tcCcQ6zJ4agc42qtoLh5y6IYbfQABxkoPQ6gUhiXIQbinkPxu3EUHaCwNfM1URkbRCC1iTOm1UO6dnPLYes5De1M8zOlL2qMkz1a5S1jiNQE/Zh+LkLolYJcwD2IkiS7kXlhLXdSUy9LVqCdTibmpp048aN0Z0w/+WHzIzz/LuDv/xuZpd3nmg/IEglXPDDwm1YVId9AhOBRe3+253FqnP+zbAvMRob49fwcXDD8wO3h9mvYeFwL3Yyio/u/7bZgclw/0xzqDfhNmApyf7+pUTQ+xAEEdmkqk3528tuReCLsGachWaX1vPVHA6V1dCTpx7SnuKz0RBXM/k65yN0F4jNjpaVS//BpJ67TvoaU//0nWSUsHSDwyrsCN3dz5sGDOh7DT93QWf1QVRAqRB5jySqt9IVQZw4zfRrRsDBzv4zwYqqzMCvvQP3d5pxQ6iz7vyZzbrq62is2O3YPruCJDVVlZmaB6UyKDjcs9beeqZ33d33udDsLs5BMY5zD4r7XgSv/XrTqm08+Nu/0qNKpQiXnTqOb8y2DbcySroiSCJOOv7OtwZu6+0euK3Y70Co+vMd7Z0Wt9bdtHMoB7SSoWIxiFXVsOED13K9w3I4biOZZ2wM09aiNjmc9L1xp72II2Ap7hobYeP1nq7a3MbKTW30ZCfhPaqs3NRG01EjUmNxWeKktvHzO4UwmELYytxDf9fPrXUE79ClQ2jnMOp4B4Y3suED1/LZDUfRme8GacGvkSyWmXWeYN1JPbd2/wMtvf0jmJ3cGQf7oGiHad//pKmZvN7TJD4DqSCIEye3xyE19qsCO5WR4XTQXphftZzag/3tFtVykIM1h8KXM6uU65vXFBQC4M8H3O0sLJRBwyJYn9ncxpOPbINed/recqy7YNL3P+4VlR1e72kSn4GySzGRKCbOyejqh48DJPP/+XfDObcPqKlLVU1mu93+MRlXazt3Ft1e7OH2ayRzk2ogivrPTulBnAYl36k4ShhTaRog+hQTbvB6T5P4DKQrgrgppLZx0uvH7FWTm2Uv7x3pYBx+T1VVKKIySLENNyH7US3BrXr3XN/csHyL7QokiR4jYWPS9z+Js2mv9zSJz0AqCJJKSHr9oFiX5osr5tBc9Z6NABigqnJ66IN6jLgJ2feavTUobtQWpRoQFRRTRuokVrLzek+T+Ayk7qMpnsh3Ge3zGqrYQ0UIlcucsHNJzKdSpM8zw0pYgTtxBgqVC+XgihomqfvoIMM6uA7PlrNs7+gOfXaRP5tu6Z1OS9d0BHh10Xm2x4ThsugmZL9HNdLso0lUWww2kjibHgwEEgQiMgJYDhwNvAbMUdW9eft8CrjTsulY4FJVXSUiPwY+yXuF7K9Q1S1B2lQO5M+K2jvfizEolUIdJsgJmEIz8ZytIIpBI0l9Mxjp7u6mtbWVDw3bz9JZYyzfvM3vf/92bO1KIsOGDaOxsZGqKvt6EfkEXREsAJ5S1WYRWZD9/GXrDqr6NDAJ+gTHy8AvLLvMU9WHA7bDP1FlkjSI1QhqDejaofUsPjiHlu7pofkkJ9HQVahNUQZQJbFvBhOtra0cdthhHH300YjY5TKJlr0dXbyxbz9dPb1UV1Zw5PBhHF5bHXezUFX27NlDa2srEyZMcHVMUEEwCzgt+/f9wFryBEEeFwM/V9WOgOf1jK2eunJ9dJkkDZJTNcysWNfPWNsou2muWgrd8Fj7wBKNJkji0jwpbUpKOwYr+/fvj0UI2A34AG17O+nN2qC6enpp25t5L+MWBiLCyJEj2bXLffWzQMZiEWlX1brs3wLszX122H8N8C1V/T/Zzz8GPgocAJ4CFqjqAYdjrwauBhg/fvyUP//5z67b6WRg2nTo9ZkiL/kUyt2TAHKqEKfcPq299VxS+6PgBsoSXC2lDF5+//vfc9xxx0V6zr0dXf0GfIAKEUSgp3fg2FldWcGxY94XZRMdsesvJ2Nx0YAyEfmliDxv82+WdT/NSBRHqSIiY8gUsV9t2XwjGZvBVGAEBVYTqrpEVZtUtWnUqFFOu9ni5E8+zCEgKoqCM0HIBeiMFRsffmCs7Amujsglq9v3V0DfWy3FVDN3AAmp5ZsyuHlj3/5+QgCgV9VWCEBmZVCKFBUEqnqGqp5o8+9R4I3sAJ8b6N8s8FNzgJ+pap9lU1Vf1wwHgP8EBtZPNICjN0fvSPsDDBecWbW5jWnNa5iw4HGmNa8JHNmai2Z9U+wF4v7a0cHVEYXKJMZN0oVUinFWbW5j5779bG1t5w+vv83ejq7iBxnAbmB/4/U2vvRPl3P+x6dw3rTJ3H7zArq7Mu2prnxvSN2xYwcXX3xx0XOce+65tLe3+2rfokWL+I//+A9fx1oJmmKiBZib/Xsu8GiBfS8DHrRusAgRAWYDoehjnLw2llZ/xj6Vg8HcPWGlOZg9uYHRF95q2/7acwIM1rmZtlMyvH2txgWbZ5IspGIk9vsSErl36GBvf318MWFgoj+sAztkDLE3fO6znHHOeTy+7llafr2Rjnff5TuLv06FSJ/94ODBg4wdO5aHHy7uB/PEE09QV1fnuW0mCSoImoEzReQl4IzsZ0SkSUSW5nYSkaOBccCv8o7/qYhsA7YB9cA3ArbHFqdcJ5POuzr03D2h5kZxylXkt/39Ztr2dNSMDj1/T1FirCGdVKLIqxQXdu9Qrypv7NvveIyp/jhy+DAqLMbp363/NUOHDeXz13yOhsNrqKmuYt7N32TVip/yfx+6n7mXXszpp5/OjBkzeO211zju+BP4w+tv89uXdnDOzAv40HHHccEFF3DqqaeSC4o9+uij2b17d2b/447jc5/7HCeccAKf/vSn6ezMTHh+9KMfMXXqVD784Q9z0UUX0dFh1t8mkNeQqu4BZths3whcZfn8GjBAV6GqkYRbFvbmCDeVQ+hBRiZTUfz8ywNn2lYqq+nqfIftFZewozrrqto7PfoUunHXkE4gSUxtbAqnd6WQPt5Uf+Q8gHJeQ6+99CKnNDX1bT+8thrGvI8JRx1FdQU8++yzbN26lREjRvDc7//IwV6lq6eXFfffyyHvG87KXz7Dvh2v8MmPTrU930svvcSDDz7Ij370I+bMmcPKlSv5zGc+w4UXXsjnPvc5AG666Sbuvfderr32WtfXUYyyiSz26k9uKi1CyQQZbV1hn/oaUIV9chiH9XZSx99A+ruqtvROjzZ61il9d4TpuJOWEz/qqOYor9/pHcpX21gx2R+H11b3Dfxrhw+j461Kx33PPPNMRowYAcDuvx0g55W5ecMzXH7lP9OrysjGDzBx4kTb4ydMmMCkSZMAmDJlCq+99hoAzz//PDfddBPt7e288847nHXWWZ6voxBpGmobTC6zTabgDZUC+vU2reed3qFU6sF+22uli/lDMgbaSAWbaZWYR5KohokytXHU12/3Dln18Tn2dnTxh9ffZmtrO6MOG2r7W0H74/jjj2fTpk39tr399tv85S9/YciQIRxyyCF927t6vHsWDR36XrsrKys5eDDzzl1xxRV897vfZdu2bdx8883s3++sFvNDKghsMKnXz3n4XHHo71hXfR2vDPv/2HTo9ZlgtiThoF9XhcUH5xR0VY1FsE2ck4n1WNSe+T/C+IYk5sSPcsIR9fXn3qEhFRldfXVlBQ2H1/QL3Mr5++cG2c98ZDxDh/Qf3kz0x4wZM+jo6OCBBx4AoKenh3/7t3/jiiuuoLa2tt++1ZXv2RYmTT2VXzz2MwD+8qc/sm3bNk/n/dvf/saYMWPo7u7mpz/9aaBrsCMVBDaYXmbPrlzPIrmHxordVKCZILakuTs66Nf3cmhG9aP1tt+/KfVll/kxicnlvBbICUIc1z97cgOjhw9jYmMdx45534Do3Xx//9M+dARf+NQHOOKwoUb7Q0T42c9+xkMPPcQxxxzD3/3d3zFs2DBuvfXWAfvWHza0Lwp6zmevZO9be7jg9I9wz7du5YQTTmD48OGuz/v1r3+dU089lWnTpnHssccGugY70jTUNhhPJ+zkjpmkCOb/86+w8T6sMYEdWs2C7qto6Z0+IJ0FkNHLx1ghLS7KPd10XNdfKLJ4a2u743ETG+vCaZALcukpOru6qejtYfyRdbz1+l8544wzePHFF6muDi8dhdHI4nLE+DI76e6OW1fAc/+FVQj0KjzU84m+ouwtvdNZqFfTUTOGJJTJjJOSsfuERBKv38lwXMigHAWH11Zz7Jj38cER1Vz5D+dy2kencsEFF/D9738/VCHglbLxGvKC8eRhSXd3tAnQqhCYUbGFm7OfG+pqmH7W56md/M3o25cwyj25nOnrN+GBdOTwYbY5gfINynFx2GGHkeRiWmUjCLw+bEbTFyfA3bEgDiuTsbIHAIFQlvxJc8EsRhLam4Q2gLn3w015Tzfk+/u7SQud1DTScVAWgsDUw+abnPokqZk8HVYsOzSTi6mu1l1xCy/Efk88koT2JqENpjEZCGf19y9GflbRJKWRjoOysBH4dnczmeEyRnfHosxYOCBnUYdWs/hgpo35/gQmcrgk0QWzEElobxLaYJq4PLCcsooWSlsxmCkLQeDrYSunDJfZAK3W3np6VWjtre/zFgLYZymFaSqYKIkumIVIQnuT0AbTRBkIZ8UpqKtU00gHpSwEga+HrdwyXE6cwyW1P+L9B37K9K67+4QA9O8nU7PSuAYAvyShvUlog2ni8kAK4mX0xhtvcPnll/P+97+fKVOm8NGPfpSf/exnppsYKWUhCHw9bEl3+QwBN/1kalaaRBfEQiShvUlog2lCDYQroNrNzyoK7ryMVJXZs2fziU98gldeeYVNmzaxbNkyWltLe1woC2OxL3e3pLt8hsDsyQ1s/PNbPPjbv9KjSqUIF03p7x1iKoleqblgJqG9cbchLI8lox56OXKqXYd65H68jADWrFlDdXU1//zP/9y37aijjuLaa6+lp6eHBQsWsHbtWg4cOMAXvvAFrrnmGtauXcuiRYuor6/n+eefZ8qUKfzkJz+JvPZyIcpCEICPhy3pLp8hsGpzGys3tdGTNaL1qLJyUxtNR43o67t5Z33Itv6zn1lpKANAiLhpb9junXH1Wcl5LBVS7WYdNbx4GeXYvn07J598su139957L8OHD2fDhg0cOHCAadOm8elPfxqAzZs3s337dsaOHcu0adNYv34906dPt/2dOCgL1ZAvYs5wGQdu9P+zJzdw0ZQGKrOzGbtVQ7mSxKykpig5j6WIVLtf+MIX+PCHP8zUqVP5xS9+wQMPPMCkSZM49dRT2bNnDy+99BIAp5xyCo2NjVRUVDBp0qS+9NJJIdCKQET+AVgEHAecki1IY7ff2cC3gUpgqarmKplNAJYBI4FNwP9Q1XCKkW5d4d2P32TRF59EGUDkRv/vZtUQB0kItCrH4jCJ9VgKSbV7wgknsHLlyr7P3/ve99i9ezdNTU2MHz+e73znOwNqBaxdu9YxvXRSCLoieB64EPi10w4iUgl8DzgHOB64TESOz359O3Cnqn4Q2AtcGbA99pSoK2jUM0w3XilJnBkmZSZecoOlB0rOY8kmNsaEavf0009n//79/OAHP+jblisbedZZZ/GDH/yA7u6Mu/Uf//hH3n333UDni4pAgkBVf6+qxUaAU4CXVfWV7Gx/GTArW7D+dCBX3fl+MgXszVOirqBRD7pReg2ZJCnCqeQGSw+UnMdSSKpdEWHVqlX86le/YsKECZxyyinMnTuX22+/nauuuorjjz+ek08+mRNPPJFrrrkmcTN/J6IwFjcA1jVaK3AqGXVQu2pf2atWbOoa5xCRq4GrAcaPH++tBSXqChr1oOvGKyWJpTeTIpxMGtKTRtweS74ISbU7ZswYli1bZvvdrbfeOqA2wWmnncZpp53W9/m73/2u8TYFpaggEJFfAqNtvvqqqj5qvkn2qOoSYAlk6hF4OrhEXUHjGHSLeaUkcbBLinAqycHSA6Xm5ZXinqKCQFXPCHiONmCc5XNjdtseoE5EhmRXBbnt5ilRV9AkDrpJHOyS1E/pYJlSikShGtoAHJP1EGoDLgUuV1UVkaeBi8nYDeYC4awwkp7904EkDrq5dsXdBitJ7acU86hqogKxkorXypOBSlWKyAXAd4BRQDuwRVXPEpGxZNxEz83udy5wFxn30ftU9ZvZ7e8nIwRGAJuBz6jqgWLnDbtUZUpKSvJ49dVXOeywwxg5cmQqDAqgquzZs4e//e1vTJgwod93TqUq05rFKSkpJUF3dzetra3s31+eqaK9MGzYMBobG6mq6l9LxEkQlE2KiZSUlNKmqqpqwAw3xQxpiomUlJSUMicVBCkpKSllTioIUlJSUsqckjQWi8gu4M8+D68HdhtsTqlQjtddjtcM5Xnd6TW74yhVHZW/sSQFQRBEZKOd1XywU47XXY7XDOV53ek1ByNVDaWkpKSUOakgSElJSSlzylEQLIm7ATFRjtddjtcM5Xnd6TUHoOxsBCkpKSkp/SnHFUFKSkpKioVUEKSkpKSUOYNWEIjI2SLyooi8LCILbL4fKiLLs9//VkSOjqGZRnFxzf8qIi+IyFYReUpEjoqjnaYpdt2W/S4SERWRknczdHPNIjIne7+3i8h/Rd3GMHDxjI8XkadFZHP2OT83jnaaRETuE5E3ReR5h+9FRO7O9slWETnZ80lUddD9I5Pu+k/A+4Fq4Dng+Lx9Pg/8MPv3pcDyuNsdwTV/CqjN/v0vpX7Nbq87u99hwK+BZ4CmuNsdwb0+hkxq98Ozn4+Iu90RXfcS4F+yfx8PvBZ3uw1c9yeAk4HnHb4/F/g5IMBHgN96PcdgXRGcArysqq+oaheZmgez8vaZBdyf/fthYIaUdpLzotesqk+rakf24zNkqsKVOm7uNcDXgduBwZDD2M01fw74nqruBVDVNyNuYxi4uW4F3pf9eziwI8L2hYKq/hp4q8Aus4AHNMMzZCo/jvFyjsEqCBoAa5Hi1uw22300UypzHzAyktaFg5trtnIlmVlEqVP0urNL5XGq+niUDQsRN/f674C/E5H1IvKMiJwdWevCw811LwI+IyKtwBPAtdE0LVa8vvsDSOsRlCEi8hmgCfhk3G0JGxGpAL4FXBFzU6JmCBn10GlkVn6/FpGTVLU9zkZFwGXAj1X1f4nIR4H/LSInqmpv3A1LMoN1RdAGjLN8bsxus91HRIaQWUbuiaR14eDmmhGRM4CvAjPVRVnQEqDYdR8GnAisFZHXyOhQW0rcYOzmXrcCLararaqvAn8kIxhKGTfXfSWwAkBVfwMMI5OcbTDj6t0vxGAVBBuAY0RkgohUkzEGt+Tt0wLMzf59MbBGs5aXEqXoNYvIZOAeMkJgMOiMoch1q+o+Va1X1aNV9WgytpGZqlrKtU7dPN+ryKwGEJF6MqqiVyJsYxi4ue6/ADMAROQ4MoJgV6StjJ4W4LNZ76GPAPtU9XUvPzAoVUOqelBEvgisJuNpcJ+qbheRW4CNqtoC3Etm2fgyGUPMpfG1ODgur/kO4FDgoaxd/C+qOjO2RhvA5XUPKlxe82rg0yLyAtADzFPVUl7xur3ufwN+JCI3kDEcX1HiEzxE5EEyQr0+a/u4GagCUNUfkrGFnAu8DHQA/+j5HCXeRykpKSkpARmsqqGUlJSUFJekgiAlJSWlzEkFQUpKSkqZkwqClJSUlDInFQQpKSkpZU4qCFJSUlLKnFQQpKSkpJQ5/z9lsQeDpGR8qwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "original_samples = to_numpy(X_u_train[:N, :])\n",
    "generated_samples = to_numpy(X_gen)\n",
    "x = original_samples[:, 0:1]\n",
    "t = original_samples[:, 1:2]\n",
    "x_g = generated_samples[:, 0:1]\n",
    "t_g = generated_samples[:, 1:2]\n",
    "plt.scatter(t, x, label='Original')\n",
    "plt.scatter(t_g, x_g, label='Gen')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}