202003_models/HybridSN_GRSL2020.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "HybridSN_GRSL2020.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyPDsO6eQBt1IpNez93Op6KT",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/OUCTheoryGroup/colab_demo/blob/master/HybridSN_GRSL2020.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NTXgOJlhyK9q",
        "colab_type": "text"
      },
      "source": [
        "# HybridSN 高光谱分类\n",
        "\n",
        "S. K. Roy, G. Krishna, S. R. Dubey, B. B. Chaudhuri\n",
        "HybridSN: Exploring 3-D–2-D CNN Feature Hierarchy for Hyperspectral Image Classification, *IEEE GRSL* 2020\n",
        "\n",
        "这篇论文构建了一个 混合网络 解决高光谱图像分类问题，首先用 3D卷积，然后使用 2D卷积，代码相对简单，下面是代码的解析。\n",
        "\n",
        "首先取得数据，并引入基本函数库。"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MEhgw9ON0DxY",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "! wget http://www.ehu.eus/ccwintco/uploads/6/67/Indian_pines_corrected.mat\n",
        "! wget http://www.ehu.eus/ccwintco/uploads/c/c4/Indian_pines_gt.mat\n",
        "! pip install spectral"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "g6fvtNAMdk3J",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import scipy.io as sio\n",
        "from sklearn.decomposition import PCA\n",
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.metrics import confusion_matrix, accuracy_score, classification_report, cohen_kappa_score\n",
        "import spectral\n",
        "import torch\n",
        "import torchvision\n",
        "import torch.nn as nn\n",
        "import torch.nn.functional as F\n",
        "import torch.optim as optim"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0cIoYQIei9OW",
        "colab_type": "text"
      },
      "source": [
        "## 1. 定义 HybridSN 类\n",
        "\n",
        "模型的网络结构为如下图所示：\n",
        "\n",
        "![HybridSN](https://gaopursuit.oss-cn-beijing.aliyuncs.com/202003/20200304203710.jpg)\n",
        "\n",
        "三维卷积部分：\n",
        "- conv1：（1, 30, 25, 25）， 8个 7x3x3 的卷积核 ==>（8, 24, 23, 23）\n",
        "- conv2：（8, 24, 23, 23）， 16个 5x3x3 的卷积核 ==>（16, 20, 21, 21）\n",
        "- conv3：（16, 20, 21, 21），32个 3x3x3 的卷积核 ==>（32, 18, 19, 19）\n",
        "\n",
        "接下来要进行二维卷积，因此把前面的 32*18 reshape 一下，得到 （576, 19, 19）\n",
        "\n",
        "二维卷积：（576, 19, 19） 64个 3x3 的卷积核，得到 （64, 17, 17）\n",
        "\n",
        "接下来是一个 flatten 操作，变为 18496 维的向量，\n",
        "\n",
        "接下来依次为256，128节点的全连接层，都使用比例为0.4的 Dropout，\n",
        "\n",
        "最后输出为 16 个节点，是最终的分类类别数。 \n",
        "\n",
        "下面是 HybridSN 类的代码："
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Iy7oOuVoQhKS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class_num = 16\n",
        "\n",
        "class HybridSN(nn.Module):\n",
        "''' your code here '''\n",
        "\n",
        "# 随机输入，测试网络结构是否通\n",
        "# x = torch.randn(1, 1, 30, 25, 25)\n",
        "# net = HybridSN()\n",
        "# y = net(x)\n",
        "# print(y.shape)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AZzxfudHh0ea",
        "colab_type": "text"
      },
      "source": [
        "## 2. 创建数据集\n",
        "\n",
        "首先对高光谱数据实施PCA降维；然后创建 keras 方便处理的数据格式；然后随机抽取 10% 数据做为训练集，剩余的做为测试集。\n",
        "\n",
        "首先定义基本函数："
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bUVru6XtPZdn",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# 对高光谱数据 X 应用 PCA 变换\n",
        "def applyPCA(X, numComponents):\n",
        "    newX = np.reshape(X, (-1, X.shape[2]))\n",
        "    pca = PCA(n_components=numComponents, whiten=True)\n",
        "    newX = pca.fit_transform(newX)\n",
        "    newX = np.reshape(newX, (X.shape[0], X.shape[1], numComponents))\n",
        "    return newX\n",
        "\n",
        "# 对单个像素周围提取 patch 时，边缘像素就无法取了，因此，给这部分像素进行 padding 操作\n",
        "def padWithZeros(X, margin=2):\n",
        "    newX = np.zeros((X.shape[0] + 2 * margin, X.shape[1] + 2* margin, X.shape[2]))\n",
        "    x_offset = margin\n",
        "    y_offset = margin\n",
        "    newX[x_offset:X.shape[0] + x_offset, y_offset:X.shape[1] + y_offset, :] = X\n",
        "    return newX\n",
        "\n",
        "# 在每个像素周围提取 patch ，然后创建成符合 keras 处理的格式\n",
        "def createImageCubes(X, y, windowSize=5, removeZeroLabels = True):\n",
        "    # 给 X 做 padding\n",
        "    margin = int((windowSize - 1) / 2)\n",
        "    zeroPaddedX = padWithZeros(X, margin=margin)\n",
        "    # split patches\n",
        "    patchesData = np.zeros((X.shape[0] * X.shape[1], windowSize, windowSize, X.shape[2]))\n",
        "    patchesLabels = np.zeros((X.shape[0] * X.shape[1]))\n",
        "    patchIndex = 0\n",
        "    for r in range(margin, zeroPaddedX.shape[0] - margin):\n",
        "        for c in range(margin, zeroPaddedX.shape[1] - margin):\n",
        "            patch = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]   \n",
        "            patchesData[patchIndex, :, :, :] = patch\n",
        "            patchesLabels[patchIndex] = y[r-margin, c-margin]\n",
        "            patchIndex = patchIndex + 1\n",
        "    if removeZeroLabels:\n",
        "        patchesData = patchesData[patchesLabels>0,:,:,:]\n",
        "        patchesLabels = patchesLabels[patchesLabels>0]\n",
        "        patchesLabels -= 1\n",
        "    return patchesData, patchesLabels\n",
        "\n",
        "def splitTrainTestSet(X, y, testRatio, randomState=345):\n",
        "    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=testRatio, random_state=randomState, stratify=y)\n",
        "    return X_train, X_test, y_train, y_test"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "g2HjOjuwhfLD",
        "colab_type": "text"
      },
      "source": [
        "下面读取并创建数据集："
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7b-811cVew0O",
        "colab_type": "code",
        "outputId": "691130ff-af07-4c7e-80ae-2fe17d0fa87b",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 311
        }
      },
      "source": [
        "# 地物类别\n",
        "class_num = 16\n",
        "X = sio.loadmat('Indian_pines_corrected.mat')['indian_pines_corrected']\n",
        "y = sio.loadmat('Indian_pines_gt.mat')['indian_pines_gt']\n",
        "\n",
        "# 用于测试样本的比例\n",
        "test_ratio = 0.90\n",
        "# 每个像素周围提取 patch 的尺寸\n",
        "patch_size = 25\n",
        "# 使用 PCA 降维，得到主成分的数量\n",
        "pca_components = 30\n",
        "\n",
        "print('Hyperspectral data shape: ', X.shape)\n",
        "print('Label shape: ', y.shape)\n",
        "\n",
        "print('\\n... ... PCA tranformation ... ...')\n",
        "X_pca = applyPCA(X, numComponents=pca_components)\n",
        "print('Data shape after PCA: ', X_pca.shape)\n",
        "\n",
        "print('\\n... ... create data cubes ... ...')\n",
        "X_pca, y = createImageCubes(X_pca, y, windowSize=patch_size)\n",
        "print('Data cube X shape: ', X_pca.shape)\n",
        "print('Data cube y shape: ', y.shape)\n",
        "\n",
        "print('\\n... ... create train & test data ... ...')\n",
        "Xtrain, Xtest, ytrain, ytest = splitTrainTestSet(X_pca, y, test_ratio)\n",
        "print('Xtrain shape: ', Xtrain.shape)\n",
        "print('Xtest  shape: ', Xtest.shape)\n",
        "\n",
        "# 改变 Xtrain, Ytrain 的形状，以符合 keras 的要求\n",
        "Xtrain = Xtrain.reshape(-1, patch_size, patch_size, pca_components, 1)\n",
        "Xtest  = Xtest.reshape(-1, patch_size, patch_size, pca_components, 1)\n",
        "print('before transpose: Xtrain shape: ', Xtrain.shape) \n",
        "print('before transpose: Xtest  shape: ', Xtest.shape) \n",
        "\n",
        "# 为了适应 pytorch 结构，数据要做 transpose\n",
        "Xtrain = Xtrain.transpose(0, 4, 3, 1, 2)\n",
        "Xtest  = Xtest.transpose(0, 4, 3, 1, 2)\n",
        "print('after transpose: Xtrain shape: ', Xtrain.shape) \n",
        "print('after transpose: Xtest  shape: ', Xtest.shape) \n",
        "\n",
        "\n",
        "\"\"\" Training dataset\"\"\"\n",
        "class TrainDS(torch.utils.data.Dataset): \n",
        "    def __init__(self):\n",
        "        self.len = Xtrain.shape[0]\n",
        "        self.x_data = torch.FloatTensor(Xtrain)\n",
        "        self.y_data = torch.LongTensor(ytrain)        \n",
        "    def __getitem__(self, index):\n",
        "        # 根据索引返回数据和对应的标签\n",
        "        return self.x_data[index], self.y_data[index]\n",
        "    def __len__(self): \n",
        "        # 返回文件数据的数目\n",
        "        return self.len\n",
        "\n",
        "\"\"\" Testing dataset\"\"\"\n",
        "class TestDS(torch.utils.data.Dataset): \n",
        "    def __init__(self):\n",
        "        self.len = Xtest.shape[0]\n",
        "        self.x_data = torch.FloatTensor(Xtest)\n",
        "        self.y_data = torch.LongTensor(ytest)\n",
        "    def __getitem__(self, index):\n",
        "        # 根据索引返回数据和对应的标签\n",
        "        return self.x_data[index], self.y_data[index]\n",
        "    def __len__(self): \n",
        "        # 返回文件数据的数目\n",
        "        return self.len\n",
        "\n",
        "# 创建 trainloader 和 testloader\n",
        "trainset = TrainDS()\n",
        "testset  = TestDS()\n",
        "train_loader = torch.utils.data.DataLoader(dataset=trainset, batch_size=128, shuffle=True, num_workers=2)\n",
        "test_loader  = torch.utils.data.DataLoader(dataset=testset,  batch_size=128, shuffle=False, num_workers=2)"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Hyperspectral data shape:  (145, 145, 200)\n",
            "Label shape:  (145, 145)\n",
            "\n",
            "... ... PCA tranformation ... ...\n",
            "Data shape after PCA:  (145, 145, 30)\n",
            "\n",
            "... ... create data cubes ... ...\n",
            "Data cube X shape:  (10249, 25, 25, 30)\n",
            "Data cube y shape:  (10249,)\n",
            "\n",
            "... ... create train & test data ... ...\n",
            "Xtrain shape:  (1024, 25, 25, 30)\n",
            "Xtest  shape:  (9225, 25, 25, 30)\n",
            "before transpose: Xtrain shape:  (1024, 25, 25, 30, 1)\n",
            "before transpose: Xtest  shape:  (9225, 25, 25, 30, 1)\n",
            "after transpose: Xtrain shape:  (1024, 1, 30, 25, 25)\n",
            "after transpose: Xtest  shape:  (9225, 1, 30, 25, 25)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "V1Jjx8ysqYXA",
        "colab_type": "text"
      },
      "source": [
        "## 3. 开始训练"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "y2higQPzadIO",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "outputId": "d02bf75a-bff0-4cee-89c1-8d2317692131"
      },
      "source": [
        "# 使用GPU训练，可以在菜单 \"代码执行工具\" -> \"更改运行时类型\" 里进行设置\n",
        "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
        "\n",
        "# 网络放到GPU上\n",
        "net = HybridSN().to(device)\n",
        "criterion = nn.CrossEntropyLoss()\n",
        "optimizer = optim.Adam(net.parameters(), lr=0.001)\n",
        "\n",
        "# 开始训练\n",
        "total_loss = 0\n",
        "for epoch in range(100):\n",
        "    for i, (inputs, labels) in enumerate(train_loader):\n",
        "        inputs = inputs.to(device)\n",
        "        labels = labels.to(device)\n",
        "        # 优化器梯度归零\n",
        "        optimizer.zero_grad()\n",
        "        # 正向传播 +　反向传播 + 优化 \n",
        "        outputs = net(inputs)\n",
        "        loss = criterion(outputs, labels)\n",
        "        loss.backward()\n",
        "        optimizer.step()\n",
        "        total_loss += loss.item()\n",
        "    print('[Epoch: %d]   [loss avg: %.4f]   [current loss: %.4f]' %(epoch + 1, total_loss/(epoch+1), loss.item()))\n",
        "\n",
        "print('Finished Training')"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[Epoch: 1]   [loss avg: 22.0578]   [current loss: 2.6951]\n",
            "[Epoch: 2]   [loss avg: 21.1697]   [current loss: 2.4139]\n",
            "[Epoch: 3]   [loss avg: 20.5813]   [current loss: 2.5084]\n",
            "[Epoch: 4]   [loss avg: 20.1962]   [current loss: 2.4382]\n",
            "[Epoch: 5]   [loss avg: 19.9130]   [current loss: 2.3063]\n",
            "[Epoch: 6]   [loss avg: 19.6433]   [current loss: 2.2822]\n",
            "[Epoch: 7]   [loss avg: 19.3210]   [current loss: 2.2242]\n",
            "[Epoch: 8]   [loss avg: 18.8895]   [current loss: 1.8448]\n",
            "[Epoch: 9]   [loss avg: 18.4002]   [current loss: 1.8071]\n",
            "[Epoch: 10]   [loss avg: 17.8621]   [current loss: 1.6316]\n",
            "[Epoch: 11]   [loss avg: 17.3014]   [current loss: 1.4054]\n",
            "[Epoch: 12]   [loss avg: 16.7531]   [current loss: 1.0829]\n",
            "[Epoch: 13]   [loss avg: 16.2118]   [current loss: 1.1922]\n",
            "[Epoch: 14]   [loss avg: 15.6703]   [current loss: 1.0161]\n",
            "[Epoch: 15]   [loss avg: 15.1338]   [current loss: 0.8618]\n",
            "[Epoch: 16]   [loss avg: 14.5846]   [current loss: 0.6342]\n",
            "[Epoch: 17]   [loss avg: 14.0319]   [current loss: 0.6196]\n",
            "[Epoch: 18]   [loss avg: 13.5041]   [current loss: 0.5754]\n",
            "[Epoch: 19]   [loss avg: 12.9558]   [current loss: 0.3481]\n",
            "[Epoch: 20]   [loss avg: 12.4367]   [current loss: 0.3802]\n",
            "[Epoch: 21]   [loss avg: 11.9417]   [current loss: 0.2809]\n",
            "[Epoch: 22]   [loss avg: 11.4674]   [current loss: 0.2940]\n",
            "[Epoch: 23]   [loss avg: 11.0233]   [current loss: 0.1445]\n",
            "[Epoch: 24]   [loss avg: 10.6067]   [current loss: 0.0417]\n",
            "[Epoch: 25]   [loss avg: 10.2273]   [current loss: 0.1158]\n",
            "[Epoch: 26]   [loss avg: 9.8660]   [current loss: 0.0978]\n",
            "[Epoch: 27]   [loss avg: 9.5320]   [current loss: 0.0638]\n",
            "[Epoch: 28]   [loss avg: 9.2183]   [current loss: 0.1600]\n",
            "[Epoch: 29]   [loss avg: 8.9181]   [current loss: 0.0551]\n",
            "[Epoch: 30]   [loss avg: 8.6377]   [current loss: 0.1212]\n",
            "[Epoch: 31]   [loss avg: 8.3814]   [current loss: 0.0531]\n",
            "[Epoch: 32]   [loss avg: 8.1361]   [current loss: 0.0632]\n",
            "[Epoch: 33]   [loss avg: 7.9017]   [current loss: 0.0711]\n",
            "[Epoch: 34]   [loss avg: 7.6787]   [current loss: 0.0405]\n",
            "[Epoch: 35]   [loss avg: 7.4688]   [current loss: 0.0547]\n",
            "[Epoch: 36]   [loss avg: 7.2757]   [current loss: 0.0477]\n",
            "[Epoch: 37]   [loss avg: 7.0911]   [current loss: 0.0191]\n",
            "[Epoch: 38]   [loss avg: 6.9156]   [current loss: 0.1832]\n",
            "[Epoch: 39]   [loss avg: 6.7468]   [current loss: 0.0815]\n",
            "[Epoch: 40]   [loss avg: 6.5903]   [current loss: 0.0150]\n",
            "[Epoch: 41]   [loss avg: 6.4381]   [current loss: 0.0347]\n",
            "[Epoch: 42]   [loss avg: 6.2904]   [current loss: 0.0292]\n",
            "[Epoch: 43]   [loss avg: 6.1473]   [current loss: 0.0067]\n",
            "[Epoch: 44]   [loss avg: 6.0115]   [current loss: 0.0122]\n",
            "[Epoch: 45]   [loss avg: 5.8820]   [current loss: 0.0044]\n",
            "[Epoch: 46]   [loss avg: 5.7571]   [current loss: 0.0157]\n",
            "[Epoch: 47]   [loss avg: 5.6375]   [current loss: 0.0384]\n",
            "[Epoch: 48]   [loss avg: 5.5240]   [current loss: 0.0018]\n",
            "[Epoch: 49]   [loss avg: 5.4141]   [current loss: 0.0067]\n",
            "[Epoch: 50]   [loss avg: 5.3088]   [current loss: 0.0278]\n",
            "[Epoch: 51]   [loss avg: 5.2079]   [current loss: 0.0088]\n",
            "[Epoch: 52]   [loss avg: 5.1092]   [current loss: 0.0115]\n",
            "[Epoch: 53]   [loss avg: 5.0156]   [current loss: 0.0261]\n",
            "[Epoch: 54]   [loss avg: 4.9321]   [current loss: 0.0229]\n",
            "[Epoch: 55]   [loss avg: 4.8470]   [current loss: 0.0198]\n",
            "[Epoch: 56]   [loss avg: 4.7631]   [current loss: 0.0116]\n",
            "[Epoch: 57]   [loss avg: 4.6818]   [current loss: 0.0235]\n",
            "[Epoch: 58]   [loss avg: 4.6045]   [current loss: 0.0164]\n",
            "[Epoch: 59]   [loss avg: 4.5325]   [current loss: 0.0681]\n",
            "[Epoch: 60]   [loss avg: 4.4630]   [current loss: 0.0059]\n",
            "[Epoch: 61]   [loss avg: 4.3946]   [current loss: 0.0497]\n",
            "[Epoch: 62]   [loss avg: 4.3327]   [current loss: 0.0958]\n",
            "[Epoch: 63]   [loss avg: 4.2662]   [current loss: 0.0127]\n",
            "[Epoch: 64]   [loss avg: 4.2014]   [current loss: 0.0102]\n",
            "[Epoch: 65]   [loss avg: 4.1389]   [current loss: 0.0122]\n",
            "[Epoch: 66]   [loss avg: 4.0796]   [current loss: 0.0259]\n",
            "[Epoch: 67]   [loss avg: 4.0216]   [current loss: 0.0177]\n",
            "[Epoch: 68]   [loss avg: 3.9684]   [current loss: 0.0060]\n",
            "[Epoch: 69]   [loss avg: 3.9142]   [current loss: 0.0285]\n",
            "[Epoch: 70]   [loss avg: 3.8611]   [current loss: 0.0087]\n",
            "[Epoch: 71]   [loss avg: 3.8088]   [current loss: 0.0023]\n",
            "[Epoch: 72]   [loss avg: 3.7571]   [current loss: 0.0115]\n",
            "[Epoch: 73]   [loss avg: 3.7080]   [current loss: 0.0018]\n",
            "[Epoch: 74]   [loss avg: 3.6591]   [current loss: 0.0415]\n",
            "[Epoch: 75]   [loss avg: 3.6129]   [current loss: 0.0043]\n",
            "[Epoch: 76]   [loss avg: 3.5671]   [current loss: 0.0745]\n",
            "[Epoch: 77]   [loss avg: 3.5228]   [current loss: 0.0284]\n",
            "[Epoch: 78]   [loss avg: 3.4795]   [current loss: 0.0200]\n",
            "[Epoch: 79]   [loss avg: 3.4367]   [current loss: 0.0017]\n",
            "[Epoch: 80]   [loss avg: 3.3959]   [current loss: 0.0405]\n",
            "[Epoch: 81]   [loss avg: 3.3556]   [current loss: 0.0115]\n",
            "[Epoch: 82]   [loss avg: 3.3150]   [current loss: 0.0003]\n",
            "[Epoch: 83]   [loss avg: 3.2780]   [current loss: 0.0104]\n",
            "[Epoch: 84]   [loss avg: 3.2410]   [current loss: 0.0113]\n",
            "[Epoch: 85]   [loss avg: 3.2046]   [current loss: 0.0045]\n",
            "[Epoch: 86]   [loss avg: 3.1682]   [current loss: 0.0094]\n",
            "[Epoch: 87]   [loss avg: 3.1324]   [current loss: 0.0016]\n",
            "[Epoch: 88]   [loss avg: 3.0970]   [current loss: 0.0011]\n",
            "[Epoch: 89]   [loss avg: 3.0637]   [current loss: 0.0091]\n",
            "[Epoch: 90]   [loss avg: 3.0309]   [current loss: 0.0027]\n",
            "[Epoch: 91]   [loss avg: 2.9992]   [current loss: 0.0014]\n",
            "[Epoch: 92]   [loss avg: 2.9677]   [current loss: 0.0148]\n",
            "[Epoch: 93]   [loss avg: 2.9366]   [current loss: 0.0113]\n",
            "[Epoch: 94]   [loss avg: 2.9063]   [current loss: 0.0042]\n",
            "[Epoch: 95]   [loss avg: 2.8805]   [current loss: 0.1005]\n",
            "[Epoch: 96]   [loss avg: 2.8511]   [current loss: 0.0133]\n",
            "[Epoch: 97]   [loss avg: 2.8231]   [current loss: 0.0075]\n",
            "[Epoch: 98]   [loss avg: 2.7946]   [current loss: 0.0003]\n",
            "[Epoch: 99]   [loss avg: 2.7677]   [current loss: 0.0170]\n",
            "[Epoch: 100]   [loss avg: 2.7406]   [current loss: 0.0065]\n",
            "Finished Training\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ixPqqlcJrVuu",
        "colab_type": "text"
      },
      "source": [
        "## 4. 模型测试"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "u3FgXj-N3rPv",
        "colab_type": "code",
        "outputId": "b4da39cd-22db-4d4c-d004-cd5afc51a4d3",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 415
        }
      },
      "source": [
        "count = 0\n",
        "# 模型测试\n",
        "for inputs, _ in test_loader:\n",
        "    inputs = inputs.to(device)\n",
        "    outputs = net(inputs)\n",
        "    outputs = np.argmax(outputs.detach().cpu().numpy(), axis=1)\n",
        "    if count == 0:\n",
        "        y_pred_test =  outputs\n",
        "        count = 1\n",
        "    else:\n",
        "        y_pred_test = np.concatenate( (y_pred_test, outputs) )\n",
        "\n",
        "# 生成分类报告\n",
        "classification = classification_report(ytest, y_pred_test, digits=4)\n",
        "print(classification)"
      ],
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "              precision    recall  f1-score   support\n",
            "\n",
            "         0.0     1.0000    0.9268    0.9620        41\n",
            "         1.0     0.9835    0.9253    0.9535      1285\n",
            "         2.0     0.9293    0.9853    0.9565       747\n",
            "         3.0     0.9534    0.8638    0.9064       213\n",
            "         4.0     0.9929    0.9609    0.9766       435\n",
            "         5.0     0.9569    0.9467    0.9518       657\n",
            "         6.0     1.0000    1.0000    1.0000        25\n",
            "         7.0     0.9795    1.0000    0.9896       430\n",
            "         8.0     0.7500    0.6667    0.7059        18\n",
            "         9.0     0.9454    0.9897    0.9671       875\n",
            "        10.0     0.9753    0.9824    0.9788      2210\n",
            "        11.0     0.9380    0.9625    0.9501       534\n",
            "        12.0     0.9944    0.9514    0.9724       185\n",
            "        13.0     0.9659    0.9956    0.9805      1139\n",
            "        14.0     0.9815    0.9193    0.9494       347\n",
            "        15.0     0.9342    0.8452    0.8875        84\n",
            "\n",
            "    accuracy                         0.9653      9225\n",
            "   macro avg     0.9550    0.9326    0.9430      9225\n",
            "weighted avg     0.9657    0.9653    0.9651      9225\n",
            "\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RD2dCvA3y04n",
        "colab_type": "text"
      },
      "source": [
        "可以看出，准确率为 96.53%，性能还可以了~~~\n",
        "\n",
        "## 5. 备用函数\n",
        "\n",
        "下面是用于计算各个类准确率，显示结果的备用函数，以供参考\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CPgJh16kzm0X",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from operator import truediv\n",
        "\n",
        "def AA_andEachClassAccuracy(confusion_matrix):\n",
        "    counter = confusion_matrix.shape[0]\n",
        "    list_diag = np.diag(confusion_matrix)\n",
        "    list_raw_sum = np.sum(confusion_matrix, axis=1)\n",
        "    each_acc = np.nan_to_num(truediv(list_diag, list_raw_sum))\n",
        "    average_acc = np.mean(each_acc)\n",
        "    return each_acc, average_acc\n",
        "\n",
        "\n",
        "def reports (test_loader, y_test, name):\n",
        "    count = 0\n",
        "    # 模型测试\n",
        "    for inputs, _ in test_loader:\n",
        "        inputs = inputs.to(device)\n",
        "        outputs = net(inputs)\n",
        "        outputs = np.argmax(outputs.detach().cpu().numpy(), axis=1)\n",
        "        if count == 0:\n",
        "            y_pred =  outputs\n",
        "            count = 1\n",
        "        else:\n",
        "            y_pred = np.concatenate( (y_pred, outputs) )\n",
        "\n",
        "    if name == 'IP':\n",
        "        target_names = ['Alfalfa', 'Corn-notill', 'Corn-mintill', 'Corn'\n",
        "                        ,'Grass-pasture', 'Grass-trees', 'Grass-pasture-mowed', \n",
        "                        'Hay-windrowed', 'Oats', 'Soybean-notill', 'Soybean-mintill',\n",
        "                        'Soybean-clean', 'Wheat', 'Woods', 'Buildings-Grass-Trees-Drives',\n",
        "                        'Stone-Steel-Towers']\n",
        "    elif name == 'SA':\n",
        "        target_names = ['Brocoli_green_weeds_1','Brocoli_green_weeds_2','Fallow','Fallow_rough_plow','Fallow_smooth',\n",
        "                        'Stubble','Celery','Grapes_untrained','Soil_vinyard_develop','Corn_senesced_green_weeds',\n",
        "                        'Lettuce_romaine_4wk','Lettuce_romaine_5wk','Lettuce_romaine_6wk','Lettuce_romaine_7wk',\n",
        "                        'Vinyard_untrained','Vinyard_vertical_trellis']\n",
        "    elif name == 'PU':\n",
        "        target_names = ['Asphalt','Meadows','Gravel','Trees', 'Painted metal sheets','Bare Soil','Bitumen',\n",
        "                        'Self-Blocking Bricks','Shadows']\n",
        "    \n",
        "    classification = classification_report(y_test, y_pred, target_names=target_names)\n",
        "    oa = accuracy_score(y_test, y_pred)\n",
        "    confusion = confusion_matrix(y_test, y_pred)\n",
        "    each_acc, aa = AA_andEachClassAccuracy(confusion)\n",
        "    kappa = cohen_kappa_score(y_test, y_pred)\n",
        "    \n",
        "    return classification, confusion, oa*100, each_acc*100, aa*100, kappa*100"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mDsF8FLA0aOz",
        "colab_type": "text"
      },
      "source": [
        "检测结果写在文件里："
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Xx5Iyy2Cz5oH",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "classification, confusion, oa, each_acc, aa, kappa = reports(test_loader, ytest, 'IP')\n",
        "classification = str(classification)\n",
        "confusion = str(confusion)\n",
        "file_name = \"classification_report.txt\"\n",
        "\n",
        "with open(file_name, 'w') as x_file:\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{} Kappa accuracy (%)'.format(kappa))\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{} Overall accuracy (%)'.format(oa))\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{} Average accuracy (%)'.format(aa))\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{}'.format(classification))\n",
        "    x_file.write('\\n')\n",
        "    x_file.write('{}'.format(confusion))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zkrwmJS00j30",
        "colab_type": "text"
      },
      "source": [
        "下面代码用于显示分类结果："
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JSWV5XxX0rLp",
        "colab_type": "code",
        "outputId": "04176ec6-bc41-48d7-d4dd-d85e8aba89cf",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 155
        }
      },
      "source": [
        "# load the original image\n",
        "X = sio.loadmat('Indian_pines_corrected.mat')['indian_pines_corrected']\n",
        "y = sio.loadmat('Indian_pines_gt.mat')['indian_pines_gt']\n",
        "\n",
        "height = y.shape[0]\n",
        "width = y.shape[1]\n",
        "\n",
        "X = applyPCA(X, numComponents= pca_components)\n",
        "X = padWithZeros(X, patch_size//2)\n",
        "\n",
        "# 逐像素预测类别\n",
        "outputs = np.zeros((height,width))\n",
        "for i in range(height):\n",
        "    for j in range(width):\n",
        "        if int(y[i,j]) == 0:\n",
        "            continue\n",
        "        else :\n",
        "            image_patch = X[i:i+patch_size, j:j+patch_size, :]\n",
        "            image_patch = image_patch.reshape(1,image_patch.shape[0],image_patch.shape[1], image_patch.shape[2], 1)\n",
        "            X_test_image = torch.FloatTensor(image_patch.transpose(0, 4, 3, 1, 2)).to(device)                                   \n",
        "            prediction = net(X_test_image)\n",
        "            prediction = np.argmax(prediction.detach().cpu().numpy(), axis=1)\n",
        "            outputs[i][j] = prediction+1\n",
        "    if i % 20 == 0:\n",
        "        print('... ... row ', i, ' handling ... ...')"
      ],
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "... ... row  0  handling ... ...\n",
            "... ... row  20  handling ... ...\n",
            "... ... row  40  handling ... ...\n",
            "... ... row  60  handling ... ...\n",
            "... ... row  80  handling ... ...\n",
            "... ... row  100  handling ... ...\n",
            "... ... row  120  handling ... ...\n",
            "... ... row  140  handling ... ...\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mV6KtTVw4h1i",
        "colab_type": "code",
        "outputId": "a5b7830e-29d6-40ba-82e1-f8d1b2a00023",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 323
        }
      },
      "source": [
        "predict_image = spectral.imshow(classes = outputs.astype(int),figsize =(5,5))"
      ],
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAATkAAAEyCAYAAACBJqcyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAcoklEQVR4nO3dfYxl9X3f8fenEHBxUh5MQ9e7tGwb\naotaSU0QO5YtyzKJjZmVoZJlgawGJ8yusus4JLVkFlvUqtxIWIniEDW7zu7gGFcITIkT0I4Tl2Ai\n/9OZeBcnPBO2+IFdLV4s16Spqya03/5xzpm5c+c+nnPuPef+7ueFRnPvufee8+PM7G8+5/d0FBGY\nmaXqHzRdADOzSXIlZ2ZJcyVnZklzJWdmSXMlZ2ZJcyVnZkmbWCUn6VpJz0s6IenApI5jZjaIJjFO\nTtJZwF8DPw+cBL4B3BQRz9R+MDOzASaV5K4GTkTEixHxd8D9wPUTOpaZWV9nT2i/24GXOp6fBHb1\ne7N0ccBlEypKFceBnx3wmk1Odt7PP/+VZg5//qvZ91fPH/kjr7763U3PL35znQWq1/efa7oEtfp+\nRPzjfi9OqpIbStJeYG/27J8Cx5oqygCif7k0zYLMoey8v+Mdh5o5/OJK9n1lceSPrKzs3/T8hi/U\nWJ6aLS80XYJafWfQi5Oq5E4Bl3Y835FvWxcRh4HDANJVLZtA6wqsaYuLDVVulpxJtcl9A7hc0k5J\n5wA3Ag9P6FhmZn1NJMlFxGuSfgX4KnAW8PmIeHoSx6pXd4KrEDCXVrduS+waIWljXKZau02sTS4i\nvgJ8ZVL7NzMbRWMdD+1SY4IbKNvvKnsAWFhayjY74c2eomOiU49N1rzkKrn1CoQjI7x7uh0MRdm2\n6L60Xd7FvHZ+LC4ebLoIo+l5Oetaro08d9XMkjZzSW5YUquW4Gq8TO1xCbply/Ly5udFoltaA7rT\nnS9pW6fEWDqbPic5M0taS5LccTanq/6JarSk1k/9CS66Pqv8GKssjbyPBfJEV6S1XsNPem3r/IxN\nnxPcTHCSM7OktSTJdduchqqlt439bdWy2WRleNCx2UBOcmaWtJYmucx6W9VE2uESNoPtdzMzPs5m\njpOcmSWt1UlugxjefjZOYkugLa6MLTMr2pPsutdi61YkvZVtV2bPTz8+8TLNk8WD9fUUr+xv18wP\nJzkzS9qMJDmop21tThNcArqTXmdWcHueDeIkZ2ZJm6EkV0XXEkeVx93NoRb32I7anmfzyUnOzJI2\nJ0kua89bmECiW8v3tatN6bDOdNUvwQ16vQXprlO/pFc54XkVkpkwJ5XcZnVerhb7amOXRiz1XiRA\n3Us81a3Fl7adBl3mDq0Ae60MbK3ky1UzS9qcJbliGEqNuWtpLfu+vGvsj27c42HCyapLkfD2LGff\nl7cu5zn8MrWMGUl4MEpnxuL6ZeodK9nA5HnPdnfs3wa07zw4yZlZ0uYsyU1AnkJaswzA8kJHYhqc\nEI8sZa8fYSNRNvL/0ZnwWpjqesmmLmWZpW3JpSmfPng6ezA4BE+dk5yZJW1Ok9zkF9FcGJKigIm1\nxUXRPrg03v41iXa4cRVlmJFEZ+3nJGdmSZudJBebvg1tO1pVvnR6jJFm1Guvo6W7KgNLh/XkjW2p\nz02sh9nUnmdVDQqjPs3T4yRnZkmbnSSXG7X3b6wEV1iPiZ15cXL9jcUCkLXYdDvDcm19sbTUUO/q\nMiyPfgvHFPRLeU549Sud5CRdKukxSc9IelrSrfn2iyQ9IumF/PuF9RXXzGw8VZLca8DHIuJxST8B\nHJf0CPBh4NGIuFPSAeAAcFv1ok5RdOQZjdoKOL5iCe+mx1kVc1n7zXWduDlLcYMM61R20htf6SQX\nEacj4vH88f8EngW2A9cD9+Rvuwe4oWohzczKqqVNTtJlwFuBNeCSiMiHPvMycEmfz+wF9tZx/Ikq\nUp26e1lbM8ehtJ4J7kgeJfb0jgzFGDwVc3ZtqnolvTrS3Tg3n6nzpjfTULmSk/TjwB8CvxYRf6OO\nYRgREdKW2qF47TBwON9HO1YqKv5hH5mhgaibOhzGU1Rumyq7PpVbwZVb+0z7EreoEGelsqs0hETS\nj5FVcPdGxJfzzd+TtC1/fRtwploRzczKK53klEW2u4FnI+K3O156GLgZuDP//lClErZFcdlaJJ2u\nURrdw0Fm8b6gwh0AKZr3GXJVLlffDvxb4ElJf5lv+wRZ5faApFuA7wAfrFZEM7PyFNF8c9hobXLF\nvK6WNPh3FyNv+BgnwR1d2ZfvqmtnFRpRosRA4PXl0Ns0PmHe44eN43hEXNXvRU/rMrOkzc60rmEJ\nrrtXcFI9pOvH2bz/Mm1wuxcPZQ/KjAbuk3REg4N6+xkyLMWpzSbJSc7MkjY7SW6YYqHIutvstuyu\n4dQxQurpvuVg48luyNg7s0lykjOzpKWT5OpIcGpJD27NbVSNT8DvNqyNzqxGTnJmlrTZSXJTCVfa\n9G0Uoyx7vnI0W958cffBnttbM9e/aoLsHmc3S3OALVlOcmaWtBma8WBm1tPAGQ+zc7lqc211wOIB\nw+5x2/3Zzvf32+9I983tMitLD/Uyznpys8aXq2aWtISS3LSveLt6C2bggjv25ItkHil3N6+B2tJ5\nMoJBqbCKWVtMcl44yZlZ0hJKcjbMRBKcWcs5yZlZ0lqd5EYZ3bJx35zuRqEZaCSzmdPZnlemB9am\nz0nOzJLW6iQ3iiLtacud7juTnVOd1a9IdXOX6FaH9E4vtOt8OMmZWdJaneSKLDZKDos+yyNlbXZu\nr7PxTGos3UwbluAGva/BdOckZ2ZJa3WSq5K3es/579Ve16tdz8w2GTXFlf38BJOek5yZJa3VSa5I\nVv3a2wYZ5TPFGLutY+2c6MyA6gmu6nFqSHitruQKUpSq6IbZuEzdcsTOd9V+XLPWm1blNgW+XDWz\npM1EkmvWiMNP2nKnrznk4R41SijBFSonOUlnSfqmpKP5852S1iSdkPQlSedUL6aZWTl1JLlbgWeB\nf5Q//wzw2Yi4X9LngFuAQ1UOMIn2uPKGlMWJbu6lvJT4LKqU5CTtABYhm7wnScC7gQfzt9wD3FDl\nGGZmVVRNcr8DfBz4ifz5G4AfRsRr+fOTwPbyu8/S0Obez/LDSiYqPPzERhSweChbIn099Q37tcl/\nvTrbHz998PSmfRTLrt+xfxsAy6v9h19UvcXuLCmd5CTtBs5ExPGSn98r6ZikY2XLYGY2TJUk93bg\n/ZKuA15H1iZ3F3CBpLPzNLcDONXrwxFxGDgM4953tb4Et3V5pio72/Jg/Y/zXMwWK/H/OPDGOi0L\n6nVb2ZcnuP2D37eq/ss5FYltpev5oAQ3j0onuYi4PSJ2RMRlwI3A1yLiQ8BjwAfyt90MPFS5lGZm\nJU1inNxtwP2S/iPwTeDuCRyjlGlPwPdsscGSu7FOwz9fJ7jeaqnkIuLPgT/PH78IXF3Hfs3Mqpqr\nGQ/dPbLas5o9WB7wJ3hpLX/Pruwp2fNlurqnuu+6s3VCbP92Jic8s4nx3FUzS9pcJblucSRLY1re\nPPaus+2ueM9G2ssS3dDWj85k1yvVdfLq7GYTM9eVXKH7MrbnQOPubeMMcRjlUrbsvqtwZVqr9eEe\nUb0DoNjXOJYWss+4A2IzX66aWdKc5MZVxwDicZPdpPgyuVZ1JLjufZVJdLaZk5yZJc1JblxF21yd\n4auNyc6prhXW0+EIvxJui+vNSc7MkpZMktudL11zdF8aCxYuFeOUd+Ubph3u3F5niXCSM7OkJZPk\nqiS4YnrX+sDfFilmlbG6efvUFz10srMZ5SRnZkmb2SRXbwdkMb2rzn1O1lLbkt0onP6sAU5yZpa0\ndie57rjSoO72uqaGsvXTeapae5MSLx7auOVdNS75PyOc5Mwsae1Oci2yvsDmurbGpRa0141L4Ztx\nT8nSWnaeZ6j5uTJXcnOgdZXelvpM83T11Kjl9fu2zk8158tVM0uak9wcal2yS0Sdi2ZOwvKugLU9\nTRdj6pzkzCxpTnI2G8NPbKhhy58vrYnl5eKH3c60OQlOcmaWNCc528TtdbNr2KKZy8ursJS/Z36C\nnJOcmaWt3UlufcVIOtYcsmlyskvH0tJCxzi5+eEkZ2ZJa3eS6xwGP0dtCG3mZDdd62PvRvgHsLyQ\nTb5fWu09fcTj5EqQdIGkByU9J+lZSW+TdJGkRyS9kH+/sK7CmpmNq2qSuwv404j4gKRzgPOATwCP\nRsSdkg4AB4Dbyu2+4TV5ltYcVYbwGLvJGueWhP3arTd+LvM5Qbh0kpN0PvBO4G6AiPi7iPghcD1w\nT/62e4AbqhbSzKysKkluJ/AK8AeSfgY4DtwKXBIRp/P3vAxc0uvDkvYCeyscf/I6e3dtKLfXNaz4\nfc07UDdmN4zxg1hIr/G7SiV3NnAl8NGIWJN0F9ml6bqICEk9rzkj4jBwGKDfe9olb9QluyRYbvF6\ncm3hSq8Z65XbkbzSS6/eGkuVjoeTwMmIKBoCHiSr9L4naRtA/v1MtSKamZVXOslFxMuSXpL0poh4\nHrgGeCb/uhm4M//+UC0lbViR4Ky8KsmumLJUTEJvs2GLHFe5bhk8lCQ/cPGWOU9whaq9qx8F7s17\nVl8EfpEsHT4g6RbgO8AHKx7DzKw0RTTfHNa/Ta75shV/HZe6bmG/pU2uzvMobflfrzKrrUhMLbr5\n2RZl2+tWWzBNaZSBukMN+/WZz9EfozoeEVf1e9HTuswsaS2f1jVv+v+5Xh8dkGjTYNn2us4U1YZU\nBxuhrFL4Ku4Od8Rd0lU5yZlZ0mYvyRVRxgN1rUt321jTyW6sRLflTXmCcw9pZU5yZpa02UtyU0xw\nS6x1DDnatb5tcob/7U+9ba6ffj3Dg9ru2pLsammjs9Kc5MwsabOT5BpoiyvSG8Di4qHswcqVEzzi\n6H/ruxOdE97GbIh+49aaTnbTSnQxYNCd5jBPOsmZWdJmJ8k13Ju6srIPaL5Nrp95S3C9rN+Sr+Vj\n7AK3z03T7FRySSt+5cefGjavl6l1m/alrDsjpseXq2aWNCe5VvDf9bZpupPC6uMkZ2ZJa/lSSzOk\nlvOoTd/abNLLNm0Z5DvG0lPTWGa9M9l1p75xfhMm8aOewyEkXmrJzOaX2+TqoiT/QrZH1+ltet56\nrwHHffPT0vrtsyofd3VA02CCN9qqhZOcmSWtVUluUFvCuMqPPBu23zESW95Ot7S2+TPF9KP1watd\nlnfVO+4tiVsBjvKD7LuKfo0pe8Cu+r0UXQlOnXHM8WvinOTMLGmtSnJ1akV3bdFO1720d1eCK24G\nvLSURS7PXiipSGyT7KzvtethtyCcSEFsVE5yZpa0ZJPcLCkS3Cxoz20NB8wSqbMNbhLcDjdVTnJm\nljQnuSka1qtq42h5WrPWcCU3Ra7czKbPl6tmljQnuTnSnk4Ds+mplOQk/bqkpyU9Jek+Sa+TtFPS\nmqQTkr4k6Zy6CmtmNq7SSU7SduBXgSsi4n9LegC4EbgO+GxE3C/pc8AtwKFaSms25zz6ZHxV2+TO\nBv6hpLOB84DTwLuBB/PX7wFuqHgMM7PSSldyEXEK+C3gu2SV26vAceCHEfFa/raTwPaqhRyFVveg\n1T3TOJRZK0X+n21WupKTdCFwPbATeCPweuDaMT6/V9IxScfKlsHMbJgqvas/B3wrIl4BkPRl4O3A\nBZLOztPcDuBUrw9HxGHgcP7Zyn9+YuFI1V2YzbR+y4DNe7qr0ib3XWBB0nmSBFwDPAM8Bnwgf8/N\nwEPVimhmVl6VNrk1sg6Gx4En830dBm4D/p2kE8AbgLtrKKdZfVTzV8up4795VGkwcER8CvhU1+YX\ngaur7NfMrC6e1jUFywuJLEFuNoNcyZlZ0jx3dUTz2p5hNutcyQ3hys1stvly1cyS5iTXg9ObWTqc\n5MwsaU5yHZzgzNLjJGdmSWtVknOSMrO6OcmZWdJcyZlZ0lzJmVnSXMmZWdJcyZlZ0lzJmVnSXMmZ\nWdJcyZlZ0lo1GNisL48Tt5Kc5MwsaU5ylWX3tFxaddQwayMnOTNLmpNcRU5w9Vje1X2X967zGjXe\nBV7+mc0TJzkzS5qTnI1k8veNdbqyyXCSM7OkuZIzs6S5kjOzpLlNzsZTZy/nIO4BtZoMTXKSPi/p\njKSnOrZdJOkRSS/k3y/Mt0vS70o6IekJSVdOsvBmZsOMkuS+APwn4Isd2w4Aj0bEnZIO5M9vA94H\nXJ5/7QIO5d/N2qdEKl3VHgA+ffD0+raV/Su1FcnqNzTJRcTXgR90bb4euCd/fA9wQ8f2L0ZmFbhA\n0ra6CmtmNq6yHQ+XRETxp+xl4JL88XbgpY73ncy3bSFpr6Rjko6VLINN0PKu6DELwRbiCAtxpOli\n2Bgq965GRFDMUh/vc4cj4qqIuKpqGax+S2tiaS2xxn8p+4qo3IGysu9oTYWySStbyX2vuAzNv5/J\nt58CLu143458m5lZI8pWcg8DN+ePbwYe6tj+C3kv6wLwasdlrVk7FImu4jCVlX1HnehmwNDeVUn3\nAe8CLpZ0EvgUcCfwgKRbgO8AH8zf/hXgOuAE8CPgFydQZjOzkQ2t5CLipj4vXdPjvQF8pGqhzKam\nO81Na7CzTY2ndZlZ0lzJmXWqoa3O2sWVnJklzZWcWS9OdMlwJWdmSfNSS5Y+J7K55iRnZklzkrPZ\n56RmAzjJmVnSnORsoI3llpyWbDa5krOBiuWWlhsuR1++VLUhfLlqZklzkrPZNKUEt3ho96bnxdJK\nndt9h4d2c5Izs6Q5ydlAbel4iPzwuxcPAc2lp84Ed8f+bY2WxUbjJGdmSXOSs5mg9UC5v5HjF6mt\n0HnfVWs3JzkzS5qTnLXK4u4+bW4NjYdb2Z+VZEt59rstblY4yZlZ0pzkrFVWju7LHjTU9mbpcZIz\ns6QpWnALNknNF8LMZtXxiLiq34tOcmaWNLfJzZHFxYMAHF3Zt2n7ntWNnsuNGQ5dvNqHzShfrpZ0\ncHGx6SL0tG/lKIcWN08q37+SDXYIBp9mRccQjqP7ul5Ms5JbPNjOn2PTiqEzM8KXq2Y2v3y5mpju\nFDdIMdl92x3ZcI1QsGe1eLxv0/48oGO+FAl3xhJdT0OTnKTPSzoj6amObb8p6TlJT0j6I0kXdLx2\nu6QTkp6X9N5JFdzMbBRD2+QkvRP4W+CLEfGWfNt7gK9FxGuSPgMQEbdJugK4D7gaeCPwZ8C/jIj/\nO+QYM9Mm19a2uEFGbZMbRIne48FtcqNpeaKr1iYXEV8HftC17b9GxGv501VgR/74euD+iPg/EfEt\n4ARZhWdm1og6Oh5+CfiT/PF24KWO107m22bewcXFmUxxZnWY5cRbqeNB0ieB14B7S3x2L7C3yvHN\nzIYpXclJ+jCwG7gmNhr2TgGXdrxtR75ti4g4DBzO99XaNjmnN7PMrPa4lrpclXQt8HHg/RHxo46X\nHgZulHSupJ3A5cBfVC+mmVk5Q5OcpPuAdwEXSzoJfAq4HTgXeETZSPjViPjliHha0gPAM2SXsR8Z\n1rPaZk5xZlvNWqIbWslFxE09Nt894P2/AfxGlUKZmdXFMx56cIKzQnEDm2nfuKap47bl+HXy3FUz\nS5qTXIciwe1fOdqxbfS5oOMqjjPJY1g1TSWZphNU08evk5damiPjTOvafTCreI/uzyriFKd1zfIA\n1zZoUceDl1oys/nly1WAChPX22yVPQAssDz2Z4sEZzbrnOTMLGlOcglb4Ej+KEtyKbarmQ3jJGdm\nSXMlZ2ZJcyVnZklzJWdmSXMlZ2ZJcyVnZklzJWdmSXMlZ2ZJcyVnZknzjAebO159pJwWrToyFldy\nZtbXrFZsnXy5amZJc5Jrk+iYQO91RM1q4SRnZklzkgNoyxJE6vvEzEpykjOzpLmSM7OkuZIzs6S5\nTW5CVlmqZT/F/S+LO5oPex8AK/lg18VsjFMKY53MynKSM7OkOcm13LAE1/v9jwPwac9eMhue5CR9\nXtIZSU/1eO1jkkLSxflzSfpdSSckPSHpykkU2sxsVIoYPLJe0juBvwW+GBFv6dh+Kdm97t4M/GxE\nfF/SdcBHgeuAXcBdEbFraCGU3vD+utrk6lDm5tJmM+R4RFzV78WhSS4ivg78oMdLnwU+zubbz19P\nVhlGRKwCF0ga73rLzKxGpToeJF0PnIqIv+p6aTvwUsfzk/m2XvvYK+mYpGNlymBmNoqxOx4knQd8\nAnhPlQNHxGHgcL7P5C5XzawdyvSu/gtgJ/BXkgB2AI9Luho4BVza8d4d+TYzs0aMXclFxJPATxbP\nJX0buCrveHgY+BVJ95N1PLwaEad778ncIWA2eaMMIbkP+G/AmySdlHTLgLd/BXgROAEcAfbXUkoz\ns5KGJrmIuGnI65d1PA7gI9WLZWZWD0/rMrOkDR0MPJVCuHfVzMqrNhjYzGyWtWWC/veB/5V/b5uL\naV+52lgmaGe52lgmaGe52lgmGF6ufzbow624XAWQdGxQ5GxKG8vVxjJBO8vVxjJBO8vVxjJB9XL5\nctXMkuZKzsyS1qZK7nDTBeijjeVqY5mgneVqY5mgneVqY5mgYrla0yZnZjYJbUpyZma1a0UlJ+la\nSc/ny6YfaKgMl0p6TNIzkp6WdGu+/SJJj0h6If9+YQNlO0vSNyUdzZ/vlLSWn68vSTqngTJdIOlB\nSc9JelbS25o+V5J+Pf/ZPSXpPkmva+Jc9bplQL9zM81bBvQp12/mP8MnJP2RpAs6Xrs9L9fzkt47\nrTJ1vFbL7RUar+QknQX8HvA+4ArgJklXNFCU14CPRcQVwALwkbwcB4BHI+Jy4NH8+bTdCjzb8fwz\nwGcj4qeA/wEMWjRhUu4C/jQi3gz8TF6+xs6VpO3Ar5KtiPMW4CzgRpo5V18Aru3a1u/cvA+4PP/a\nCxyacrkeAd4SET8N/DVwO0D+u38j8K/yzxzM/61Oo0zF7RXeA3y3Y3O5cxURjX4BbwO+2vH8duD2\nFpTrIeDngeeBbfm2bcDzUy7HDrJ/FO8GjgIiGxh5dq/zN6UynQ98i7xNt2N7Y+eKjVWpLyIb5H4U\neG9T5wq4DHhq2LkBfh+4qdf7plGurtf+DXBv/njTv0Pgq8DbplUm4EGyP57fBi6ucq4aT3KMsWT6\ntEi6DHgrsAZcEhtr4r0MXDLl4vwO2b00/l/+/A3ADyPitfx5E+drJ/AK8Af5ZfSypNfT4LmKiFPA\nb5H95T8NvAocp/lzVeh3btr0+/9LwJ/kjxsrVx23V+jUhkquVST9OPCHwK9FxN90vhbZn4+pdUdL\n2g2ciYjj0zrmiM4GrgQORcRbyabkbbo0beBcXUh2I6WdwBuB19PjMqgNpn1uRiHpk2RNNvc2XI7i\n9gr/vq59tqGSa82S6ZJ+jKyCuzcivpxv/l5xx7H8+5kpFuntwPvz1ZfvJ7tkvYvsLmjFvOMmztdJ\n4GRErOXPHySr9Jo8Vz8HfCsiXomIvwe+THb+mj5XhX7npvHff0kfBnYDH8or4CbL1Xl7hW+zcXuF\nf1K2TG2o5L4BXJ73gp1D1tj58LQLIUnA3cCzEfHbHS89DNycP76ZrK1uKiLi9ojYEdnCpDcCX4uI\nDwGPAR9ookx5uV4GXpL0pnzTNcAzNHiuyC5TFySdl/8sizI1eq469Ds3DwO/kPccLjDlWwZIupas\nOeT9EfGjrvLeKOlcSTvJGvv/YtLliYgnI+InI+Ky/Pf+JHBl/jtX7lxNqoFzzIbH68h6dv478MmG\nyvAOskuIJ4C/zL+uI2sDexR4Afgz4KKGyvcu4Gj++J+T/cKdAP4LcG4D5fnXwLH8fP0xcGHT5wr4\nD8BzwFPAfwbObeJcAfeRtQv+ff6P9JZ+54asI+n38t/9J8l6h6dZrhNk7VzF7/znOt7/ybxczwPv\nm1aZul7/NhsdD6XOlWc8mFnS2nC5amY2Ma7kzCxpruTMLGmu5Mwsaa7kzCxpruTMLGmu5Mwsaa7k\nzCxp/x+NkbioTOxXeAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 360x360 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    }
  ]
}