Vanilla Test Set - Sequential (best).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "65d9e98e-ca31-47c5-811d-13ba75bcb40c",
   "metadata": {},
   "source": [
    "Initiliaze and prepare all models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5cc6ee28-d34f-401b-abaf-8eae39d417da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load model from checkpoint\n",
    "import torch\n",
    "import yaml\n",
    "from models.decoder import architectures\n",
    "\n",
    "with open(\"models/decoder/sequential_decoder_6_params.yaml\") as file:\n",
    "    params = yaml.safe_load(file)\n",
    "\n",
    "#ckpt_file_path = \"models/decoder/sequential_decoder_6.ckpt\"\n",
    "#model = architectures.LightningDecoderNet.load_from_checkpoint(ckpt_file_path, \n",
    "#                                                               inner_dims=params[\"inner_dims\"], \n",
    "#                                                               seq_len=params[\"seq_len\"], \n",
    "#                                                               learning_rate=params[\"learning_rate\"])\n",
    "#torch.save(model.state_dict(), \"models/decoder/sequential_decoder_6.pth\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "831d662a-6e4a-491d-9b5e-cbaec4d6a2a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ~~~~~~~~~~~~~~~~\n",
    "# PyTorch Model\n",
    "\n",
    "import torch\n",
    "import yaml\n",
    "from models.decoder import architectures\n",
    "\n",
    "with open(\"models/decoder/sequential_decoder_length6_params.yaml\") as file:\n",
    "    params = yaml.safe_load(file)\n",
    "\n",
    "decoder = architectures.LightningDecoderNet(params[\"inner_dims\"], params[\"seq_len\"], params[\"learning_rate\"])\n",
    "decoder.load_state_dict(torch.load(\"models/decoder/sequential_decoder_length6.pth\"))\n",
    "decoder.eval();\n",
    "\n",
    "# ~~~~~~~~~~~~~~~~\n",
    "# ESM\n",
    "\n",
    "# To Do"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "863e1f89-0001-42b7-823e-353bb0526e53",
   "metadata": {},
   "source": [
    "### Data & Metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cff9d927-c602-44a6-8131-a620048f3279",
   "metadata": {},
   "source": [
    "Specify data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e09be215-572c-464e-8294-d8dbf96406cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.utils import data\n",
    "from data.gdm_dataset import GasDataSet\n",
    "\n",
    "dataset = GasDataSet(\"data/30x25/test.pt\", seq_len=params[\"seq_len\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7a4bad58-3fa1-46d3-ba8c-5a1651ed8f51",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "\n",
    "import pickle\n",
    "\n",
    "X_path = \"models/esm/X.pkl\"\n",
    "X_all = pickle.load(open(X_path, \"rb\"))\n",
    "\n",
    "y_path = \"models/esm/y.pkl\"\n",
    "y_all = pickle.load(open(y_path, \"rb\" ) )\n",
    "\n",
    "\n",
    "# Load ESM predictions\n",
    "\n",
    "from scipy.io import loadmat\n",
    "\n",
    "y_esm_path = \"models/esm/preds.mat\"\n",
    "y_esm_all = loadmat(y_esm_path, verify_compressed_data_integrity=False)\n",
    "y_esm_all = y_esm_all[\"export\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc88ca43-a1c3-40ba-aca1-214a076826fb",
   "metadata": {},
   "source": [
    "Specify metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "31bc2295-1248-4bae-9d9f-626583828ab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gdm_metrics\n",
    "\n",
    "rmse = {\"decoder\": 0,\n",
    "        \"esm\": 0,\n",
    "        }\n",
    "\n",
    "kld = {\"decoder\": 0,\n",
    "       \"esm\": 0,\n",
    "       }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbd06c75-e4ea-4946-8641-b35abc2abb66",
   "metadata": {},
   "source": [
    "## Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "63ec30ac-f22a-439e-b306-41decf063d06",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████| 31800/31800 [03:19<00:00, 159.26it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "########\n",
      "# RMSE #\n",
      "########\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'X_all' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[6], line 20\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m########\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m# RMSE #\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m########\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     19\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m elem \u001b[38;5;129;01min\u001b[39;00m rmse:\n\u001b[0;32m---> 20\u001b[0m     \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00melem\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m: \u001b[39m\u001b[38;5;130;01m\\t\u001b[39;00m\u001b[38;5;124m \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mrmse[elem]\u001b[38;5;241m/\u001b[39m\u001b[38;5;28mlen\u001b[39m(X_all)\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     22\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m########\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m# KL D #\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m########\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     23\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m elem \u001b[38;5;129;01min\u001b[39;00m kld:\n",
      "\u001b[0;31mNameError\u001b[0m: name 'X_all' is not defined"
     ]
    }
   ],
   "source": [
    "from tqdm import tqdm\n",
    "\n",
    "rmse = {\"decoder\": 0,\n",
    "        }\n",
    "\n",
    "kld = {\"decoder\": 0,\n",
    "       }\n",
    "\n",
    "\n",
    "dataloader = iter(data.DataLoader(dataset, shuffle=False, drop_last=True))\n",
    "\n",
    "for X, y in tqdm(dataloader):\n",
    "    with torch.no_grad(): \n",
    "        y_decoder = decoder(X.squeeze(2))\n",
    "        rmse[\"decoder\"] += gdm_metrics.rmse(y_decoder, y)\n",
    "        kld[\"decoder\"] += gdm_metrics.kld(y_decoder, y)\n",
    "            \n",
    "print(f\"########\\n# RMSE #\\n########\")\n",
    "for elem in rmse:\n",
    "    print(f\"{elem}: \\t {rmse[elem]/len(dataloader)}\")\n",
    "    \n",
    "print(f\"########\\n# KL D #\\n########\")\n",
    "for elem in kld:\n",
    "    print(f\"{elem}: \\t {kld[elem]/len(dataloader)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "846a55ba-1f73-4f8e-8500-082bdb29f775",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11.845172577123666"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "0.0664244294166565 * np.sqrt(len(dataloader)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ca2b46dd-c9c4-4714-a75b-51d078a80717",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "########\n",
      "# RMSE #\n",
      "########\n",
      "decoder: \t 0.0664244294166565\n",
      "########\n",
      "# KL D #\n",
      "########\n",
      "decoder: \t 0.0029262625612318516\n"
     ]
    }
   ],
   "source": [
    "print(f\"########\\n# RMSE #\\n########\")\n",
    "for elem in rmse:\n",
    "    print(f\"{elem}: \\t {rmse[elem]/len(dataloader)}\")\n",
    "    \n",
    "print(f\"########\\n# KL D #\\n########\")\n",
    "for elem in kld:\n",
    "    print(f\"{elem}: \\t {kld[elem]/len(dataloader)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3a8fda30-a5dc-42f2-b4b5-b2cab2c9991f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████| 30120/30120 [00:49<00:00, 605.02it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "########\n",
      "# RMSE #\n",
      "########\n",
      "decoder: \t 0.07328736037015915\n",
      "esm: \t 0.1408499926328659\n",
      "########\n",
      "# KL D #\n",
      "########\n",
      "decoder: \t 0.003549593733623624\n",
      "esm: \t 0.013593791984021664\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "rmse = {\"decoder\": 0,\n",
    "        \"esm\": 0,\n",
    "        }\n",
    "\n",
    "kld = {\"decoder\": 0,\n",
    "       \"esm\": 0,\n",
    "       }\n",
    "\n",
    "for i in tqdm(range(len(X_all))):\n",
    "    X = X_all[i].unsqueeze(0)\n",
    "    y = y_all[i].unsqueeze(0).unsqueeze(0)\n",
    "    \n",
    "    with torch.no_grad(): \n",
    "        y_decoder = decoder(X)\n",
    "        rmse[\"decoder\"] += gdm_metrics.rmse(y_decoder, y)\n",
    "        kld[\"decoder\"] += gdm_metrics.kld(y_decoder, y) \n",
    "    \n",
    "    y_esm = torch.tensor(np.reshape(y_esm_all[i],[30,25], order=\"F\"))\n",
    "    y_esm = y_esm.unsqueeze(0).unsqueeze(0)\n",
    "    rmse[\"esm\"] += gdm_metrics.rmse(y_esm, y)\n",
    "    kld[\"esm\"] += gdm_metrics.kld(y_esm, y)\n",
    "    \n",
    "print(f\"########\\n# RMSE #\\n########\")\n",
    "for elem in rmse:\n",
    "    print(f\"{elem}: \\t {rmse[elem]/len(X_all)}\")\n",
    "    \n",
    "print(f\"########\\n# KL D #\\n########\")\n",
    "for elem in kld:\n",
    "    print(f\"{elem}: \\t {kld[elem]/len(X_all)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "43f461fd-16bf-422b-8b9c-c1488e6df0fb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Decoder RMSE: 12.719105320877341\n",
      "ESM RMSE: 24.44467752275151\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Correct RMSE calculation: \n",
    "print(f\"Decoder RMSE: {0.07328736037015915 * 30120/np.sqrt(30120)}\")\n",
    "print(f\"ESM RMSE: {0.1408499926328659 * 30120/np.sqrt(30120)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "caa85295-2fa1-4581-b111-4336e8a98588",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 1000\n",
    "\n",
    "X = X_all[i].unsqueeze(0)\n",
    "y = y_all[i].unsqueeze(0).unsqueeze(0)\n",
    "y_esm = torch.tensor(np.reshape(y_esm_all[i],[30,25], order=\"F\"))\n",
    "y_esm = y_esm.unsqueeze(0).unsqueeze(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ddd2cb8b-ddb6-406b-8e92-39845c38e16b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fe8d6f78220>"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAGdCAYAAADHWhKIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhHElEQVR4nO3dbWxU173v8d8eezwew+CEEnvGwfFxK6L2Bi5XDSkEJcFEjRWfW9SEtCKJVIHURkl5kJATRaWoitsXOEoVhI5oqBodUVBDw5uERgKFuCI2jSgVQUTh0AiRGxPcA64Vl/jZ44dZ90XKnDo8zX+YsRfj70caCc/8995reY1/sxnP/jtwzjkBAKZUaKoHAAAgjAHAC4QxAHiAMAYADxDGAOABwhgAPEAYA4AHCGMA8EDxVA/gy1KplM6fP69YLKYgCKZ6OACQNeec+vr6VFVVpVDo2ue+3oXx+fPnVV1dPdXDAICc6ejo0Ny5c69Z410Yx2IxSdKyW59UcVCS0TZueNh8nCASsdXPKDMfw8oNDNrqk0lTvXXOUv7nbZ2zlP95s9b5MSlrXVpqqy+PmeolKRXN/Bhj40kdPv0f6Vy7lryF8SuvvKJf/vKXunDhgu666y5t27ZN999//3W3u/TWRHFQouJQhmEcpMzjCzIM+nR9yP7ktnLBmLHe1lbEOmcp//O2zvmLbfI7b9Y6PyZlrTPMjP+pt885VZTFC10Gb7nm5Rd4e/fu1caNG7V582adOHFC999/vxoaGnTu3Ll8HA4Abnp5CeOtW7fqhz/8oX70ox/pG9/4hrZt26bq6mrt2LEjH4cDgJtezsN4ZGREx48fV319/YT76+vrdeTIkcvqk8mkent7J9wAYLrJeRh/9tlnGh8fV2Vl5YT7Kysr1dnZeVl9c3OzysvL0zc+SQFgOsrbRR9ffsPaOXfFN7E3bdqknp6e9K2joyNfQwIAb+X80xRz5sxRUVHRZWfBXV1dl50tS1IkElEki4/hAEAhyfmZcUlJie6++261tLRMuL+lpUVLly7N9eEAoCDk5XPGjY2N+sEPfqBFixbp3nvv1W9+8xudO3dOzzzzTD4OBwA3vbyE8apVq9Td3a1f/OIXunDhgubPn68DBw6opqYm43244eGML+ZIDdqv7AmljBeKlITNx7ByQ0Om+pTxykPznCUFEfvFAxbWOUv5n3e+5yz5udb5fo5Pylob9x+U2a86DAx/w9lSm7cr8NauXau1a9fma/cAUFBooQkAHiCMAcADhDEAeIAwBgAPEMYA4AHCGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHggb70pblQQiWT8F26zaoBTavzz7eH8NwqScUzmpijG/UuSivP7FMlmTHmfd57nLNnHNBlrnffn+GSstbU3enGR8QiSy+AvPWdTy5kxAHiAMAYADxDGAOABwhgAPEAYA4AHCGMA8ABhDAAeIIwBwAOEMQB4gDAGAA8QxgDgAX97U5RFFYQyvM68xH5NvfU6fFeaWZ+MGxFopm0D43X42fQeyPe8zXOW8j5v1jo/JmWtS2xzcOEsItByCuvys1sAQJ4QxgDgAcIYADxAGAOABwhjAPAAYQwAHiCMAcADhDEAeIAwBgAPEMYA4AHCGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHiAMAYADxDGAOABwhgAPEAYA4AHCGMA8ABhDAAeIIwBwAOEMQB4oHiqB3A1bnBILhjPrHZoyH6A0oipPNBM+zGMXH+/rX44aTuAcc5S/udtnbOU/3mz1vkxKWsdjZrKg+Ii2/4luZJw5sWpzEs5MwYAD+Q8jJuamhQEwYRbPB7P9WEAoKDk5W2Ku+66S3/84x/TXxcV2f8rAADTSV7CuLi4mLNhADDIy3vGZ86cUVVVlWpra/X444/rk08+uWptMplUb2/vhBsATDc5D+PFixdr9+7dOnjwoF599VV1dnZq6dKl6u7uvmJ9c3OzysvL07fq6upcDwkAvJfzMG5oaNBjjz2mBQsW6Nvf/rb2798vSdq1a9cV6zdt2qSenp70raOjI9dDAgDv5f1zxjNmzNCCBQt05syZKz4eiUQUidg/EwkAhSTvnzNOJpP66KOPlEgk8n0oALhp5TyMn3vuObW1tam9vV1/+ctf9L3vfU+9vb1avXp1rg8FAAUj529T/O1vf9MTTzyhzz77TLfddpuWLFmio0ePqqamJteHAoCCkfMwfv3113OyH5dMygUuo9rU8LB5/+b/EkzC+9rW6/Ct887mv0FBaWkWW2XO3HtA+Z93vucs+bnW+X6OT8pah2wzD8Yy638zYRuXWS5Za+lNAQAeIIwBwAOEMQB4gDAGAA8QxgDgAcIYADxAGAOABwhjAPAAYQwAHiCMAcADhDEAeCDv/YyzFcycqSBUklFtUXEW0zBehx+UGq/bT2V+TXpaNGoqt/dcyKL3QDbfW4NsxpT3eed5zpJ9TJOx1kE4bN7GZDLW2tpfo9j+x5JdEOSlljNjAPAAYQwAHiCMAcADhDEAeIAwBgAPEMYA4AHCGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHjA20ZBqdvnKFWUYdMPQzOONGdr5BOMpWz1w6Omekkyz6LE1tglm0YwrjSzZk3ZCjTTvpG1yZNx3vmes2Sfd1BaajtAFs2OfFxr87ytax3OIgItp7CGmOHMGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHiAMAYADxDGAOABwhgAPEAYA4AHCGMA8IC3vSl67pypopLMrksfLbP3pggP2npTlHWNmeojfx8w1UtSMJw0b5NvQcr2fXKhLPqEAODMGAB8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHiAMAYADxDGAOABwhgAPEAYA4AHCGMA8IC3vSkufiNQqDSzPgejs219IyQp3G2duq2+pNv+OucGBk31qb5+U30omlmvj38VzJxpqzfu3/Xb5iBJztrDozRiKg9km3M2rPO2zjkwzlnK/7wnY62DaNRWX1xkqpckVxLOvDiVeSlnxgDgAcIYADxgDuPDhw9rxYoVqqqqUhAE2rdv34THnXNqampSVVWVotGo6urqdOrUqVyNFwAKkjmMBwYGtHDhQm3fvv2Kj7/00kvaunWrtm/frmPHjikej+uhhx5SX1/fDQ8WAAqV+Rd4DQ0NamhouOJjzjlt27ZNmzdv1sqVKyVJu3btUmVlpfbs2aOnn376xkYLAAUqp+8Zt7e3q7OzU/X19en7IpGIli1bpiNHjlxxm2Qyqd7e3gk3AJhuchrGnZ2dkqTKysoJ91dWVqYf+7Lm5maVl5enb9XV1bkcEgDcFPLyaYogmPhpU+fcZfddsmnTJvX09KRvHR0d+RgSAHgtpxd9xONxSV+cIScSifT9XV1dl50tXxKJRBSJ2D+gDgCFJKdnxrW1tYrH42ppaUnfNzIyora2Ni1dujSXhwKAgmI+M+7v79fHH3+c/rq9vV0ffPCBZs+erTvuuEMbN27Uli1bNG/ePM2bN09btmxRWVmZnnzyyZwOHAAKiTmM33//fS1fvjz9dWNjoyRp9erV+u1vf6vnn39eQ0NDWrt2rS5evKjFixfrnXfeUSwWy92oAaDAmMO4rq5OzrmrPh4EgZqamtTU1HQj49Lo7DGFopk1ALo1bv843EXNso2nzPitusovLK9pZNRU7pK2JioumzHl+f18c9MfSanhYVO99b24oNTeUMnKOu98z1lSYax1yDbzYGzcVC9JwTXy70Zq6U0BAB4gjAHAA4QxAHiAMAYADxDGAOABwhgAPEAYA4AHCGMA8ABhDAAeIIwBwAOEMQB4IKf9jHMpfLFYoaHMhncxsPWZkKTwP2xTLx4yHiCLPhBBbKapvihkPEa4xFYvKSjNc6/p8RnmTULhsKk+iBjnbdx/NoJo1FRv769hX7cg3/POYkzmeVv7axQXGY9g6/FiqeXMGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHiAMAYADxDGAOABwhgAPEAYA4AHCGMA8IC3vSlmnZGKMmwpMDbDfk198YAz1ZdeHDfVj5fav7WpmttsGwQVxgPY5ixJoVHbvDVuO0Zwi60fhyTJ2Y5hnbVL2eqDkVHjESRr5xJzf41i+/PPldp7l1gEsq91UFpq28DYX8OFs4hAyyms4cnHmTEAeIAwBgAPEMYA4AHCGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHiAMAYADxDGAOABwhgAPOBto6CZ/z2i4uLMXitSYftrSmBsmhMYu80MVdqbrgxUFpnqR2O2/Yf7bfWSVPZ3W6Og8KCty854xL5242Fbm53iYduYIv+wNf4pvmgqlyQFQ0n7RihonBkDgAcIYwDwAGEMAB4gjAHAA4QxAHiAMAYADxDGAOABwhgAPEAYA4AHCGMA8ABhDAAe8LY3RfHgmIqLxzKqdSFbrwJJUpFtm9GZtm9V3+22PhOS1PfNYVP9HYl/mOo7/n6rqV6SRk9FTfWRf9he30dm2ddurMxWX9JjG1O5bfcq7hkybiG5Ids2rn/AVB+URkz1khRopnkbC9dvb47ihm09PIKo7fkaFNt/Tl1JOPNiQ1sUzowBwAPmMD58+LBWrFihqqoqBUGgffv2TXh8zZo1CoJgwm3JkiW5Gi8AFCRzGA8MDGjhwoXavn37VWsefvhhXbhwIX07cODADQ0SAAqd+T3jhoYGNTQ0XLMmEokoHo9nPSgAmG7y8p5xa2urKioqdOedd+qpp55SV1dXPg4DAAUj55+maGho0Pe//33V1NSovb1dP/vZz/Tggw/q+PHjikQu/w1vMplUMvk/vzHt7e3N9ZAAwHs5D+NVq1al/z1//nwtWrRINTU12r9/v1auXHlZfXNzs37+85/nehgAcFPJ+0fbEomEampqdObMmSs+vmnTJvX09KRvHR0d+R4SAHgn7xd9dHd3q6OjQ4lE4oqPRyKRK759AQDTiTmM+/v79fHHH6e/bm9v1wcffKDZs2dr9uzZampq0mOPPaZEIqGzZ8/qpz/9qebMmaNHH300pwMHgEJiDuP3339fy5cvT3/d2NgoSVq9erV27NihkydPavfu3fr888+VSCS0fPly7d27V7GY8e/KA8A0Yg7juro6Oeeu+vjBgwdvaECXBCmnIHX149woF9h6IqTCxl4Ws0zlkqSaqm5T/b8n/stU/07R/zLVS9Kn/z3XVB8atX2fRm6xr/FozLqN7VcjY6X57xLgkiOm+tSwrW9JVjPI89uF1j4TUhbzDtlmHoyNm+olKbhG/t1ILb0pAMADhDEAeIAwBgAPEMYA4AHCGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHiAMAYADxDGAOCBvPczzlaQcgrGM2yykceGQpeEkrZjhPvsxzjXOdtU/3Zwl6n+079/xVQvSZEeW+Of4kHb92m8xLZ/SZJxk+IhW33RqPH5ZGw6JUlBpMRUHxotte2/1N70JwiHzduYZDEm69liYG12VFxkPIKtyZilljNjAPAAYQwAHiCMAcADhDEAeIAwBgAPEMYA4AHCGAA8QBgDgAcIYwDwAGEMAB4gjAHAA/72phgZV1A0llltFvt3o7bXoXDIdpTY3+yvcy4UNdWf/39zTfXRflO5JKmsy9anITD2CUmFs1k9a78M494z7YnyTy5kX+vQjDLbBkW2HgrZ9JlwpbZ+GVaBZtq3KbX15JBx3i6cRQRaltvwVOLMGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHiAMAYADxDGAOABwhgAPEAYA4AHCGMA8IC3vSlC/UMKhVKZFY+Pm/dvvSY9cLa+EaXd9te5IMPpXjJWauvREMqs1ceXtrH1aRj6im3eg1W2/UvS6K22iYx229a6pNfWByIasdVLkoJsenIYOPv31dpXxBn7teDaODMGAA8QxgDgAcIYADxAGAOABwhjAPAAYQwAHiCMAcADhDEAeIAwBgAPEMYA4AHCGAA8QBgDgAe8bRTkevrkQsnMaoeGzfsPymyNf4IiWzOY0IwSU70khUZsjVqKrC+l9t4x5mZE/dW2/X/1W+dsG0haEf/QVL/vwv8x1Z/vs01i1jl7oyDrczb1eY+pPhSJmOolKdBMW72x2ZHr7zfVS5IbziwDLgmixp/r4izWriScebGh+RdnxgDgAVMYNzc365577lEsFlNFRYUeeeQRnT59ekKNc05NTU2qqqpSNBpVXV2dTp06ldNBA0ChMYVxW1ub1q1bp6NHj6qlpUVjY2Oqr6/XwMBAuuall17S1q1btX37dh07dkzxeFwPPfSQ+vr6cj54ACgUpveM33777Qlf79y5UxUVFTp+/LgeeOABOee0bds2bd68WStXrpQk7dq1S5WVldqzZ4+efvrp3I0cAArIDb1n3NPzxS8VZs+eLUlqb29XZ2en6uvr0zWRSETLli3TkSNHrriPZDKp3t7eCTcAmG6yDmPnnBobG3Xfffdp/vz5kqTOzk5JUmVl5YTaysrK9GNf1tzcrPLy8vStutr463gAKABZh/H69ev14Ycf6ve///1lj335Iy/Ouat+DGbTpk3q6elJ3zo6OrIdEgDctLL6nPGGDRv01ltv6fDhw5o7d276/ng8LumLM+REIpG+v6ur67Kz5UsikYgiWXwmEgAKienM2Dmn9evX64033tChQ4dUW1s74fHa2lrF43G1tLSk7xsZGVFbW5uWLl2amxEDQAEynRmvW7dOe/bs0R/+8AfFYrH0+8Dl5eWKRqMKgkAbN27Uli1bNG/ePM2bN09btmxRWVmZnnzyybxMAAAKgSmMd+zYIUmqq6ubcP/OnTu1Zs0aSdLzzz+voaEhrV27VhcvXtTixYv1zjvvKBaL5WTAAFCITGHs3PWbGwRBoKamJjU1NWU7pi+ONTwsF2R2YXdqcNC8/5DxuvpgfNx8DDPbkORCxg2yMF5iO8bILYaL8SX938r/MtVL0rpbbL/kTaYMvQQk7bhl7vWL/kXK3CRE0tiYqdwlbT0abKvwhVC0NIutMmftMyFJqWFbD49QyLYWwZj95zrIIAezqaU3BQB4gDAGAA8QxgDgAcIYADxAGAOABwhjAPAAYQwAHiCMAcADhDEAeIAwBgAPEMYA4IGs+hlPhiASURCUZFQbStmvxA9KjT2Uzb0s7GMKjdi2sbZEyLDVxwQusB0k0l1kqt93YaGpXpJGne0Yb53/36b6ks+Na23oP5DepqzMVF+UMh6jxNaPQ8riZ8I6pmjUVi/72WJg7Y1ebHsuSZIzZIGlljNjAPAAYQwAHiCMAcADhDEAeIAwBgAPEMYA4AHCGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHjA30ZBZVEFoQybfmTTFKUksyZEl7iQ7XUrSI6b6iWppGfEVO/68/9aWjxga6Tiimxr0Tky11QvSTvKbzfVl/TYGv/M+tTWUSlVZNu/JI189TbjFrb6YMzevCiUHLPVDyRN9fbvksw/20HYVu/CWUSg5cfOsAycGQOABwhjAPAAYQwAHiCMAcADhDEAeIAwBgAPEMYA4AHCGAA8QBgDgAcIYwDwAGEMAB7wtjdF3jnjtfspW7+C0KDtun1JUr/tGOY5GPtrSFJRqe1a/+J+W/3MDvtTMBU29gkxfp+svSaGv2LvjdJXbZvD6CzbHKz9OCRp1llbP5WZn9r2XzQ4bNtgmuHMGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHiAMAYADxDGAOABwhgAPEAYA4AHCGMA8IC3vSnc4JBckNm18m5oyH6A0oip3Hylf2DvDeAGB231w7b+F4FxzpIUmjnDVt9rm7cbyqJfwcioqTwwzmHkq7eZ6q19JiQpuO+iqX5Z1VlT/Z/P/5upXpJ6daupvqzTFh8h4/NVklxPr22DaNRUHoSziEDDNsF45j1FODMGAA8QxgDgAVMYNzc365577lEsFlNFRYUeeeQRnT59ekLNmjVrFATBhNuSJUtyOmgAKDSmMG5ra9O6det09OhRtbS0aGxsTPX19RoYGJhQ9/DDD+vChQvp24EDB3I6aAAoNKZ3r99+++0JX+/cuVMVFRU6fvy4HnjggfT9kUhE8Xg8NyMEgGnght4z7unpkSTNnj17wv2tra2qqKjQnXfeqaeeekpdXV1X3UcymVRvb++EGwBMN1mHsXNOjY2Nuu+++zR//vz0/Q0NDXrttdd06NAhvfzyyzp27JgefPBBJZNX/lhLc3OzysvL07fq6upshwQAN62sP2e8fv16ffjhh3rvvfcm3L9q1ar0v+fPn69FixappqZG+/fv18qVKy/bz6ZNm9TY2Jj+ure3l0AGMO1kFcYbNmzQW2+9pcOHD2vu3LnXrE0kEqqpqdGZM2eu+HgkElEkYr8YAQAKiSmMnXPasGGD3nzzTbW2tqq2tva623R3d6ujo0OJRCLrQQJAoTO9Z7xu3Tr97ne/0549exSLxdTZ2anOzk4N/fNy5P7+fj333HP685//rLNnz6q1tVUrVqzQnDlz9Oijj+ZlAgBQCExnxjt27JAk1dXVTbh/586dWrNmjYqKinTy5Ent3r1bn3/+uRKJhJYvX669e/cqFovlbNAAUGjMb1NcSzQa1cGDB29oQOljJZNyQWZNNlLD9mYz5o+RTML72m7A1ijIOu/QeGaNlyYozm8vqdTnPeZt3FU+mXM1ReYj2BoFjc7KvBnMJdbGPz++rdV4hDpjvdRWfoup3hUbm2EZGzxJWTzHjfu3NpHKJ3pTAIAHCGMA8ABhDAAeIIwBwAOEMQB4gDAGAA8QxgDgAcIYADxAGAOABwhjAPAAYQwAHshv44EbEJSWKgiVZFSbzStKYOw1EZRkNpa06/TxuKJo1FQeCtlmbp2zJAXhsHkbi1C01LyNC2w9EYKIce3GbWtX0mPs0SDp6Pl/M9UXBcvyun8pi3kYv0/mdZAUKiszHsP4HC+2dy6xPP8stZwZA4AHCGMA8ABhDAAeIIwBwAOEMQB4gDAGAA8QxgDgAcIYADxAGAOABwhjAPAAYQwAHvC3N0V5TEEos+vMA+P165LM16S7cP6/VYHxGMHYuO0A2VyHn+d5W+cs2eftIrb+GqER2/5j54zrIKnnz7ea6v94yy2m+pLP7f0yyo3zsH6fXGyGqV6SQtbeKNaf61J7vwwVWb639KYAgJsKYQwAHiCMAcADhDEAeIAwBgAPEMYA4AHCGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHjA20ZBqWipUkUZNgpyzrx/FxgbqUzCy5YrsTVFyWbeVs7UFCULk9CAybrWwXjKVF/WOWKql6TwgG3e4xHbE7AoaZuDJIX7xkz11u+Ti2TRlMf4M2H+uc7m+W05hqGWM2MA8ABhDAAeIIwBwAOEMQB4gDAGAA8QxgDgAcIYADxAGAOABwhjAPAAYQwAHvDucmj3z0t8x8aTGW8zKZdD5//KY8l4BeukXA6t/F4OHYxPwhyMa+3Gi0z142P2S4/HRo2XQ4ds501u1D6mYMx2ObQbH7Xtf3zcVC/Zn+Pmn+tsnt+GY1zKMZfBPLwL476+PknS4dP/McUjAYDc6OvrU3l5+TVrApdJZE+iVCql8+fPKxaLKfjSK1Bvb6+qq6vV0dGhWbNmTdEIJ9d0nLM0Pec9HecsFfa8nXPq6+tTVVWVQtf53413Z8ahUEhz5869Zs2sWbMKbtGuZzrOWZqe856Oc5YKd97XOyO+hF/gAYAHCGMA8MBNFcaRSEQvvPCCIpHMms4Xguk4Z2l6zns6zlmavvP+Mu9+gQcA09FNdWYMAIWKMAYADxDGAOABwhgAPHDThPErr7yi2tpalZaW6u6779af/vSnqR5SXjU1NSkIggm3eDw+1cPKucOHD2vFihWqqqpSEATat2/fhMedc2pqalJVVZWi0ajq6up06tSpqRlsjlxvzmvWrLls7ZcsWTI1g82R5uZm3XPPPYrFYqqoqNAjjzyi06dPT6gpxLW2uCnCeO/evdq4caM2b96sEydO6P7771dDQ4POnTs31UPLq7vuuksXLlxI306ePDnVQ8q5gYEBLVy4UNu3b7/i4y+99JK2bt2q7du369ixY4rH43rooYfSPUxuRtebsyQ9/PDDE9b+wIEDkzjC3Gtra9O6det09OhRtbS0aGxsTPX19RoYGEjXFOJam7ibwLe+9S33zDPPTLjv61//uvvJT34yRSPKvxdeeMEtXLhwqocxqSS5N998M/11KpVy8Xjcvfjii+n7hoeHXXl5ufv1r389BSPMvS/P2TnnVq9e7b773e9OyXgmS1dXl5Pk2tranHPTY62vx/sz45GRER0/flz19fUT7q+vr9eRI0emaFST48yZM6qqqlJtba0ef/xxffLJJ1M9pEnV3t6uzs7OCWsfiUS0bNmygl/71tZWVVRU6M4779RTTz2lrq6uqR5STvX09EiSZs+eLWl6r/Ul3ofxZ599pvHxcVVWVk64v7KyUp2dnVM0qvxbvHixdu/erYMHD+rVV19VZ2enli5dqu7u7qke2qS5tL7Tbe0bGhr02muv6dChQ3r55Zd17NgxPfjgg0omM+/x7TPnnBobG3Xfffdp/vz5kqbvWv8r77q2Xc2X22k65y67r5A0NDSk/71gwQLde++9+trXvqZdu3apsbFxCkc2+abb2q9atSr97/nz52vRokWqqanR/v37tXLlyikcWW6sX79eH374od57773LHptua/2vvD8znjNnjoqKii57dezq6rrsVbSQzZgxQwsWLNCZM2emeiiT5tKnR6b72icSCdXU1BTE2m/YsEFvvfWW3n333QmtclnrmyCMS0pKdPfdd6ulpWXC/S0tLVq6dOkUjWryJZNJffTRR0okElM9lElTW1ureDw+Ye1HRkbU1tY2rda+u7tbHR0dN/XaO+e0fv16vfHGGzp06JBqa2snPM5a6+b4NMXrr7/uwuGw+8///E/317/+1W3cuNHNmDHDnT17dqqHljfPPvusa21tdZ988ok7evSo+853vuNisVjBzbmvr8+dOHHCnThxwklyW7dudSdOnHCffvqpc865F1980ZWXl7s33njDnTx50j3xxBMukUi43t7eKR559q41576+Pvfss8+6I0eOuPb2dvfuu++6e++9191+++039Zx//OMfu/Lyctfa2uouXLiQvg0ODqZrCnGtLW6KMHbOuV/96leupqbGlZSUuG9+85vpj8QUqlWrVrlEIuHC4bCrqqpyK1eudKdOnZrqYeXcu+++6/TFn3udcFu9erVz7ouPPL3wwgsuHo+7SCTiHnjgAXfy5MmpHfQNutacBwcHXX19vbvttttcOBx2d9xxh1u9erU7d+7cVA/7hlxpvpLczp070zWFuNYWtNAEAA94/54xAEwHhDEAeIAwBgAPEMYA4AHCGAA8QBgDgAcIYwDwAGEMAB4gjAHAA4QxAHiAMAYADxDGAOCB/w8JSMbzYq+5KwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.imshow(y_esm.squeeze())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "9bb91357-9533-4796-82b4-0eb54947f63e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plot_gdm(data, ax=None, vmin=None, vmax=None):\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots();\n",
    "        \n",
    "    ax.axis(\"off\");\n",
    "    return ax.imshow(data.squeeze().detach().numpy(), vmin=vmin, vmax=vmax)\n",
    "\n",
    "\n",
    "def plot_dict(input_dict, figsize=[5,5], cbar_base_col=0):\n",
    "    \"\"\" Plots the input of a dict. The dict should contain GDM data.\n",
    "    \n",
    "    Input:\n",
    "    input_dict: ...\n",
    "    figsize: ... (automatically in the future?)\n",
    "    cbar_base_col: Which col should be used as baseline vmin and vmax of plt.imshow.\n",
    "    \"\"\"\n",
    "    \n",
    "    rows = len(next(iter(input_dict.values())))\n",
    "    cols = len(input_dict)\n",
    "\n",
    "    fig, axes = plt.subplots(rows, len(input_dict), constrained_layout=True, figsize=figsize)  \n",
    "    \n",
    "    # Loop over rows (amount of samples)\n",
    "    for row in range(rows):\n",
    "        vmin=None\n",
    "        vmax=None\n",
    "\n",
    "        vmin = input_dict[list(input_dict.keys())[cbar_base_col]][row].min()\n",
    "        vmax = input_dict[list(input_dict.keys())[cbar_base_col]][row].max()   \n",
    "        \n",
    "        # Loop over cols (amount of models)\n",
    "        for col, i in zip(input_dict, range(cols)):\n",
    "            try:\n",
    "                # Try for multiple rows\n",
    "                plot_gdm(input_dict[col][row].squeeze(), axes[row][i], vmin=vmin, vmax=vmax);\n",
    "                if row == 0:\n",
    "                    axes[row][i].set_title(col)\n",
    "        \n",
    "            except:\n",
    "                # otherwise only single row\n",
    "                plot_gdm(input_dict[col][row].squeeze(), axes[i], vmin=vmin, vmax=vmax);\n",
    "                if row == 0:\n",
    "                    axes[i].set_title(col)\n",
    "    \n",
    "    return fig, axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "52c60f95-d4cb-4131-8419-3a0eb1382b03",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([1, 20, 6, 5])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_list[0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "b549a267-13b5-4809-a973-e6c61023f8b6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAGbCAYAAAAoSIKLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqn0lEQVR4nO3de5Bcdfnn8eecvnfPZO6Z3K+EJEBC1AhLRAIighr+EKG2SkUWBPlxcasWVqpw/yAUbEFAqkCUYmvdaLGYWssii6BukJuKJCng9xNQViOXJORKLnOf6bl093f/kMwm5nk6czo9mfQ371cVf/DMmfM9ffp0P3Pg8/2ewDnnBAAAeCWc6AMAAADVR4MHAMBDNHgAADxEgwcAwEM0eAAAPESDBwDAQzR4AAA8RIMHAMBDNHgAADxEgz/Mxo0bZfXq1dLV1XVC9z1nzhxZtWpV1cdE7Vi9erUEQTDRh6E6mY8NgI0Gf5iNGzfK3XffPW4Nfrz2DQDAP6PBVyifz0/0IQA1aWBgYKIPATgl0OA/tnr1avnud78rIiJz586VIAgkCAL53e9+N/qf0NevXy+f+MQnJJ1Oy9133y3btm2TIAjkpz/96VH7C4JAVq9efcx9H27Dhg3yyU9+UjKZjCxatEjWrl07ni8ZE+TXv/61LFu2TFKplMydO1e+//3vH7WNc04ee+wxWbZsmWQyGWlqapIrr7xSPvjgg6O23bBhg1x88cXS0NAg2WxWFi9eLPfdd98R2zzzzDNy3nnnSTablfr6ernkkktk06ZNFR1blOO78MIL5ayzzpI//OEPsmLFCslms3LddddFOV2oIe+++6587Wtfk8mTJ0sqlZLFixfLj370o9Gfl0oluffee2XhwoWSyWSksbFRli5dKo888sjoNof+l9Dbb78tV111lTQ0NEhzc7PcdtttUigUZMuWLXLZZZdJfX29zJkzRx544IGJeKm1wcE559yOHTvcd77zHScibv369W7Tpk1u06ZNrru7282ePdtNnTrVzZs3z61du9a9/PLL7rXXXnNbt251IuJ+8pOfHLU/EXF33XXXMfftnHOzZ892M2bMcGeccYZ74okn3HPPPeeuuuoqJyLu97///Qk8CxhvL7zwgovFYu78889369evd7/4xS/cpz/9aTdr1ix3+MfxhhtucIlEwt1+++1uw4YNbt26dW7RokWuvb3d7d27d3S7H//4xy4IAnfhhRe6devWuRdeeME99thj7uabbx7d5mc/+5kTEfeFL3zBPf300+7nP/+5+9SnPuWSyaR75ZVXIh9blONbuXKla25udjNnznSPPvqoe/nll7mmPfXOO++4hoYGt2TJEvfEE0+43/72t+722293YRi61atXO+ecu++++1wsFnN33XWXe/HFF92GDRvcww8/PPpz55y76667nIi4hQsXunvuucc9//zz7o477nAi4m699Va3aNEi94Mf/MA9//zz7tprr3Ui4p566qmJetknNRr8YR588EEnIm7r1q1H1GfPnu1isZjbsmXLEfWxNvhy+z60/3Q67bZv3z5ay+fzrrm52d14443H85Jwkjn33HPdtGnTXD6fH6319PS45ubm0Sa6adMmJyLuoYceOuJ3d+zY4TKZjLvjjjucc8719va6SZMmufPPP9+VSiV1vGKx6KZNm+aWLFniisXiaL23t9dNnjzZrVixItKxRTk+5/7R4EXEvfjii2M+R6hNl156qZsxY8bojcsht956q0un066jo8OtWrXKLVu2rOx+DjX4f76+li1bNnqTdMjIyIhra2tzV1xxRfVeiEf4T/RjtHTpUjn99NPHbf/Lli2TWbNmjf57Op2W008/XbZv3z5uY+LE6u/vl9dff12uuOIKSafTo/X6+nq5/PLLR//9V7/6lQRBIN/4xjekUCiM/jNlyhQ5++yzR//XzsaNG6Wnp0duvvlmM+W+ZcsW2b17t1x99dUShv//415XVydf/epXZfPmzTIwMDDmY4tyfIc0NTXJ5z73uUpPG2rA4OCgvPjii/KVr3xFstnsEdfFl770JRkcHJTNmzfLOeecI2+99ZbcfPPN8txzz0lPT4+5z3+eWbR48WIJgkC++MUvjtbi8bicdtppfE8a4hN9ALVi6tSp47r/lpaWo2qpVIown0c6OzulVCrJlClTjvrZ4bWPPvpInHPS3t6u7mfevHkiIrJ//34REZkxY4Y55sGDB0VEv36nTZsmpVJJOjs7xTk3pmOLcnyHjPdnBxPv4MGDUigU5NFHH5VHH31U3ebAgQNy5513Si6XkyeffFIef/xxicVicsEFF8iaNWtk+fLlR2zf3Nx8xL8nk0nJZrNH/AF6qF7uD4VTGQ1+jLQ7pEMX2tDQ0BH1Q1+qwOGampokCALZu3fvUT87vNba2ipBEMgrr7wiqVTqqG0P1dra2kREZOfOneaYh/5w3LNnz1E/2717t4RhKE1NTeKcG9OxRTm+Q5hD77+mpiaJxWJy9dVXyy233KJuM3fuXInH43LbbbfJbbfdJl1dXfLCCy/I9773Pbn00ktlx44dks1mT/CR+43/RH+YQ19MY71rbm9vl3Q6LW+//fYR9V/+8pfHvW/4J5fLyTnnnCPr16+XwcHB0Xpvb688++yzo/++atUqcc7Jrl27ZPny5Uf9s2TJEhERWbFihTQ0NMjjjz8uzjl1zIULF8r06dNl3bp1R2zT398vTz311GiyfqzHFuX4cOrIZrNy0UUXyZ/+9CdZunSpel3883+lbGxslCuvvFJuueUW6ejokG3btk3MwXuMO/jDHPpieuSRR+Saa66RRCIhCxcuNLc/9P8h165dK/Pnz5ezzz5bXnvtNVm3bt2Y911fXz8+LwYnpXvuuUcuu+wyueSSS+T222+XYrEoa9askVwuJx0dHSIi8pnPfEa+/e1vy7XXXitvvPGGXHDBBZLL5WTPnj3yxz/+UZYsWSI33XST1NXVyUMPPSTXX3+9fP7zn5cbbrhB2tvb5b333pO33npLfvjDH0oYhvLAAw/I17/+dVm1apXceOONMjQ0JA8++KB0dXXJ/fffH+nYohwfTi2PPPKInH/++fLZz35WbrrpJpkzZ4709vbKe++9J88++6y89NJLcvnll8tZZ50ly5cvl7a2Ntm+fbs8/PDDMnv2bFmwYMFEvwT/TFy+7+R05513umnTprkwDJ2IuJdfftnNnj3bffnLX1a37+7udtdff71rb293uVzOXX755W7btm1HpeitfTvnzP2vXLnSrVy5ssqvEBPtmWeecUuXLnXJZNLNmjXL3X///aPJ4cOtXbvWnXvuuS6Xy7lMJuPmz5/vvvnNb7o33njjiO1+85vfuJUrV7pcLuey2aw744wz3Jo1a47Y5umnn3bnnnuuS6fTLpfLuYsvvti9+uqrFR/bWI9v5cqV7swzzzye04UasnXrVnfddde56dOnu0Qi4dra2tyKFSvcvffe65xz7qGHHnIrVqxwra2to9fYt771Lbdt27bRfRy63vbv33/Evq+55hqXy+WOGpNrzBY4Z/y3PQAAULP4f/AAAHiIBg8AgIdo8AAAeIgGDwCAh2jwAAB4iAYPAICHaPAAAHhozCvZ/a/3lqv1ncNHPyRFROTJDz6t1lvXZNR68OqbYz0UeOj50i8q/t2LXrpdrZ/bsk2tD5X0y/5/b9av2Wkv6eNO+vMBtV58Tx9XRERKRftnOClVem3OefI+tb509i61/vb26Wq99aWj1/sXEWnduE+tu51HP3dARKRUbpnsQL/XC2Ix+3e0sYvG9e1KXo5b1bET+veS9SyH5wb+5zGPjTt4AAA8RIMHAMBDNHgAADxEgwcAwENjDtn1FvVw3G3NH6j19kS3Wr/v/H+v1qe/OtYjAY7Uku5X68ty29V6VzGn78jIyqQPjug/6O7V61a4BqcUK0z3L9N/p9Yfk4vU+rbmeWrdJRPRDqjcc8WccfGHesDL3o9x7ZtjW9vXyLjVHLtkbB+r4Jg+xh08AAAeosEDAOAhGjwAAB6iwQMA4CEaPAAAHhpzin7Nm5eq9eKy59X6Sx2L1Hq6o0ySE6jA63+Zr9YHF+sp40JJ/7u2fqu+TGXigD4jpHSwQz+gcmllnDKspWettPxftk9T663Gd2YwrM/uMK8+Y8nTihjLrUoQMTlu7ifivedEjVvNsaPOWBgD7uABAPAQDR4AAA/R4AEA8BANHgAAD9HgAQDw0JhT9NOeTKr1B0p6uj62Na3W5/z3TWMdEhiThnf0y3j7Fn0N72S3nnqdvvmAPsDuj9RyqVA49sHhlNX6UkqtW2vLW2n5xi0Daj3o6VPrrqSnt4OYPkukrIi/E4ixvbXGe5US5RM1blXHtoSV34dzBw8AgIdo8AAAeIgGDwCAh2jwAAB4iAYPAICHxpyiT/36dbW+8E19/eTCrt2VHREQ0dQ/dOo/KOqp5LCzR60X9uhpeSkVKzksnOJaN+5T6y6pPyPBWlveSsuX+vV0vSsaa6CXU8ka7Djp8a4CAOAhGjwAAB6iwQMA4CEaPAAAHqLBAwDgoTGn6C2k5THRgp3GWvHdRlqeNeRxAride6Jtb9WNteWttLwrGrM+XJl0faCPbq6zbog8do2PW9WxQ6MdG+//WHAHDwCAh2jwAAB4iAYPAICHaPAAAHiIBg8AgIeOO0UPTLRiZ7f+g3KpYWCclfJ5/QfOyMsHgV6ORUx2W9e9Na6IiFi/ox9T9cau8XGrOXbJ2D5WwTF9jDt4AAA8RIMHAMBDNHgAADxEgwcAwEM0eAAAPESKHrWvZKwFDUykwLh/clW6Xq39G2udm+ntsvuKeA8YdexaH7eaY4eVp+Ut3MEDAOAhGjwAAB6iwQMA4CEaPAAAHqLBAwDgIRo8AAAeYpocAIwD8yExUadDRXzYTCDG9uUepFKlqVuRx67xcas6tiWs/D6cO3gAADxEgwcAwEM0eAAAPESDBwDAQzR4AAA8RIoeAE4GlTzoBCiDKwoAAA/R4AEA8BANHgAAD9HgAQDwEA0eAAAPkaIHgHHgikXjByW9Hji9bK11Xq1xJ3LsGh+3qmOHRjsulXnfjoE7eAAAPESDBwDAQzR4AAA8RIMHAMBDNHgAADxEih4AxoOVonZ6ilrE2j4Y53EncuwaH7eaY5eM7WMVHNPHuIMHAMBDNHgAADxEgwcAwEM0eAAAPESDBwDAQ6ToAWA8BNb9k7VGubG9uR9r3Ijp7Ykcu9bHrebYYeVpeQt38AAAeIgGDwCAh2jwAAB4iAYPAICHaPAAAHiIFD0AjIMgFtN/YK13XqV0dSARx53IsWt83KqObQkrvw/nDh4AAA/R4AEA8BANHgAAD9HgAQDwEA0eAAAP0eABAPAQDR4AAA/R4AEA8BANHgAAD9HgAQDwEA0eAAAPsRY9AIwDVywaPyjp9cDpZWut82qNO5Fj1/i4VR07NNpxqcz7dgzcwQMA4CEaPAAAHqLBAwDgIRo8AAAeosEDAOAhUvQAMB6sFLXTU9Qi1vbBOI87kWPX+LjVHLtkbB+r4Jg+xh08AAAeosEDAOAhGjwAAB6iwQMA4CEaPAAAHiJFDwDjIbDun6w1yo3tzf1Y40ZMb0/k2LU+bjXHDitPy1u4gwcAwEM0eAAAPESDBwDAQzR4AAA8RIMHAMBDpOgBYBwEsZj+A2u98yqlqwOJOO5Ejl3j41Z1bEtY+X04d/AAAHiIBg8AgIdo8AAAeIgGDwCAh2jwAAB4iAYPAICHaPAAAHiIBg8AgIdo8AAAeIgGDwCAh2jwAAB4iLXoAWAcuGLR+EFJrwdOL1trnVdr3Ikcu8bHrerYodGOS2Xet2PgDh4AAA/R4AEA8BANHgAAD9HgAQDwEA0eAAAPkaIHgPFgpaidnqIWsbYPxnnciRy7xset5tglY/tYBcf0Me7gAQDwEA0eAAAP0eABAPAQDR4AAA/R4AEA8BApegAYD4F1/2StUW5sb+7HGjdiensix671cas5dlh5Wt7CHTwAAB6iwQMA4CEaPAAAHqLBAwDgIRo8AAAeosEDAOAhpskBwDgIYjH9B9YDTao0fSqQiONO5Ng1Pm5Vx7aEld+HcwcPAICHaPAAAHiIBg8AgIdo8AAAeIgGDwCAh2jwAAB4iAYPAICHaPAAAHiIBg8AgIdo8AAAeIgGDwCAh1iLHgDGgSsWjR+U9Hrg9LK11nm1xp3IsWt83KqOHRrtuFTmfTsG7uABAPAQDR4AAA/R4AEA8BANHgAAD9HgAQDwECl6ABgPVora6SlqEWv7YJzHncixa3zcao5dMraPVXBMH+MOHgAAD9HgAQDwEA0eAAAP0eABAPAQDR4AAA+RokftC4yUadnUMPBPrOuo4v1Z90/WGuXG9uZ+rHEjprcncuxaH7eaY4dVvv6EO3gAALxEgwcAwEM0eAAAPESDBwDAQzR4AAA8RIoetY+0PKqhytdRkDC+Xq01x6uUog7CiONO5Ng1Pm41xw6sWRxh5ffh3MEDAOAhGjwAAB6iwQMA4CEaPAAAHqLBAwDgIVL0AE68cuu+W2n2qM8cmOBnFJip6FjEBHfUFHXJWAM96rgTOXatjDvRYx9rl1XfIwAAmHA0eAAAPESDBwDAQzR4AAA8RIMHAMBDgXMs5A0AgG+4gwcAwEM0eAAAPESDBwDAQzR4AAA8RIMHAMBDNHgAADxEgwcAwEM0eAAAPDTmx8VeEl41nsdxSuq6+jy1vu9zw2o9vi+p72hmXi0Xh2Pm2K1tPWq9o6tO39egvq+65gG1PikzaI6t2fSFNZG2P1zkazM0zoszHvt4ItaCqtKjTWOtLfpupk9W64X6lFrvmZdR6/3T9OPsn1XQDyiln9NEdkTfXkRG8gn9B0P6/UhQsM6dXk706PuJ9+v7+et//U/6jo7h0vr/oNbdsP7ag4T+dRyk9PfI4oaG9PqI8R6JSJDUz3mQzRr1tD6G8RoiPwrV+DyUsvq5KDQa9bT+WY8N6ddlvE//7i2nmNPPXTFljD1YVOuJLv07M+jX6xu2HPs7kzt4AAA8RIMHAMBDNHgAADxEgwcAwENjDtnh2Nx5Z6v1Dy/LqfViWk8BWXmr7MIutd4QMdAmItI/rAdDWpt61Xos1EMpyZgeGBku6gGTloweyjuhrDDdiWAF/CxOP7/WRRKk9fBToU4PIfXN0rfvmafvf7BdD2o1Te/W61k9AFqf0INgIiLdw/oxHejTP0eWwbweSi2UjIBYxLfmWKwwnRsxglzGdRnEo31Nu2F9/65gh+wsQSbiZ8X68rLqlph+7+kSet0K0xWyxj2scThhmWByYIQ2rbGL6Wj3z/G4ESKNeu4Owx08AAAeosEDAOAhGjwAAB6iwQMA4CEaPAAAHiJFX4H4nFlq/aMz9WUdU2d3qvW6tJ4kbsv0q/WuIX350Oa0vn2hTCx4zqSDan2goCePp6T1dP17va1qPRnq6e+SqzwRaoq6zGu1lp6tIN0apvU0u7UcaWBs75ob1HrfXL3euVD/qA9M0c9FbFafWp/ZqF8HCxv3qfVJcT1FXxcvk6Iv6Nf5e4k2tZ4v6DNCRrL69d+b02edVPvatJaeNdPySWMp6pjxOQ6NmRTW0rbWfioZ+ziS3WNifEaDon7urKVnrbS8tX1YKDNroGR8VoaiTb+IDenfjYE19nF8X3EHDwCAh2jwAAB4iAYPAICHaPAAAHiIBg8AgIdI0Zfx0X9codZ7FuopyFiTnhg+u2W/Ws/G9TWj9+Xr1fr8SQfUes5IJIfW4ski0l/Qk7aNGX1N8QNDdfr2Kf01xwM9ETouKXpL1HR95P0ba0cb6WYRkbBhklq3UvFDbfr664MtenK8Y5Ge6M3P06+1nHHNntaiX2uzch1qvSWhz+SYkdS3T4f6Ou0iIv0lIwVubW9cy0Ml/Vz0ZfTtO4f0WTCVMmdGWGvLG4n1IBUt4W4m4ssJjXu9uJXgr9L2JSM5biTWAyOBHh/Qr6dwWB/XSsuHefu6tFhN1ErXByNGit6qGzMHxoI7eAAAPESDBwDAQzR4AAA8RIMHAMBDNHgAADxEil5EOq85T60PfVZfd3vGJH2dbmt9d0vMSLmfVq8nmFNG8jgVFtT6SJm16EMz5a7/zdeY1NPWluKJTMtbIqblrXSzlUo21/xuaTTHGJrZpNZ7Z+r7yrfq53FED+PL0Fx9nfWZU/TnIbRm9Gt5fp1+DU5P6fuZEtdnXzTH9P1PieufLRGR3QVjPf2kPqMgkdLTx9b1//f+drXePayvgX/CWLMvrPXgjXS9Sxv1uH0/FxSNz4qVcjdS8aWMPrujlNBfg5VmD4aMNPt4r4FfCeN9s74CT+Qr4A4eAAAP0eABAPAQDR4AAA/R4AEA8BANHgAAD9Vsij48a5H5s92XNKv1ul16YnPoK11qvS2jJ5LnN+gJY0uzsU53nbWGvOiJ1r6inrS20sKJQE8XlzNipOgLxrrezUn9tVnH1DuSjnxMJ4qZls/p6W03tUWt52fozxIQEemar6eM+2bp73mxSV9DPp7VZ04saNevTes5Bta1aa0h3xzXU/G5UL+Wk8Y1OFhmhkfW2FerkbyPGTNCrDXtUzH93CXC6J+XctyQ/jrcsP6emmvXG9ellZYvNOqzAYpp++s+HDaerzGonysrkT8yyTimjL59bFB/7xI9xnMejJkxhZz+uSqmo41bSUOMPHbeWqPeWB9/OPr6+KO/W/FvAgCAkxYNHgAAD9HgAQDwEA0eAAAP0eABAPBQzabo3/+Gvqa3iEjz0n1qffeHerr+c1N2qPV8UU9Hxo3UbiZmJJ5Dffv9w3raesRIrOeM1L11nNa68iIie/L6YuZ1CX2MkrGwsjV2x1DWHPuEMdatDoy1va0Uc2nGZLXeu0B//zoX2uc9P1NPxDa26wnx9nq9nozpqecZ2S613pbU99MU11P0ReNv/66iPqNgxOlfJTuK+kwDKykvYj+joTE2oNa3D7eq9Y6CfqxDRf1YXZWfn+BGjAR6Qa+ba85b+zeS7FZavmgk2UVErK+KoKS/Fy5mzLYxxrDqlvigvr0zlsa3EusFo24Jh6Lf80Ye23gNCetZAcex/j538AAAeIgGDwCAh2jwAAB4iAYPAICHaPAAAHjopE/RH/j2efoP5ujpXxGR/LCe7D5twR61nokZyeaEntq1UvGFkv73krUmfNFI7XYb67UfNJLpoZE63tXXoNZFRDp79X1Nbeoxf0fTktbfh65BfT3s+pSdnq6YsT51ENcvb6suk/XEd/ciPS3fcab+/pXm6deNiMi8yfoa79NzXWp9kjFzoiT62NNS+n7qYvpzFSwf5NvUeirUE+BWffeQfQ1a6mL6a54U11/DQElf/7y3oH+O+ov69sPG7JVKBUn9e0hCY3aHsea8hMa67EX9urfWlS8zqUZi1jroRt0Z6fq4tc66kRyPDek/CIaNuvFZj+WjJdatca1z948f6u+bta698UgRc2wpWPsxdjQG3MEDAOAhGjwAAB6iwQMA4CEaPAAAHqLBAwDgIRo8AAAeOmmmyVnT4TqX6FMHmnL2tJ+pxgM6mlP6lK6WZJ9atx76Yj0kxnrYzNYeffpVpzHtrWhMt7Om1VmPIjjQqR+niEg6ox/rYEG/JAaG9Ck8yVCfVtKYzqv1vmH9gS7HxXoYQ6Cfx7BFf+hQ/3z9AUbWdDhZoF9PCyYf0LcXkakZfRpiKqZPM7Omw1nXZndBn544YEwN22NMY3u/W3+AS8J4yI01XXK3MVUzP2J/9STj+hgNKf0zbz0gKW2dU+NzZL0HlQon2Z8//ReM+62Eca5K+ndjbFB/HdaDY0Ts6XBBXp9CHBpTt0JjbOvBOIExNSwYjvZehMZxOmuKoXHughF7mpy1r4QxXTGWML5/rCmJQ/prCIrG9Lkx4A4eAAAP0eABAPAQDR4AAA/R4AEA8BANHgAAD53wFH1s8QK1bqXl20/TE8kpI80rItKe1lP0ltB4KoD1EAvrwRpvHpyh1keMVLz1kJjd2/QEc5A2HiJRMlLeZcKX/Z16mr0vVhdp7Pyg/kCNTFpPhFozBI5HmNGT40FGf9hIqWWSWu9coL+Wwhx9RsBCIy0/r+6gWi9nqKR/FEPjKR1WErzHeMDKR3n9Ne/t15Pe+/br2ycz+vvak9XH7e7T6yMDxoNVRCRM6Ndaf73+Oy05/eE+9Uk9XT9c1GcgHBzImcdUCdekn0Mz2W09VMRIfLuE/jpKRt0ZqW4R++ExVlo+GNBnNAQF43vZeA32DBijbh1P1AeyWPuPl3ngUEq//qz3wTzb4/BQGQt38AAAeIgGDwCAh2jwAAB4iAYPAICHaPAAAHho3FL0sUZ9Deq//me9nm3U14OfXtet1luMdeVFRGamO9X6PmMN+d3Getz5op6q3pfX99Od1xPDxc36Guf9s/Q0fuqAnsqc+qqevtzxBf1tDKba6/W7Tv21JfqNtc9P19euL+7R19Mv7dH/dnTn6u/n8Sjl9ZR7rF6fEZCfob9/A1P1FOuCafvUupWWn5nuUOsiIiNOf293DOrXSMnp59Fao75zWH8/dvXq1/iBfXrSO7FPvz5GjCTxgZx+7ce69WvTCLiLiEgxrf9O74h+7qxnMVizVHoG9WPt7NHPXaUKjfrsjlLSSF2ba5Qba+ob6e2RBv09KqatMyUSHzDO7ZCeio/1GDN6uvVnLZSG9Dc8iBvfXeloz6xwg/r+XcFYlz+pn6Oyzw/I6MfkUtFmM4TGevfOuL7FmLU1FtzBAwDgIRo8AAAeosEDAOAhGjwAAB6iwQMA4KHjTtFbafntN5+pD5jR142e2dSl1qcY68pPSdlp7GyoJ76t9b735PXXsKOrUa13d+lp28zf9HTujPs3qnXr3BW7oiXNG+acp9Y7U3YSteldvd76pn6+P3R6ujRtLLveP0NPMJeG9HT28QhiRvq0pVEt903Tr4PiFP26mZnr0nef1Gd+pEJ9vXYRkYTTE7SFkv4arGu2Z0S/1vpH9HSwlRCPHdTfj9RBPXEdFPV6MaXfKyT0UySB/SgJGanXxxge1o+1N6an1a0U/aCxn9JwmXXIK1Ay0tWFtF6PhdYzJfTX4eL6ObfS8iMZ+37OfD+MMax10620vDPqUtQHNj/TBmv/VoreZByPiJjr17uYMdMlGe3cGY8oOC7cwQMA4CEaPAAAHqLBAwDgIRo8AAAeosEDAOCh407Rd3x5sVp3n9TXJJ7TpCfEW9P62vKNCT11H4qeRBQR+WhEX187ZqRqt3fq64Dn/2+jWl+4Tl/rvvSXfzOPSRM1LW9p/W+b9HpV9v4Pkyefo9Z3XGmkTnv1Syv2oZ54Hg/FOn0WwcBUPQ3b0KRfg+0p/Vq2Zmv0FfWEu4i9Fn0qpqd9t3RNVusHe3NqfbBLHzuzXU+O1+3UPxOpbv19TXXpMwQCI+kdFPR60UiSi4gMN+jXzkCbfj/SV9Tf564B4+stbnx3GDMEKhUO62vLx4w0trkWvbV2uZHGNteVLxMQj+WNNLsxtsVaW95My1vbR0zRW2vLm9tb4xrvjYiIFI33ZzjiOTLOaVA0PkPGZ2ssuIMHAMBDNHgAADxEgwcAwEM0eAAAPESDBwDAQ2NO0fd87d+p9c7FRiI5o68NvLhhr1qfYyxq3hzXF7PeOdys1kVEduUb1fqHvXpavnennrqf+4Keki795W/m2L7Kva/PHEjs0rP6iR79upi+Rl+XX+6o6LBERCTWrifNe6br668Ptukp1pl1+rVWHxtU61Yifsegfp2JiKRCPS1/YEhPxX/UoV+bxf16Wj63W/+bvenv+muu26o/eyDs1OuuT59p4Ab1z7uVbk7U6a9XRCSV1V9beor+PIREn56iHzLS+AVjaGPZ/4qFQ8Y66EYqOizoKe1gUN9PaKTogyEj1W2tKy9lkt0Dxhry1thp/b0w15a36lbK3WBl3820fCL6mx0M6zNIQn2SjYTGGvXWuZMR63rRr4ux4A4eAAAP0eABAPAQDR4AAA/R4AEA8BANHgAAD405Snjw8rxaXz7rQ7U+PdOl1pdkd6j1pLFQ8t5Cg1r/W98UtS4i8ud9U9V6z0d1ar19k57BjL/0r+YYp5riX99V63P+i14/kVyjnq7ub9f/fo216tdye1aPwxadvp+tA/oMgr8Z68eLiHT26cn+/D69Xv++/hFNH9STuA1b9deW2Kun4uWg8VyFAX0/zkj6uhF9xomVPA6G9e1FRIJB/VwkjX019uif62JOX3+/kNPPaSlZ3fudUkJPiDtjHCsrba3KHuT1hHusR5/pYKa3y7F+x1hb3pohEFlY3ecCHMU4TlcmsR5YM0WMunnurPXuzZkGlV+X3MEDAOAhGjwAAB6iwQMA4CEaPAAAHqLBAwDgoTGn6Nua9BTueY0fqPUz0jvVupVI3jHSotZ3Denreu/sa1TrIiLuFf13Fv+Pd/Rj6uo294WT30iznroebtTTqi2N+przzckBtT5Q0tdT39anPw9h9w79WhYRyWzXk93T/q6nd+s/0JP9YY+ecpf9HWq51GusLW+kiQMjxeys9LTB3H6ozPraRvrYDejvT9ijv7aYsQ5+0lgvvZTL2MdUgZEGffxCVv8OjA/o5yQ01qK31i533fo1Uxoy0t5ir9lurS1vptCNMazZF4GxJnzUDP2JGNdZnwkjRe8KxtjWuc7q32NBSr+OxoI7eAAAPESDBwDAQzR4AAA8RIMHAMBDNHgAADw05hT9yvb31Pp3mrar9a0jelL5z8P6Ot0bu0/Ttz+oryvf3WcnXpv36GlU0vJ+cnE93TpSpyd9V07Vr+WS0/ez6cBctf7Bjja13vBnPSkvItLyjp64Tb+/Xz+mvfv0HRkJcSspbSV6La5MyD2Skp6iL7f/qMcq+UG1HFhrextpaHP7ChUy+v3TSMYYv2RkuOPGfZgx28C8Bsqk6K215aOek6jPKrBYSfMJHdeYEGIm+K3r2Jq5krYGqHx9f+7gAQDwEA0eAAAP0eABAPAQDR4AAA/R4AEA8NCYo4pdBX2d3N8O6Inh1waWq/V/65qp1t/Zo6flc7/LqfW5r+vrT4uIuDc2mz+Df+J9I2o92Z1W6292zFDrQ0X94/Dhu+1qvfGvesK4fbO+FriISPDXrWq90N9v/o66HyOha60t7zUjku8KRlTfWuu+WsfzsfiAkUw3nsdhrUUfjERb/99Mgpd5jkBgzMqQiCl6a433yNvXyLgiZc5d1LFDY+wg6sr8h+2y4t8EAAAnLRo8AAAeosEDAOAhGjwAAB6iwQMA4KExRw//z78uVesvtyxQ68O79PT7pHf1vymaO/UE6aR1m9T6KZgVhiHctletT9msX967YrPUeiyv73/mB3piPff+Af0X9h3U62KvEx5V5PXafXYca3WPp9Q+44Iy1pa30vLBQLRrJkin9Hq5hLj1MyuRb63nbx2TtR9r3KT9PIcJG9daQ75KYwcx436bFD0AADgcDR4AAA/R4AEA8BANHgAAD9HgAQDw0JhT9Kf/y2tq3Z13tloPNrEePE6M4v79aj006rONteVd0VjTfGBAHzc/qNbLpZVPybXiT1Gx/V3jO4A1k6Ka15iRlj8l1eC54A4eAAAP0eABAPAQDR4AAA/R4AEA8BANHgAAD405RW8JNr1VjeMAKmet1WysUV7YZ6wh7/QUvcnYvyvpa4rj1FLq7lHrbnhYrVtrl1try1tpeWc878CN2M8vCBLG2OZv6KKOXevjVnNsa93/IF75rAju4AEA8BANHgAAD9HgAQDwEA0eAAAP0eABAPDQcafogQlnpNlNpNxxArhBK12tp+ilqF+X5Z5toO9fT2+b45ZhJfvHe+xaGVfEnhXhrGcFWGNb30tRv98Owx08AAAeosEDAOAhGjwAAB6iwQMA4CEaPAAAHiJFDwDjwFxz3HjmQZBM6ttHTNGb41byO+M8dq2PK1LmfQv0+2dz7NAY23rWxhhwBw8AgIdo8AAAeIgGDwCAh2jwAAB4iAYPAICHaPAAAHiIaXIAMA6CVEqvWw80saZoWduH+vQpa1JV2QepWGMnE/bvVGPsGh+3mmMHMeN+m2lyAADgcDR4AAA8RIMHAMBDNHgAADxEgwcAwEM0eAAAPESDBwDAQzR4AAA8RIMHAMBDNHgAADxEgwcAwEOsRQ8A48ANDen14WG1HiSTen2cx53IsWt93GqObT1zIIi7Co7qH7iDBwDAQzR4AAA8RIMHAMBDNHgAADxEgwcAwEOk6FH7AiP76ipPn+IUZF1HFXIjBb1e0OuWwEhXV3XcQL/XG/exa3xcETstH3nsUtEYgBQ9AAA4DA0eAAAP0eABAPAQDR4AAA/R4AEA8BApegCnFiMtH8Ri1R0mmYi2fcL4Oo54XOZ+TsKxa31ckTJry1tpeeu6CI2xj2N2B3fwAAB4iAYPAICHaPAAAHiIBg8AgIdo8AAAeIgUPWofa87XnnLJ4Al6P13RWAu8QkE2q9fTxjhGijqIRbwPM9ZTN9c6n8ixa33cao4dJ0UPAADGgAYPAICHaPAAAHiIBg8AgIdo8AAAeIgUPVBtJ2FC3Dwm63iibh/5eMrdW5SqM/ZxpI+rIUgZa5RHPecRX0cQN/Zf7vxN1Ng1Pu4JGTskRQ8AAA5DgwcAwEM0eAAAPESDBwDAQzR4AAA8FDjHQt4AAPiGO3gAADxEgwcAwEM0eAAAPESDBwDAQzR4AAA8RIMHAMBDNHgAADxEgwcAwEM0eAAAPPT/AGlK9u5JQutyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 500x400 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "min_idx = 237\n",
    "max_idx = 11277\n",
    "\n",
    "X_list = [X_all[min_idx].unsqueeze(0), X_all[max_idx].unsqueeze(0)]\n",
    "y_list = [y_all[min_idx].unsqueeze(0).unsqueeze(0), y_all[max_idx].unsqueeze(0).unsqueeze(0)]\n",
    "y_esm_list = [torch.tensor(np.reshape(y_esm_all[min_idx],[30,25], order=\"F\")),\n",
    "              torch.tensor(np.reshape(y_esm_all[max_idx],[30,25], order=\"F\"))]\n",
    "#plot_gdm(y_list[1])\n",
    "\n",
    "\n",
    "data = {\n",
    "        \"truth\": y_list,\n",
    "        \"decoder\": [decoder(x) for x in X_list],\n",
    "        \"esm\": y_esm_list,\n",
    "        #\"decoderdenoise\": [model_denoise(x.squeeze(1)) for x in X_list],\n",
    "        #\"kdm\": [kdm.calculate(x.squeeze()) for x in X_list],\n",
    "        #\"gmrf\": [gmrf.calculate(x.squeeze()) for x in X_list],\n",
    "        }\n",
    "\n",
    "fig, axes = plot_dict(data, [5,4]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "783735d2-2a7d-4b3e-b0b4-f3046b3c6cca",
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import plot_dict\n",
    "\n",
    "X_list = [data_iter.next(), data_iter.next(), data_iter.next()]\n",
    "\n",
    "data = {\n",
    "        \"input\": X_list,\n",
    "        \"decoder\": [model(x.squeeze(1)) for x in X_list],\n",
    "        \"decoderdenoise\": [model_denoise(x.squeeze(1)) for x in X_list],\n",
    "        \"kdm\": [kdm.calculate(x.squeeze()) for x in X_list],\n",
    "        \"gmrf\": [gmrf.calculate(x.squeeze()) for x in X_list],\n",
    "        }\n",
    "\n",
    "fig, axes = plot_dict(data, [5,5]);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}