examples.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dp_model.model_files.sfcn import SFCN\n",
    "from dp_model import dp_loss as dpl\n",
    "from dp_model import dp_utils as dpu\n",
    "import torch\n",
    "import torch.nn.functional as F\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example for Brain Age Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Label shape: torch.Size([1, 40])\n",
      "Input data shape: torch.Size([1, 1, 160, 192, 160])\n",
      "dtype: torch.float32\n",
      "Output shape: torch.Size([1, 40])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAUeUlEQVR4nO3df5Bd5X3f8ffHIkprAcYxMrElsDBWQ9XW2MwW3EJxqAODTBqZSTMW45jYCaPSQfWQxpMqnUzqjJsWZ+zW45ZYVQgpdgYTx2NN1aAYMukfboLxaKkpIAyOLAtrEUYL+EdwXIPMt3/co/ay3NUeoV3dq8fv18zOvef5cc73Prvz2bNn749UFZKkdr1s3AVIkpaWQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXj/UklyVZH+SZ5K8ucf4n0wy03Pf70ny5y+xrpc8V5rLoNcJL8nFSe5O8u0kTyf5iyR/v+f0DwObq+rkqvpSkn1Jfmop65WOt5PGXYB0LJKcCvwx8M+BTwPLgX8EfL/nLl4H7F6a6qTJ4Bm9TnR/C6CqPlVVP6iq71XVXVV1P0CSlyX59SSPJjmY5BNJXpHkR5M8AywD/neSryb5JHAW8N+7Szm/utDBk2zp5v5VkoeSXPXiIflP3V8bDyd521DHK5L8XpLHkzyW5N8mWbZoKyN1DHqd6L4C/CDJrUnWJ3nlnP73dF+XAq8HTgb+c1V9v6pO7sacV1XnVNW7ga8D/6S7lPPbPY7/VQZ/QbwC+E3gD5K8Zqj/QmAvcDrwb4DPJvmxru9W4BDwBuDNwOXAtf0futSPQa8TWlV9B7gYKOB3gdkkO5Kc0Q15F/AfqmpvVT0D/BqwMcmiXLasqj+qqgNV9XxV/SHwl8AFQ0MOAh+tque6/keAK7v61gM3VNV3q+og8B+BjYtRlzTMoNcJr6q+XFXvqarVwN8FXgt8tOt+LfDo0PBHGfxv6gwWQZJrktyX5FtJvtUd//ShIY/VC9858NGuptcBPwI8PjT3vwCvXoy6pGEGvZpSVQ8D/5VB4AIcYBCqh53F4HLJE/Ptou+xkryOwV8Rm4FXVdVpwINAhoatSjK8fVZX034G/zA+vapO675Oraq/0/f4Ul8GvU5oSc5N8itJVnfbZwJXA/d0Qz4F/HKSs5OcDPw74A+r6tA8u3yCwbX8PlYw+MUw2x37vfz/XzCHvRp4X5IfSfJzwN8GdlbV48BdwEeSnNr90/icJG/teWypN4NeJ7q/YvAPzy8m+S6DgH8Q+JWu/xbgk8Dnga8B/wf4F0fY378Hfr27nPL+Ix24qh4CPgJ8gcEviL8H/MWcYV8E1gJPAr8F/NOqeqrru4bB00EfAr4JfAZ4DdIiix88Iklt84xekhpn0EtS4wx6SWqcQS9JjZvINzU7/fTTa82aNeMuQ5JOGPfee++TVbVyVN9EBv2aNWuYnp4edxmSdMJI8uh8fV66kaTGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxk3kK2MlnVjWbLljZPu+G688zpVoFM/oJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDWuV9AnuSLJI0n2JNkyon9DkvuT3JdkOsnFQ337kjxwuG8xi5ckLWzB97pJsgy4CbgMmAF2JdlRVQ8NDfszYEdVVZI3Ap8Gzh3qv7SqnlzEuiVJPfU5o78A2FNVe6vqWeB2YMPwgKp6pqqq21wBFJKkidAn6FcB+4e2Z7q2F0hyVZKHgTuAXxzqKuCuJPcm2TTfQZJs6i77TM/OzvarXpK0oD5BnxFtLzpjr6rtVXUu8A7gg0NdF1XV+cB64Pokl4w6SFVtq6qpqppauXJlj7IkSX30CfoZ4Myh7dXAgfkGV9XngXOSnN5tH+huDwLbGVwKkiQdJ32CfhewNsnZSZYDG4EdwwOSvCFJuvvnA8uBp5KsSHJK174CuBx4cDEfgCTpyBZ81k1VHUqyGbgTWAbcUlW7k1zX9W8Ffha4JslzwPeAd3bPwDkD2N79DjgJuK2qPrdEj0WSNEKvjxKsqp3AzjltW4fufwj40Ih5e4HzjrFGSdIx8JWxktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIa1yvok1yR5JEke5JsGdG/Icn9Se5LMp3k4r5zJUlLa8GgT7IMuAlYD6wDrk6ybs6wPwPOq6o3Ab8I3HwUcyVJS6jPGf0FwJ6q2ltVzwK3AxuGB1TVM1VV3eYKoPrOlSQtrT5BvwrYP7Q907W9QJKrkjwM3MHgrL73XEnS0ukT9BnRVi9qqNpeVecC7wA+eDRzAZJs6q7vT8/OzvYoS5LUR5+gnwHOHNpeDRyYb3BVfR44J8npRzO3qrZV1VRVTa1cubJHWZKkPvoE/S5gbZKzkywHNgI7hgckeUOSdPfPB5YDT/WZK0laWictNKCqDiXZDNwJLANuqardSa7r+rcCPwtck+Q54HvAO7t/zo6cu0SPRZI0woJBD1BVO4Gdc9q2Dt3/EPChvnMlScePr4yVpMYZ9JLUOINekhrX6xq9pB9ua7bcMbJ9341XHudK9FJ4Ri9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjegV9kiuSPJJkT5ItI/rfleT+7uvuJOcN9e1L8kCS+5JML2bxkqSFLfjBI0mWATcBlwEzwK4kO6rqoaFhXwPeWlXfTLIe2AZcONR/aVU9uYh1S5J66nNGfwGwp6r2VtWzwO3AhuEBVXV3VX2z27wHWL24ZUqSXqo+Qb8K2D+0PdO1zeeXgD8Z2i7griT3Jtk036Qkm5JMJ5menZ3tUZYkqY8+nxmbEW01cmByKYOgv3io+aKqOpDk1cCfJnm4qj7/oh1WbWNwyYepqamR+5ckHb0+Z/QzwJlD26uBA3MHJXkjcDOwoaqeOtxeVQe624PAdgaXgiRJx0mfoN8FrE1ydpLlwEZgx/CAJGcBnwXeXVVfGWpfkeSUw/eBy4EHF6t4SdLCFrx0U1WHkmwG7gSWAbdU1e4k13X9W4HfAF4F/E4SgENVNQWcAWzv2k4Cbquqzy3JI5EkjdTnGj1VtRPYOadt69D9a4FrR8zbC5w3t12SdPz4ylhJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY3rFfRJrkjySJI9SbaM6H9Xkvu7r7uTnNd3riRpaS0Y9EmWATcB64F1wNVJ1s0Z9jXgrVX1RuCDwLajmCtJWkJ9zugvAPZU1d6qeha4HdgwPKCq7q6qb3ab9wCr+86VJC2tPkG/Ctg/tD3Ttc3nl4A/eYlzJUmL7KQeYzKirUYOTC5lEPQXv4S5m4BNAGeddVaPsiRJffQJ+hngzKHt1cCBuYOSvBG4GVhfVU8dzVyAqtpGd21/ampq5C8DSSeeNVvuGNm+78Yrj3MlP7z6XLrZBaxNcnaS5cBGYMfwgCRnAZ8F3l1VXzmauZKkpbXgGX1VHUqyGbgTWAbcUlW7k1zX9W8FfgN4FfA7SQAOVdXUfHOX6LFIkkboc+mGqtoJ7JzTtnXo/rXAtX3nSpKOH18ZK0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDWuV9AnuSLJI0n2JNkyov/cJF9I8v0k75/Tty/JA0nuSzK9WIVLkvpZ8MPBkywDbgIuA2aAXUl2VNVDQ8OeBt4HvGOe3VxaVU8ea7GSpKPX54z+AmBPVe2tqmeB24ENwwOq6mBV7QKeW4IaJUnHoE/QrwL2D23PdG19FXBXknuTbJpvUJJNSaaTTM/Ozh7F7iVJR9In6DOirY7iGBdV1fnAeuD6JJeMGlRV26pqqqqmVq5ceRS7lyQdSZ+gnwHOHNpeDRzoe4CqOtDdHgS2M7gUJEk6TvoE/S5gbZKzkywHNgI7+uw8yYokpxy+D1wOPPhSi5UkHb0Fn3VTVYeSbAbuBJYBt1TV7iTXdf1bk/w4MA2cCjyf5AZgHXA6sD3J4WPdVlWfW5qHIkkaZcGgB6iqncDOOW1bh+5/g8Elnbm+A5x3LAVKko6Nr4yVpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGtcr6JNckeSRJHuSbBnRf26SLyT5fpL3H81cSdLSWjDokywDbgLWA+uAq5OsmzPsaeB9wIdfwlxJ0hLqc0Z/AbCnqvZW1bPA7cCG4QFVdbCqdgHPHe1cSdLS6hP0q4D9Q9szXVsfvecm2ZRkOsn07Oxsz91LkhbSJ+gzoq167r/33KraVlVTVTW1cuXKnruXJC2kT9DPAGcOba8GDvTc/7HMlSQtgj5BvwtYm+TsJMuBjcCOnvs/lrmSpEVw0kIDqupQks3AncAy4Jaq2p3kuq5/a5IfB6aBU4Hnk9wArKuq74yau1QPRpL0YgsGPUBV7QR2zmnbOnT/Gwwuy/SaK0k6fnxlrCQ1rtcZvaS2rdlyx8j2fTdeeZwr0VLwjF6SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1LheQZ/kiiSPJNmTZMuI/iT5WNd/f5Lzh/r2JXkgyX1JphezeEnSwhb8hKkky4CbgMuAGWBXkh1V9dDQsPXA2u7rQuDj3e1hl1bVk4tWtSSptz5n9BcAe6pqb1U9C9wObJgzZgPwiRq4BzgtyWsWuVZJ0kvQJ+hXAfuHtme6tr5jCrgryb1JNs13kCSbkkwnmZ6dne1RliSpjz5BnxFtdRRjLqqq8xlc3rk+ySWjDlJV26pqqqqmVq5c2aMsSVIffYJ+BjhzaHs1cKDvmKo6fHsQ2M7gUpAk6TjpE/S7gLVJzk6yHNgI7JgzZgdwTffsm7cA366qx5OsSHIKQJIVwOXAg4tYvyRpAQs+66aqDiXZDNwJLANuqardSa7r+rcCO4G3A3uAvwbe200/A9ie5PCxbquqzy36o5AkzWvBoAeoqp0Mwny4bevQ/QKuHzFvL3DeMdYoqWFrttwxsn3fjVce50ra5StjJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN6/XBI5JOfH7Axw8vz+glqXEGvSQ1rlfQJ7kiySNJ9iTZMqI/ST7W9d+f5Py+cyVJS2vBa/RJlgE3AZcBM8CuJDuq6qGhYeuBtd3XhcDHgQt7zpWkefm/hWPX54z+AmBPVe2tqmeB24ENc8ZsAD5RA/cApyV5Tc+5kqQl1OdZN6uA/UPbMwzO2hcas6rnXACSbAI2dZvPJHmkR23H0+nAk+MuYgTrOnqTWttY6sqHFuybt64jze2572Ppn9TvI4zne/m6+Tr6BH1GtFXPMX3mDhqrtgHbetQzFkmmq2pq3HXMZV1Hb1Jrs66jM6l1weTV1ifoZ4Azh7ZXAwd6jlneY64kaQn1uUa/C1ib5Owky4GNwI45Y3YA13TPvnkL8O2qerznXEnSElrwjL6qDiXZDNwJLANuqardSa7r+rcCO4G3A3uAvwbee6S5S/JIlt6kXlayrqM3qbVZ19GZ1LpgwmpL1chL5pKkRvjKWElqnEEvSY0z6OeRZFmSLyX54277A0keS3Jf9/X2MdW1L8kDXQ3TXduPJfnTJH/Z3b5yQuoa+5olOS3JZ5I8nOTLSf7BhKzXqLomYb1+Yuj49yX5TpIbxr1mR6hrEtbsl5PsTvJgkk8l+RvjXq8X1eg1+tGS/EtgCji1qn46yQeAZ6rqw2Ouax8wVVVPDrX9NvB0Vd3YvZ/QK6vqX01AXR9gzGuW5Fbgf1bVzd0zv14O/GvGv16j6rqBCfgZO6x7C5PHGLzI8XrGvGbz1PVexrhmSVYBfw6sq6rvJfk0gyenrGNC1gs8ox8pyWrgSuDmcdfS0wbg1u7+rcA7xljLxEhyKnAJ8HsAVfVsVX2LMa/XEeqaNG8DvlpVjzJZP2PDdU2Ck4C/meQkBr+wDzBZ62XQz+OjwK8Cz89p35zBu3PeMsY/xQq4K8m93dtGAJzRvW6B7vbVE1IXjHfNXg/MAr/fXYa7OckKxr9e89UFk/EzdthG4FPd/XGv2bDhumCMa1ZVjwEfBr4OPM7gNUR3MVnrZdDPleSngYNVde+cro8D5wBvYvAN/cjxrq1zUVWdz+AdQ69PcsmY6phrVF3jXrOTgPOBj1fVm4HvApPwVtnz1TXu9fp/ustJPwP80bhqGGVEXWNds+4XywbgbOC1wIokP388a+jDoH+xi4Cf6a453w784yR/UFVPVNUPqup54HcZvDPncVdVB7rbg8D2ro4nMni3ULrbg5NQ1wSs2QwwU1Vf7LY/wyBgx71eI+uagPUath74X1X1RLc97jUbWdcErNlPAV+rqtmqeg74LPAPmZz1Agz6F6mqX6uq1VW1hsGfiP+jqn7+8DetcxXw4PGuLcmKJKccvg9c3tWxA/iFbtgvAP9tEuoa95pV1TeA/Ul+omt6G/AQY16v+eoa93rNcTUvvDwy1jUb8oK6JmDNvg68JcnLk4TB9/LLTM56AT7r5oiS/CTw/u5ZN59k8OdhAfuAf3b4GtxxrOf1DM6WYfDn/21V9VtJXgV8GjiLwQ/ez1XV0xNQ1ySs2ZsY/FN9ObCXwbM0XsYY1+sIdX2MMa9XV9vLGby9+Our6ttd21h/xo5Q1yT8jP0m8E7gEPAl4FrgZMa8Xi+o0aCXpLZ56UaSGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMb9XyrNTNq12lbaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEXCAYAAACjyo8UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfZRcVZnv8e+PhMi7EWkxrxA0A0RHMLYhDiM4vJlEJTqDaxKVIHc0chdRQBxFr17R0VFZOiNxMDEDQVAwF1GvkQlGvMgVxguTIBhIME4bAukkQCMkvAQnBJ77x95tTipV3aeT7q7Kye+z1llddfY+p56zu/upXfucs0sRgZmZVdc+zQ7AzMwGlhO9mVnFOdGbmVWcE72ZWcU50ZuZVZwTvZlZxTnRW9NJOlJSSBqan98s6Zxd2M9YSc9IGtL/UZrtuZzorRRJayU9lxPpo5KulnTQQLxWREyNiGtKxnRaYbuHI+KgiHhhIOJqNZKOknSTpKclPS7psrz+JZKukvRQLrtH0tRe9nWRpEckbZa0UNJLCmXHSro1l3VIetdAH5v1Lyd664t3RMRBwETgjcCnayso8d/VAJM0DLgFuBV4JTAa+G4uHgqsA04GXgp8BrhB0pEN9vVW4BLgVOBI4Cjgc7lsKPBj4CbgUGA28F1Jf9b/R2UDxf+Q1mcRsR64GXgtgKTbJH1R0r8DW4CjJL009yo3Slov6QvdQyqShkj6au6FrgHeVtx/3t8HCs8/KOmB3DtdJWmipO8AY4Gf5E8ZH68zBDRS0mJJT+Se6AcL+7xU0g2Srs37XSmpvWwbSLpc0jpJT0m6W9KbC2X7S7pG0pM57o9L6iyUj5T0A0ldkh6U9JE+/QKS9wMbIuKfIuLZiPhjRKwAyM8vjYi1EfFiRNwEPAi8ocG+zgGuioiVEfEk8A95/wDHACOBf46IFyLiVuDfgbN3IWZrEid66zNJY4BpwD2F1WeTensHAw8B1wDbgFcDrwfOALqT9weBt+f17cBZPbzWu4FLgVnAIcCZwB8i4mzgYfKnjIi4rM7m3wM6SYnqLOAfJZ1aKD8TWAQMBxYD/1J43W9K+mYPzbAMOJ7Uy70e+L6k/XLZZ9neMz4deF9hv/sAPwF+A4wi9aIvzL1qJL1H0qYelrF5V5OBtfl8xuP5zfHPG7Th4cCfASsbHMtrcjzdfgMcLunlgOrtkvwmb3uIiPDipdcFWAs8A2wiJfJvAvvnstuAzxfqHg78V3d5XjcT+EV+fCtwXqHsDCCAoYX9fSA/Xgpc0ENMpxWeH9m9H2AM8AJwcKH8S8C38+NLgZ8XyiYAz+1G+zwJHJcfrwHeWij7ANCZH58APFyz7SeBq/v4ej8DngemAsOAv8+vO6ym3r7Az4Fv9bCv3wNTaraJ3J775v1+PD8+A9gKLG3236SX8svQXt8JzLZ7Z0T8vEHZusLjI0hJYaP0pw7hPoU6I2vqP9TDa44hJaK+Ggk8ERFP17xOcXjmkcLjLcB+koZGxLbedi7pYlICH0lKiocAhxVeu3h8tW0zUtKmwrohwO29vWaN54A7IuLmHM9XSedMjiX3zvOnh++QEvOcHvb1TI6/W/fjpyPieUnvBL4BfAJYDtxAeiO3PYQTvfWX4jSo60iJ4LAGSXMjKYF3G1unTnFfryrxmrU2AIdKOriQ7McC63vYppQ8Hv8J0rDLyoh4UdKTbB/m2Eg6OboqPy8e6zrgwYgY32Df7wW+1cPLT4iIh4EVwIk9xCjgKtKnq2kR8XwP+1wJHEdK4OTHj0bEHwAijf2fXNj3r0hDc7aH8Bi99buI2EgaWviapEMk7SPpVZK6k8UNwEckjZb0MtIVH41cCXxM0hvyFT2vlnRELnuUNA5eL4Z1wK+AL0naT9LrgL8DruuHQzyYdP6hCxgq6X+yY4/4BuCTkl4maRQ79qb/A3hK0ifySdshkl4r6Y057usinXNotDyc9/NdYLKk0/JJ7guBx4EHcvk8Uu/+HRHxXC/Hcy3wd5Im5N/Hp4FvdxdKel1uwwMkfQwYUSy31udEbwNlFmnseBVp/PpGUoIA+FfS2PtvgF8DP2y0k4j4PvBF0gnPp4H/TToBCmnM/dP5JOXH6mw+kzTOvAH4EfDZiLilTPCS5kua36B4Kemqo9+RhoP+yI7DM58nnQR+kDQ+fiN5qCPSNf7vIJ3IfZCUnK8kXQZZWkSsJp3knU9q3+nAmRGxNb8Rfii/xiP5qqRn8qeF4o1lY/O+fgpcBvwiH89DpBPK3c4mfUp5jPQp5vSI8NDNHkQR/uIRs4Ek6b8DMyLi5F4rmw0A9+jN+pmkEZJOzENWRwMXkz5RmDWFT8aa9b9hpBOq40iXoy4iXY5q1hQeujEzqzgP3ZiZVZwT/V5CNTM9mtnew4neBoWkc/LkX09J6pR0WffkYyW2CxUmOcvr607RWyifkScUe1bS77snHcvXii/PE449KennkiYUtntJvrTyUaXJ0H6Sr4WvjevkHNcXata3Sbo+X/L5pKSdrtuXdKjShGZ31Kx/h6T786WPv6qJq8f2K1xC2b28IOkbhfIDlObveVxpuuFfFspurtl2q6T7CuW/yPE+Jek3kqYXyv5K0n35eP8g6UcN2qvuMdvgcKK3wXIA6aaew0jzvZwK1Lv2/U/yzTufpGYyLvU8RS+STge+ApxLurnpJNJ8LZCuqT+LdC3+YaTJzBYVdn8B8CbgdaSpDDaRbv8vvv6+wOXAXXXC/iFpaoUjgFcAX61T5ytsv7Gpe5/jSTdznUeaZO0nwOJCMu+x/Yo3VZHuhn0O+H7hJRbkYz42/7yosO3Umu1/VbPtBcCIiDiE7dMUd98TsYo0r8/w3F7/SbpZq9djtsHjRL+XyT3Wr0vakJevK3/JhKTDci95U+7N3q48t7zSnZzrcw96tXacBbJXETEvIm6PiK2Rpjm+jh5u4c++BMwl3VRU9H4aTNGbfY40ydqdkabpXZ9fk4jYFGn63iBNWfACaYbNbuNIE3Y9GhF/JL0JvKbm9S8m3fn72+JKSWeQpjv4+4jYHBHPR8Q9NXXeRJr58eqafb4VuD0i7sjTRnyFNLvlyTnuvrTfWaSbm27Pr3k0aabO2RHRFWm64bvrbag0Z/2bSXPkkF97RWEqiyDNYzQmlz0aERsKu6htz56O2QaJE/3e53+Qprg9njSnySS2f4HIxaQ7OttIvcJPAZETxRzgjRFxMCkprYU+Tatb6yQaT5uLpEmkCcjq3Z3acIpepekA2oE2pTnoOyX9i6T9a/a/iXRH6zeAfywUXQWcqDRn/AHAe0l3wXZvdwTw30h3v9aLazVwTR7GWKbt0z50x3YFqS1rL3cTO04J3P280XTAPbXfOcC1sf2SuhNId7t+LrfXfZL+psG2s0hvOA/uEFzqAPyR9CnmNtLkZt1lY3N7Pkf6lHFZoaynY7ZB4kS/93kvqbf7WER0kXq/3V8i8TxpmoIjcm/09pwsXgBeAkyQtG/uEf8eICKuj4jhPSwP1wYg6VxSMq43rNGdHL4JfDgiXqxTZTQwg9TbHwn8G/DjPKRzOKnHeRapZ3o8ad77Hb4NKw81vJSUgIq97t+R5rlfDzxFGuooJvW5wGci4pkGcZ1BmkrglcDXclzds1p+BLirQW/6FuBkSW/Jx/Ep0vX4B9RW7Kn98hvryew46dho0hvGZlJ7zSG9GR1bJ45Z1JnHJiLeThoGm0b6xPNioezh3J6Hkdq5+Emnp2O2wdIfcx17af2FPHc7qdf1msL6Y4Ct+fHBpOS0Ji+XFOq9B7iDNK/KImDkLsbxTtJkZH/eQ50PAwsLz28jz0+fn/+YPLd9fi5SEjsOeBmp53hOofxvgHsavNY+wB+AV+Tn15HuYj2U9Ob2GVKigjRHza2Fbb8NfKHw/HLSzJTF/d9HmodmJGlum0Pz+veTphku1j0LuD/Hc3l+fHZf2o+UaP9vzbqLSFMVDy2s+wk18/wDf0masvigXn6HPyXNq1Ov7JU5vqFljtnL4Czu0e99NpBOFHYbm9cREU9HxMURcRQpqX20eyw+Us/9L/O2QRpDRtJ761zxUVz+NHQjaQppQrN3RMR9NHYq8C6lL6t+BPgL0kyY3d8AtYIGwwCRvgqvs1F5HfuQes3dV4ocR/pykiciTdz1DWBS7pWfCrQX4vpb0rdD/bi3uEhDZCOAVXnby/N+H8mfYIiIGyPitRHxctKkYkeQvskKKN1+s9h5CuEV9SrWcQ7ww6j/aaVoKI2njh5KOgl9CCWO2QZJs99pvAzOwvYe/RdIV1W0kT5q30HulZK+3u/VpB7yGNKMhW8BjgZOIfVwhwELyd/U1IfXP4XUUz2pRN3hpJ5h9/Ir4KPAS3P50aQvCjmN9KUdF5G+nGRYLv88KUG+gtTDvx34h1x2OmkoZwgpGc0lvdHtl8uvBn5AGtbZlzSEsj6XHVwT1/8C/pntPdZDSZ94zsn7Pwt4IrfzS2q2vYA03v3KwnG/IW/Xlvd9fV/aj/SG+CyFb9XK6/cFOkifToaSTuI+DRxTqLM/6QqjU2q2PYb0LVb75/28j/TpYGIu/+v8+9gnx30D8Otc1usxexmk//9mB+BlkH7R2xP9fjm5bczL3EKSuyjXe5bUK/5MXv860jzqT+fEdRN9HLohjVtvIw0NdC83F8pvBj7VYNvbKAzd5HV/nZPXU7m8OBy1L2mMfxPpUsfiMb6bNIb8DGk++SXA6wrbvpw0fPNY3v4OYFKDuL5NYegmr3szabjmGdIJyzc32Pb97Dx0c0ehjb8FHFi2/XKdbwHfafB6rwH+X/7drgLeVVM+k3TCVjXrj83J+encHsuK25KG2R7M+32ENKx3RNlj9jI4i+e6MTOrOI/Rm5lVXKlEL2lKvkmmQ9JOX/umZG4uXyFpYl5/tKR7C8tTki7s74MwM7PGeh26yWfHf0c6idVJGqObGRGrCnWmkcbqppFuzrg8Ik6os5/1wAkR8VB/HoSZmTVWpkc/CeiIiDURsZV0smV6TZ3p5DvxIuJOYHhhLoxupwK/d5I3MxtcZb5hahQ7fvFxJ6nX3ludUaSrOrrNAL7X6EUkzSZNmMSBBx74hmOOOaZEaGZmBnD33Xc/HhFt9crKJHrVWVdvno6GdfIt3WeSZiKsKyIWkGbYo729PZYvX96oqpmZ1ZDUcLSkzNBNJ3mmumw0+U7KPtSZSrqJ4tESr2dmZv2oTKJfBoyXNC73zGeQ5vAuWgzMylffTAY2R0Rx2GYmPQzbmJnZwOl16CYitkmaAywl3Z69MCJWSjovl88n3V04jXSn4hbSFz4A6ZttSFfsfKj/wzczs96UGaMnIpaQknlx3fzC4wDOb7DtFtJt5WZm1gS+M9bMrOKc6M3MKs6J3sys4pzozcwqzonezKziSl11Y1YlR17ybw3L1n75bYMYidngcI/ezKzinOjNzCrOid7MrOKc6M3MKs6J3sys4pzozcwqzpdXmtVodPmlL720PZUTvVWSk7XZdh66MTOrOCd6M7OKc6I3M6s4J3ozs4pzojczqzgnejOzinOiNzOrOCd6M7OKK5XoJU2RtFpSh6RL6pRL0txcvkLSxELZcEk3SvqtpAckvak/D8DMzHrWa6KXNAS4ApgKTABmSppQU20qMD4vs4F5hbLLgZ9GxDHAccAD/RC3mZmVVKZHPwnoiIg1EbEVWARMr6kzHbg2kjuB4ZJGSDoEOAm4CiAitkbEpn6M38zMelEm0Y8C1hWed+Z1ZeocBXQBV0u6R9KVkg7cjXjNzKyPyiR61VkXJesMBSYC8yLi9cCzwE5j/ACSZktaLml5V1dXibDMzKyMMrNXdgJjCs9HAxtK1gmgMyLuyutvpEGij4gFwAKA9vb22jcSs5bhmTFtT1OmR78MGC9pnKRhwAxgcU2dxcCsfPXNZGBzRGyMiEeAdZKOzvVOBVb1V/BmZta7Xnv0EbFN0hxgKTAEWBgRKyWdl8vnA0uAaUAHsAU4t7CLDwPX5TeJNTVlZmY2wEp98UhELCEl8+K6+YXHAZzfYNt7gfbdiNHMzHaD74w1M6s4f5Wg7ZF8QtSsPPfozcwqzonezKzinOjNzCrOid7MrOKc6M3MKs6J3sys4pzozcwqzonezKzifMOUWT9qdCMX+GYuax736M3MKs6J3sys4pzozcwqzonezKzinOjNzCrOid7MrOKc6M3MKs7X0ZsNIn9hijWDe/RmZhXnHr1ZC3GP3waCe/RmZhVXKtFLmiJptaQOSZfUKZekubl8haSJhbK1ku6TdK+k5f0ZvJmZ9a7XoRtJQ4ArgNOBTmCZpMURsapQbSowPi8nAPPyz25/FRGP91vUZmZWWpkx+klAR0SsAZC0CJgOFBP9dODaiAjgTknDJY2IiI39HrHtNTxebdY/yiT6UcC6wvNOduytN6ozCtgIBPAzSQF8KyIW7Hq4Znsvv/HZriqT6FVnXfShzokRsUHSK4BbJP02In6504tIs4HZAGPHji0RlpmZlVHmZGwnMKbwfDSwoWydiOj++RjwI9JQ0E4iYkFEtEdEe1tbW7nozcysV2US/TJgvKRxkoYBM4DFNXUWA7Py1TeTgc0RsVHSgZIOBpB0IHAGcH8/xm9mZr3odegmIrZJmgMsBYYACyNipaTzcvl8YAkwDegAtgDn5s0PB34kqfu1ro+In/b7UZiZWUOl7oyNiCWkZF5cN7/wOIDz62y3BjhuN2M0M7Pd4DtjzcwqzonezKzinOjNzCrOid7MrOKc6M3MKs6J3sys4vzFI2YV4blwrBH36M3MKs6J3sys4pzozcwqzonezKzinOjNzCrOid7MrOKc6M3MKs6J3sys4pzozcwqznfGWtP4Tk6zweEevZlZxTnRm5lVnBO9mVnFOdGbmVWcE72ZWcWVSvSSpkhaLalD0iV1yiVpbi5fIWliTfkQSfdIuqm/Ajczs3J6vbxS0hDgCuB0oBNYJmlxRKwqVJsKjM/LCcC8/LPbBcADwCH9FLeZ9ZEvZ917lenRTwI6ImJNRGwFFgHTa+pMB66N5E5guKQRAJJGA28DruzHuM3MrKQyiX4UsK7wvDOvK1vn68DHgRd3MUYzM9sNZRK96qyLMnUkvR14LCLu7vVFpNmSlkta3tXVVSIsMzMro0yi7wTGFJ6PBjaUrHMicKaktaQhn1Mkfbfei0TEgohoj4j2tra2kuGbmVlvyiT6ZcB4SeMkDQNmAItr6iwGZuWrbyYDmyNiY0R8MiJGR8SRebtbI+J9/XkAZmbWs16vuomIbZLmAEuBIcDCiFgp6bxcPh9YAkwDOoAtwLkDF7KZmfVFqdkrI2IJKZkX180vPA7g/F72cRtwW58jNDOz3eI7Y83MKs6J3sys4pzozcwqzonezKzinOjNzCrOid7MrOKc6M3MKq7UdfRmu8LT4pq1BvfozcwqzonezKzinOjNzCrOid7MrOJ8MtbMfOK84tyjNzOrOCd6M7OKc6I3M6s4J3ozs4pzojczqzgnejOzinOiNzOrOCd6M7OKc6I3M6u4Uole0hRJqyV1SLqkTrkkzc3lKyRNzOv3k/Qfkn4jaaWkz/X3AZiZWc96TfSShgBXAFOBCcBMSRNqqk0FxudlNjAvr/8v4JSIOA44HpgiaXI/xW5mZiWU6dFPAjoiYk1EbAUWAdNr6kwHro3kTmC4pBH5+TO5zr55if4K3szMelcm0Y8C1hWed+Z1pepIGiLpXuAx4JaIuGvXwzUzs74qk+hVZ11tr7xhnYh4ISKOB0YDkyS9tu6LSLMlLZe0vKurq0RYZmZWRplE3wmMKTwfDWzoa52I2ATcBkyp9yIRsSAi2iOiva2trURYZmZWRpn56JcB4yWNA9YDM4D31NRZDMyRtAg4AdgcERsltQHPR8QmSfsDpwFf6b/wrdk8j7lZ6+s10UfENklzgKXAEGBhRKyUdF4unw8sAaYBHcAW4Ny8+Qjgmnzlzj7ADRFxU/8fhpmZNVLqG6YiYgkpmRfXzS88DuD8OtutAF6/mzGamdlu8J2xZmYV50RvZlZxTvRmZhXnRG9mVnFO9GZmFVfqqhsz27v5fok9m3v0ZmYV50RvZlZxTvRmZhXnRG9mVnFO9GZmFedEb2ZWcU70ZmYV50RvZlZxTvRmZhXnRG9mVnFO9GZmFedEb2ZWcZ7UzHrkyazM9nzu0ZuZVZwTvZlZxTnRm5lVnBO9mVnFlUr0kqZIWi2pQ9IldcolaW4uXyFpYl4/RtIvJD0gaaWkC/r7AMzMrGe9JnpJQ4ArgKnABGCmpAk11aYC4/MyG5iX128DLo6IY4HJwPl1tjUzswFUpkc/CeiIiDURsRVYBEyvqTMduDaSO4HhkkZExMaI+DVARDwNPACM6sf4zcysF2US/ShgXeF5Jzsn617rSDoSeD1wV70XkTRb0nJJy7u6ukqEZWZmZZRJ9KqzLvpSR9JBwA+ACyPiqXovEhELIqI9Itrb2tpKhGVmZmWUSfSdwJjC89HAhrJ1JO1LSvLXRcQPdz1UMzPbFWWmQFgGjJc0DlgPzADeU1NnMTBH0iLgBGBzRGyUJOAq4IGI+Kd+jNvMWoinymhtvSb6iNgmaQ6wFBgCLIyIlZLOy+XzgSXANKAD2AKcmzc/ETgbuE/SvXndpyJiSf8ehpmZNVJqUrOcmJfUrJtfeBzA+XW2u4P64/dmZjZIfGesmVnFOdGbmVWcE72ZWcU50ZuZVZwTvZlZxTnRm5lVnL8zdi/nG13Mqs89ejOzinOiNzOrOCd6M7OKc6I3M6s4J3ozs4pzojczqzgnejOzinOiNzOrOCd6M7OK852xZjagfPd187lHb2ZWcU70ZmYV50RvZlZxTvRmZhXnk7F7AZ8MM9u7lerRS5oiabWkDkmX1CmXpLm5fIWkiYWyhZIek3R/fwZuZmbl9JroJQ0BrgCmAhOAmZIm1FSbCozPy2xgXqHs28CU/gjWzMz6rkyPfhLQERFrImIrsAiYXlNnOnBtJHcCwyWNAIiIXwJP9GfQZmZWXplEPwpYV3jemdf1tU6PJM2WtFzS8q6urr5samZmPSiT6FVnXexCnR5FxIKIaI+I9ra2tr5samZmPSiT6DuBMYXno4ENu1DHzMyaoEyiXwaMlzRO0jBgBrC4ps5iYFa++mYysDkiNvZzrGZmtgt6vY4+IrZJmgMsBYYACyNipaTzcvl8YAkwDegAtgDndm8v6XvAW4DDJHUCn42Iq/r7QMxsz+T7PAZeqRumImIJKZkX180vPA7g/AbbztydAM3MbPd4CgQzs4rzFAgV4I++ZtYT9+jNzCrOid7MrOKc6M3MKs6J3sys4nwy1sxami822H17VaJv9AcD/qMxs+qqXKLfnXf/Vu05tGpcZlZfb53Kwf6f9hi9mVnFVa5Hb2Y2GPakT9pO9H2wJ/1izfYG/p8sx4m+RfgP1swGisfozcwqzj16M6usKl6Ftyuc6PtJq11OZWY925vuq/HQjZlZxTnRm5lVnBO9mVnFOdGbmVWcE72ZWcU50ZuZVVypRC9piqTVkjokXVKnXJLm5vIVkiaW3dbMzAZWr4le0hDgCmAqMAGYKWlCTbWpwPi8zAbm9WFbMzMbQGV69JOAjohYExFbgUXA9Jo604FrI7kTGC5pRMltzcxsACkieq4gnQVMiYgP5OdnAydExJxCnZuAL0fEHfn5/wE+ARzZ27aFfcwmfRoAOBpYvXuH1u8OAx5vdhB1OK6+a9XYHFfftGpc0JzYjoiItnoFZaZAUJ11te8OjeqU2TatjFgALCgRT1NIWh4R7c2Oo5bj6rtWjc1x9U2rxgWtF1uZRN8JjCk8Hw1sKFlnWIltzcxsAJUZo18GjJc0TtIwYAawuKbOYmBWvvpmMrA5IjaW3NbMzAZQrz36iNgmaQ6wFBgCLIyIlZLOy+XzgSXANKAD2AKc29O2A3IkA69Vh5UcV9+1amyOq29aNS5osdh6PRlrZmZ7Nt8Za2ZWcU70ZmYV50TfgKQhku7J9wgg6VJJ6yXdm5dpTYprraT7cgzL87pDJd0i6T/zz5e1SFxNbzNJwyXdKOm3kh6Q9KYWaa96cbVCex1deP17JT0l6cJmt1kPcbVCm10kaaWk+yV9T9J+zW6vnWL0GH19kj4KtAOHRMTbJV0KPBMRX21yXGuB9oh4vLDuMuCJiPhynk/oZRHxiRaI61Ka3GaSrgFuj4gr85VfBwCfovntVS+uC2mBv7FueQqT9cAJwPk0uc0axHUuTWwzSaOAO4AJEfGcpBtIF6dMoEXaC9yjr0vSaOBtwJXNjqWk6cA1+fE1wDubGEvLkHQIcBJwFUBEbI2ITTS5vXqIq9WcCvw+Ih6itf7GinG1gqHA/pKGkt6wN9Ba7eVE38DXgY8DL9asn6M0O+fCJn4UC+Bnku7O00YAHJ7vWyD/fEWLxAXNbbOjgC7g6jwMd6WkA2l+ezWKC1rjb6zbDOB7+XGz26yoGBc0sc0iYj3wVeBhYCPpHqKf0Vrt5URfS9Lbgcci4u6aonnAq4DjSb/Qrw12bNmJETGRNCPo+ZJOalIcterF1ew2GwpMBOZFxOuBZ4FWmCq7UVzNbq8/ycNJZwLfb1YM9dSJq6ltlt9YpgPjgJHAgZLeN5gxlOFEv7MTgTPzmPMi4BRJ342IRyPihYh4EfhX0sycgy4iNuSfjwE/ynE8qjRbKPnnY60QVwu0WSfQGRF35ec3khJss9urblwt0F5FU4FfR8Sj+Xmz26xuXC3QZqcBD0ZEV0Q8D/wQ+Atap70AJ/qdRMQnI2J0RBxJ+oh4a0S8r/uXlr0LuH+wY5N0oKSDux8DZ+Q4FgPn5GrnAD9uhbia3WYR8QiwTtLRedWpwCqa3F6N4mp2e9WYyY7DI01ts4Id4mqBNnsYmCzpAEki/S4foHXaC/BVNz2S9BbgY/mqm++QPh4GsBb4UPcY3CDGcxSptwzp4//1EfFFSS8HbgDGkv7w3h0RT7RAXK3QZseTTqoPA9aQrtLYhya2Vw9xzaXJ7ZVjOwBYBxwVEZvzuqb+jfUQVyv8jX0O+FtgG3AP8AHgIJrcXjvE6ERvZlZtHroxM6s4J3ozs4pzojczqzgnejOzinOiN+0Ze7wAAAAUSURBVDOrOCd6M7OKc6I3M6u4/w8q17s5SLEleQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Example\n",
    "model = SFCN()\n",
    "model = torch.nn.DataParallel(model)\n",
    "fp_ = './brain_age/run_20190719_00_epoch_best_mae.p'\n",
    "model.load_state_dict(torch.load(fp_))\n",
    "model.cuda()\n",
    "\n",
    "# Example data: some random brain in the MNI152 1mm std space\n",
    "data = np.random.rand(182, 218, 182)\n",
    "label = np.array([71.3,]) # Assuming the random subject is 71.3-year-old.\n",
    "\n",
    "# Transforming the age to soft label (probability distribution)\n",
    "bin_range = [42,82]\n",
    "bin_step = 1\n",
    "sigma = 1\n",
    "y, bc = dpu.num2vect(label, bin_range, bin_step, sigma)\n",
    "y = torch.tensor(y, dtype=torch.float32)\n",
    "print(f'Label shape: {y.shape}')\n",
    "\n",
    "# Preprocessing\n",
    "data = data/data.mean()\n",
    "data = dpu.crop_center(data, (160, 192, 160))\n",
    "\n",
    "# Move the data from numpy to torch tensor on GPU\n",
    "sp = (1,1)+data.shape\n",
    "data = data.reshape(sp)\n",
    "input_data = torch.tensor(data, dtype=torch.float32).cuda()\n",
    "print(f'Input data shape: {input_data.shape}')\n",
    "print(f'dtype: {input_data.dtype}')\n",
    "\n",
    "# Evaluation\n",
    "model.eval() # Don't forget this. BatchNorm will be affected if not in eval mode.\n",
    "with torch.no_grad():\n",
    "    output = model(input_data)\n",
    "\n",
    "# Output, loss, visualisation\n",
    "x = output[0].cpu().reshape([1, -1])\n",
    "print(f'Output shape: {x.shape}')\n",
    "loss = dpl.my_KLDivLoss(x, y).numpy()\n",
    "\n",
    "# Prediction, Visualisation and Summary\n",
    "x = x.numpy().reshape(-1)\n",
    "y = y.numpy().reshape(-1)\n",
    "\n",
    "plt.bar(bc, y)\n",
    "plt.title('Soft label')\n",
    "plt.show()\n",
    "\n",
    "prob = np.exp(x)\n",
    "pred = prob@bc\n",
    "plt.bar(bc, prob)\n",
    "plt.title(f'Prediction: age={pred:.2f}\\nloss={loss}')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CPU version\n",
    "- Get rid of model.cuda()\n",
    "- Get rid of data.cuda()\n",
    "- Deal with [DataParallel issue with CPU](https://discuss.pytorch.org/t/missing-keys-unexpected-keys-in-state-dict-when-loading-self-trained-model/22379)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Label shape: torch.Size([1, 40])\n",
      "Input data shape: torch.Size([1, 1, 160, 192, 160])\n",
      "dtype: torch.float32\n",
      "Output shape: torch.Size([1, 40])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAUeUlEQVR4nO3df5Bd5X3f8ffHIkprAcYxMrElsDBWQ9XW2MwW3EJxqAODTBqZSTMW45jYCaPSQfWQxpMqnUzqjJsWZ+zW45ZYVQgpdgYTx2NN1aAYMukfboLxaKkpIAyOLAtrEUYL+EdwXIPMt3/co/ay3NUeoV3dq8fv18zOvef5cc73Prvz2bNn749UFZKkdr1s3AVIkpaWQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXj/UklyVZH+SZ5K8ucf4n0wy03Pf70ny5y+xrpc8V5rLoNcJL8nFSe5O8u0kTyf5iyR/v+f0DwObq+rkqvpSkn1Jfmop65WOt5PGXYB0LJKcCvwx8M+BTwPLgX8EfL/nLl4H7F6a6qTJ4Bm9TnR/C6CqPlVVP6iq71XVXVV1P0CSlyX59SSPJjmY5BNJXpHkR5M8AywD/neSryb5JHAW8N+7Szm/utDBk2zp5v5VkoeSXPXiIflP3V8bDyd521DHK5L8XpLHkzyW5N8mWbZoKyN1DHqd6L4C/CDJrUnWJ3nlnP73dF+XAq8HTgb+c1V9v6pO7sacV1XnVNW7ga8D/6S7lPPbPY7/VQZ/QbwC+E3gD5K8Zqj/QmAvcDrwb4DPJvmxru9W4BDwBuDNwOXAtf0futSPQa8TWlV9B7gYKOB3gdkkO5Kc0Q15F/AfqmpvVT0D/BqwMcmiXLasqj+qqgNV9XxV/SHwl8AFQ0MOAh+tque6/keAK7v61gM3VNV3q+og8B+BjYtRlzTMoNcJr6q+XFXvqarVwN8FXgt8tOt+LfDo0PBHGfxv6gwWQZJrktyX5FtJvtUd//ShIY/VC9858NGuptcBPwI8PjT3vwCvXoy6pGEGvZpSVQ8D/5VB4AIcYBCqh53F4HLJE/Ptou+xkryOwV8Rm4FXVdVpwINAhoatSjK8fVZX034G/zA+vapO675Oraq/0/f4Ul8GvU5oSc5N8itJVnfbZwJXA/d0Qz4F/HKSs5OcDPw74A+r6tA8u3yCwbX8PlYw+MUw2x37vfz/XzCHvRp4X5IfSfJzwN8GdlbV48BdwEeSnNr90/icJG/teWypN4NeJ7q/YvAPzy8m+S6DgH8Q+JWu/xbgk8Dnga8B/wf4F0fY378Hfr27nPL+Ix24qh4CPgJ8gcEviL8H/MWcYV8E1gJPAr8F/NOqeqrru4bB00EfAr4JfAZ4DdIiix88Iklt84xekhpn0EtS4wx6SWqcQS9JjZvINzU7/fTTa82aNeMuQ5JOGPfee++TVbVyVN9EBv2aNWuYnp4edxmSdMJI8uh8fV66kaTGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxk3kK2MlnVjWbLljZPu+G688zpVoFM/oJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDWuV9AnuSLJI0n2JNkyon9DkvuT3JdkOsnFQ337kjxwuG8xi5ckLWzB97pJsgy4CbgMmAF2JdlRVQ8NDfszYEdVVZI3Ap8Gzh3qv7SqnlzEuiVJPfU5o78A2FNVe6vqWeB2YMPwgKp6pqqq21wBFJKkidAn6FcB+4e2Z7q2F0hyVZKHgTuAXxzqKuCuJPcm2TTfQZJs6i77TM/OzvarXpK0oD5BnxFtLzpjr6rtVXUu8A7gg0NdF1XV+cB64Pokl4w6SFVtq6qpqppauXJlj7IkSX30CfoZ4Myh7dXAgfkGV9XngXOSnN5tH+huDwLbGVwKkiQdJ32CfhewNsnZSZYDG4EdwwOSvCFJuvvnA8uBp5KsSHJK174CuBx4cDEfgCTpyBZ81k1VHUqyGbgTWAbcUlW7k1zX9W8Ffha4JslzwPeAd3bPwDkD2N79DjgJuK2qPrdEj0WSNEKvjxKsqp3AzjltW4fufwj40Ih5e4HzjrFGSdIx8JWxktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIa1yvok1yR5JEke5JsGdG/Icn9Se5LMp3k4r5zJUlLa8GgT7IMuAlYD6wDrk6ybs6wPwPOq6o3Ab8I3HwUcyVJS6jPGf0FwJ6q2ltVzwK3AxuGB1TVM1VV3eYKoPrOlSQtrT5BvwrYP7Q907W9QJKrkjwM3MHgrL73XEnS0ukT9BnRVi9qqNpeVecC7wA+eDRzAZJs6q7vT8/OzvYoS5LUR5+gnwHOHNpeDRyYb3BVfR44J8npRzO3qrZV1VRVTa1cubJHWZKkPvoE/S5gbZKzkywHNgI7hgckeUOSdPfPB5YDT/WZK0laWictNKCqDiXZDNwJLANuqardSa7r+rcCPwtck+Q54HvAO7t/zo6cu0SPRZI0woJBD1BVO4Gdc9q2Dt3/EPChvnMlScePr4yVpMYZ9JLUOINekhrX6xq9pB9ua7bcMbJ9341XHudK9FJ4Ri9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjegV9kiuSPJJkT5ItI/rfleT+7uvuJOcN9e1L8kCS+5JML2bxkqSFLfjBI0mWATcBlwEzwK4kO6rqoaFhXwPeWlXfTLIe2AZcONR/aVU9uYh1S5J66nNGfwGwp6r2VtWzwO3AhuEBVXV3VX2z27wHWL24ZUqSXqo+Qb8K2D+0PdO1zeeXgD8Z2i7griT3Jtk036Qkm5JMJ5menZ3tUZYkqY8+nxmbEW01cmByKYOgv3io+aKqOpDk1cCfJnm4qj7/oh1WbWNwyYepqamR+5ckHb0+Z/QzwJlD26uBA3MHJXkjcDOwoaqeOtxeVQe624PAdgaXgiRJx0mfoN8FrE1ydpLlwEZgx/CAJGcBnwXeXVVfGWpfkeSUw/eBy4EHF6t4SdLCFrx0U1WHkmwG7gSWAbdU1e4k13X9W4HfAF4F/E4SgENVNQWcAWzv2k4Cbquqzy3JI5EkjdTnGj1VtRPYOadt69D9a4FrR8zbC5w3t12SdPz4ylhJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY3rFfRJrkjySJI9SbaM6H9Xkvu7r7uTnNd3riRpaS0Y9EmWATcB64F1wNVJ1s0Z9jXgrVX1RuCDwLajmCtJWkJ9zugvAPZU1d6qeha4HdgwPKCq7q6qb3ab9wCr+86VJC2tPkG/Ctg/tD3Ttc3nl4A/eYlzJUmL7KQeYzKirUYOTC5lEPQXv4S5m4BNAGeddVaPsiRJffQJ+hngzKHt1cCBuYOSvBG4GVhfVU8dzVyAqtpGd21/ampq5C8DSSeeNVvuGNm+78Yrj3MlP7z6XLrZBaxNcnaS5cBGYMfwgCRnAZ8F3l1VXzmauZKkpbXgGX1VHUqyGbgTWAbcUlW7k1zX9W8FfgN4FfA7SQAOVdXUfHOX6LFIkkboc+mGqtoJ7JzTtnXo/rXAtX3nSpKOH18ZK0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDWuV9AnuSLJI0n2JNkyov/cJF9I8v0k75/Tty/JA0nuSzK9WIVLkvpZ8MPBkywDbgIuA2aAXUl2VNVDQ8OeBt4HvGOe3VxaVU8ea7GSpKPX54z+AmBPVe2tqmeB24ENwwOq6mBV7QKeW4IaJUnHoE/QrwL2D23PdG19FXBXknuTbJpvUJJNSaaTTM/Ozh7F7iVJR9In6DOirY7iGBdV1fnAeuD6JJeMGlRV26pqqqqmVq5ceRS7lyQdSZ+gnwHOHNpeDRzoe4CqOtDdHgS2M7gUJEk6TvoE/S5gbZKzkywHNgI7+uw8yYokpxy+D1wOPPhSi5UkHb0Fn3VTVYeSbAbuBJYBt1TV7iTXdf1bk/w4MA2cCjyf5AZgHXA6sD3J4WPdVlWfW5qHIkkaZcGgB6iqncDOOW1bh+5/g8Elnbm+A5x3LAVKko6Nr4yVpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGtcr6JNckeSRJHuSbBnRf26SLyT5fpL3H81cSdLSWjDokywDbgLWA+uAq5OsmzPsaeB9wIdfwlxJ0hLqc0Z/AbCnqvZW1bPA7cCG4QFVdbCqdgHPHe1cSdLS6hP0q4D9Q9szXVsfvecm2ZRkOsn07Oxsz91LkhbSJ+gzoq167r/33KraVlVTVTW1cuXKnruXJC2kT9DPAGcOba8GDvTc/7HMlSQtgj5BvwtYm+TsJMuBjcCOnvs/lrmSpEVw0kIDqupQks3AncAy4Jaq2p3kuq5/a5IfB6aBU4Hnk9wArKuq74yau1QPRpL0YgsGPUBV7QR2zmnbOnT/Gwwuy/SaK0k6fnxlrCQ1rtcZvaS2rdlyx8j2fTdeeZwr0VLwjF6SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1LheQZ/kiiSPJNmTZMuI/iT5WNd/f5Lzh/r2JXkgyX1JphezeEnSwhb8hKkky4CbgMuAGWBXkh1V9dDQsPXA2u7rQuDj3e1hl1bVk4tWtSSptz5n9BcAe6pqb1U9C9wObJgzZgPwiRq4BzgtyWsWuVZJ0kvQJ+hXAfuHtme6tr5jCrgryb1JNs13kCSbkkwnmZ6dne1RliSpjz5BnxFtdRRjLqqq8xlc3rk+ySWjDlJV26pqqqqmVq5c2aMsSVIffYJ+BjhzaHs1cKDvmKo6fHsQ2M7gUpAk6TjpE/S7gLVJzk6yHNgI7JgzZgdwTffsm7cA366qx5OsSHIKQJIVwOXAg4tYvyRpAQs+66aqDiXZDNwJLANuqardSa7r+rcCO4G3A3uAvwbe200/A9ie5PCxbquqzy36o5AkzWvBoAeoqp0Mwny4bevQ/QKuHzFvL3DeMdYoqWFrttwxsn3fjVce50ra5StjJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN6/XBI5JOfH7Axw8vz+glqXEGvSQ1rlfQJ7kiySNJ9iTZMqI/ST7W9d+f5Py+cyVJS2vBa/RJlgE3AZcBM8CuJDuq6qGhYeuBtd3XhcDHgQt7zpWkefm/hWPX54z+AmBPVe2tqmeB24ENc8ZsAD5RA/cApyV5Tc+5kqQl1OdZN6uA/UPbMwzO2hcas6rnXACSbAI2dZvPJHmkR23H0+nAk+MuYgTrOnqTWttY6sqHFuybt64jze2572Ppn9TvI4zne/m6+Tr6BH1GtFXPMX3mDhqrtgHbetQzFkmmq2pq3HXMZV1Hb1Jrs66jM6l1weTV1ifoZ4Azh7ZXAwd6jlneY64kaQn1uUa/C1ib5Owky4GNwI45Y3YA13TPvnkL8O2qerznXEnSElrwjL6qDiXZDNwJLANuqardSa7r+rcCO4G3A3uAvwbee6S5S/JIlt6kXlayrqM3qbVZ19GZ1LpgwmpL1chL5pKkRvjKWElqnEEvSY0z6OeRZFmSLyX54277A0keS3Jf9/X2MdW1L8kDXQ3TXduPJfnTJH/Z3b5yQuoa+5olOS3JZ5I8nOTLSf7BhKzXqLomYb1+Yuj49yX5TpIbxr1mR6hrEtbsl5PsTvJgkk8l+RvjXq8X1eg1+tGS/EtgCji1qn46yQeAZ6rqw2Ouax8wVVVPDrX9NvB0Vd3YvZ/QK6vqX01AXR9gzGuW5Fbgf1bVzd0zv14O/GvGv16j6rqBCfgZO6x7C5PHGLzI8XrGvGbz1PVexrhmSVYBfw6sq6rvJfk0gyenrGNC1gs8ox8pyWrgSuDmcdfS0wbg1u7+rcA7xljLxEhyKnAJ8HsAVfVsVX2LMa/XEeqaNG8DvlpVjzJZP2PDdU2Ck4C/meQkBr+wDzBZ62XQz+OjwK8Cz89p35zBu3PeMsY/xQq4K8m93dtGAJzRvW6B7vbVE1IXjHfNXg/MAr/fXYa7OckKxr9e89UFk/EzdthG4FPd/XGv2bDhumCMa1ZVjwEfBr4OPM7gNUR3MVnrZdDPleSngYNVde+cro8D5wBvYvAN/cjxrq1zUVWdz+AdQ69PcsmY6phrVF3jXrOTgPOBj1fVm4HvApPwVtnz1TXu9fp/ustJPwP80bhqGGVEXWNds+4XywbgbOC1wIokP388a+jDoH+xi4Cf6a453w784yR/UFVPVNUPqup54HcZvDPncVdVB7rbg8D2ro4nMni3ULrbg5NQ1wSs2QwwU1Vf7LY/wyBgx71eI+uagPUath74X1X1RLc97jUbWdcErNlPAV+rqtmqeg74LPAPmZz1Agz6F6mqX6uq1VW1hsGfiP+jqn7+8DetcxXw4PGuLcmKJKccvg9c3tWxA/iFbtgvAP9tEuoa95pV1TeA/Ul+omt6G/AQY16v+eoa93rNcTUvvDwy1jUb8oK6JmDNvg68JcnLk4TB9/LLTM56AT7r5oiS/CTw/u5ZN59k8OdhAfuAf3b4GtxxrOf1DM6WYfDn/21V9VtJXgV8GjiLwQ/ez1XV0xNQ1ySs2ZsY/FN9ObCXwbM0XsYY1+sIdX2MMa9XV9vLGby9+Our6ttd21h/xo5Q1yT8jP0m8E7gEPAl4FrgZMa8Xi+o0aCXpLZ56UaSGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMb9XyrNTNq12lbaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEXCAYAAACjyo8UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAfD0lEQVR4nO3df5wdVX3/8debjZFfwaCsml8Y0EiNfgXTGGKx6FeoJohEi9pEBYpipBIVxWr80Yq//VLar1IhKUKACJIiapvyDaJWqFDFkoBEY4xdQyBLAiwiIRAkBD/fP865ZTLcuzub7ObeTN7Px2Mee++cMzOfObv7uXPPzJxRRGBmZvW1V7sDMDOz4eVEb2ZWc070ZmY150RvZlZzTvRmZjXnRG9mVnNO9NZ2kiZKCkkj8vtrJZ2yA+s5WNLDkrqGPkqz3ZcTvVUiaZ2kR3MivVfSJZL2H45tRcTMiLisYkzHFpa7KyL2j4gnhiOuTiPpUEnXSNos6X5J5+T5T5d0saQ7c9ltkmYOsK4PSrpH0iZJiyQ9vVA2T9JySY9JunSYd8uGgRO9DcYbImJ/YArwcuCT5QpK/Hc1zCSNBL4P/BB4LjAeuDwXjwDWA68CngH8DXCVpIkt1vU6YD5wDDAROBT4dKHKBuBzwKKh3QvbVfwPaYMWEXcD1wIvAZB0g6TPS/pPYAtwqKRn5KPKjZLulvS5RpeKpC5J5+aj0LXA64vrz+s7rfD+3ZJW56PTX0qaIunrwMHAv+VvGR9p0gU0VtJSSQ9I6pH07sI6z5Z0laTFeb2rJE2t2gaSviJpvaSHJK2Q9KeFsn0kXSbpdznuj0jqLZSPlfQtSX2S7pD0/kH9ApK/BDZExD9ExCMR8fuIWAmQ358dEesi4g8RcQ1wB/DHLdZ1CnBxRKyKiN8Bn83rJ6/v2xHxL8BvdyBO6wBO9DZokiYAxwG3FWafBMwFRgF3ApcB24AXAC8DXgs0kve7gePz/KnAm/vZ1luAs4GTgQOAE4DfRsRJwF3kbxkRcU6Txa8EeoGxeRtfkHRMofwEYAkwGlgKfLWw3QskXdBPM9wCHAE8E/gG8E1Je+eyT/HkkfGfAe8orHcv4N+A24FxpKPoM/NRNZLeJunBfqaD86qmA+vy+Yz784fj/2rRhs8BXgisarEvL87xNNwOPEfSs/rZf9udRIQnTwNOwDrgYeBBUiK/ANgnl90AfKZQ9znAY43yPG8OcH1+/UPg9ELZa4EARhTWd1p+fR3wgX5iOrbwfmJjPcAE4AlgVKH8i8Cl+fXZwA8KZZOBR3eifX4HHJ5frwVeVyg7DejNr48E7iot+zHgkkFu73vA48BMYCTw13m7I0v1ngb8APinftb1G2BGaZkAJpbqfa7Rfp52r2lEtY8DMwDeGBE/aFG2vvD6eaRksVFSY95ehTpjS/Xv7GebE0iJaLDGAg9ExObSdordM/cUXm8B9pY0IiK2DbRySWeREvhYUlI8ADiosO3i/pXbZqykBwvzuoAbB9pmyaPATRFxbY7nXNI5kxeRj87zt4evA1uBef2s6+Ecf0Pj9eYmdW035K4bGyrFYVDXk47oD4qI0Xk6ICJenMs3khJ4w8G0th54foVtlm0AnilpVGk7d/ezTCW5P/6jwFuBAyNiNLAJaHyqbSSdHG0o7ut64I5Cu4yOiFERcVxe99vzOYdWU6OtVtLP/it9wl5M+nZ1YkQ83s8urQIOL7w/HLg3ItwnXxNO9DbkImIjqWvh7yUdIGkvSc+X9Kpc5Srg/ZLGSzqQdMVHKxcBH5b0x/mKnhdIel4uu5fUD94shvXAj4EvStpb0kuBdwFXDMEujiKdf+gDRkj6W7Y/Ir4K+JikAyWNY/uj6f8CHpL00XzStkvSSyS9PMd9RaRzDq2mu/J6LgemSzo2n+Q+E7gfWJ3LF5CO7t8QEY8OsD+LgXdJmpx/H58ELm0UShqRzz90AV25Pd0bsBtxorfhcjKp7/iXpP7rq4ExuexrpL7324FbgW+3WklEfBP4POmE52bgX0gnQCH1uX8yn6T8cJPF55D67TcA3wE+FRHfrxK8pIWSFrYovo501dGvSd1Bv2f77pnPkE4C30HqH7+a9A2HSNf4v4F0IvcOUnK+iHQZZGURsYZ0knchqX1nASdExNb8QfievI17Ct8G3p73rXFj2cF5Xd8FzgGuz/tzJ+mEcsMnSV1F8/M2H6XJpbXWuRThB4+YDSdJfwXMjohXDVjZbBj4iN5siEkaI+mo3GV1GHAW6RuFWVu4n81s6I0E/gk4hHQ56hLS5ahmbeGuGzOzmnPXjZlZzTnR7yFUGunRzPYcTvS2S0g6JQ/+9ZCkXknnVLkWOy8X2n6Qs4Wlm4gek7S5UH650mBqD0n6dXHZ0ro/ldd9bGHe6Dwg2X15Oru0zERJ10vaIulXpWUl6ROS7srbXiLpgEL5pZK2lmJvDPR2kKT/lPTbfLnoTyQd1SLuH6oweFupbJKk30u6vDCvMdhbcbt/02TZkXmfekvzr1cagO0hSbdLmtUirkvydl5Qmn+spFslPaI0ENxbmy1vw8eJ3naVfUk39RxEGu/lGKDZte//I9+88zFKg3FFxOnFm4hIg5d9s1Dli6RxWhqDoH1O0nYjN0p6Pmmgs42lzf7fHOtEYBpwkqRTC+VXkgZzexbwCeBqSd257GTS4G5HkYZB2Af4x9L6zyndANUYO/9h4J1AN3Ag8H9II3Nul8zztfD9fUCeTxpwrZnRhe1+tkn5XwP3NZn/AWBMbs+5wOWSxhQrSHolTe5gljSZdA/EJ0j3ChwBrOgnfhsGTvR7GKWHUnxZ0oY8fVn5IRP5qPKafET5gKQblceWV7qT826lIX3XaPtRIAcUEQsi4saI2BppmOMrSAmxP18EziPdVNRqf/YDTiSNltnY1qqIeKzxNk/lJPRV0jAGW0vz30BKxlsiYh1pGIF35m29kDQW/6ci4tGI+Bbw87z9xrIXR8T6iHiYlKz/QtK+A+wnkYYZXhMRfyANpfAEKeE3bg5D0jNINzJ9pEVbzCZd5fPvA22vybKHkG6G+mKT2FYWxv8J0jhGEwrLjiB9oDUbT+eTpAHVro2IbRHx24jYkbGLbCc40e95PkEa4vYI0pgm03jyLsezSHd0dpPGSPk4EPla8HnAyyNiFPA60siRgxlWt+xoWg+bi6RppAHIWt2d2nAiaSiCH5WWv0DSFuBXpKP2ZYWytwBbI2IZzan0+iX59YuBtaWB0m7P8xt1y8s+HZhUmPfe/CG6QtKJlEhaSbrTdilwUUQUj7C/QBra4J4myx1AuiP3rBb7BHCnUrfZJZIOKpX9I+n33XS4hHwA8Hvgp6TRRZcXij8I/CjyePgl0/PyP8/daZdLemaTejaMnOj3PG8nDSl8X0T0kZ4kdFIue5w0TMHzIuLxfAQepKPLpwOTJT0t0gMtfgMQEd8oDdBVnu4qB5C7QqYC5zYLMPdbXwC8Lx/h9ucUYHGUrhOOiPeSxqT5U9IQC4/lde9PSphntljfd4H5kkblvuZ3krpyAPYnDV5WtClvB9KwCKflPvFnkL4xUFj+PFLSfzbpqU+XlvvhI+KlpHFz3gbc1Jiv9FCUo3hqV1DDZ8nfJpqU3U96ItjzSA8fGUVhzB9JbyINEd3ypq6IOD4vdxxwXeP3ovRsgvcAf9ti0fGkv68T8743686yYeZEv+cZy/bDAt+Z5wH8HdADfE/SWknzASKih5QYzwbuyycZx7IDJL0R+BIwMyJadcm8F1gZET8ZYF0TSI/LW9ysPCKeiIibSMnmr/LsTwNfj4g7Wqz2/aSj2v8G/pUnH14CTx3Ol/y+cYS/KNe/gfRt5fo8vzfHc2vuutiWv01cAfx5k7h/HxFXkj5wDs/dZxeQxuV/yhDKko4AjiWdX2jWDg9HxPK83XtJ385eqzTg3H6kcW7e16I9iut5PA+L/DpJJ+TZXyYdOJQ/ABseJY21/+vcnfUF0oeF7UJO9HueDaQju4aD8zwiYnNEnBURh5L6mz/U6IvPR+6vzMsGqf95MMPqImkGaUCzN0TEz/uJ8RjgTUoPq74H+BPSSJhfLdU7GfhxRKwdYJ9H8GQf/TGkkTMb655Aep7qR/N+PhARb4+I50YaVnkv0oiTkJL3odp+6OPD83wiPbbvUxExMSLG5/l303po5GD7rp6yp5FG5zyA9A3on3PMjZOtvUpDJr+adPL4rlz+YeBESbf2s13ytiflZW/My34bGJPbZ2KL5cvt+XeF9gT4iaS35df9Dqdsu0h0wNNPPA3/RH4aE+kpQT8m9cMfROoe+Fyuczzp0X8iJcCNpCRyGPAaUvfNSNKR66WD3P5rSM8cPbpC3dGkB143ph8DHwKeUaq3Bnhnad6zgdmkbpYu0vmER4BZufxZpXWvB94C7J/Ln5/rdJGe3nQ/8OLC+m8mdTntDbyJdPKzO5c9My8v0hOrfgHMLSz75hzXXqSnam0GXp3LpgOvzO27D6nbZzPp25ZKMb+clDzH5fr7lsrPJY2Y2YjryPw73Cvv2z/z5NO+RpSW/XPSB/9zcxv8UW6HfUgfPO8gncCeUmjv4vKR96Xx9LF3kkbpPDTHeRXpG1Xb/yf2pKntAXjaRb/oJxP93qS+4o15Og/YO9f5YK73CKm74W/y/JeSjmo3Aw8A1wBjB7n960ljuD9cmK4tlF8LfLzFsjeQHy1YmPeKHOeo0vxu4D9yAn6IdFXMuwdql8L7t+ZEtwX4GYVHAubyiTmeR0kfNMVlX5jnbSF1iX2otOyNpD79h0gncWcXyl6V5zXa+D9o8aFI4ZGJLcrPBi4vvJ+Tk+0j+Xe+GHhui2VfTX7sYX7/ItIJ2M25TW8B3tRPewbwgtK8T5NOmPeRnnh1YLv/H/a0yWPdmJnVnPvozcxqzonezKzmnOjNzGrOid7MrOY68glTBx10UEycOLHdYZiZ7TZWrFhxf0R0NyvryEQ/ceJEli9fPnBFMzMDQNKdrcrcdWNmVnNO9GZmNedEb2ZWc070ZmY150RvZlZzTvRmZjXnRG9mVnNO9GZmNedEb2ZWcx15Z6zZcJo4//+1LFv3pdfvwkjMdg0f0ZuZ1ZwTvZlZzTnRm5nVnBO9mVnNOdGbmdWcE72ZWc358kqrpVaXUFa5fHJnljXrRD6iNzOrOSd6M7Oac6I3M6s5J3ozs5pzojczqzknejOzmnOiNzOrOSd6M7Oac6I3M6u5Sole0gxJayT1SJrfpFySzsvlKyVNKZSNlnS1pF9JWi3pFUO5A2Zm1r8BE72kLuB8YCYwGZgjaXKp2kxgUp7mAgsKZV8BvhsRfwQcDqwegrjNzKyiKkf004CeiFgbEVuBJcCsUp1ZwOJIbgZGSxoj6QDgaOBigIjYGhEPDmH8ZmY2gCqJfhywvvC+N8+rUudQoA+4RNJtki6StF+zjUiaK2m5pOV9fX2Vd8DMzPpXJdGrybyoWGcEMAVYEBEvAx4BntLHDxARF0bE1IiY2t3dXSEsMzOrokqi7wUmFN6PBzZUrNML9EbET/P8q0mJ38zMdpEqif4WYJKkQySNBGYDS0t1lgIn56tvpgObImJjRNwDrJd0WK53DPDLoQrezMwGNuCDRyJim6R5wHVAF7AoIlZJOj2XLwSWAccBPcAW4NTCKt4HXJE/JNaWysx2O34wie1uKj1hKiKWkZJ5cd7CwusAzmix7M+AqTsRo5mZ7QQ/StB2Sz6qNqvOQyCYmdWcE72ZWc050ZuZ1ZwTvZlZzTnRm5nVnBO9mVnNOdGbmdWcE72ZWc050ZuZ1ZwTvZlZzTnRm5nVnMe6MRtCrcbgAY/DY+3jI3ozs5pzojczqzknejOzmnOiNzOrOSd6M7Oac6I3M6s5J3ozs5rzdfRmu5CfdWvt4CN6M7Oaq3REL2kG8BWgC7goIr5UKlcuPw7YAvxlRNyay9YBm4EngG0RMXXIorda89Gv2dAYMNFL6gLOB/4M6AVukbQ0In5ZqDYTmJSnI4EF+WfD/46I+4csarOa8oebDYcqXTfTgJ6IWBsRW4ElwKxSnVnA4khuBkZLGjPEsZqZ2Q6okujHAesL73vzvKp1AviepBWS5rbaiKS5kpZLWt7X11chLDMzq6JKoleTeTGIOkdFxBRS984Zko5utpGIuDAipkbE1O7u7gphmZlZFVUSfS8wofB+PLChap2IaPy8D/gOqSvIzMx2kSpX3dwCTJJ0CHA3MBt4W6nOUmCepCWkk7CbImKjpP2AvSJic379WuAzQxe+2Z7DJ2ptRw2Y6CNim6R5wHWkyysXRcQqSafn8oXAMtKllT2kyytPzYs/B/hOuvqSEcA3IuK7Q74XZmbWUqXr6CNiGSmZF+ctLLwO4Iwmy60FDt/JGM3MbCf4zlgzs5pzojczqzknejOzmnOiNzOrOSd6M7Oac6I3M6s5J3ozs5pzojczqzk/StCsJjxEgrXiI3ozs5pzojczqzknejOzmnOiNzOrOSd6M7Oa81U31ja+SsRs1/ARvZlZzTnRm5nVnBO9mVnNOdGbmdWcE72ZWc050ZuZ1ZwTvZlZzVVK9JJmSFojqUfS/CblknReLl8paUqpvEvSbZKuGarAzcysmgETvaQu4HxgJjAZmCNpcqnaTGBSnuYCC0rlHwBW73S0ZmY2aFWO6KcBPRGxNiK2AkuAWaU6s4DFkdwMjJY0BkDSeOD1wEVDGLeZmVVUZQiEccD6wvte4MgKdcYBG4EvAx8BRu14mGa2szzkxJ6ryhG9msyLKnUkHQ/cFxErBtyINFfScknL+/r6KoRlZmZVVEn0vcCEwvvxwIaKdY4CTpC0jtTl8xpJlzfbSERcGBFTI2Jqd3d3xfDNzGwgVRL9LcAkSYdIGgnMBpaW6iwFTs5X30wHNkXExoj4WESMj4iJebkfRsQ7hnIHzMysfwP20UfENknzgOuALmBRRKySdHouXwgsA44DeoAtwKnDF7KZmQ1GpfHoI2IZKZkX5y0svA7gjAHWcQNww6AjNDOzneI7Y83Mas6J3sys5pzozcxqzonezKzmnOjNzGqu0lU3ZjvCt9ybdQYf0ZuZ1ZwTvZlZzTnRm5nVnBO9mVnNOdGbmdWcE72ZWc050ZuZ1ZwTvZlZzfmGKTPzzW015yN6M7Oac6I3M6s5J3ozs5pzojczqzknejOzmnOiNzOrOSd6M7Oaq5ToJc2QtEZSj6T5Tcol6bxcvlLSlDx/b0n/Jel2SaskfXqod8DMzPo3YKKX1AWcD8wEJgNzJE0uVZsJTMrTXGBBnv8Y8JqIOBw4ApghafoQxW5mZhVUOaKfBvRExNqI2AosAWaV6swCFkdyMzBa0pj8/uFc52l5iqEK3szMBlZlCIRxwPrC+17gyAp1xgEb8zeCFcALgPMj4qc7Hq51Gt86b9b5qhzRq8m88lF5yzoR8UREHAGMB6ZJeknTjUhzJS2XtLyvr69CWGZmVkWVRN8LTCi8Hw9sGGydiHgQuAGY0WwjEXFhREyNiKnd3d0VwjIzsyqqJPpbgEmSDpE0EpgNLC3VWQqcnK++mQ5sioiNkroljQaQtA9wLPCrIYzfzMwGMGAffURskzQPuA7oAhZFxCpJp+fyhcAy4DigB9gCnJoXHwNclvvp9wKuiohrhn43zMyslUrj0UfEMlIyL85bWHgdwBlNllsJvGwnYzQzs53gO2PNzGrOid7MrOac6M3Mas6J3sys5pzozcxqzonezKzmKl1eaWZ7No9ptHvzEb2ZWc050ZuZ1ZwTvZlZzTnRm5nVnBO9mVnNOdGbmdWcE72ZWc050ZuZ1ZwTvZlZzfnOWOuX74g02/35iN7MrOac6M3Mas6J3sys5pzozcxqzonezKzmKiV6STMkrZHUI2l+k3JJOi+Xr5Q0Jc+fIOl6SaslrZL0gaHeATMz69+AiV5SF3A+MBOYDMyRNLlUbSYwKU9zgQV5/jbgrIh4ETAdOKPJsmZmNoyqHNFPA3oiYm1EbAWWALNKdWYBiyO5GRgtaUxEbIyIWwEiYjOwGhg3hPGbmdkAqiT6ccD6wvtenpqsB6wjaSLwMuCngw3SzMx2XJVErybzYjB1JO0PfAs4MyIearoRaa6k5ZKW9/X1VQjLzMyqqJLoe4EJhffjgQ1V60h6GinJXxER3261kYi4MCKmRsTU7u7uKrGbmVkFVRL9LcAkSYdIGgnMBpaW6iwFTs5X30wHNkXERkkCLgZWR8Q/DGnkZmZWyYCDmkXENknzgOuALmBRRKySdHouXwgsA44DeoAtwKl58aOAk4CfS/pZnvfxiFg2tLthZmatVBq9MifmZaV5CwuvAzijyXI30bz/3sxqxKOcdjbfGWtmVnNO9GZmNedEb2ZWc070ZmY150RvZlZzfmbsHs5XS5jVn4/ozcxqzonezKzmnOjNzGrOid7MrOac6M3Mas6J3sys5pzozcxqzonezKzmnOjNzGrOid7MrOac6M3Mas5j3ZjZsPJ4Su3nI3ozs5pzojczqzknejOzmnMf/R7AfaRmezYf0ZuZ1VylRC9phqQ1knokzW9SLknn5fKVkqYUyhZJuk/SL4YycDMzq2bARC+pCzgfmAlMBuZImlyqNhOYlKe5wIJC2aXAjKEI1szMBq/KEf00oCci1kbEVmAJMKtUZxawOJKbgdGSxgBExI+AB4YyaDMzq65Koh8HrC+8783zBlunX5LmSlouaXlfX99gFjUzs35USfRqMi92oE6/IuLCiJgaEVO7u7sHs6iZmfWjSqLvBSYU3o8HNuxAHTMza4Mq19HfAkySdAhwNzAbeFupzlJgnqQlwJHApojYOKSRmlkt+T6P4TfgEX1EbAPmAdcBq4GrImKVpNMlnZ6rLQPWAj3A14D3NpaXdCXwE+AwSb2S3jXE+2BmZv2odGdsRCwjJfPivIWF1wGc0WLZOTsToJmZ7RwPgVAD/uprZv3xEAhmZjXnRG9mVnNO9GZmNedEb2ZWc070ZmY156tuzKyj+aqynVe7RN/fH0WrsoHK2/0H1alxmdnuoXaJ3sys3TrtoNJ99GZmNecjejOzHbA7dak60Q/C7vSLNdsT+H+yGif6DuE/WDMbLu6jNzOrOR/RD5FOO8tuZjv3TblO/7NO9Ga2Rxro4KxO3HVjZlZzTvRmZjXnRG9mVnNO9GZmNedEb2ZWc070ZmY1VynRS5ohaY2kHknzm5RL0nm5fKWkKVWXNTOz4TVgopfUBZwPzAQmA3MkTS5VmwlMytNcYMEgljUzs2FU5Yh+GtATEWsjYiuwBJhVqjMLWBzJzcBoSWMqLmtmZsNIEdF/BenNwIyIOC2/Pwk4MiLmFepcA3wpIm7K7/8d+CgwcaBlC+uYS/o2AHAYsGbndm3IHQTc3+4gmnBcg9epsTmuwenUuKA9sT0vIrqbFVQZAkFN5pU/HVrVqbJsmhlxIXBhhXjaQtLyiJja7jjKHNfgdWpsjmtwOjUu6LzYqiT6XmBC4f14YEPFOiMrLGtmZsOoSh/9LcAkSYdIGgnMBpaW6iwFTs5X30wHNkXExorLmpnZMBrwiD4itkmaB1wHdAGLImKVpNNz+UJgGXAc0ANsAU7tb9lh2ZPh16ndSo5r8Do1Nsc1OJ0aF3RYbAOejDUzs92b74w1M6s5J3ozs5pzom9BUpek2/I9Akg6W9Ldkn6Wp+PaFNc6ST/PMSzP854p6fuS/jv/PLBD4mp7m0kaLelqSb+StFrSKzqkvZrF1QntdVhh+z+T9JCkM9vdZv3E1Qlt9kFJqyT9QtKVkvZud3s9JUb30Tcn6UPAVOCAiDhe0tnAwxFxbpvjWgdMjYj7C/POAR6IiC/l8YQOjIiPdkBcZ9PmNpN0GXBjRFyUr/zaF/g47W+vZnGdSQf8jTXkIUzuBo4EzqDNbdYirlNpY5tJGgfcBEyOiEclXUW6OGUyHdJe4CP6piSNB14PXNTuWCqaBVyWX18GvLGNsXQMSQcARwMXA0TE1oh4kDa3Vz9xdZpjgN9ExJ101t9YMa5OMALYR9II0gf2BjqrvZzoW/gy8BHgD6X585RG51zUxq9iAXxP0oo8bATAc/J9C+Sfz+6QuKC9bXYo0AdckrvhLpK0H+1vr1ZxQWf8jTXMBq7Mr9vdZkXFuKCNbRYRdwPnAncBG0n3EH2PzmovJ/oySccD90XEilLRAuD5wBGkX+jf7+rYsqMiYgppRNAzJB3dpjjKmsXV7jYbAUwBFkTEy4BHgE4YKrtVXO1ur/+Ru5NOAL7ZrhiaaRJXW9ssf7DMAg4BxgL7SXrHroyhCif6pzoKOCH3OS8BXiPp8oi4NyKeiIg/AF8jjcy5y0XEhvzzPuA7OY57lUYLJf+8rxPi6oA26wV6I+Kn+f3VpATb7vZqGlcHtFfRTODWiLg3v293mzWNqwPa7Fjgjojoi4jHgW8Df0LntBfgRP8UEfGxiBgfERNJXxF/GBHvaPzSsjcBv9jVsUnaT9KoxmvgtTmOpcApudopwL92QlztbrOIuAdYL+mwPOsY4Je0ub1axdXu9iqZw/bdI21ts4Lt4uqANrsLmC5pX0ki/S5X0zntBfiqm35JejXw4XzVzddJXw8DWAe8p9EHtwvjOZR0tAzp6/83IuLzkp4FXAUcTPrDe0tEPNABcXVCmx1BOqk+ElhLukpjL9rYXv3EdR5tbq8c277AeuDQiNiU57X1b6yfuDrhb+zTwF8A24DbgNOA/Wlze20XoxO9mVm9uevGzKzmnOjNzGrOid7MrOac6M3Mas6J3sys5pzozcxqzonezKzm/j94bO3z/rsnMAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Example\n",
    "model = SFCN()\n",
    "model = torch.nn.DataParallel(model)\n",
    "fp_ = './brain_age/run_20190719_00_epoch_best_mae.p'\n",
    "model.load_state_dict(torch.load(fp_))\n",
    "\n",
    "# Example data: some random brain in the MNI152 1mm std space\n",
    "data = np.random.rand(182, 218, 182)\n",
    "label = np.array([71.3,]) # Assuming the random subject is 71.3-year-old.\n",
    "\n",
    "# Transforming the age to soft label (probability distribution)\n",
    "bin_range = [42,82]\n",
    "bin_step = 1\n",
    "sigma = 1\n",
    "y, bc = dpu.num2vect(label, bin_range, bin_step, sigma)\n",
    "y = torch.tensor(y, dtype=torch.float32)\n",
    "print(f'Label shape: {y.shape}')\n",
    "\n",
    "# Preprocessing\n",
    "data = data/data.mean()\n",
    "data = dpu.crop_center(data, (160, 192, 160))\n",
    "\n",
    "# Move the data from numpy to torch tensor on GPU\n",
    "sp = (1,1)+data.shape\n",
    "data = data.reshape(sp)\n",
    "input_data = torch.tensor(data, dtype=torch.float32)\n",
    "print(f'Input data shape: {input_data.shape}')\n",
    "print(f'dtype: {input_data.dtype}')\n",
    "\n",
    "# Evaluation\n",
    "model.eval() # Don't forget this. BatchNorm will be affected if not in eval mode.\n",
    "with torch.no_grad():\n",
    "    output = model.module(input_data)\n",
    "\n",
    "# Output, loss, visualisation\n",
    "x = output[0].reshape([1, -1])\n",
    "print(f'Output shape: {x.shape}')\n",
    "loss = dpl.my_KLDivLoss(x, y).numpy()\n",
    "\n",
    "# Prediction, Visualisation and Summary\n",
    "x = x.numpy().reshape(-1)\n",
    "y = y.numpy().reshape(-1)\n",
    "\n",
    "plt.bar(bc, y)\n",
    "plt.title('Soft label')\n",
    "plt.show()\n",
    "\n",
    "prob = np.exp(x)\n",
    "pred = prob@bc\n",
    "plt.bar(bc, prob)\n",
    "plt.title(f'Prediction: age={pred:.2f}\\nloss={loss}')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example for Sex Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Input data shape: torch.Size([1, 1, 160, 192, 160])\n",
      "dtype: torch.float32\n",
      "\n",
      "Predicted probability: \n",
      "Female\t35.79%,\n",
      "Male\t64.21%\n"
     ]
    }
   ],
   "source": [
    "# Example\n",
    "model = SFCN(output_dim=2, channel_number=[28, 58, 128, 256, 256, 64])\n",
    "model = torch.nn.DataParallel(model)\n",
    "fp_ = './sex_prediction/run_20191008_00_epoch_last.p'\n",
    "model.load_state_dict(torch.load(fp_))\n",
    "model.cuda()\n",
    "\n",
    "# Example data: some random brain in the MNI152 1mm std space\n",
    "data = np.random.rand(182, 218, 182)\n",
    "y = torch.tensor([1]) # Assuming Sex is Male (0=Female, 1=Male)\n",
    "\n",
    "# Preprocessing\n",
    "data = data/data.mean()\n",
    "data = dpu.crop_center(data, (160, 192, 160))\n",
    "\n",
    "# Move the data from numpy to torch tensor on GPU\n",
    "sp = (1,1)+data.shape\n",
    "data = data.reshape(sp)\n",
    "input_data = torch.tensor(data, dtype=torch.float32).cuda()\n",
    "print(f'Input data shape: {input_data.shape}')\n",
    "print(f'dtype: {input_data.dtype}')\n",
    "\n",
    "# Evaluation\n",
    "model.eval() # Don't forget this. BatchNorm will be affected if not in eval mode.\n",
    "with torch.no_grad():\n",
    "    output = model(input_data)\n",
    "\n",
    "# Output, loss, visualisation\n",
    "x = output[0].cpu().reshape([1, -1])\n",
    "loss = F.nll_loss(x, y)\n",
    "\n",
    "# Prediction, Visualisation and Summary\n",
    "x = np.exp(x.numpy().reshape(-1))\n",
    "\n",
    "print('\\nPredicted probability: \\nFemale\\t%.2f%%,\\nMale\\t%.2f%%'%(x[0]*100, x[1]*100))"
   ]
  }
 ],
 "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.6.10"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}