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

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[MLENS] backend: threading\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[pyKeOps]: Warning, no cuda detected. Switching to cpu only.\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 = \"/Users/pongpisit/Desktop/research/pinn/Solving-Differential-Equations-with-Neural-Networks/SymbolicMathematics/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, prob_activation=torch.sigmoid, 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.prob_activation = prob_activation\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 = self.prob_activation(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], prob_activation=F.hardsigmoid, 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.586485374602489e-06\n",
      "Test MSE: 8.9e-04\n",
      "Epoch 100:  2.962250107430009e-07\n",
      "Test MSE: 1.0e-03\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": "ab1ddb4e1c8b4b74967436a7c0b7b58c",
       "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.59E-07\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuKklEQVR4nO3deXxU5dn/8c812TcSIGFNQthlDxATFgU3UEGhValFrcUNta3W2sdan1qXPk+f2p/dXFtttS51LWIVwUrrAojsIKvsBEiCEALZyJ5cvz9moCGEkISZnMzM9X695pWZM+ecue5Mkm/OOffct6gqxhhjgpfL6QKMMcY4y4LAGGOCnAWBMcYEOQsCY4wJchYExhgT5CwIjDEmyIU6XUBLJSYmalpamtNlGGOMX1mzZs1hVU1q7Dm/C4K0tDRWr17tdBnGGONXRGTv6Z6zU0PGGBPkLAiMMSbIWRAYY0yQ87trBMaYlquuriYnJ4eKigqnSzE+FhkZSXJyMmFhYc3exoLAmCCQk5NDXFwcaWlpiIjT5RgfUVUKCgrIycmhd+/ezd7OTg0ZEwQqKiro3LmzhUCAExE6d+7c4iO/oAmCiupa5q7NwYbdNsHKQiA4tOZ9DpogeO/LXO59ez0Ltxx0uhRj2j9VWL4c3n3X/dVH/0D94Q9/oKyszCf7bq7CwkKeffbZNnu9tLQ0Dh8+DMC4ceNavZ+XXnqJvLw8r9QUNEFw9ahk+nWJ5VcLvqKqps7pcoxpvxYsgNRUmDQJZs1yf01NdS/3skAJgpqamlZt98UXX7T6Nf0iCEQkUkRWish6EdksIo82ss4sEckXkS89t1t9VU9oiIufTR1EdkEZryzL9tXLGOPfFiyAa66BnBwoLYXiYvfXnBz38laGwbFjx5g6dSojRoxg6NChvPXWWzz55JPk5eVx4YUXcuGFFwKwcOFCxo4dy6hRo5gxYwalpaUArFmzhokTJzJ69GguvfRSDhw4AMAFF1zAD3/4Q9LT0xk6dCgrV6488Xo333wzmZmZjBw5kvfeew+AzZs3k5mZSXp6OsOHD2fHjh389Kc/ZdeuXaSnp3PfffedUvv//M//MHDgQM477zxmzpzJb37zmxOvfc8995CRkcETTzzBvHnzyMrKYuTIkVxyySUcPOg++1BQUMDkyZMZMmQIt95660mnp2NjY0/cf/zxxzn33HMZPnw4Dz/8MADZ2dkMGjSI2267jSFDhjB58mTKy8uZM2cOq1ev5vrrryc9PZ3y8vJWvS8nqKpPboAAsZ77YcAKYEyDdWYBT7dkv6NHj9bWqqur0xv+slyHPfxPPVJa2er9GONvtmzZcuaV6upUe/ZUdZ8IavyWnOxer4XmzJmjt95664nHhYWFqqraq1cvzc/PV1XV/Px8Pf/887W0tFRVVR977DF99NFHtaqqSseOHauHDh1SVdU333xTb7rpJlVVnThx4on9Llq0SIcMGaKqqg888IC++uqrqqp69OhR7d+/v5aWluoPfvAD/dvf/qaqqpWVlVpWVqZ79uw5sV1DK1eu1BEjRmh5ebkWFxdrv3799PHHHz/x2nfeeeeJdY8cOaJ1nu/Nn//8Z7333ntVVfWuu+7SRx99VFVVP/jgAwVOtDkmJkZVVT/66CO97bbbtK6uTmtra3Xq1Km6aNEi3bNnj4aEhOi6detUVXXGjBkn2jVx4kRdtWpVo3U39n4Dq/U0f1d91n3U88KlnodhnpujV2pFhAenDubyJxbzxMc7eGTaECfLMaZ9WbECioqaXqewEFauhKysFu162LBh/PjHP+b+++/niiuu4Pzzzz9lneXLl7NlyxbGjx8PQFVVFWPHjmXbtm1s2rSJSZMmAVBbW0v37t1PbDdz5kwAJkyYQHFxMYWFhSxcuJD333//xH/vFRUV7Nu3j7Fjx/LLX/6SnJwcrrrqKvr3799k3UuXLmX69OlERkYSGRnJlVdeedLz11577Yn7OTk5XHvttRw4cICqqqoT3TcXL17M3LlzAZg6dSodO3Y85XUWLlzIwoULGTlyJAClpaXs2LGD1NRUevfuTXp6OgCjR48mOzu7yZpbw6efIxCREGAN0A94RlVXNLLa1SIyAdgO/EhV9/uypoHd4piZmcqry/dyw5he9OsSe+aNjAkGBw6A6wxni10uaMV56QEDBrB27VoWLFjAgw8+yMUXX8xDDz100jqqyqRJk3jjjTdOWr5x40aGDBnCsmXLGt13w14yIoKq8s477zBw4MCTnhs0aBBZWVnMnz+fKVOm8Nxzz9GnT58Wt+e4mJiYE/fvuusu7r33XqZNm8Znn33GI4880uz9qCoPPPAAt99++0nLs7OziYiIOPE4JCTk7E8DNcKnF4tVtVZV04FkIFNEhjZYZR6QpqrDgX8BLze2HxGZLSKrRWR1fn7+Wdf1o0kDiA4L4VcLvjrrfRkTMLp3h7ozdKSoq4MePVq867y8PKKjo7nhhhu47777WLt2LQBxcXGUlJQAMGbMGJYuXcrOnTsB93n+7du3M3DgQPLz808EQXV1NZs3bz6x77feeguAzz//nPj4eOLj47n00kt56qmnTpyPX7duHQC7d++mT58+3H333UyfPp0NGzacVEND48ePZ968eVRUVFBaWsoHH3xw2jYWFRXRs2dPAF5++T9/yiZMmMDrr78OwIcffsjRo0dP2fbSSy/lxRdfPHFNJDc3l0OHDjX5PW2q7pZqk15DqloIfApc1mB5gapWeh7+BRh9mu2fV9UMVc1ISmp0OO0WSYyN4PsX9ePjrYf4fMfhs96fMQEhKwvi45teJyEBMjNbvOuNGzeeuEj76KOP8uCDDwIwe/ZsLrvsMi688EKSkpJ46aWXmDlzJsOHD2fs2LFs3bqV8PBw5syZw/3338+IESNIT08/qbdNZGQkI0eO5I477uCFF14A4Oc//znV1dUMHz6cIUOG8POf/xyAt99+m6FDh5Kens6mTZu48cYb6dy5M+PHj2fo0KGnXCw+99xzmTZtGsOHD+fyyy9n2LBhxJ/me/TII48wY8YMRo8eTWJi4onlDz/8MIsXL2bIkCHMnTuX1NTUU7adPHky1113HWPHjmXYsGFcc801Z/wjP2vWLO644452f7E4CUjw3I8ClgBXNFine7373wSWn2m/Z3OxuL7yqhod/9jHOuWJxVpb2/KLX8b4k2ZdLFZVnT9fNSqq8QvFUVHu59uRpi6YektJSYmqqh47dkxHjx6ta9as8enreUNLLxb78oigO/CpiGwAVgH/UtUPROQXIjLNs87dnq6l64G7cfciahORYSH86JIBbM4r5p+bv26rlzWmfZsyBebMgeRkiI2FDh3cX5OT3cunTHG6wjY3e/Zs0tPTGTVqFFdffTWjRo1yuiSvE/WzIRcyMjLUWzOU1dYpl/5hMarKR/dMIDQkaD5fZ4LMV199xaBBg5q/gaq7d1BenvuaQGYm2BAVfqOx91tE1qhqRmPrB/VfvhCXcO+kAezKP2ZHBcbUJ+K+ZvDNb7q/WggEtKAOAoBLh3SjT2IMzy3abQPSmYBmP9/BoTXvc9AHQYhLuG1CHzbmFrFsV4HT5RjjE5GRkRQUFFgYBDj1zEcQGRnZou1sYhrgmyN78tuF2/nT4t2M65d45g2M8TPJycnk5OTgjc/hmPbt+AxlLWFBgLsH0U3j03j8o21syStmcI8OTpdkjFeFhYW1aMYqE1yC/tTQcTdk9SImPITnF+9yuhRjjGlTFgQe8dFhzMxMZd6GAxwo8v5YHsYY015ZENTz3XFp1Kny5kqfjntnjDHtigVBPSmdojm/fxJvrdpPTa3NYmaMCQ4WBA1cl5nK18UVfLrNelcYY4KDBUEDFw/qQtcOETadpTEmaFgQNBAW4uLGsWks2XGYbV97Z6xvY4xpzywIGnFdZiqRYS5e/HyP06UYY4zPWRA0omNMONeMTubddbkUlFaeeQNjjPFjFgSnMWtcGlW1dby9OsfpUowxxqcsCE6jX5c4snp34vWVe6mrs4G6jDGBy4KgCTeM6cX+I+Us3mFdSY0xgcuCoAmXDulGYmw4f1u+z+lSjDHGZywImhAe6uJbGSl8svUguYU2/pAxJjBZEJzBzMxUFHhzpR0VGGMCkwXBGaR0iubCgV14c9V+qm38IWNMALIgaIYbxqSSX1LJv7YcdLoUY4zxOguCZpg4oAs9E6L42/K9TpdijDFeZ0HQDCEu4bqsVL7YVcCu/FKnyzHGGK+yIGimb2WkEOaCRS/Pg3ffheXLQe2DZsYY/2eT1zdT0pKPWfHczUSUlqCRYUhdHSQkwHPPwZQpTpdnjDGt5rMjAhGJFJGVIrJeRDaLyKONrBMhIm+JyE4RWSEiab6q56wsWADXXEOnI4eIqSpHiouhtBRycuCaa9zPG2OMn/LlqaFK4CJVHQGkA5eJyJgG69wCHFXVfsDvgV/7sJ7WUYXZs6H8NB8oKy+H22+300TGGL/lsyBQt+NXVsM8t4Z/LacDL3vuzwEuFhHxVU2tsmIFFBU1vU5hIaxc2SblGGOMt/n0YrGIhIjIl8Ah4F+quqLBKj2B/QCqWgMUAZ0b2c9sEVktIqvz89t4ALgDB8B1hm+TywV5eW1TjzHGeJlPg0BVa1U1HUgGMkVkaCv387yqZqhqRlJSkldrPKPu3aHuDJ8orquDHj3aph5jjPGyNuk+qqqFwKfAZQ2eygVSAEQkFIgHCtqipmbLyoL4+KbXSUiAzMw2KccYY7zNl72GkkQkwXM/CpgEbG2w2vvAdz33rwE+UW1nV11F4PnnISqq8eejotxdSNvZpQ1jjGkuXx4RdAc+FZENwCrc1wg+EJFfiMg0zzovAJ1FZCdwL/BTH9bTelOmwJw5kJwMsbHQoQMVkdF83SGRyjfess8RGGP8ms8+UKaqG4CRjSx/qN79CmCGr2rwqilTYN8+d++gvDx2agxXrKji8Z4j/KQBxhjTOPtkcUuIuK8ZAENU6btnEa+v3MeMjBSHCzPGmNazsYZaSUSYmZnKun2FfHWg2OlyjDGm1SwIzsI1o5MJD3Xx+gqbvcwY478sCM5CQnQ4U4d15x/rcimrqnG6HGOMaRULgrN0XVYqJZU1zFtvnyw2xvgnC4KzlNGrI/27xNrpIWOM37IgOEsi7tnL1ucUsSn3DIPTGWNMO2RB4AVXjUwmMszF6yvtqMAY438sCLwgPjqMK4b34L11uZRW2kVjY4x/sSDwkuuyUjlWVct7X+Y6XYoxxrSIBYGXjExJ4Jxucby+Yh/tbdw8Y4xpigWBl4gI12elsjmvmA05dtHYGOM/LAi8aPrInkSFhVhXUmOMX7Eg8KIOkWFMG9GD99fnUVxR7XQ5xhjTLBYEXnZdVirl1bW8t84uGhtj/IMFgZcNT45nSI8OvGYXjY0xfsKCwMuOf9J469clrNtf6HQ5xhhzRhYEPjA9vScx4XbR2BjjHywIfCA2IpTpI3syb30eRWV20dgY075ZEPjIdZmpVNbUMXddjtOlGGNMkywIfGRoz3hGJMfbJ42NMe2eBYEPXZeVyo5Dpazee9TpUowx5rQsCHzoyhE9iIsItYvGxph2zYLAh6LDQ/nGyJ7M33iAo8eqnC7HGGMaZUHgY9dlpVJVU8c7a+2isTGmffJZEIhIioh8KiJbRGSziPywkXUuEJEiEfnSc3vIV/U4ZVD3DoxMTeD1lXbR2BjTPvnyiKAG+LGqDgbGAN8XkcGNrLdEVdM9t1/4sB7HXJ/Vi935x1i++4jTpRhjzCl8FgSqekBV13rulwBfAT199Xrt2RXDuxMfFcary7OdLsUYY07RJtcIRCQNGAmsaOTpsSKyXkQ+FJEhbVFPW4sMC2FmZir/3PQ1OUfLnC7HGGNO4vMgEJFY4B3gHlUtbvD0WqCXqo4AngL+cZp9zBaR1SKyOj8/36f1+sqNY3shIrz8RbbTpRhjzEl8GgQiEoY7BF5T1bkNn1fVYlUt9dxfAISJSGIj6z2vqhmqmpGUlOTLkn2mR0IUlw/txpur9lNaWeN0OcYYc4Ivew0J8ALwlar+7jTrdPOsh4hkeuop8FVNTrv5vN6UVNTwzhrrSmqMaT9Cfbjv8cB3gI0i8qVn2X8DqQCq+ifgGuBOEakByoFvawD3sRyV2pGRqQn8dekevjOmFy6XOF2SMcb4LghU9XOgyb90qvo08LSvamiPbh7fm7veWMcnWw9xyeCuTpdjjDH2yeK2dtnQbnSPj+TFpXucLsUYYwALgjYXFuLiu+PS+GJXAVvyGnaiMsaYtmdB4IBvn5tCVFgIf7WjAmNMO2BB4ICE6HCuHt2T977M41BJhdPlGGOCnAWBQ245rw/VdXW88sVep0sxxgQ5CwKH9E6MYfLgrry6fC/H7ANmxhgHWRA4aPaEvhSVV/P31fudLsUYE8QsCBw0upf7A2YvL9tLXV3Afo7OGNPOWRA4bNa4NPYcPsai7f45mJ4xxv9ZEDjs8qHd6RIXwV9tVFJjjEMsCBwWHuriO2N6sXh7PjsPlTpdjjEmCDUrCEQkRkRcnvsDRGSaZ4hp4wUzs1IJD3HxyrJsp0sxxgSh5h4RLAYiRaQnsBD3qKIv+aqoYJMYG8GVI3owZ00OxRXVTpdjjAkyzQ0CUdUy4CrgWVWdAQTktJJOuWl8GmVVtby5cp/TpRhjgkyzg0BExgLXA/M9y0J8U1JwGtoznvH9OvP84j1UVNc6XY4xJog0NwjuAR4A3lXVzSLSB/jUZ1UFqbsv6s/h0kresKMCY0wbalYQqOoiVZ2mqr/2XDQ+rKp3+7i2oJPVpzOZvTvx3KLdVNbYUYExpm00t9fQ6yLSQURigE3AFhG5z7elBacfXtyfr4sr+Ptqm9fYGNM2mntqaLCqFgPfAD4EeuPuOWS8bFzfzoxKTeCPn+2iqqbO6XKMMUGguUEQ5vncwDeA91W1GrDBcXxARLj74v7kFpYzd60dFRhjfK+5QfAckA3EAItFpBdg8yz6yMQBSYxISeCpT3baUYExxueae7H4SVXtqapT1G0vcKGPawtaIsI9l7iPCuassaMCY4xvNfdicbyI/E5EVntuv8V9dGB85IIBSaSnJPDMp3ZUYIzxreaeGnoRKAG+5bkVA3/1VVHGfVTwo0kDyC0s522buMYY40PNDYK+qvqwqu723B4F+viyMAMT+icyKtV9VGCfKzDG+Epzg6BcRM47/kBExgPlvinJHHf8qOBAUQVvr7KjAmOMbzQ3CO4AnhGRbBHJBp4Gbm9qAxFJEZFPRWSLiGwWkR82so6IyJMislNENojIqBa3IMCd1y+RjF4deebTXTYGkTHGJ5rba2i9qo4AhgPDVXUkcNEZNqsBfqyqg4ExwPdFZHCDdS4H+ntus4E/tqT4YHD8qODr4gresqMCY4wPtGiGMlUt9nzCGODeM6x7QFXXeu6XAF8BPRusNh14xdMldTmQICLdW1JTMBjX1z0G0TOf7rSjAmOM153NVJXS7BVF0oCRwIoGT/UE6v+bm8OpYYGIzD7edTU/P/gmeRcR/mvyQA6VVNosZsYYrzubIGjWEBMiEgu8A9xT72iiZS+k+ryqZqhqRlJSUmt24fcye3di4oAknv1sFyU2i5kxxouaDAIRKRGR4kZuJUCPM+3cMz7RO8Brqjq3kVVygZR6j5M9y0wj/mvyQArLqvnzkj1Ol2KMCSBNBoGqxqlqh0Zucaoa2tS2IiLAC8BXqvq706z2PnCjp/fQGKBIVQ+0qiVBYFhyPFOGdeOFJbspKK10uhxjTIA4m1NDZzIe91DVF4nIl57bFBG5Q0Tu8KyzANgN7AT+DHzPh/UEhHsnDaC8upZnP9vldCnGmADR5H/1Z0NVP+cMF5RVVYHv+6qGQNSvSxxXj0rm1WV7mTUujZRO0U6XZIzxc748IjA+8uPJAwlxCb/68CunSzHGBAALAj/ULT6SOy/oy4KNX7Nid4HT5Rhj/JwFgZ+67fw+9IiP5BcfbKG2ziaLM8a0ngWBn4oKD+H+y89hc14x79iUlsaYs2BB4MemjejByNQEHv9oG6WVNU6XY4zxUxYEfkxEeOiKweSXVPIn605qjGklCwI/NzK1I99I78HzS3aTc7TM6XKMMX7IgiAA/OSyc3AJPPbhVqdLMcb4IQuCANAjIYrbJ/Tlgw0HWJ19xOlyjDF+xoIgQNw+sQ/dOri7k9ZZd1JjTAtYEASI6PBQ7r98IBtyinh3nQ3gaoxpPguCADJ9RE9GJMfz/z7aSlmVdSc1xjSPBUEAcbmEh64czMHiSp749w6nyzHG+AkLggAzulcnrs1I4S+f72FTbpHT5Rhj/IAFQQD67ymD6Bgdzk/mbKCyxia7N8Y0zYIgAMVHh/Grq4ax5UAxv/5wm9PlGGPaOQuCADVpcFe+O7YXLy7dw2fbDjldjjGmHbMgCGAPTBlE/y6x/PSdjRRXVDtdjjGmnbIgCGCRYSH8ZsYIDpVU8MsPbDYzY0zjLAgC3IiUBG6f2Je3Vu+3U0TGmEZZEASBH17cn/5dYnlgrp0iMsacyoIgCESGhfD4jBEcLK7g/+bbKSJjzMksCIJEekoCsyf05c1V+1m0Pd/pcowx7YgFQRC555L+9OsSy0/mrKegtNLpcowx7YQFQRCJDAvhiW+nc7Ssmh//fb0NV22MAXwYBCLyoogcEpFNp3n+AhEpEpEvPbeHfFWL+Y8hPeL5+dRBfLYtn798vtvpcowx7YAvjwheAi47wzpLVDXdc/uFD2sx9dwwpheXD+3G//vnNtbuO+p0OcYYh/ksCFR1MWDzJrZDIsJjVw+nW3wkd72+jqIy61JqTDBz+hrBWBFZLyIfisgQh2sJKvFRYTx93SgOFlfwk3fWo2rXC4wJVk4GwVqgl6qOAJ4C/nG6FUVktoisFpHV+fnW9dFb0lMS+Onl5/DR5oO8smyv0+UYYxziWBCoarGqlnruLwDCRCTxNOs+r6oZqpqRlJTUpnUGulvO683F53Thl/O/solsjAlSjgWBiHQTEfHcz/TUUuBUPcFKRPjNjBF0jg3nztfWUFhW5XRJxpg25svuo28Ay4CBIpIjIreIyB0icodnlWuATSKyHngS+LbaiWpHdIwJ59nrR3GwqJIfvL6Omto6p0syxrQh8be/vRkZGbp69WqnywhIb6/ez0/mbOCW83rz8ysGO12OMcaLRGSNqmY09lxoWxdj2q9vZaSwJa+YFz7fw+DuHbh6dLLTJRlj2oDT3UdNO/OzqYMY26czD8zdyOps+xiIMcHAgsCcJCzExbPXj6JHQiS3v7qG/UfKnC7JGONjFgTmFB1jwnlh1rlU19Zx80urbDIbYwKcBYFpVN+kWP54w2j2HD7G9/62lsqaWqdLMsb4iAWBOa3x/RL51VXD+HznYb7/2jqqrVupMQHJgsA0aUZGCr+YPoR/f3WQe99eT63NYWBMwLHuo+aMbhybRllVLY99uJWoMBePXTUcl0ucLssY4yUWBKZZ7pjYl7KqWp78eAdRYSE8Mm0InhFCjDF+zoLANNuPLulPeVUNf16yh6jwUO6/bKCFgTEBwILANJuI8N9TBlFeXcufFu0iPNTFjy7pb2FgjJ+zIDAtIiL8YtpQKqvrePLjHZRV1vCzqYMsDIzxYxYEpsVcLuHXVw8nJiKUv3y+h5KKGv7vqmGE2AVkY/ySBYFpFZdLePjKwXSICuPJj3dQWlnD769NJzzUeiQb428sCEyriQj3ThpAh8hQ/nf+V5RU1vDs9aOIjbAfK2P8if37Zs7aref34ddXD2PpzsNc9exSsg8fc7okY0wLWBAYr7j23FReuTmTQyWVTHv6cxZtz3e6JGNMM1kQGK8Z3y+ReT84jx4JUdz015U8t2gX/jYDnjHByILAeFVKp2jmfm8clw/tzq8+3Mr3XlvL0WNVTpdljGmCBYHxuujwUJ6+biQPXH4O//7qIJf+YTGL7VSRMe2WBYHxCRHh9ol9efd744mPCuPGF1fy8HubKK+yeQ2MaW8sCIxPDe0Zz7y7zuPm8b15edlernhqCRtzipwuyxi/syr7iM/+kbIgMD4XGRbCQ1cO5m+3ZHGsspZvPruU3y7cZkcHxjTTptwivvPCCv53/haf7N+CwLSZ8/on8tE9E7hyRA+e+mQnF//2M+atz7OeRcY0YeehUm5/dQ0do8O555IBPnkN8bdfwoyMDF29erXTZZiztHLPER55fzNbDhST2bsTj1w5hME9OjhdljFedfRYFSEhQlxE6ImBGUsra9iSV8z+I2XUqhITHkpMRAhRYSG4XMKRY1Vsyi1iQ04Ruw+Xsv9IOZFhLt6+fSzDkxNaXYuIrFHVjEafsyAwTqmtU95atZ/HP9pKUXk1MzNT+fHkgXSKCXe6NGPOyqLt+Tz58Q7W7D0KQGxEKMkdoyirqiXnaBlnmvE1xCUM7BpH78QYMtI6cvnQ7nSLjzyrmhwJAhF5EbgCOKSqQxt5XoAngClAGTBLVdeeab8WBIGnqKya3/97O68u30tkqItbzuvNLef3IT4qzOnSjDmJqjY55HpdnfKHf2/nyU92ktIpim+NTiEyLISco2XkFpYTFR5Kn8QY0lMT6NUpmrAQF2VVtZRW1lBZXUutKvFRYfRJivX6mF1OBcEEoBR45TRBMAW4C3cQZAFPqGrWmfZrQRC4dh4q4ff/2sH8jQfoEBnK9y7sx03j04gIDXG6NBPEjlXWEB0ewnOLd/P7f22na4dIMnp1ZOLAJOKjwggPcREe6sLlEn67cBtLdxYwY3Qy//vNoe3qZ9exU0MikgZ8cJogeA74TFXf8DzeBlygqgea2qcFQeDbnFfEbz7axqfb8knpFMXN43szIyPFRjU1bW7Jjny+88JKEmPDOVxaxYQBScRGhPDFrgIKy6pPWT8qLISHrxzMteemtLvJmpoKAid/s3oC++s9zvEsOyUIRGQ2MBsgNTW1TYozzhnSI56/3pTJ4u35/O5f23l03hZ+t3A7MzJSuGFMKn2SYp0u0QSJ5xfvJjE2gqzenRjcowN3TuyLyyXU1Nax/WApFTW1VNfUUVlTx7HKGoYlx5PcMdrpslvML/7FUtXngefBfUTgcDmmjUwYkMSEAUms23eUvy7N5pVl2by4dA/j+3XmhqxeXDK4K2Eh1gPa+MaOgyUs2XGY/5o8gB9c1P+k50JDXAHVy83JIMgFUuo9TvYsM+YkI1M7MjK1Iw9eMYi3V+3n9RX7uPO1tXTtEMG156YyMzOF7vFRTpdpAoiq8uxnuwgPdTEzM/DPQjj579T7wI3iNgYoOtP1ARPcusRF8oOL+rPk/ov4840ZnNOtA099soPxj33CDX9ZwZw1OZRUnHre1piWqK6t4+43v+TddbncND6NzrERTpfkc77sNfQGcAGQCBwEHgbCAFT1T57uo08Dl+HuPnqTqp7xKrBdLDb17Sso4+9r9vOPL3PZf6SciFAXkwZ3ZdLgrozrm0hSXOD/Ehvvmr/hAN9/fS33ThrAXRf1a3cXfVvLPlBmAp6qsnZfIf9Yl8v8jQc44pkD4ZxucVw8qAtTh/VgUPe4gPmlNr7zwNyNfLA+j3UPTSI0gK5BtddeQ8Z4jYgwuldHRvfqyCPThrApt4jPdx5myY58/vjZLp75dBe9E2OYOqw7lw3txuDuHXC5LBTMqZbtOkxWn04BFQJnYkFgAk6ISxiRksCIlAS+f2E/DpdW8tHmr1mw8QDPfraTpz/dScfoMMb27cy4vomM69uZ3okxdrRgyC0sJ7ugjBvHpjldSpuyIDABLzE2guuzenF9Vi8Ol1ayaFs+X+wq4Itdh1mw8WsAkuIiyEzrxLlpHTm3dyfO6daBEDtiCDpf7DwMwLh+nR2upG1ZEJigkhgbwdWjk7l6dDKqSnZBGV/sOszKPUdYtecI8ze6O651iAxlTJ/OjOvbmbF9E+nfJdZOJQWBL3YV0DkmnIFd45wupU1ZEJigJSL0Toyhd2IM12f1AiDnaBmrso+wfNcRlu46zMItBwH36JEjUuJJT0lgZIr7WkRHGyU1oJRUVLNw89dcPqx70J0mtCAwpp7kjtEkd4zmmyOTAXf31JXZR/hy/1G+3F/InxbtptYzhvDArnFk9u7E8OR4hvSIp3/XWPuksx+buzaXY1W13DCml9OltDkLAmOakNo5mtTO0Vwz2h0M5VW1bMwtYlX2EVbsOcLctTm8unwvAOEhLgZ2i2NgN/c48mmdY0hLjCatcwwxNmBeu6aqvLIsmxHJ7qO+YGM/nca0QFR4CJm9O5HZuxPfv9A9uU52wTE25RaxJa+YTXlFLN6ez5w1OSdt1yUu4j/BkBhD784xpCXG0Ccppl0NVRysXv4im135x/jNjBFOl+IICwJjzkKIS+ibFEvfpFimp/c8sfxYZQ3ZBcfIPlzm+XqM7IJjfLI1n8OlOSdt3zsxhoFd3UcS/brEktopmpRO0TYxj49V19Yxd20Om/OKeWXZXi4Z1JVvpPdwuixHWBAY4wMxEaEM6eG+dtBQaWUN2YePsfvwMXYcLGHr1yVszC060WPpuA6Roe5TU55gSOnovh8fFUaIS0jtHE2HSAuL1nrm05384d87CHUJlw/txh++nR5UHyKrz4LAmDYWGxHK0J7xDO15ckgcP4rYf6SMfUfK2H+knH1Hyth6oIR/bzlEVW3dKfvqEhdBWmIMCVFhxEaGEh8VRvf4SHokRNEjIYqeCVEkxUYERNfXg8UVRIS6CA1xnfgeHS6tpKyylrKqWsJDXfRIiKSypo7IsBAGdo2je0LkSRPHH7fjYAnPfLqTK0f04Ilr0wPi+3M2LAiMaSeaOoqoq1MOllSwt6CMsqoaKqvryC4oY1d+KXsLjrHvSBklFTUUlVdTWllz0rZhIUKXuEg6RIXRJS6C5I5RJHeMpktcBAnRYSREh5MQHUbH6PATRxvtyYacQp7+ZOeJrrwtFRHqIjE2gqS4CBJjIwgPFZZsP0xMRCgPXzk46EMALAiM8Qsul9A9PqpZ8y4UV1STV1hOXmE5uYUV5B4t52BxBSUV1RwsrmR9TmGj0ywe1yEylPjoMCJDQ4gMCyE2IpTuCZH09BxlHD/S6JEQSXR4439CjlXWsCGniP1HyhCB3YePkXu0nONDXDYc7LLh0JcdIsOoq1M25Bbx1YFi4iJDufuifnSICqO6VknpFEVqp2i6xEUSExFCdHgo5dW1HCgsJzIshNLKGrYfLOFQcSX5pZUcLnF/zTnqDsxLh3Zj1rg0EoNgiOnmsNFHjQlCJRXVFJRWUVhezdGyKgrLqigsq+ZoWTVFZVUUlVdTUV1HZU0tJRU1HCiq4EBROXUN/lx0jA6jR0IU3TpEEh0RSnlVDTlHy9l+sOSkdUNdQs+OUYTUP0XT4B/x4w9V3WGmCv26xDJlWHeuGtWTOLseclZs9FFjzEniIsNa/Ie1praOgyWV7iONo+Xkeo468grLOVBUQVlVDdHhoXSLj+TSId1IT02gX1IsqtA1PsK6ybZjFgTGmGYJDXHR03Na6Nw0p6sx3hScfaWMMcacYEFgjDFBzoLAGGOCnAWBMcYEOQsCY4wJchYExhgT5CwIjDEmyFkQGGNMkPO7ISZEpATY1srN44GiVq7T2PKGy+o/PtP9ROBws6pufo3NWccb7ai/zJft8FYb6j+296JlNTbn+eZ+7xt73LAd/vxe1L/fHtvRS1WTGt2bqvrVDVh9Fts+39p1GlvecFn9x2e678/taLDMZ+3wVhuaqN3ei7N8L1ryvW9OO/z5vfCXdjR2C7ZTQ/POYp3GljdcNq+F91vL6XZ4ow3N2Y+32lD/sb0Xza+luc8393vf2GNvtsPp96K5NZyJr9txCn88NbRaTzOCnj+xdrQfgdAGCIx2BEIbwP/a4Y9HBM87XYCXWDvaj0BoAwRGOwKhDeBn7fC7IwJjjDHe5Y9HBMYYY7zIgsAYY4KcBYExxgS5gAoCETlfRP4kIn8RkS+crqe1RMQlIr8UkadE5LtO19MaInKBiCzxvB8XOF3P2RCRGBFZLSJXOF1La4jIIM/7MEdE7nS6ntYSkW+IyJ9F5C0Rmex0Pa0lIn1E5AURmeN0Lce1myAQkRdF5JCIbGqw/DIR2SYiO0Xkp03tQ1WXqOodwAfAy76s93S80Q5gOpAMVAM5vqr1dLzUBgVKgUgcaAN4rR0A9wNv+6bKpnnp9+Irz+/Ft4Dxvqz3dLzUjn+o6m3AHcC1vqz3dLzUjt2qeotvK22h1n76zds3YAIwCthUb1kIsAvoA4QD64HBwDDcf+zr37rU2+5tIM5f2wH8FLjds+0cP22Dy7NdV+A1P34vJgHfBmYBV/hjGzzbTAM+BK7z1/ei3na/BUYFQDva/Hf7dLd2M3m9qi4WkbQGizOBnaq6G0BE3gSmq+qvgEYP00UkFShS1RJf1ns63miHiOQAVZ6HtT4st1Heei88jgIRPin0DLz0XlwAxOD+xS4XkQWqWufLuuvz1nuhqu8D74vIfOB1H5bcKC+9FwI8Bnyoqmt9XHKjvPy70W60myA4jZ7A/nqPc4CsM2xzC/BXn1XUOi1tx1zgKRE5H1jsy8JaoEVtEJGrgEuBBOBpn1bWMi1qh6r+DEBEZgGH2zIEmtDS9+IC4CrcgbzAl4W1UEt/L+4CLgHiRaSfqv7Jl8W1QEvfj87AL4GRIvKAJzAc1d6DoMVU9WGnazhbqlqGO9D8lqrOxR1oAUFVX3K6htZS1c+Azxwu46yp6pPAk07XcbZUtQD3dY52o91cLD6NXCCl3uNkzzJ/EwjtCIQ2QGC0IxDaANaOdqO9B8EqoL+I9BaRcNwX7d53uKbWCIR2BEIbIDDaEQhtAGtH++H01ep6V9DfAA7wny6Tt3iWTwG2474q/zOn6wyGdgRCGwKlHYHQBmtH+7/ZoHPGGBPk2vupIWOMMT5mQWCMMUHOgsAYY4KcBYExxgQ5CwJjjAlyFgTGGBPkLAhMwBCR0jZ+vTad80JEEkTke235miY4WBAYcxoi0uRYXKo6ro1fMwGwIDBeZ0FgApqI9BWRf4rIGs+Maed4ll8pIitEZJ2I/FtEunqWPyIir4rIUuBVz+MXReQzEdktInfX23ep5+sFnufniMhWEXnNM2QyIjLFs2yNiDwpIh80UuMsEXlfRD4BPhaRWBH5WETWishGEZnuWfUxoK+IfCkij3u2vU9EVonIBhF51JffSxPAnP5os93s5q0bUNrIso+B/p77WcAnnvsd4cQn628Ffuu5/wiwBoiq9/gL3EM4JwIFQFj91wMuAIpwDzbmApYB5+GenW0/0Nuz3hvAB43UOAv3cAWdPI9DgQ6e+4nATkCANE6eEGUy8LznORfuiU8mOP0+2M3/bgE3DLUxx4lILDAO+LvnH3T4zyQ5ycBbItId96xSe+pt+r6qltd7PF9VK4FKETmEe9a1htNvrlTVHM/rfon7j3YpsFtVj+/7DWD2acr9l6oeOV468H8iMgGowz3efddGtpnsua3zPI4F+tN+5rAwfsKCwAQyF1CoqumNPPcU8DtVfd8zccsj9Z471mDdynr3a2n896Y56zSl/mteDyQBo1W1WkSycR9dNCTAr1T1uRa+ljEnsWsEJmCpajGwR0RmgHuqQxEZ4Xk6nv+MGf9dH5WwDehTb2rD5k64Hg8c8oTAhUAvz/ISIK7eeh8BN3uOfBCRniLS5ezLNsHGjghMIIn2zPd83O9w/3f9RxF5EAgD3sQ9ufgjuE8ZHQU+AXp7uxhVLfd09/yniBzDPW59c7wGzBORjcBqYKtnfwUislRENuGet/c+ERkELPOc+ioFbgAOebstJrDZMNTG+JCIxKpqqacX0TPADlX9vdN1GVOfnRoyxrdu81w83oz7lI+dzzftjh0RGGNMkLMjAmOMCXIWBMYYE+QsCIwxJshZEBhjTJCzIDDGmCBnQWCMMUHu/wO54BXsQmDVlAAAAABJRU5ErkJggg==\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:23<00:00, 12.95it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -5.713312149047852\n",
      "Semi-supervised solver loss @Epoch 0:  4.001071453094482\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 1/300 [00:00<00:31,  9.57it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:23<00:00, 12.51it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.6514633297920227\n",
      "Semi-supervised solver loss @Epoch 100:  0.92142653465271\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:23, 12.64it/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.44it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.532879114151001\n",
      "Semi-supervised solver loss @Epoch 200:  0.8275467157363892\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:25, 11.64it/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.17it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.5176467895507812\n",
      "Semi-supervised solver loss @Epoch 300:  0.7745926380157471\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:20, 14.40it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:23<00:00, 12.97it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.5450326204299927\n",
      "Semi-supervised solver loss @Epoch 400:  0.7810180187225342\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 1/300 [00:00<00:34,  8.57it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:21<00:00, 13.67it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.5630310773849487\n",
      "Semi-supervised solver loss @Epoch 500:  0.7555843591690063\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:22, 13.50it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:23<00:00, 12.70it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.5669344067573547\n",
      "Semi-supervised solver loss @Epoch 600:  0.7515137195587158\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:21, 13.80it/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.33it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.5778838396072388\n",
      "Semi-supervised solver loss @Epoch 700:  0.7290953993797302\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:26, 11.27it/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.48it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.5912109613418579\n",
      "Semi-supervised solver loss @Epoch 800:  0.7037064433097839\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  1%|          | 2/300 [00:00<00:20, 14.87it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the generator for \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 300/300 [00:23<00:00, 12.66it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best generator loss: -0.5944154858589172\n",
      "Semi-supervised solver loss @Epoch 900:  0.7234848737716675\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:  6.062612101231935e-07\n",
      "Epoch 100:  2.1892678603308013e-07\n",
      "Epoch 200:  1.067091304207679e-07\n",
      "Epoch 300:  7.798048073937025e-08\n",
      "Epoch 400:  5.308918460400491e-08\n",
      "Epoch 500:  4.020244759317393e-08\n",
      "Epoch 600:  3.242866952746226e-08\n",
      "Epoch 700:  2.7727995899340385e-08\n",
      "Epoch 800:  2.3977570151600958e-08\n",
      "Epoch 900:  1.831545937136525e-08\n",
      "Epoch 1000:  1.724212594922392e-08\n",
      "Epoch 1100:  1.6634308153129496e-08\n",
      "Epoch 1200:  1.5379793438796696e-08\n",
      "Epoch 1300:  1.442038755072872e-08\n",
      "Epoch 1400:  1.3292114076079997e-08\n",
      "Epoch 1500:  1.2243337010886535e-08\n",
      "Epoch 1600:  1.1420932644057302e-08\n",
      "Epoch 1700:  1.0419674900674636e-08\n",
      "Epoch 1800:  9.117528954050158e-09\n",
      "Epoch 1900:  8.308806975776406e-09\n",
      "Epoch 2000:  7.700505122443246e-09\n",
      "Epoch 2100:  7.53306128586928e-09\n",
      "Epoch 2200:  6.878121183717667e-09\n",
      "Epoch 2300:  6.716982081655942e-09\n",
      "Epoch 2400:  6.451289280562378e-09\n",
      "Epoch 2500:  5.886655163322985e-09\n",
      "Epoch 2600:  5.821041870746058e-09\n",
      "Epoch 2700:  9.286266866581627e-09\n",
      "Epoch 2800:  5.062712471470832e-09\n",
      "Epoch 2900:  4.798273334216674e-09\n",
      "Epoch 3000:  4.103773321162407e-09\n",
      "Epoch 3100:  4.0638648002300215e-09\n",
      "Epoch 3200:  3.887029365046146e-09\n",
      "Epoch 3300:  3.82348819272238e-09\n",
      "Epoch 3400:  3.6556802029963364e-09\n",
      "Epoch 3500:  3.587372621183249e-09\n",
      "Epoch 3600:  3.53371576444772e-09\n",
      "Epoch 3700:  3.4927678527196804e-09\n",
      "Epoch 3800:  3.358192612878952e-09\n",
      "Epoch 3900:  3.1973135250495943e-09\n",
      "Epoch 4000:  2.9703191017915742e-09\n",
      "Epoch 4100:  2.9368556475617424e-09\n",
      "Epoch 4200:  6.657388894382166e-06\n",
      "Epoch 4300:  3.4518097891123034e-06\n",
      "Epoch 4400:  2.428248080832418e-06\n",
      "Epoch 4500:  1.953161699930206e-06\n",
      "Epoch 4600:  1.7415312640878255e-06\n",
      "Epoch 4700:  1.5388001202154555e-06\n",
      "Epoch 4800:  1.322260573033418e-06\n",
      "Epoch 4900:  1.2499633612605976e-06\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: 8.5e-03\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.41916308, 0.36172065, 0.5795206 , 0.45167705, 0.5093097 ],\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.3134551048278809\n",
      "[4 3 1 2 0]\n",
      "Epoch 5:  1.151564598083496\n",
      "[4 3 1 2 0]\n",
      "Epoch 10:  1.1834001541137695\n",
      "[4 3 1 2 0]\n",
      "Epoch 15:  1.192171573638916\n",
      "[4 3 1 2 0]\n",
      "Epoch 20:  1.1328505277633667\n",
      "[4 3 1 2 0]\n",
      "Epoch 25:  1.1666581630706787\n",
      "[4 3 1 2 0]\n",
      "Epoch 30:  1.165968656539917\n",
      "[4 3 1 2 0]\n",
      "Epoch 35:  1.1812328100204468\n",
      "[4 3 1 2 0]\n",
      "Epoch 40:  1.150465488433838\n",
      "[4 3 1 2 0]\n",
      "Epoch 45:  1.1654245853424072\n",
      "[4 3 1 2 0]\n",
      "Epoch 50:  1.148111343383789\n",
      "[4 3 1 2 0]\n",
      "Epoch 55:  1.1458015441894531\n",
      "[4 3 1 2 0]\n",
      "Epoch 60:  1.147080659866333\n",
      "[4 3 1 2 0]\n",
      "Epoch 65:  1.1605441570281982\n",
      "[4 3 1 2 0]\n",
      "Epoch 70:  1.1624443531036377\n",
      "[4 3 1 2 0]\n",
      "Epoch 75:  1.1684749126434326\n",
      "[4 3 1 2 0]\n",
      "Epoch 80:  1.182016134262085\n",
      "[4 3 1 2 0]\n",
      "Epoch 85:  1.1607873439788818\n",
      "[4 3 1 2 0]\n",
      "Epoch 90:  1.1495277881622314\n",
      "[4 3 1 2 0]\n",
      "Epoch 95:  1.1685261726379395\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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Feature importance ranking ---\n",
      "uf 0.32843083\n",
      "u_xx 0.31017172\n",
      "u_x 0.16392581\n",
      "u_xxx 0.12693694\n",
      "x 0.07053461\n"
     ]
    }
   ],
   "source": [
    "feature_importance = (semisup_model.selector.latest_weighted_features).detach().numpy()\n",
    "feature_importance = feature_importance/feature_importance.sum()\n",
    "print(\"--- Feature importance ranking ---\")\n",
    "for idx in np.argsort(feature_importance)[::-1]:\n",
    "    print(feature_names[idx], feature_importance[idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}