hw9_anomaly_detection.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU",
    "gpuClass": "standard"
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 599
        },
        "id": "GVrLsEg3Jvww",
        "outputId": "0bb61e60-992c-4d8c-c3fc-d14d898fd1d7"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0th epoch starting.\n",
            "1th epoch starting.\n",
            "2th epoch starting.\n",
            "3th epoch starting.\n",
            "4th epoch starting.\n",
            "5th epoch starting.\n",
            "6th epoch starting.\n",
            "7th epoch starting.\n",
            "8th epoch starting.\n",
            "9th epoch starting.\n",
            "Time ellapsed in training is: 63.605716705322266\n",
            "---------------------\n",
            "10.422374729761481 6.172500033003813\n",
            "---------------------\n",
            "threshold:  28.939874828772922\n",
            "Type 1 Error: 0.8900 %\n",
            "Type 2 Error: 4.2800 % \n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZgU1b3/8fcXGBgERDY3UBaDIioOMqCIcnG5CGrQeDXRy43yM+4aF/SHaG6EaJLHa8z9GbdHcU+CisZoiKKIV2RzgQGRKwKCCDqIgIAssgzL9/fHqWZ6pmugZ5juHmY+r+fpp6tOVXWfGpr+dJ1zqsrcHRERkfLq5boCIiJSMykgREQklgJCRERiKSBERCSWAkJERGIpIEREJJYCQqQCZvammV1W3etWsg79zKy4ul9XJB0Ncl0BkepkZhuTZvcDtgI7ovmr3X10uq/l7gMzsa7IvkIBIbWKuzdNTJvZEuAKd3+n/Hpm1sDdt2ezbiL7GjUxSZ2QaKoxs9vN7FvgGTNrYWavm9kqM1sbTbdL2uY9M7simh5iZlPN7P5o3S/NbGAV1+1oZpPNbIOZvWNmj5jZX9Pcj6Oj9/rezOaa2aCkZWeb2WfR6y4zs9ui8tbRvn1vZmvMbIqZ6f++7JE+JFKXHAy0BNoDVxE+/89E84cDm4GHd7P9icACoDVwH/CUmVkV1n0emA60AkYCP0+n8maWB/wTeBs4EPglMNrMjopWeYrQjNYMOBZ4Nyq/FSgG2gAHAXcCusaO7JECQuqSncAId9/q7pvdfbW7v+Lum9x9A/A74F92s/1Sd3/C3XcAzwGHEL5w017XzA4HegJ3uXuJu08FxqZZ/5OApsC90bbvAq8Dl0TLtwFdzWx/d1/r7rOSyg8B2rv7Nnef4roIm6RBASF1ySp335KYMbP9zOxxM1tqZuuBycABZla/gu2/TUy4+6Zosmkl1z0UWJNUBvB1mvU/FPja3XcmlS0F2kbT/wacDSw1s0lm1jsq/wOwCHjbzBab2fA030/qOAWE1CXlfzXfChwFnOju+wN9o/KKmo2qw3KgpZntl1R2WJrbfgMcVq7/4HBgGYC7z3D38wjNT68BL0XlG9z9VnfvBAwChprZGXu5H1IHKCCkLmtG6Hf43sxaAiMy/YbuvhQoAkaaWcPoV/6P09z8I2ATMMzM8sysX7Tti9FrDTaz5u6+DVhPaFLDzM41sx9FfSDrCMN+d8a/hUgpBYTUZQ8AjYHvgA+Bt7L0voOB3sBq4LfAGML5Grvl7iWEQBhIqPOjwKXuPj9a5efAkqi57JrofQA6A+8AG4EPgEfdfWK17Y3UWqa+KpHcMrMxwHx3z/gRjEhl6AhCJMvMrKeZHWFm9cxsAHAeoc9ApEbRmdQi2Xcw8HfCeRDFwLXu/nFuqySSSk1MIiISS01MIiISq9Y0MbVu3do7dOiQ62qIiOxTZs6c+Z27t4lbVmsCokOHDhQVFeW6GiIi+xQzW1rRMjUxiYhILAWEiIjEUkCIiEisWtMHISI117Zt2yguLmbLli17XlkyIj8/n3bt2pGXl5f2NgoIEcm44uJimjVrRocOHaj4HkuSKe7O6tWrKS4upmPHjmlvpyYmEcm4LVu20KpVK4VDjpgZrVq1qvQRnAJCRLJC4ZBbVfn7KyA2bIARI2D69FzXRESkRlFAbN0Kd98NH32U65qISIasXr2agoICCgoKOPjgg2nbtu2u+ZKSkt1uW1RUxI033rjH9zj55JOrpa7vvfce5557brW81t5SJ/V+0Z0fN23a/Xoiss9q1aoVs2fPBmDkyJE0bdqU2267bdfy7du306BB/NdhYWEhhYWFe3yP999/v3oqW4PoCCI/Pzxv3pzbeohIVg0ZMoRrrrmGE088kWHDhjF9+nR69+5N9+7dOfnkk1mwYAFQ9hf9yJEjufzyy+nXrx+dOnXiwQcf3PV6TZs23bV+v379uPDCC+nSpQuDBw8mcdXscePG0aVLF3r06MGNN964xyOFNWvWcP7559OtWzdOOukk5syZA8CkSZN2HQF1796dDRs2sHz5cvr27UtBQQHHHnssU6ZM2eu/kY4g6tWDRo0UECLZcvPNEP2arzYFBfDAA5XerLi4mPfff5/69euzfv16pkyZQoMGDXjnnXe48847eeWVV1K2mT9/PhMnTmTDhg0cddRRXHvttSnnFnz88cfMnTuXQw89lD59+jBt2jQKCwu5+uqrmTx5Mh07duSSSy7ZY/1GjBhB9+7dee2113j33Xe59NJLmT17Nvfffz+PPPIIffr0YePGjeTn5zNq1CjOOussfvWrX7Fjxw42VUOriAICQjOTmphE6pyLLrqI+vXrA7Bu3Touu+wyFi5ciJmxbdu22G3OOeccGjVqRKNGjTjwwANZsWIF7dq1K7NOr169dpUVFBSwZMkSmjZtSqdOnXadh3DJJZcwatSo3dZv6tSpu0Lq9NNPZ/Xq1axfv54+ffowdOhQBg8ezAUXXEC7du3o2bMnl19+Odu2beP888+noKBgr/42oIAI8vNBZ3iKZEcVfulnSpMmTXZN//rXv+a0007j1VdfZcmSJfTr1y92m0aNGu2arl+/Ptu3b6/SOntj+PDhnHPOOYwbN44+ffowfvx4+vbty+TJk3njjTcYMmQIQ4cO5dJLL92r91EfBECDBlDN/4Aism9Zt24dbdu2BeDZZ5+t9tc/6qijWLx4MUuWLAFgzJgxe9zm1FNPZfTo0UDo22jdujX7778/X3zxBccddxy33347PXv2ZP78+SxdupSDDjqIK6+8kiuuuIJZs2btdZ0VEKCAEBGGDRvGHXfcQffu3av9Fz9A48aNefTRRxkwYAA9evSgWbNmNG/efLfbjBw5kpkzZ9KtWzeGDx/Oc889B8ADDzzAscceS7du3cjLy2PgwIG89957HH/88XTv3p0xY8Zw00037XWda809qQsLC73KNww68kjo0QNeeKF6KyUiAMybN4+jjz4619XIuY0bN9K0aVPcneuvv57OnTtzyy23ZO394/4dzGymu8eO49URBOgIQkSy4oknnqCgoIBjjjmGdevWcfXVV+e6SrulTmqAvDwFhIhk3C233JLVI4a9pSMI0BGESBbUlubsfVVV/v4KCFBAiGRYfn4+q1evVkjkSOJ+EPmJK0ekSU1MoIAQybB27dpRXFzMqlWrcl2VOitxR7nKUECAAkIkw/Ly8ip1JzOpGdTEBAoIEZEYCggIAVHBdVdEROoqBQSEK7ru3JnrWoiI1CgKCAAz0OgKEZEyFBAQjiAUECIiZSggIBxBqIlJRKQMBQSoiUlEJIYCAtTEJCISQwEBamISEYmhgAA1MYmIxMhoQJjZADNbYGaLzGx4zPKhZvaZmc0xs/8xs/ZJyy4zs4XR47JM1lNNTCIiqTIWEGZWH3gEGAh0BS4xs67lVvsYKHT3bsDfgPuibVsCI4ATgV7ACDNrkam6qolJRCRVJo8gegGL3H2xu5cALwLnJa/g7hPdfVM0+yGQuNTgWcAEd1/j7muBCcCAjNVUTUwiIikyGRBtga+T5oujsor8AnizMtua2VVmVmRmRXt1GWE1MYmIpKgRndRm9h9AIfCHymzn7qPcvdDdC9u0abM3FVATk4hIOZkMiGXAYUnz7aKyMszsTOBXwCB331qZbauNmphERFJkMiBmAJ3NrKOZNQQuBsYmr2Bm3YHHCeGwMmnReKC/mbWIOqf7R2WZoYAQEUmRsTvKuft2M7uB8MVeH3ja3eea2d1AkbuPJTQpNQVeNjOAr9x9kLuvMbN7CCEDcLe7r8lUXXW5bxGRVBm95ai7jwPGlSu7K2n6zN1s+zTwdOZql0RHECIiKWpEJ3XOKSBERFIoIEDDXEVEYiggQMNcRURiKCBATUwiIjEUEKAmJhGRGAoIUBOTiEgMBQSoiUlEJIYCAtTEJCISQwEBamISEYmhgAA1MYmIxFBAgJqYRERiKCBATUwiIjEUEKAmJhGRGAoIUBOTiEgMBQSoiUlEJIYCAtTEJCISQwEBamISEYmhgAA1MYmIxFBAgJqYRERiKCBAASEiEkMBAQoIEZEYCghQQIiIxFBAQBjFpE5qEZEyFBCgIwgRkRgKCFBAiIjEUEBACAgRESlDAQGlAaGjCBGRXRQQoIAQEYmhgIAwigkUECIiSRQQUHoEoaGuIiK7KCBATUwiIjEUEKCAEBGJoYAABYSISAwFBCggRERiKCBAASEiEkMBARrmKiISI6MBYWYDzGyBmS0ys+Exy/ua2Swz225mF5ZbtsPMZkePsZmsp4a5ioikapCpFzaz+sAjwL8CxcAMMxvr7p8lrfYVMAS4LeYlNrt7QabqV4aamEREUmQsIIBewCJ3XwxgZi8C5wG7AsLdl0TLcvvTXQEhIpIik01MbYGvk+aLo7J05ZtZkZl9aGbnx61gZldF6xStWrWq6jVVQIiIpKjJndTt3b0Q+HfgATM7ovwK7j7K3QvdvbBNmzZVfycFhIhIikwGxDLgsKT5dlFZWtx9WfS8GHgP6F6dlStDo5hERFJkMiBmAJ3NrKOZNQQuBtIajWRmLcysUTTdGuhDUt9FtdMoJhGRFBkLCHffDtwAjAfmAS+5+1wzu9vMBgGYWU8zKwYuAh43s7nR5kcDRWb2CTARuLfc6KfqpSYmEZEUmRzFhLuPA8aVK7sraXoGoemp/HbvA8dlsm5lKCBERFLU5E7q7FFAiIikUECAAkJEJIYCAhQQIiIxFBCgYa4iIjEUEKBhriIiMRQQoCYmEZEYCghQQIiIxFBAgAJCRCSGAgIUECIiMRQQoIAQEYmhgIDSYa4axSQisosCAnQEISISI62AMLMmZlYvmj7SzAaZWV5mq5ZFCggRkRTpHkFMJtwCtC3wNvBz4NlMVSrrFBAiIinSDQhz903ABcCj7n4RcEzmqpVlCggRkRRpB4SZ9QYGA29EZfUzU6UcUECIiKRINyBuBu4AXo3uCteJcKe32kEX6xMRSZHWHeXcfRIwCSDqrP7O3W/MZMWyShfrExFJke4opufNbH8zawJ8CnxmZv83s1XLIjUxiYikSLeJqau7rwfOB94EOhJGMtUOCggRkRTpBkRedN7D+cBYd98G1J5vUwWEiEiKdAPicWAJ0ASYbGbtgfWZqlTWKSBERFKk20n9IPBgUtFSMzstM1XKAQWEiEiKdDupm5vZf5tZUfT4I+FoonbQMFcRkRTpNjE9DWwAfho91gPPZKpSWadhriIiKdJqYgKOcPd/S5r/jZnNzkSFckJNTCIiKdI9gthsZqckZsysD7A5M1XKAQWEiEiKdI8grgH+bGbNo/m1wGWZqVIOKCBERFKkO4rpE+B4M9s/ml9vZjcDczJZuaxRQIiIpKjUHeXcfX10RjXA0AzUJzc0iklEJMXe3HLUqq0WuaZRTCIiKfYmIGrPz201MYmIpNhtH4SZbSA+CAxonJEa5YICQkQkxW4Dwt2bZasiOaWAEBFJsTdNTLWHAkJEJIUCAhQQIiIxMhoQZjbAzBaY2SIzGx6zvK+ZzTKz7WZ2Yblll5nZwuiR2ZPyNMxVRCRFxgLCzOoDjwADga7AJWbWtdxqXwFDgOfLbdsSGAGcCPQCRphZi0zVVcNcRURSZfIIohewyN0Xu3sJ8CJwXvIK7r7E3ecA5b+ZzwImuPsad18LTAAGZKymamISEUmRyYBoC3ydNF8clVXbtmZ2VeIeFatWrapyRRUQIiKp9ulOancf5e6F7l7Ypk2bqr+QAkJEJEUmA2IZcFjSfLuoLNPbVp4CQkQkRSYDYgbQ2cw6mllD4GJgbJrbjgf6m1mLqHO6f1SWGQoIEZEUGQsId98O3ED4Yp8HvOTuc83sbjMbBGBmPc2sGLgIeNzM5kbbrgHuIYTMDODuqCwzEsNcNYpJRGSXdG8YVCXuPg4YV67srqTpGYTmo7htnybcCzvzdAQhIpJin+6krjYKCBGRFAoIUECIiMRQQIACQkQkhgICFBAiIjEUEKCL9YmIxFBAgC7WJyISQwEBamISEYmhgAAFhIhIDAUEKCBERGIoIECX2hARiaGAAMjLC8/bt+e2HiIiNYgCAqBBdEmqwYNhxYrc1kVEpIZQQEBpQAAce2zu6iEiUoMoIKBsQHz3Xe7qISJSgyggoGxAAAwcqBFNIlLnKSCgtJM64a23YOvW3NRFRKSGUEBA6hEEQOPG2a+HiEgNooCA+IAAnRchInWaAgIqDoiFC7NbDxGRGkQBARUHRJcu2a2HiEgNooCA0msxiYjILgqIhJKS+PLx47NbDxGRGkIBkVB+qGvC5MnZrYeISA2hgNiT3/8+1zUQEckJBUSykpLwOPPMsuWNGsGTT+amTiIiOaKASJaXFx4TJpQtLymBK6/MTZ1ERHJEAVGRG2/MdQ1ERHJKAVGRCy9MLdMF/ESkDlFAVOTww1PL6tWDjz/WnedEpE5QQFSkfXuYNi21/IQTKh4SKyJSiyggdufkk+Ghh3JdCxGRnFBA7Ml11+W6BiIiOaGA2JN69eAXv8h1LUREsk4BkY7hw1PL/uM/1FktIrWaAiIdnTrB+eeXLRs9Gs44Izf1ERHJAgVEOurVg1dfhbFjy5brQn4iUotlNCDMbICZLTCzRWaW0k5jZo3MbEy0/CMz6xCVdzCzzWY2O3o8lsl6pu3HP4YxY8qWbd2am7qIiGRYxgLCzOoDjwADga7AJWbWtdxqvwDWuvuPgP8H/FfSsi/cvSB6XJOpelbaT39adj4/H+67T2dZi0itk8kjiF7AIndf7O4lwIvAeeXWOQ94Lpr+G3CG2T54e7fbb4cPPsh1LUREqlUmA6It8HXSfHFUFruOu28H1gGtomUdzexjM5tkZqfGvYGZXWVmRWZWtGrVquqt/e58/31qmUY0iUgtU1M7qZcDh7t7d2Ao8LyZ7V9+JXcf5e6F7l7Ypk2b7NWuefPUsilTsvf+IiJZkMmAWAYcljTfLiqLXcfMGgDNgdXuvtXdVwO4+0zgC+DIDNa18ho2LDv/n/+Zm3qIiGRIJgNiBtDZzDqaWUPgYqDcOFHGApdF0xcC77q7m1mbqJMbM+sEdAYWZ7CulbdmDbRuXbbskktyUxcRkQzIWEBEfQo3AOOBecBL7j7XzO42s0HRak8BrcxsEaEpKTEUti8wx8xmEzqvr3H3NZmqa5U0aQKrVkG3bqVlL76Yu/qIiFQz81oyPLOwsNCLioqy/8YlJeGe1cn++lcYPDj7dRERqSQzm+nuhXHLamon9b6jfF8EwJ13Zr8eIiLVTAFRHcrfnvSrr2D+/NzURUSkmiggqsNLL6WWHX00LFqU/bqIiFQTBUR1qOjk786dYc6cML15c3iIiOwjFBDVpaQkPiiOPz6U77cftGgB48enNkmJiNRADXJdgVojLw82bYLGjSteZ+tWGDAgTH//PRxwQHbqJiJSBTqCqE75+dC9e3rrtmgBCxbA4pp1/p+ISIICorpNmwbXXQfz5u153S5d4IgjQpPTa6+lt42ISJYoIKpb48bwyCPhy79Ll/S2eeUV+MlPoGtX+OUvQ5/FBRfA6tWZrauIyG4oIDLpgw/gs88qt83DD4fnV18N13oqKYEbb4THH6/++omI7IY6qTPpgAPC4x//gPPK3yspTcmX8bj66rLLdu6EH36AZs2qXkcRkQroCCIbBg2CHTugbVt4/fVwKY7p0yv/OmbhMXVqOMeidWvYf/9w1vbEiWXXLSqCd9+tnvqLSJ2ki/XlmjvUq+acdi89J6OW/PuKSGboYn01mRl880143H9/9b1mwpdfwocfhs5zszDK6ocfwr0rvv22et5PRGolBURNcMgh4XHrreEX/6mxt+CumsGDoXdv2LIlzJ9yCjRtGu5dccghoQO8uLh0/e3by86LSJ2lgKiJJk+GSZPg4ov3/rU++GD3yx96CA47LBxdHH986Nc47LDQh3HQQaH89tvD87JlsH596E/Zvj0cmYhIraU+iH3F2rXh6OLww0MTUU1y0EHhPI7LLw8d8YceGsonTgzNZn//e+pNlUSkRlAfRG3QogW0bBn6Kt56q7TfYscO+N3vclu3FSvgscegV68QEPffDxs2wOmnw7hxcO65cM898PXXsHAhfPppOL8DQtPXqlWlr/XRR2XnRSRndARRG7iHs65XrYIOHeDzz8PjtNPg5ZfDpT8AOnYMndb7gksvDeeP7NwJv/0tHHUUDBsWQvGJJ2D48BA0mzaF+4PXr5/rGovsk3Z3BKGAqAsefRSuvz4ESOvWoaxXL5gxA0aOhDPOqN6O8Wzr1q30vhsAf/oT3HQT/POf8OMfh33+8svQOR9n27ZwNd7def31cMn200+PXz51Khx8MPzoR1XbB5EcUUDUdYl/4+Thrzt3hkeD6GT6J54IJ+999RWsWxeaeiCMrPrjH7Nb30xo0iRc52rlSnj66RCQZ58dAhJCoCxZEsq//RbOOSeEZl5eeCSfVzJ1arjI4rx58OabIXw7dixdvjdKSsLRUOKIqKQk/r7nItVEASGV4x76Ctq2DV+sAGvWhBFRZ54JV14JEyaEL8lp00q3Gz8+lBUVVc8IrH3RkUeG5r1OncJgghUrwuivyZNDsHz0Ufj7uId7gmzdGgKoZcuwvVkIp9dfLw2lb74JIdG4cWhKbNQoHBUlzqxP2L69NPDTkbh/SUV3RJQ6QQEhmbVjR+iUTr4B0vjxcPLJ4Zf5M8+UlnfuHMIn4eijoXlzDZmFcI7K1KlhujJHbj16QLt2oc/mlVfCUcfPfhYGBTRsGAYKvPMOLF8OQ4aEy8svXBjuXTJkSNl/H6lzdhcQuHutePTo0cOlhpoxw332bPcdO8L899+7L1/uvmSJ+8aNoWzkSPff/z6st3Wr+zffuP/lL+5ffOE+erR7+M2tRyYep58enps0Cc/Tp7uvXOnepo17797u55xTuu7LL4fnN95wnzXL/Z573CdNCmULFrhPneq+YYP72rXuJSXhtd55J5QlvPlmWP+OO1I/K1OmhGW33176nDB1anjPODt2uD/5ZPhclVdS4r5iRdU+u3UAUOQVfK9m5cs7Gw8FRC03ZUr4YjnkkPBFMXNmxeuuXOk+b577c8+VfrEdd1wIoHS/NJcsyf0Xtx7h8fOfV7ysZ8/Usttuc9+5071ePfe77nLv3DmUn3Za+CGSWO+880J4LFjg/uij7kuXhs/PqlVhm9/+1v2f/3T//PPwegmvveZ+6qnu994b3v/Xvw6fzWefdf/ww9IQW726dHr58tIfSOXt2OG+fbv7ww+7b9lSuf8XTz3l/swzldumHAWESEJRkfvzz6eW79zp/vrr7tu2lZY1blz6q/qxx0LZ22+XfsE8+aT7zTe7//Wv4Vdyrr9I9dg3Hg8/vPvlQ4aUTp90Uul0w4Zl1xs6tHT6/fer/F9idwGhPgiRylq4MAyb7d+/bPmyZeHeHO+/H87b+OSTMMw2MSLp4othzBi4+27YuBGOOSZ0Ol9xRThT/qyzwqXbd+6EP/8ZDjwQ2rcPQ3hPOy2MlnrlFfjNb8I5IEccAffem/39l5qpit/l6qQW2Rds2hRGPSWGzKZj6lQ46SS46y645ZYwUGDr1nDOx9Kl4YKM9eqFUWiNGoVraeXnh7L8/FC2bVsIqrfeCud63HFH2D5xN8S//CV0cL/6atlre82YEd5v2LCwTYcO4aZW8+end3/19u1DHePUrx8GP0h6+vcPA0OqQJ3UIlJ5H34Y+nOSbd3qvmiRe3Fx5V9vy5b4ppDXXnPfvNn9hhvcR4wou2zjRvfrrgt9Sl98EfoHli1z/+EH92nTQjPLffe5z50b1h81qmwzzBNPuL/3nvsDD7jPmRPKbr3V/dxz3f/whzA/aFDo13rwwVCHAQNKt+/bt3S6f3/3Qw8NfVtnnFFafuedpdPXX5+dZqoLLyw7n9xHUkmoiUlEJMaOHfGXaZk2LVyEMt0z4zduDOeTJM4bgjDcOPkkS4AXXoBjj4XjjgtNiWbhqG706HBi5hFHhKO4detCU2ZhYbi22YABMGtWOGp76KHwWlu3hvdt1myvTqZUE5OIiMTS1VxFRKTSFBAiIhJLASEiIrEUECIiEksBISIisRQQIiISSwEhIiKxFBAiIhKr1pwoZ2argAou7JKW1sB31VSdfUFd21/QPtcV2ufKae/ubeIW1JqA2FtmVlTR2YS1UV3bX9A+1xXa5+qjJiYREYmlgBARkVgKiFKjcl2BLKtr+wva57pC+1xN1AchIiKxdAQhIiKxFBAiIhKrzgeEmQ0wswVmtsjMhue6PnvDzJ42s5Vm9mlSWUszm2BmC6PnFlG5mdmD0X7PMbMTkra5LFp/oZldlot9SZeZHWZmE83sMzOba2Y3ReW1cr/NLN/MppvZJ9H+/iYq72hmH0X7NcbMGkbljaL5RdHyDkmvdUdUvsDMzsrNHqXPzOqb2cdm9no0X6v32cyWmNn/mtlsMyuKyrL7ua7oXqR14QHUB74AOgENgU+Arrmu117sT1/gBODTpLL7gOHR9HDgv6Lps4E3AQNOAj6KylsCi6PnFtF0i1zv2272+RDghGi6GfA50LW27ndU76bRdB7wUbQfLwEXR+WPAddG09cBj0XTFwNjoumu0ee9EdAx+n9QP9f7t4d9Hwo8D7wezdfqfQaWAK3LlWX1c13XjyB6AYvcfbG7lwAvAufluE5V5u6TgTXlis8DnoumnwPOTyr/swcfAgeY2SHAWcAEd1/j7muBCcCAzNe+atx9ubvPiqY3APOAttTS/Y7qvTGazYseDpwO/C0qL7+/ib/D34AzzMyi8hfdfau7fwksIvx/qJHMrB1wDvBkNG/U8n2uQFY/13U9INoCXyfNF0dltclB7r48mv4WOCiarmjf99m/SdSU0J3wq7rW7nfU1DIbWEn4D/8F8L27b49WSa77rv2Klq8DWrEP7W/kAWAYsDOab0Xt32cH3jazmWZ2VVSW1c91g6rUWvZN7u5mVivHNZtZU+AV4GZ3Xx9+MAa1bb/dfQdQYGYHAK8CXXJcpYwys3OBle4+08z65bo+WXSKuy8zswOBCWY2P3lhNj7Xdf0IYhlwWNJ8u6isNlkRHdbnq70AAAPNSURBVGoSPa+Myiva933ub2JmeYRwGO3uf4+Ka/1+u/v3wESgN6FJIfGDL7nuu/YrWt4cWM2+tb99gEFmtoTQDHw68Cdq9z7j7sui55WEHwK9yPLnuq4HxAygczQaoiGhQ2tsjutU3cYCiZELlwH/SCq/NBr9cBKwLjp0HQ/0N7MW0QiJ/lFZjRS1LT8FzHP3/05aVCv328zaREcOmFlj4F8J/S4TgQuj1crvb+LvcCHwrofey7HAxdGIn45AZ2B6dvaictz9Dndv5+4dCP9H33X3wdTifTazJmbWLDFN+Dx+SrY/17nuqc/1g9D7/zmhHfdXua7PXu7LC8ByYBuhrfEXhLbX/wEWAu8ALaN1DXgk2u//BQqTXudyQgfeIuD/5Hq/9rDPpxDaaucAs6PH2bV1v4FuwMfR/n4K3BWVdyJ82S0CXgYaReX50fyiaHmnpNf6VfR3WAAMzPW+pbn//SgdxVRr9znat0+ix9zEd1O2P9e61IaIiMSq601MIiJSAQWEiIjEUkCIiEgsBYSIiMRSQIiISCwFhEjEzDZGzx3M7N+r+bXvLDf/fnW+vkgmKCBEUnUAKhUQSWf0VqRMQLj7yZWsk0jWKSBEUt0LnBpdh/+W6OJ4fzCzGdG19q8GMLN+ZjbFzMYCn0Vlr0UXV5ubuMCamd0LNI5eb3RUljhasei1P42u/f+zpNd+z8z+ZmbzzWx0dNY4ZnavhftfzDGz+7P+15E6QxfrE0k1HLjN3c8FiL7o17l7TzNrBEwzs7ejdU8AjvVw+WiAy919TXQZjBlm9oq7DzezG9y9IOa9LgAKgOOB1tE2k6Nl3YFjgG+AaUAfM5sH/ATo4u6euOyGSCboCEJkz/oTrnMzm3Ap8VaE6/gATE8KB4AbzewT4EPCRdI6s3unAC+4+w53XwFMAnomvXaxu+8kXEKkA+HS1VuAp8zsAmDTXu+dSAUUECJ7ZsAv3b0genR098QRxA+7VgqXoj4T6O3uxxOumZS/F++7NWl6B9DAw/0NehFuhHMu8NZevL7IbikgRFJtINy+NGE8cG10WXHM7MjoCpvlNQfWuvsmM+tCuPVjwrbE9uVMAX4W9XO0Idw2tsIrjEb3vWju7uOAWwhNUyIZoT4IkVRzgB1RU9GzhHsPdABmRR3Fqyi91WOyt4Bron6CBYRmpoRRwBwzm+XhUtUJrxLu5/AJ4aq0w9z92yhg4jQD/mFm+YQjm6FV20WRPdPVXEVEJJaamEREJJYCQkREYikgREQklgJCRERiKSBERCSWAkJERGIpIEREJNb/B1wUCPme/JelAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ],
      "source": [
        "import warnings\n",
        "\n",
        "import torch\n",
        "import torch.nn as nn\n",
        "import torch.nn.functional as F\n",
        "from torch.utils.data import DataLoader\n",
        "from torchvision import datasets\n",
        "import torch.optim as optim\n",
        "from torchvision.transforms import transforms\n",
        "from torchvision.utils import save_image\n",
        "\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "warnings.filterwarnings(action='ignore')\n",
        "\n",
        "lr = 0.001\n",
        "batch_size = 100\n",
        "epochs = 10\n",
        "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
        "\n",
        "'''\n",
        "Step 1:\n",
        "'''\n",
        "\n",
        "# MNIST dataset\n",
        "dataset = datasets.MNIST(root='./mnist_data/',\n",
        "                               train=True, \n",
        "                               transform=transforms.ToTensor(),\n",
        "                               download=True)\n",
        "\n",
        "train_dataset, validation_dataset = torch.utils.data.random_split(dataset, [50000, 10000])\n",
        "\n",
        "test_dataset = datasets.MNIST(root='./mnist_data/',\n",
        "                              train=False, \n",
        "                              transform=transforms.ToTensor())\n",
        "\n",
        "# KMNIST dataset, only need test dataset\n",
        "anomaly_dataset = datasets.KMNIST(root='./kmnist_data/',\n",
        "                              train=False, \n",
        "                              transform=transforms.ToTensor(),\n",
        "                              download=True)\n",
        "\n",
        "# print(len(train_dataset))  # 50000\n",
        "# print(len(validation_dataset))  # 10000\n",
        "# print(len(test_dataset))  # 10000\n",
        "# print(len(anomaly_dataset))  # 10000\n",
        "\n",
        "'''\n",
        "Step 2: AutoEncoder\n",
        "'''\n",
        "# Define Encoder\n",
        "class Encoder(nn.Module):\n",
        "    def __init__(self):\n",
        "        super(Encoder, self).__init__()\n",
        "        self.fc1 = nn.Linear(784, 256)\n",
        "        self.fc2 = nn.Linear(256, 128)\n",
        "        self.fc3 = nn.Linear(128, 32)\n",
        "    def forward(self, x):\n",
        "        x = x.view(x.size(0), -1)\n",
        "        x = F.relu(self.fc1(x))\n",
        "        x = F.relu(self.fc2(x))\n",
        "        z = F.relu(self.fc3(x))\n",
        "        return z\n",
        "\n",
        "# Define Decoder\n",
        "class Decoder(nn.Module):\n",
        "    def __init__(self):\n",
        "        super(Decoder, self).__init__()\n",
        "        self.fc1 = nn.Linear(32, 128)\n",
        "        self.fc2 = nn.Linear(128, 256)\n",
        "        self.fc3 = nn.Linear(256, 784)\n",
        "    def forward(self, z):\n",
        "        z = F.relu(self.fc1(z))\n",
        "        z = F.relu(self.fc2(z))\n",
        "        x = F.sigmoid(self.fc3(z))  # to make output's pixels are 0~1\n",
        "        x = x.view(x.size(0), 1, 28, 28) \n",
        "        return x\n",
        "    \n",
        "'''\n",
        "Step 3: Instantiate model & define loss and optimizer\n",
        "'''\n",
        "enc = Encoder().to(device)\n",
        "dec = Decoder().to(device)\n",
        "loss_function = nn.MSELoss()\n",
        "optimizer = optim.Adam(list(enc.parameters()) + list(dec.parameters()), lr=lr)\n",
        "\n",
        "\n",
        "'''\n",
        "Step 4: Training\n",
        "'''\n",
        "train_loader = torch.utils.data.DataLoader(dataset=train_dataset, batch_size=batch_size, shuffle=True)\n",
        "\n",
        "\n",
        "train_loss_list = []\n",
        "\n",
        "import time\n",
        "start = time.time()\n",
        "for epoch in range(epochs) :\n",
        "    print(\"{}th epoch starting.\".format(epoch))\n",
        "    enc.train()\n",
        "    dec.train()\n",
        "    for batch, (images, _) in enumerate(train_loader) :\n",
        "        images = images.to(device)\n",
        "        z = enc(images)\n",
        "        reconstructed_images = dec(z)\n",
        "        \n",
        "        optimizer.zero_grad()\n",
        "        train_loss = loss_function(images, reconstructed_images)\n",
        "        train_loss.backward()\n",
        "        train_loss_list.append(train_loss.item())\n",
        "\n",
        "        optimizer.step()\n",
        "\n",
        "        print(f\"[Epoch {epoch:3d}] Processing batch #{batch:3d} reconstruction loss: {train_loss.item():.6f}\", end='\\r')\n",
        "end = time.time()\n",
        "print(\"Time ellapsed in training is: {}\".format(end - start))\n",
        "\n",
        "# plotting train loss\n",
        "plt.plot(range(1,len(train_loss_list)+1), train_loss_list, 'r', label='Training loss')\n",
        "plt.title('Training loss')\n",
        "plt.xlabel('Iterations')\n",
        "plt.ylabel('Loss')\n",
        "plt.legend()\n",
        "plt.savefig('loss.png')\n",
        "\n",
        "enc.eval()\n",
        "dec.eval()\n",
        "\n",
        "'''\n",
        "Step 5: Calculate standard deviation by using validation set\n",
        "'''\n",
        "validation_loader = torch.utils.data.DataLoader(dataset=validation_dataset, batch_size=batch_size)\n",
        "\n",
        "def score_function(true_image, reconstructed_image):\n",
        "\n",
        "    return torch.linalg.norm(true_image - reconstructed_image) **2\n",
        "\n",
        "\n",
        "score_list = []\n",
        "for images, _ in validation_loader:\n",
        "\n",
        "    enc.eval()\n",
        "    dec.eval()\n",
        "\n",
        "    for idx in range(batch_size):\n",
        "      \n",
        "      image = images[idx].to(device)\n",
        "      z = enc(image)\n",
        "      reconstructed_image = dec(z)\n",
        "      score = score_function(image, reconstructed_image)\n",
        "\n",
        "      score_list.append(score.item())\n",
        "    \n",
        "print(\"---------------------\")\n",
        "\n",
        "\n",
        "mean = np.mean(score_list)\n",
        "std = np.std(score_list)\n",
        "print(mean, std)\n",
        "print(\"---------------------\")\n",
        "\n",
        "threshold = mean + 3 * std\n",
        "print(\"threshold: \", threshold)\n",
        "\n",
        "\n",
        "'''\n",
        "Step 6: Anomaly detection (mnist)\n",
        "'''\n",
        "test_loader = torch.utils.data.DataLoader(dataset=test_dataset, batch_size=batch_size)\n",
        "\n",
        "cnt_all_MNIST, cnt_error_MNIST = 0, 0\n",
        "for images, _ in test_loader:\n",
        "    for idx in range(batch_size):\n",
        "      image = images[idx].to(device)\n",
        "      z = enc(image)\n",
        "      reconstructed_image = dec(z)\n",
        "      score = score_function(image.unsqueeze(0), reconstructed_image)\n",
        "\n",
        "      score = score.item()\n",
        "      cnt_all_MNIST += 1\n",
        "      if score > threshold :\n",
        "          cnt_error_MNIST += 1\n",
        "\n",
        "print(f\"Type 1 Error: {100 * cnt_error_MNIST/cnt_all_MNIST:.4f} %\")\n",
        "\n",
        "'''\n",
        "Step 7: Anomaly detection (kmnist)\n",
        "'''\n",
        "anomaly_loader = torch.utils.data.DataLoader(dataset=anomaly_dataset, batch_size=batch_size)\n",
        "\n",
        "cnt_all_KMNIST, cnt_error_KMNIST = 0, 0\n",
        "for images, _ in anomaly_loader:\n",
        "    for idx in range(batch_size):\n",
        "      image = images[idx].to(device)\n",
        "      z = enc(image)\n",
        "      reconstructed_image = dec(z)\n",
        "      score = score_function(image.unsqueeze(0), reconstructed_image)\n",
        "\n",
        "      score = score.item()\n",
        "      cnt_all_KMNIST += 1\n",
        "      if score < threshold:\n",
        "          cnt_error_KMNIST += 1\n",
        "\n",
        "print(f\"Type 2 Error: {100 * cnt_error_KMNIST / cnt_all_KMNIST:.4f} % \")\n"
      ]
    }
  ]
}