CASSL_Validation.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "CASSL_Validation.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uzxWXO1V2oh6"
      },
      "source": [
        "# Importing libraries"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sXG7u5it2WrK"
      },
      "source": [
        "import math\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import random\n",
        "from sklearn.cluster import KMeans\n",
        "from sklearn.preprocessing import LabelEncoder\n",
        "from sklearn.metrics.cluster import adjusted_rand_score, adjusted_mutual_info_score, fowlkes_mallows_score\n",
        "from sklearn.metrics import silhouette_score, davies_bouldin_score, calinski_harabasz_score\n",
        "import seaborn as sns\n",
        "from sklearn.manifold import TSNE\n",
        "from sklearn.decomposition import NMF\n",
        "import tensorflow as tf\n",
        "from sklearn.metrics.pairwise import euclidean_distances"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "azZAvGkqCVP3",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "2c1991fd-a4b7-4af2-8741-7473c0e39458"
      },
      "source": [
        "!pip install pytictoc"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Collecting pytictoc\n",
            "  Downloading pytictoc-1.5.2-py2.py3-none-any.whl (4.0 kB)\n",
            "Installing collected packages: pytictoc\n",
            "Successfully installed pytictoc-1.5.2\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "I_g0ZJR-jOi_",
        "outputId": "aa4a065e-c828-461e-8480-d232c7ebfcdf"
      },
      "source": [
        "!pip install kneed"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Collecting kneed\n",
            "  Downloading kneed-0.7.0-py2.py3-none-any.whl (9.4 kB)\n",
            "Requirement already satisfied: numpy>=1.14.2 in /usr/local/lib/python3.7/dist-packages (from kneed) (1.19.5)\n",
            "Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (from kneed) (1.4.1)\n",
            "Requirement already satisfied: matplotlib in /usr/local/lib/python3.7/dist-packages (from kneed) (3.2.2)\n",
            "Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->kneed) (1.3.2)\n",
            "Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->kneed) (2.8.2)\n",
            "Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.7/dist-packages (from matplotlib->kneed) (0.11.0)\n",
            "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->kneed) (2.4.7)\n",
            "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.7/dist-packages (from python-dateutil>=2.1->matplotlib->kneed) (1.15.0)\n",
            "Installing collected packages: kneed\n",
            "Successfully installed kneed-0.7.0\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "i9ks8QJYCe3j"
      },
      "source": [
        "from pytictoc import TicToc\n",
        "t = TicToc()"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CEhM9idP3AuQ"
      },
      "source": [
        "# Importing scRNA-seq data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "bRiS4oDCi7aX",
        "outputId": "af27eb5e-1e47-4ebe-d75c-89cf5010daff"
      },
      "source": [
        "#Mounting google drive\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "vnkqz3ZpjErq",
        "outputId": "772ccbe8-b48d-40cf-db2d-60f98b1fa59f"
      },
      "source": [
        "ls \"/content/drive/My Drive/scRNAseq Data/Preprocessed_Data_Annotations/\""
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "AMB_annotations.csv\n",
            "AMB_preprocessed_data_scprep_percentile_5_95.csv\n",
            "Baron_Human_annotations.csv\n",
            "Baron_Human_preprocessed_data_scprep_percentile_5_95.csv\n",
            "Baron_Mouse_annotations.csv\n",
            "Baron_Mouse_preprocessed_data_scprep_percentile_5_95.csv\n",
            "GSE57249_fpkm_labels.csv\n",
            "GSE57249_fpkm_preprocessed.csv\n",
            "missingAnnotations_Test.csv\n",
            "PBMC_3k_annotations.csv\n",
            "PBMC3K_preprocessed_data_minmax_normalized.csv\n",
            "PBMC_3k_preprocessed_data_scanpy.csv\n",
            "Preimplantation_annotations.csv\n",
            "Preimplantation_preprocessed_data_scprep_percentile_20_90.csv\n",
            "Segerstolpe_annotations.csv\n",
            "Segerstolpe_preprocessed_data_scprep_percentile_10_95.csv\n",
            "Usoskin_annotations.csv\n",
            "Usoskin_preprocessed_data_scprep_percentile_5_95.csv\n",
            "Zeisel_annotations.csv\n",
            "Zeisel_preprocessed_data_scprep_percentile_5_95.csv\n",
            "Zheng_PBMC_DownSampled_annotations.csv\n",
            "Zheng_PBMC_DownSampled_preprocessed_data_scprep_percentile_5_90.csv\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "A4nNDNWpi4AY"
      },
      "source": [
        "data = pd.read_csv(\"/content/drive/My Drive/scRNAseq Data/Preprocessed_Data_Annotations/Segerstolpe_preprocessed_data_scprep_percentile_10_95.csv\", header=0, index_col=0, sep=\",\")"
      ],
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "T6Dzw-jBcMoH"
      },
      "source": [
        "labels = pd.read_csv(\"/content/drive/My Drive/scRNAseq Data/Preprocessed_Data_Annotations/Segerstolpe_annotations.csv\", header=0, index_col=None, sep=\",\")"
      ],
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EPrGapzDjtH_"
      },
      "source": [
        "data.reset_index(inplace=True, drop=True)"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 253
        },
        "id": "lp27GOLalq25",
        "outputId": "71ca82e0-7693-4d7f-cdcc-75bc9fdfa81b"
      },
      "source": [
        "data.head()"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>CLIC4</th>\n",
              "      <th>TGFBR3</th>\n",
              "      <th>PRUNE</th>\n",
              "      <th>EYA3</th>\n",
              "      <th>STMN1</th>\n",
              "      <th>PTPRF</th>\n",
              "      <th>RAB3B</th>\n",
              "      <th>LEPR</th>\n",
              "      <th>DDAH1</th>\n",
              "      <th>GNAI3</th>\n",
              "      <th>RHOC</th>\n",
              "      <th>NOTCH2</th>\n",
              "      <th>SOAT1</th>\n",
              "      <th>EPRS</th>\n",
              "      <th>C1orf127</th>\n",
              "      <th>AGTRAP</th>\n",
              "      <th>SH3BGRL3</th>\n",
              "      <th>CMPK1</th>\n",
              "      <th>TTC39A</th>\n",
              "      <th>PPAP2B</th>\n",
              "      <th>PGM1</th>\n",
              "      <th>JAK1</th>\n",
              "      <th>WLS</th>\n",
              "      <th>GNG12</th>\n",
              "      <th>PRKACB</th>\n",
              "      <th>VCAM1</th>\n",
              "      <th>VAV3</th>\n",
              "      <th>KIAA1324</th>\n",
              "      <th>TSPAN2</th>\n",
              "      <th>ATP1A1</th>\n",
              "      <th>VTCN1</th>\n",
              "      <th>CTSS</th>\n",
              "      <th>IL6R</th>\n",
              "      <th>LMNA</th>\n",
              "      <th>SYT11</th>\n",
              "      <th>IFI16</th>\n",
              "      <th>RGS5</th>\n",
              "      <th>RXRG</th>\n",
              "      <th>ATP1B1</th>\n",
              "      <th>PAPPA2</th>\n",
              "      <th>...</th>\n",
              "      <th>CECR1</th>\n",
              "      <th>ADRBK2</th>\n",
              "      <th>SEZ6L</th>\n",
              "      <th>HMOX1</th>\n",
              "      <th>RPL3</th>\n",
              "      <th>PACSIN2</th>\n",
              "      <th>PLXNB2</th>\n",
              "      <th>SLC25A1</th>\n",
              "      <th>CLDN5</th>\n",
              "      <th>GGT1</th>\n",
              "      <th>TPST2</th>\n",
              "      <th>XBP1</th>\n",
              "      <th>LIF</th>\n",
              "      <th>SELM</th>\n",
              "      <th>PIK3IP1</th>\n",
              "      <th>PISD</th>\n",
              "      <th>RTCB</th>\n",
              "      <th>TIMP3</th>\n",
              "      <th>RAC2</th>\n",
              "      <th>CDC42EP1</th>\n",
              "      <th>LGALS2</th>\n",
              "      <th>LGALS1</th>\n",
              "      <th>H1F0</th>\n",
              "      <th>KDELR3</th>\n",
              "      <th>ATF4</th>\n",
              "      <th>TSPO</th>\n",
              "      <th>PIM3</th>\n",
              "      <th>ERCC_30000:mix1_7500:mix2</th>\n",
              "      <th>ERCC_7500:mix1_1875:mix2</th>\n",
              "      <th>ERCC_1875:mix1_468.75:mix2</th>\n",
              "      <th>ERCC_937.5:mix1_234.375:mix2</th>\n",
              "      <th>ERCC_468.75:mix1_117.1875:mix2</th>\n",
              "      <th>ERCC_15000:mix1_15000:mix2</th>\n",
              "      <th>ERCC_3750:mix1_3750:mix2</th>\n",
              "      <th>ERCC_937.5:mix1_937.5:mix2</th>\n",
              "      <th>ERCC_15000:mix1_22500:mix2</th>\n",
              "      <th>ERCC_3750:mix1_5625:mix2</th>\n",
              "      <th>ERCC_15000:mix1_30000:mix2</th>\n",
              "      <th>ERCC_937.5:mix1_1875:mix2</th>\n",
              "      <th>ERCC_234.375:mix1_468.75:mix2</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>0.185423</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.229955</td>\n",
              "      <td>1.637609</td>\n",
              "      <td>2.172956</td>\n",
              "      <td>1.048909</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.502127</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.03239</td>\n",
              "      <td>1.201677</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.00000</td>\n",
              "      <td>1.521521</td>\n",
              "      <td>0.889256</td>\n",
              "      <td>2.232785</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.718139</td>\n",
              "      <td>1.637609</td>\n",
              "      <td>0.355058</td>\n",
              "      <td>1.289099</td>\n",
              "      <td>0.902052</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.483096</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.918028</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.192101</td>\n",
              "      <td>1.794550</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.914669</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>2.914147</td>\n",
              "      <td>0.0</td>\n",
              "      <td>...</td>\n",
              "      <td>1.884878</td>\n",
              "      <td>0.975307</td>\n",
              "      <td>1.054358</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.364897</td>\n",
              "      <td>0.466637</td>\n",
              "      <td>1.689281</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.547660</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.793933</td>\n",
              "      <td>1.215898</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.735640</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.525283</td>\n",
              "      <td>0.302794</td>\n",
              "      <td>0.908382</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.107054</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.151397</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.670667</td>\n",
              "      <td>2.068522</td>\n",
              "      <td>1.491443</td>\n",
              "      <td>1.184334</td>\n",
              "      <td>1.481912</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.219695</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.11903</td>\n",
              "      <td>1.254060</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.00000</td>\n",
              "      <td>1.135475</td>\n",
              "      <td>1.579112</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.119030</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.196237</td>\n",
              "      <td>0.291563</td>\n",
              "      <td>1.166252</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>3.041678</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>3.255424</td>\n",
              "      <td>0.0</td>\n",
              "      <td>...</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>2.359679</td>\n",
              "      <td>0.461002</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.023511</td>\n",
              "      <td>0.119030</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.890740</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.383005</td>\n",
              "      <td>1.330797</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.201252</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.858338</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.995877</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.119030</td>\n",
              "      <td>0.119030</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.626815</td>\n",
              "      <td>1.732498</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>2.079618</td>\n",
              "      <td>1.596779</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>2.342371</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.00000</td>\n",
              "      <td>2.559642</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.00000</td>\n",
              "      <td>2.349250</td>\n",
              "      <td>1.136217</td>\n",
              "      <td>1.943230</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>4.096688</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.616865</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.032527</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>2.236645</td>\n",
              "      <td>0.0</td>\n",
              "      <td>...</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.272434</td>\n",
              "      <td>2.307667</td>\n",
              "      <td>2.515123</td>\n",
              "      <td>1.959768</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.270329</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.538955</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.901255</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.179652</td>\n",
              "      <td>0.179652</td>\n",
              "      <td>0.254066</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>0.161254</td>\n",
              "      <td>0.161254</td>\n",
              "      <td>1.206712</td>\n",
              "      <td>1.259431</td>\n",
              "      <td>1.045043</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>2.127081</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.926332</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.00000</td>\n",
              "      <td>0.161254</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.00000</td>\n",
              "      <td>1.007029</td>\n",
              "      <td>1.339473</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>2.524018</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.319918</td>\n",
              "      <td>1.007029</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.766445</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.873598</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.822235</td>\n",
              "      <td>0.702888</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.788389</td>\n",
              "      <td>2.731810</td>\n",
              "      <td>0.0</td>\n",
              "      <td>...</td>\n",
              "      <td>2.816173</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.817236</td>\n",
              "      <td>0.161254</td>\n",
              "      <td>1.195889</td>\n",
              "      <td>1.007029</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.387156</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>2.026927</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.529786</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.994034</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.161254</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.705230</td>\n",
              "      <td>1.665843</td>\n",
              "      <td>2.592007</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.931586</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.45298</td>\n",
              "      <td>1.693510</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.29344</td>\n",
              "      <td>1.457542</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>2.038605</td>\n",
              "      <td>0.115229</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.609733</td>\n",
              "      <td>1.600809</td>\n",
              "      <td>0.907312</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>3.351554</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>1.475648</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.230458</td>\n",
              "      <td>0.877556</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>2.753462</td>\n",
              "      <td>0.0</td>\n",
              "      <td>...</td>\n",
              "      <td>1.588319</td>\n",
              "      <td>2.012384</td>\n",
              "      <td>0.914600</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.896905</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.115229</td>\n",
              "      <td>0.415463</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.115229</td>\n",
              "      <td>0.620526</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.397075</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>2.411560</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.115229</td>\n",
              "      <td>0.0</td>\n",
              "      <td>1.406547</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.162958</td>\n",
              "      <td>0.000000</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0.115229</td>\n",
              "      <td>0.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "<p>5 rows × 1830 columns</p>\n",
              "</div>"
            ],
            "text/plain": [
              "      CLIC4    TGFBR3  ...  ERCC_937.5:mix1_1875:mix2  ERCC_234.375:mix1_468.75:mix2\n",
              "0  0.185423  0.000000  ...                   0.151397                            0.0\n",
              "1  0.000000  0.000000  ...                   0.000000                            0.0\n",
              "2  0.000000  0.000000  ...                   0.000000                            0.0\n",
              "3  0.161254  0.161254  ...                   0.000000                            0.0\n",
              "4  0.000000  0.000000  ...                   0.115229                            0.0\n",
              "\n",
              "[5 rows x 1830 columns]"
            ]
          },
          "metadata": {},
          "execution_count": 10
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eM5GWDv_ScjG"
      },
      "source": [
        "# Scaling data (if required)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NFsmycwG3-mO"
      },
      "source": [
        "%%script false\n",
        "#only for scanpy pre-processed data (PBMC - 10x Genomics)\n",
        "from sklearn.preprocessing import MinMaxScaler\n",
        "scaler = MinMaxScaler()\n",
        "data = pd.DataFrame(scaler.fit_transform(data))\n",
        "data.head()"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PQKgW5eASrlA"
      },
      "source": [
        "# Checking loaded data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "pvcCrZtejwjW",
        "outputId": "86032044-1e0f-4b4b-ed44-6295931e902c"
      },
      "source": [
        "data.shape"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(1630, 1830)"
            ]
          },
          "metadata": {},
          "execution_count": 12
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "it-bUYvxj7hd",
        "outputId": "e786bcce-5d0b-49fc-b5a8-d1955edfee22"
      },
      "source": [
        "labels.shape"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(1630, 2)"
            ]
          },
          "metadata": {},
          "execution_count": 13
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "xzLNQ2-yM5a9",
        "outputId": "d1be3e64-181d-42c4-8006-5fa4f32cf137"
      },
      "source": [
        "labels.nunique()"
      ],
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Unnamed: 0    1630\n",
              "x               13\n",
              "dtype: int64"
            ]
          },
          "metadata": {},
          "execution_count": 14
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "yOfWyAQI1DJs",
        "outputId": "366b6af7-52bf-424a-f477-57b0748e6866"
      },
      "source": [
        "labels = labels['x']\n",
        "labels.head()"
      ],
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0    delta\n",
              "1    alpha\n",
              "2    delta\n",
              "3    gamma\n",
              "4    alpha\n",
              "Name: x, dtype: object"
            ]
          },
          "metadata": {},
          "execution_count": 17
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Nmxn0dEhwf0s"
      },
      "source": [
        "# Finding the optimal dimensions for embedding using reconstruction error"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "v7PqelPewzCd",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "f958c067-c26e-4798-991e-3eb276ba1438"
      },
      "source": [
        "d = range(2,16)\n",
        "k=[]\n",
        "error=[]\n",
        "# try each value of k\n",
        "for i in d:\n",
        "    print(\"Applying NMF for k=%d ...\" % i)\n",
        "    # run NMF\n",
        "    model = NMF(init=\"nndsvd\", n_components=i) \n",
        "    W = model.fit_transform(data)\n",
        "    H = model.components_    \n",
        "    k.append(i)\n",
        "    error.append(model.reconstruction_err_)"
      ],
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Applying NMF for k=2 ...\n",
            "Applying NMF for k=3 ...\n",
            "Applying NMF for k=4 ...\n",
            "Applying NMF for k=5 ...\n",
            "Applying NMF for k=6 ...\n",
            "Applying NMF for k=7 ...\n",
            "Applying NMF for k=8 ...\n",
            "Applying NMF for k=9 ...\n",
            "Applying NMF for k=10 ...\n",
            "Applying NMF for k=11 ...\n",
            "Applying NMF for k=12 ...\n",
            "Applying NMF for k=13 ...\n",
            "Applying NMF for k=14 ...\n",
            "Applying NMF for k=15 ...\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "m-waFpIxxIwa",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 288
        },
        "outputId": "ce0a259e-9437-4924-e344-a0277d66bdcb"
      },
      "source": [
        "# Plotting a bar graph to compare the results \n",
        "plt.plot(k, error) \n",
        "plt.xlabel('Number of dimensions', fontsize = 20) \n",
        "plt.ylabel('Reconstruction Error', fontsize = 20) \n",
        "plt.show() "
      ],
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEPCAYAAABlZDIgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5xU1fnH8c93l6UtfZcmbemIiIobLCiCxhKjYk8swU4SSzRVjUkwahLTE8WGCUH92WswMXZBQJAioAgiS196r9Kf3x/3LgzLbBl2Zu/O7vN+veY1O+eeuffZ2dl55txzzzkyM5xzzrmKyog6AOecc9WDJxTnnHNJ4QnFOedcUnhCcc45lxSeUJxzziVFragDiEpubq7l5eVFHYZzzqWVqVOnrjGz5vG21diEkpeXx5QpU6IOwznn0oqkRSVt81NezjnnksITinPOuaTwhOKccy4pPKE455xLCk8ozjnnksITinPOuaQod0KR9CtJ30llMM4559JXIi2UXwBHpiqQdLHxq1385e05FKzaEnUozjlXpSSSUJYCjVIVSLrYvWcvw8fO57Ex86IOxTnnqpREEsqrwNcl1UtVMOkgp0EdvpXfjtemL2X5xq+iDsc556qMRBLKUGA98JqkXimKJy1cf3In9hr8c+yCqENxzrkqI5GEMgNoDXwdmCFpq6QFkuYXu1X7c0HtmtXnvKMO45lJi9mwbWfU4TjnXJWQSELJAHYBi8PbqrBcxW414lLk757SiW079/DkhBLnSXPOuRql3LMNm1leCuNIOz1aNeLUHi0Y+dFCbji5E/VqZ0YdknPORapGtCZS5fsDOrNu605emLIk6lCccy5yh5xQJDWU1E5Sjb2U+Gt5zcjv0JThH85n1569UYfjnHORSiihSKol6Q5JBcAGYCGwXlJBWF7jFuz6/oDOLN3wFf/5dFnUoTjnXKQSmXqlNvA28BsgD1gCTArv88Lyd8N6NcbA7i3o1rIBj4yex969FnU4zjkXmURaKD8CBgD/BQ43szwzOyHsrO8OvA6cHNarMTIyxPdO6cyXK7fwwZxVZT/BOeeqqUQSyuXATOB8M5sbu8HM5gEXAp8DVyQvvPRw7lGH0aZJPR716VicczVYIgmlC/A/M4vb+xyW/w/onIzA0klWZgY3nNyRyQvXM3nhuqjDcc65SCSSUHYCDcqok00w+LHGufRr7WhaP4tHR3srxTlXMyWSUD4FLpbUPN5GSbnAxQRTtNQ49WvX4uoTO/LeF6uYs2Jz1OE451ylSyShDAOaA5MkXSepk6R6kjpKugb4ONw+LBWBpoPBJ3Sgfu1Mn9reOVcjlTuhmNkLwP1AB2A4MBfYAhQA/wA6An8M65WLpBGSVkmaWaz8FklfSPpc0h9iyu8Mx7zMkXRmTPlZYVmBpDvKe/xka5pdm8v6tuffM5ZRuH5bVGE451wkEhrYaGY/B04ERgDTgPnh/Qign5kl+mE+EjgrtkDSQGAQcJSZHQH8KSzvCXwbOCJ8zsOSMiVlAg8B3wB6ApeFdSNx/ckdyRD8w6e2d87VMOUe2S5pMLDSzN4CJibj4Gb2oaS8YsXfB+43sx1hnaLBHYOA58LyBeFo/b7htgIzmx/G+VxYd1YyYkxU68b1GHR0G56bvJhbTu1CToM6UYThnHOVLpEWygiKtSZSpBtwsqSPJY2R9LWwvA3BqPwihWFZSeUHkTRE0hRJU1avXp2C0APfO6UT23ft5Qmf2t45V4MkklBWJFj/UNUCmgHHAz8FXpCkZOzYzIabWb6Z5TdvHvditaTo0qIhZ/RsyRMfLWTrjt0pO45zzlUliSSIN4GBklKdVAqBVywwCdgL5AJLgXYx9dqGZSWVR+p7Azqz8atdPDtpcdShOOdcpUgkOdwFNAT+GY45SZXXgIEAkroBtYE1wCjg25LqSOoIdCWYnHIy0DW8fLk2Qcf9qBTGVy592jfluI7N+MfYBezc7VPbO+eqv0Smm38W2AgMJvhgX0hwGqz4FLtmZqeVZ4eSniWYcDJXUiEwlKCvZkR4KfFO4CozM+BzSS8QdLbvBm4ysz3hfm4G3gIygRFm9nkCv1fKfH9AZ67+12Rem76US/Pblf0E55xLYwo+q8tRUSrv12wzsyq/Hm5+fr5NmTIlpccwM85+YBw7d+/hnR+eQkZGUrqCnHMuMpKmmll+vG2JDGzMKOetyieTyiKJ7w/ozLzVW3ln9sqow3HOuZRKZIGtwbGj0135nN2rFe2b1efh0fMob2vQOefSUVUch1Kt1MrM4Ib+nZixZAMT5/vU9s656qsqjkOpdi45ti25DWr7AlzOuWqtKo5DqXbqZmVyTb+OjPlyNZ8v2xh1OM45lxJVcRxKtXTl8R1oUKcWj46ZH3UozjmXEokklNhxKEskzZb0gaT3i93eS02o6a1xvSyuOK49//10GYvWbo06HOecS7pEEsoAoDcgoA7QHTglLC9+c3Fce1JHamVk8PhYb6U456ofH4dSiVo2qstFx7bhhSmFrN68I+pwnHMuqbyDvZIN6d+ZXXv28q/xvgCXc6568YRSyTrmZnN2r9Y8NWERm7bvijoc55xLmlITStjJPrhY2XGSflBC/VsleQdBGb53Smc279jNMx/71PbOueqjrBbKACCvWNlZwF9LqN8E6FCxkKq/I9s25qQuufxz3AK279oTdTjOOZcUfsorIt8f0JnVm3fw6rTI1wJzzrmk8IQSkRM759C7bWMeGzOPPXt90kjnXPrzhBIRSXzvlM4sXLuNN2euiDoc55yrME8oETrziFZ0zM3m0TE+tb1zLv15QolQZob4bv9OfLZ0I+ML1kYdjnPOVUh51pQfIB2wdO0AAEm/JJiG5aBtrvwu6NOGv7zzJY+MKeCkrj7npnMufZUroRA/Ufy62GMjSDB+7iYBdWplcv3JHfntG1/waeEGerdtEnVIzjl3SMpKKMWThkuBy/q258H3C3h0zDwevuLYqMNxzrlDUmpCMTNPKJWgYd0sBp/QgYdHz2P+6i10at4g6pCccy5hkXbKSxohaZWkmTFld0taKml6eDs7ZtudkgokzZF0Zkz5WWFZgaQ7Kvv3SIarT+xI7cwMhn/oM9c459JT1Fd5jSSYyqW4v5rZ0eHtDQBJPYFvA0eEz3lYUqakTOAh4BtAT+CysG5aad6wDpfmt+PlTwpZsXF71OE451zCIk0oZvYhsK6c1QcBz5nZDjNbABQAfcNbgZnNN7OdwHNh3bRzw8md2LPX+Oc4b6U459JP1C2Uktws6dPwlFjTsKwNsCSmTmFYVlL5QSQNkTRF0pTVq1enIu4KaZ9Tn0FHt+GJCYuYv3pL1OE451xCqmJCeQToDBwNLAf+nKwdm9lwM8s3s/zmzZsna7dJdec3elC3VgZ3vPIZe32OL+dcGqlyCcXMVprZHjPbCzxOcEoLYCnQLqZq27CspPK01KJRXe765uFMWrCOZyf7einOufRR5RKKpNYxDy8Aiq4AGwV8W1IdSR2BrsAkYDLQVVJHSbUJOu5HVWbMyXZpfjtO7JzD/W984R30zrm0EfVlw88CE4DukgolXQf8QdJnkj4FBgI/BDCzz4EXgFnAm8BNYUtmN3Az8BYwG3ghrJu2JPG7C49k1969/OK1mT5xpHMuLSjRDytJLYB8oCmQGa+OmT1Z8dBSKz8/36ZMmRJ1GKUa/uE8fvvGFwy7/BjO6X1Y1OE45xySpppZfrxt5ZnLq2gnWcCjwGBKbtkUzeVV5RNKOri2X0den7Gcu0d9Tr/OuTTNrh11SM45V6JyJxTgXuAaYB7wNMGlurtTEZQL1MrM4PcX9ea8YeO477+z+fOlR0UdknPOlSiRhHI58CVwjJl9laJ4XDE9D2vEd0/pxEMfzGPQ0YfRv1vVvNzZOecS6ZRvAbzhyaTy3XJqVzo1z+bnr37G1h3eKHTOVU2JJJTFQKNUBeJKVjcrk/sv7E3h+q/489tfRh2Oc87FlUhCGQl8Q1LjFMXiStG3YzOuPL49//poAdMWr486HOecO0giCeV+YBzwrqSBkry1UsluP6sHrRrV5Y6XP2Pn7r1Rh+OccwdIJKHsAi4GjgXeBdZL2hPn5if5U6Rh3SzuO78Xc1Zu5pHR86IOxznnDpDIVV5j8fXiI3fa4S0596jDGPbBXM4+shVdWzaMOiTnnAMSSChmNiCFcbgEDD23J2Pnrub2lz/lxe+dSGaGog7JOeeq3uSQrmy5Derwq3N68sniDTw1YWHU4TjnHHCICUVSlqQjJZ0sqXc4LYurRBcc04b+3Zrzh7fmULh+W9ThOOdcYglFUiNJjwIbgOnAaGAasEHSo5KaJD9EF48kfntBLwDuetVnJHbORa/cCSW8THg8MIRgDq+xBNPJjyW4AmwIMM4vJ648bZvW56dndmfMl6t5bXrarinmnKsmEmmh3AkcQbBEbwczG2Bml4Wd9R2Ah4CeYT1XSQafkMcx7Ztwz+uzWLtlR9ThOOdqsEQSyoXARDO7ycw2xG4ws41mdgvBYlkXJTNAV7rMDPH7i3qzZcdufv36rKjDcc7VYIkklA4EfSalGcOB67u7StCtZUNuGtiFUTOW8f4XK6MOxzlXQyWSULYSzDhcmuaAX3IUgRsHdKFbywbc9epMNm/fFXU4zrkaKJGEMhm4RFLXeBsldQYuDeu5Sla7Vgb3X9SbFZu284c350QdjnOuBkokofwRaABMlnSvpFMlHR5OFPlrgkTSAPhTKgJ1ZevTvilXn5jHUxMXMXnhuqjDcc7VMOVOKGb2HnAjUBf4OfAOMJNgoshfAtnAzWb2bgridOX0kzO606ZJPW5/+VO279oTdTjOuRokoYGNZvYY0A34FfAq8H54/0ugm5k9kvQIXUKy69TidxceyfzVWxn2fkHU4TjnapCEp14xs8Vm9hszu9jMTg/vf2NmixLdl6QRklZJmhln248lmaTc8LEkPSCpQNKnkvrE1L1K0tzwdlWicVQ3/bs158I+bXh0zDxmL98UdTjOuRoi6skhRwJnFS+U1A44g2DZ4SLfALqGtyEEAyyR1AwYChwH9AWGSmqa0qjTwC+/2ZPG9bK4/eVP2b3HF+NyzqVeidPXS+of/jjJzLbHPC6TmX1Y3nqS8uJs+ivwM+DfMWWDgCctmLRqoqQmkloDA4B3zGxdGPc7BEnq2fLGWx01za7N3ecdwS3PTuNf4xdyQ/9OUYfknKvmSlsPZTTBglqHA1/GPC6PzEMNSNIgYKmZzZAOWOejDbAk5nFhWFZSebx9DyFo3dC+fftDDTFtnNO7Nf+evpQ/vzOHM45oSYec7KhDcs5VY6UllHsIEsiaYo9TRlJ9givIzkjF/s1sODAcID8/v9pPzyuJe8/vxel/+ZA7X/mMp68/jmJJ2jnnkqbEhGJmd5f2OEU6Ax2BotZJW+ATSX2BpRw4rUvbsGwpwWmv2PLRlRBrWmjduB53fKMHv3htJi9OKeTSr/nMOM651Ehk+vr2ZU1NL6mhpEM+l2Rmn5lZCzPLM7M8gtNXfcxsBTAKGBxe7XU8sNHMlgNvAWdIahp2xp8RlrnQ5X3b07djM+777yxWbdoedTjOuWoqkau8FgC3llHnB2G9cpH0LMEMxd0lFUq6rpTqbwDzgQLgcYJBloSd8fcSjNSfDNxT1EHvAhkZ4v4Lj2T77r0MHfV51OE456qp0vpQilN4Sxozu6yM7XkxPxtwUwn1RgAjkhlbddOpeQNuPa0rf3xrDv+evpRBR8e9bsE55w5ZssehtCKYldhVQd/t34n8Dk35+SufMX/1lqjDcc5VM6W2UCQNLlZ0dJwyCC4Tbg9cCXyWpNhcktXKzOCBy47h7AfGctMz03j1xhOpm3XIV3g759wBFJxJKmGjtJfyXSpcdCpsG3Chmb2dhNhSKj8/36ZMmRJ1GJF4/4uVXDtyClce3577zj8y6nCcc2lE0lQzy4+3raw+lGuK9kHQR/EaB45eL7IHWAtMKL48sKt6Tu3RkiH9OzH8w/kc3ymHc3ofFnVIzrlqoNSEYmZPFP0cTrr4mpk9mfKoXMr99MzuTF64jjte/owj2zT2UfTOuQpLZD2UgZ5Mqo+szAwevOwYMjPETc98wo7dvnaKc65iEhnY2EfSryS1LGF7q3D70ckLz6VS26b1+ePFvZm5dBO/e+OLqMNxzqW5RC4b/glwPbCqhO0rgeuAH1U0KFd5zjiiFdf268jIjxby5szlUYfjnEtjiSSUE4APrITLwsLy94F+yQjMVZ47vtGDo9o25qcvfcqSdduiDsc5l6YSSSitCObWKs0yoPWhh+OiULtWBsMuDxbAvPmZT9i52xfkcs4lLpGEsg1oXkad5sCOQw/HRaVds6A/ZUbhRn7/pvenOOcSl0hCmQ4MktQg3sZwJuJBYT2Xhs7q1ZqrTujAP8ct4J1ZK6MOxzmXZhJJKMMJWiDvSOodu0HSUcDbQG5Yz6Wpn3/zcHq1acRPXpxB4XrvT3HOlV8i41CeB54EjgOmSVomabKkZcAnQF/gKTOr0Wu5p7s6tTIZdlkf9uw1bnl2Grv2eH+Kc658Eppt2MyuBr4HzCLopD82vP8cGBJud2kuLzeb+y86kmmLN/Cnt+ZEHY5zLk0ksh4KsH9d9nD99ybABjPzcyPVzDm9D2PCvLU89uF8juvUjFN7xB3P6pxz+xzyeihmts3Mlnkyqb5+eU5PerRqyI9fmMHyjV9FHY5zropL9gJbrhqpm5XJQ1f0Ycfuvfzg2Wns9v4U51wpEpnLa345b/NSGbCrXJ2bN+C3FxzJ5IXr+cs7X0YdjnOuCkukhZLB/nXlY29NgbzwVjvBfbo0cP4xbfhWfjseHj2PMV+ujjoc51wVlchlw3lm1jHOrSnQDXgTmAccnqpgXXTuPu8IurdsyI+en87KTdujDsc5VwUlpTVhZgXAhUAbYGgy9umqlnq1M3noimPYtnOP96c45+JK2ukpM9sOvANcVt7nSBohaZWkmTFl90r6VNJ0SW9LOiwsl6QHJBWE2/vEPOcqSXPD21XJ+p3cgbq0aMi95/fi4wXreOC9uVGH45yrYpLd37GbYKBjeY0EzipW9kcz621mRwP/AX4Vln8D6BrehgCPAEhqRtAqOo5gtP5QSU0P9Rdwpbv42LZc1KctD35QwLi5a6IOxzlXhSQtoUjKBS4AlpT3OWb2IbCuWNmmmIfZQNH6K4OAJy0wEWgiqTVwJvCOma0zs/UEraTiScol0b3nH0Hn5g247fnprNrs/SnOuUC5R8pL+lUJm2oB7Qg+8BsDd1Y0KEm/AQYDG4GBYXEbDkxWhWFZSeXx9juEoHVD+/btKxpmjVW/di0eurwPgx4ax23PTeep644jM0NRh+Wci1giU6/cXcb2TcB9ZvaHQw8nYGZ3AXdJuhO4mSR19BdNGwOQn58fd+VJVz7dWzXk1+cdwe0vf8aw9wu49etdow7JORexRBLKwBLK9wLrgS/MbHfFQzrA08AbBAllKUFLqEjbsGwpMKBY+egkx+HiuDS/HRPmreXv731J347NOKFzTtQhOecilMg4lDEl3Maa2cxkJRNJsV91BwFFyweOAgaHV3sdD2w0s+XAW8AZkpqGnfFnhGUuxSRx3wVHkpeTza3PTWPNFl+s07maLNGpV4Yl8+CSngUmAN0lFUq6Drhf0kxJnxIkh1vD6m8A84EC4HHgRgAzWwfcC0wOb/eEZa4SNKhTi2GX92HDV7v44fPT2bvXzyQ6V1PJrHwfAJI2Aw+a2c9TG1LlyM/PtylTpkQdRrXx9MeLuOvVmQzp34k7zupBhnfSO1ctSZpqZvnxtiXSh/I50Dk5Ibnq5vK+7fl82SaGfzifxWu38edLjyK7TsLL7Tjn0lgi41AeAM4tvp68cxD0p/zm/F788pyevD1rBRc+/BGL1/pSOc7VJIl8hSwE3gXGS3qMoL9iBfsHHu4TDlh0NYwkrjupI91aNuDmZ6Zx3kPjeOjyPvTrkht1aM65SpBIH8peguRRdHK8xCeaWWbFQ0st70NJrYVrtnLDk1OYv2Yrd519ONf0y0PyfhXn0l2y+lDuoZQk4lysvNxsXr2pHz98fjr3/GcWs5Zv4r7ze1E3q8p/13DOHaJyJxQzuzuFcbhqqEGdWjx25bH87b25PPDeXApWbeGx7xxLy0Z1ow7NOZcCiYxDaS+pURl1GkrySbLcPhkZ4kend+PRK/vw5crNnPvgOKYtXh91WM65FEjkKq8F7B9kWJIfhPWcO8BZvVrzyo0nUicrg289NpGXphZGHZJzLskSSShFa8g7d0h6tGrEqJtOIj+vKT95cQb3vD7LV350rhpJ9gJbrYCtSd6nq0aaZtfmyWv7ck2/PEaMX8BV/5rE+q07ow7LOZcEpXbKSxpcrOjoOGUAmUB74ErgsyTF5qqpWpkZDD33CA5v3YhfvDqTQQ+N5/HB+XRv1TDq0JxzFVDqOJSYsSdl7ie83wZcaGZvJyG2lPJxKFXDJ4vX892nprJ1x27+cunRnNUrkRWknXOVrbRxKGUllKuKfgRGAK8B/45TdQ+wFphgZhsqFm7l8IRSdazYuJ3v/t9UZizZwG1f78oPTu3qk0s6V0Ud8sBGM3siZidXAa+Z2ZNJjs/VcK0a1+X5Icdz16sz+du7c5m1bBN/+dbRNPDJJZ1LK4kssDXQk4lLlbpZmfzpkt788pyevDt7JRc+PJ5Fa/36DufSSSIDG5tK6impTrHyayT9W9IzkvomP0RXUxRNLvnktcexctMOzhs2nnFz10QdlnOunBK5bPi3wMexz5F0C/AP4Fzg28BoST2TGqGrcU7qmsuom/vRslEdBo/4mH+OW0B5JzF1zkUnkYTSD3jPzL6KKfsJsBToD1walv0oSbG5GqxDTjav3NiPrx/eknv/M4ufvPgp23ftiTos51wpEkkobYiZViVsibQjWBZ4nJm9BLxOkFycq7AGdWrx6JXHcutpXXn5k0IueXQCkxasizos51wJEkko9YDtMY/7EYxReTembB5B4nEuKTIyxA9P78ajVx7L8o3bufSxCVz5j4+ZusgTi3NVTSIJZSnQI+bxmcAmYEZMWVMg9pSYc0lxVq9WjP3ZQH7xzcOZvXwTFz0ygcEjJvnMxc5VIYkklA+AsyXdLOl64DzgTTOLnd2vM7AkmQE6V6Re7UyuP7kTY28fyJ3f6MFnhRu44OGPuHbkZD4r3Bh1eM7VeIkklN8BW4C/A8MJTn/dXbQxXCvlJOCj8u5Q0ghJqyTNjCn7o6QvJH0q6VVJTWK23SmpQNIcSWfGlJ8VlhVIuiOB38mlofq1a/HdUzoz9vZT+emZ3Zm6aD3nDhvH9U9MYeZSTyzORaXca8oDSGoFXBw+HGVmi2O29QG+AzxjZpPLub/+BEnqSTPrFZadAbxvZrsl/R7AzG4PLwJ4FugLHEbQd9Mt3NWXwOlAITAZuMzMZpV2bJ96pfrYvH0XI8cv5PGx89m0fTdnHdGK207vSo9Wpa4H55w7BMlaUx4zWwEMK2HbJ8AnCe7vQ0l5xcpiJ5acyP4ENgh4zsx2AAskFRAkF4ACM5sPIOm5sG6pCcVVHw3rZnHLaV0ZfGIeI8YtYMS4Bbz5+Qq+eWRrbv16V7q19FmMnasMyV4PJdmuBf4X/tyGA/tnCsOyksoPImmIpCmSpqxevToF4booNa6XxQ9P78a420/lllO7MHrOKs7824fc8uw0ClZtiTo856q9hFookrIIvv33JbiiKzNONTOz6yoamKS7gN3A0xXdVxEzG07Q/0N+fr4Pva6mGtfP4sdndOfafh15fOx8Rn60kP9+uoxBR7fhB6d1pWNudtQhOlctlTuhSDoMeIfg0uHS5hY3oEIJRdLVwDnAaba/k2cpwUDKIm3DMkopdzVY0+za/OysHlx3UkeGfzifJyYs5N/Tl3LBMW35wWld6JDjicW5ZEqkhfJn4HCCjvHHCU4z7U52QJLOAn4GnGJm22I2jQKekfQXgk75rsAkguTWVVJHgkTybeDyZMfl0ldOgzrcefbhXH9yJx4bM4+nJi7itelLuahPG245tSvtmtWPOkTnqoVyX+UlaS3wmZkNSNrBpWeBAUAusBIYCtwJ1CFYsAtgopl9L6x/F0G/ym7gNjP7X1h+NvA3glNwI8zsN2Ud26/yqrlWbdrOw6Pn8cykxezda1yS346bT+1Cmyb1og7NuSrvkFdsLLaTrcAwM7s9mcFFxROKW7FxOw+PLuC5SUswjEvz23HTwC4c5onFuRIlK6F8DCwws28nM7ioeEJxRZZt+IphHxTw4pQlCPGtr7XjxoGdad3YE4tzxSUroVwMPAnklzVoMB14QnHFFa7fxkMfFPDilEIyJC7r244bB3ahZaO6UYfmXJWRrITSH7iJYFLIvwNTgQ3x6prZh4cWauXxhOJKsmRdkFhemlpIRoa4vG97bhzQmRaeWJxLWkLZS3BJcNElwyU+0czijU+pUjyhuLIsXruNYR/M5eVPllIrQ1xxXAe+N6ATLRp6YnE1V7ISyt2UkkRimdmvyx1dRDyhuPJatHYrD75fwKvTlpKVKa48rgPfPaUzzRvWiTo05ypdUhJKdeMJxSVq4ZqtPPD+XF6btpTatTL4zvFBYslt4InF1RyeUOLwhOIO1fzVW3jw/QL+PX0pdWplMviEDgzp34kcTyyuBkh6QpF0EnAM0ATYCHxiZuMqFGUl84TiKmre6i088N5cRs1YRr2sTAafkMeQ/p1oll076tCcS5mkJRRJxwJPAd2LitjfrzIHGGxmafEp7QnFJUvBqs088F4Br3+6jPpZmVx1Yh43nNyJpp5YXDWUrE75LsAUoBEwDngfWA60Bk4lWK1xI9DXzOYmIe6U8oTikm3uys38/b25/Pez5dTPyuTqfkFiaVLfE4urPpKVUJ4ArgS+bWYvxtl+MfAc8LSZXVWBeCuFJxSXKnNWbOaBMLE0qFOLa/rlcf1JnWhcPyvq0JyrsGQllKXABDO7uJQ6LwPHm1ncBa6qEk8oLtW+WLGJv787l//NXEH92plc1KctV52YR5cWDaIOzblDlqwlgHOBL8qo8wXBOibO1Xg9WjXikSuPZfbyTfxj7KMnW0UAABkRSURBVAKen7yEpyYu4uSuuVx9Yh4Du7cgI6O0pYWcSy+JtFAKgUlmdmEpdbyF4lwJ1mzZwbMfL+b/Pl7Eyk076JBTn8En5HFJflsa1fXTYS49JOuU15MEC1ddaWbPxdl+EfA83ofiXKl27dnLmzNXMPKjhUxdtD7mdFgHurRoGHV4zpUqmVd5TQUaAB8BHxBc5dWKYJGsk4DNwNf8Ki/nyuezwo2M/Gghr89Yxs49e/10mKvykjkO5WsEU9gXjUOJnSxyDnCVmU2qQKyVxhOKq0rWbNnBc5MW89REPx3mqrZUjJQ/EegDNCYYezLNzMZXKMpK5gnFVUVFp8Oe+GghU/x0mKuCfC6vODyhuKpu5tLgdNio6X46zFUdyepDqQc0B1aY2c442+sALYFVZra9AvFWCk8oLl0UnQ77v4mLWbFpOx1y6vOd4ztwSX47Gtfz02GuciUrofwOuA1oY2br4mxvBhQCfzKzX1Ug3krhCcWlm1179vLW58HpsMkL/XSYi0ayEsp0YImZnVtKnVFAWzPrc0iRViJPKC6d7TsdNmMZO3fvpXPzbE7qkku/Lrkc3znHO/JdypSWUDIS2E8e8GUZdb4M65WLpBGSVkmaGVN2iaTPJe2VlF+s/p2SCiTNkXRmTPlZYVmBpDvKe3zn0lWvNo350yVHMeGOU/nFNw+nXbP6vDClkCFPTeWYe97hgofH85e35/Dx/LXs3L036nBdDZHI1CtZQFnvTAMSWXB7JDCM4FLkIjOBC4HHYitK6gl8GzgCOAx4V1K3cPNDwOkEp9wmSxplZrMSiMO5tJTToA7Xn9yJ60/uxM7de5m2eD3jC9YwtmANwz4o4IH3C6iXlclxnZrta8F0b9nQO/VdSiSSUOYDp5RRZwCwqLw7NLMPJeUVK5sNIB30hh8EPGdmO4AFkgqAvuG2AjObHz7vubCuJxRXo9SulcFxnXI4rlMOPzqjO5u272LivLWML1jDuII13Pff2QDkNqjNiZ1zgwTTNZc2TepFHLmrLhJJKKOAOyT9zMz+UHxjeKqpD3DQtiRpA0yMeVwYlgEsKVZ+XIpicC5tNKqbxRlHtOKMI1oBsHzjV4wv2J9gRs1YBkCn3Gz6ha2XEzrl+DT77pAlklD+BFwB/E7SpcDbwFKCD/UzgaOBxaQuoVSYpCHAEID27dtHHI1zlat143pcfGxbLj62LWbG3FVbGDt3DeML1vDyJ4U8NXERGYIj2zbhpC459OuSy7EdmlKnVmbUobs0Ue6EYmbrJQ0AngGOJ2iNxE698hHBxJHrkx1kaCnQLuZx27CMUsoPYGbDgeEQXOWVghidSwuS6NayId1aNuS6kzqyc/deZhRuYFyYYB4dM5+HPphH3awM+nbMYWD35gzo3oKOudlRh+6qsEOdeqUPQVJpAmwAJprZJ4cUQNCH8h8z61WsfDTwk6I16iUdQZDM+hJ0yr8HdCVIaF8CpxEkksnA5Wb2eWnH9cuGnSvZ5u27mLRgHWPnruHDuauZv3orAB1y6jOgW5Bcju+UQ73a3nqpaars1CuSniXoyM8FVgJDgXXAgwSj8jcA083szLD+XcC1wG7gNjP7X1h+NvA3IBMYYWa/KevYnlCcK7/Fa7cx+stVjJ6zmo/mrWH7rr3UqZXB8Z1yGOCtlxolFZNDZgPdgAZmNraC8UXCE4pzh2b7rj1MWrCO0XNWM3rOKuav8dZLTZLM6evbAn8HziVoDZiZ1Qq3nUTQP3GjmY2uaNCp5gnFueQoqfVyXKccBnRrzsAe3nqpTpI19UprYArBBJCjgBbACWaWGW7PIlhw60Uz+34yAk8lTyjOJZ+3Xqq/ZCWUR4HrgDPM7ANJQ4FfFSWUsM6rQGcz652EuFPKE4pzqVdW66V/t+Z0ys32kftppLSEksg4lLOBUWb2QSl1FgMnJxKcc676ah+uPDn4hLyDWi/3/CeYzCK7dibdWzXk8NaN9t16tGpIdp1EPp5cVZDIX6wlUNZa8bsAP1nqnDtI3axM+oetkl+d25PFa7cxYf4aZi/fzKzlmxg1YxlPf7x4X/0OOfU5vFVRkgkSTtum9eJNy+SqiEQSyjoOHEAYTzdgxaGH45yrKdrn1Kd9zv4ZK8yMpRu+Yvbyzcxevmnf7a1ZKyg6M9+wbq0wyexv0XRv1ZC6Wd4nUxUkklDGA+dJamVmByUNSV2Bs4D/S1ZwzrmaQxJtm9anbdP6nN6z5b7yrTt288WKzXyxoijJbOalqYVs3bkHgAxBx9zsfQmmZ3jfslEdb81UskQSyh8JZvEdI+k2oD7sG5PSH/grwfT2f052kM65miu7Ti2O7dCUYzs03Ve2d6+xZP02Zi/fxKywRTN9yQb+8+nyfXWa1s+iS4sGdMjJJi+nPnm52eTlZNMhpz4NfQGylEh0HMq1wCPET0S7gWvN7OkkxZZSfpWXc9XPpu27+CLmlNn8NVtZtHYrKzftOKBeboPaYaIJkk2H3Gw65mTTIbe+r3ZZhqSOlA9Pbd1IMJdXDrCRYFr5YWY2p4KxVhpPKM7VHNt27mbR2m0sWruVBWuC+4Vrt7JwzTZWbNp+QN1m2bWDFk1ONnm5QYum6OfG9TzZVOpcXpKam9nqpO40BTyhOOcAvtq5h0XrguQSm2gWrd3Kso0HJpum9bPokJNNxzDRdAxPo9WkZJOscShlHaQxcDtwM9AoWft1zrlUqlc7kx6tGtGj1cEfW9t37WHxum0sXBMmmrCVM2nBOl6bvpTY7+P7Wjbh6bO83CDx5OVm06CGjKkp128pqQNwLME4k0lmtjJmW13gh8BPgKbAthTE6Zxzla5uVua+dWOKK0o2C9Zs3ZdwFqzZykcFa3nlkwOXZMptUIeOuftPne1v2dSnfu3qk2zK/E0kPUDQZ1J0/d1OST82s4fDBbeeIFjUaifBxJG/S1GszjlXZZSWbL7auSc8dbaVBeH9wjXbGP3lalZPLTygbstGdcjL2d+ayYs5pZZu42tKTSiSriI4hbUXmB0W9wAekLQVeIxg1uHHgPvMbFkKY3XOubRQr3bmvnExxW3ZsXv/KbQ1wUUCC9du5Z1ZK1m7decBdVs3rrvvooDYy5875FTNlk1ZEV1N0PIYaGYTACT1B94B/gkUAuea2WepDNI556qLBnVq0atNY3q1aXzQtk3bd4VJZiuL125jwdqtLFq7jXdnr2TNlgOTTYuGdfaNq4kdYxPlOJuyEkpv4NWiZAJgZh9Keg24mGDciScT55xLgkZ1s+jdtgm92zY5aNvm7btYtDZozSxaG1wosGht/NNo8cbZ5OXUp0NOaq9GKyuhNAYK4pQXTRI5Ic4255xzSdawblaJLZutO2LG2azdyqLwNNr4gjW8/MnB42z6dcnlwcuOSXqMZSWUDIIru4rbBWBmXyU9IueccwnJrlOLnoc1oudhB/fZbNu5O7z0uWiczTaaZaemlVKeXp3kjnx0zjlXaerXrlXiOJtkK09CuVvS3fE2SNoTp3jfOvPOOedqjvJ88Cc6/7PPF+2cczVQqQnFzDIqKxDnnHPpLdKEIWmEpFWSZsaUNZP0jqS54X3TsFySHpBUIOlTSX1innNVWH9uOBjTOedcJYu6BTKSYJXHWHcA75lZV+C98DHAN4Cu4W0IwbosSGoGDAWOA/oCQ4uSkHPOucoTaUIxsw8J1qqPNYhgfjDC+/Njyp+0wESgiaTWwJnAO2a2zszWE4ziL56knHPOpVjULZR4WppZ0TqeK4CixaXbAEti6hWGZSWVH0TSEElTJE1ZvbrKL9ninHNppSomlH0sWP0raeNgzGy4meWbWX7z5s2TtVvnnHMkcYGtJFopqbWZLQ9Paa0Ky5cC7WLqtQ3LlgIDipWPLusgU6dOXSNpUQXizAXWVOD5UUnXuMFjj4rHHo2qGnuHkjZUxYQyCrgKuD+8/3dM+c2SniPogN8YJp23gN/GdMSfAdxZ1kHMrEJNFElTSloGsypL17jBY4+Kxx6NdIw90oQi6VmC1kWupEKCq7XuB16QdB2wCLg0rP4GcDbBZJXbgGsAzGydpHuByWG9e8yseEe/c865FIs0oZjZZSVsOi1OXQNuKmE/I4ARSQzNOedcgqp0p3wVNzzqAA5RusYNHntUPPZopF3sCr74O+eccxXjLRTnnHNJ4QnFOedcUnhCSYCkdpI+kDRL0ueSbo06pkRJypQ0TdJ/oo4lEZKaSHpJ0heSZks6IeqYykvSD8P3y0xJz0qqG3VMJUlkwtaqpoTY/xi+Zz6V9Kqkgxdrj1i8uGO2/ViSScqNIrZEeUJJzG7gx2bWEzgeuElSz4hjStStwOyogzgEfwfeNLMewFGkye8gqQ3wAyDfzHoBmcC3o42qVCMp/4StVc1IDo79HaCXmfUGvqQcY9QiMJI48w9Kakcwrm5xZQd0qDyhJMDMlpvZJ+HPmwk+1OLOG1YVSWoLfBP4R9SxJEJSY6A/8E8AM9tpZhuijSohtYB6kmoB9YFlEcdTogQnbK1S4sVuZm+b2e7w4USCmTSqlBJec4C/Aj8jjZZh94RyiCTlAccAH0cbSUL+RvAG3Rt1IAnqCKwG/hWervuHpOyogyoPM1sK/IngW+Zyghke3o42qoSVNGFrurkW+F/UQZSHpEHAUjObEXUsifCEcggkNQBeBm4zs01Rx1Meks4BVpnZ1KhjOQS1gD7AI2Z2DLCVqnva5QBhf8MggqR4GJAt6cpoozp0yZ6wtbJIuovglPXTUcdSFkn1gZ8Dv4o6lkR5QkmQpCyCZPK0mb0SdTwJ6AecJ2kh8BxwqqT/izakcisECs2sqDX4EkGCSQdfBxaY2Woz2wW8ApwYcUyJWhlO1EqxCVvTgqSrgXOAKyw9Bt51JvgCMiP8f20LfCKpVaRRlYMnlARIEsF5/Nlm9peo40mEmd1pZm3NLI+gU/h9M0uLb8pmtgJYIql7WHQaMCvCkBKxGDheUv3w/XMaaXJBQYyiCVvhwAlbqzxJZxGc5j3PzLZFHU95mNlnZtbCzPLC/9dCoE/4f1CleUJJTD/gOwTf7qeHt7OjDqqGuAV4WtKnwNHAbyOOp1zCVtVLwCfAZwT/c1V2So1wwtYJQHdJheEkrfcDp0uaS9Diuj/KGEtSQuzDgIbAO+H/66ORBhlHCXGnJZ96xTnnXFJ4C8U551xSeEJxzjmXFJ5QnHPOJYUnFOecc0nhCcU551xSeEJxVY6k0ZKq1eWHkrqGs92uCGePrdBcZJIWhoPeYsuuDvd9dUX2nU6q43slnUW6prxLnZh/ssVAdzPbHqfOQqADkBUzgZ5LMkmZwGtAF+ApgoFqB/09nEt3nlCqv/bAbVTRwWg1REegJ/C4mQ1J4XFeJZhRd3lZFauRwQQzOLsqwBNK9baeYCK/OyT9w8zWRB1QDXVYeJ/SaevNbCOwMZXHqGrMLG3WCqkJvA+letsG3As0BoaW5wmSBoTn4e8uYXup5+4lnS5prKQtklZL+lfRKnmSjpH0H0nrw+2jwmUASoqljqT7JC2QtEPSPElDJdUuoX4PSSMlLZG0U9JKSc/EzAEWW3dkGHMnSbeEK/p9JWl0OV+nYyW9HK60t0PSIkkPF02iGFPPgDHhw6HhMUt8fYs9V5JuVrDa43ZJSyUNC9eHiVc/bh9K0d9MUgNJfw1fn6/CqUjOD+vUknSXglUZt4ev9c2lxHampDckrYn52/xRcVZEjDl+dlhncficAkm3S1Kc55wn6T1Jy8O6yySNkXRjsXpx+1AkZUj6nqTJ4Xtta/jz9yUd9LkXvm6jJeVKGh5z3M8lXROnviRdJemj8H2+PXxd35L0rZJet+rOWyjV30PAzcB3JT1gZnNTeKzzCGZ1/Q/wKMGsulcDeZLuJFjtbyzBBJtHAucCnST1NrN4a7S8AHyNYC6sXQTTwN8N5Es6L3bmWAWTAL4CZAGvAwUEs7ReCHxT0sCixdGK+TtwMvBf4A1gT1m/pIKlAF4GFMa2CDgW+D4wSNJJZrYgrP5rII9gUsUxwOiwfDRl+xvBao/LCeb/KnoNjgNqAzvLsY8iWQSrFzYjmNyxNnAZ8LKkM4Abw/3+D9gBXAI8KGm1mT1f7PcfSvB3WEfwt14F9AZ+Apwt6YQ4yzpkAW8RtNb+RzCV/PkEp2LrErxORfsfAjxGsPbK68AaoEV4jGuAh8vx+z4FXA4sIVhQzoALwueeBFwR5zlNgPEEr+tLQJ3wdRghaa+ZPRFT9zcEqz8uIHifbgRaE7xfLwEOeM1qDDPzWzW8EfwDFYY/Xxw+fqVYnYVhea2YsgFh2d0l7HchsLBY2dXhc3YDp8SUZxB8iBnBh88VxZ73z3DboGLlo8PyL4GmMeV1CSbRM+A7MeVNCU7vrQF6FttXL2AL8Emx8pHhfpYCHRN4XRsAawkSz8nFtt0e7vPtYuWlvqYlHOfE8DkFQLMSXoOS/g5Xl/B3fh2oE1N+cszfZjLQJGZbJ4IP1mnF9jUwfM5HsfWLHf+vJRz/DaBeTHkLYEN4y4opn0qQ1FrEeV1y471XipVdFh7vE6BBTHk2MCXcdnmc/xcjSD6ZMeU9Cd7Xs4rVX0twcUX9smKsSTc/5VUDmNlLBB9CF0g6KYWHetbMik7vYEGr46nw4UwzK7640ZPh/dEl7O9eM1sfs7/t7F8T/NqYeoMJvl0ONbMDprU3s5nA48AxknrGOcYfbH9rojwGEXzLf97Mxhbb9meCD8/TJbVPYJ/xFJ1m+Y2Z7VsetthrkKjbzGxHzL7GEnzDbgrcbjHLKpvZfIJv670UXKVW5Afh/Q1WbBlmMxsJTCf+t3+AH5jZVzH1VxG0lhoDxU9L7iZokR3AytcPWPTeuMPMtsQ8dytB0ge4Ps7ztgE/MrM9Mc+ZRfA6HK5gYb1Yu4jToi1njNWSn/KqOX5M8K3yT8DxKTrGlDhlRR3R8VaKXBrel7TO95g4ZeMI/omPiSk7Ibw/qoS+iW7h/eEcvI7KpBKOXZKihb3eL77BzHZL+pDgFNcxBJdsH6qi45T2GiRig5nNi1O+jOAqtJL+PrWAVuz/W51A8EF6iaRL4jynNtBcUo6ZrY0p32hmBXHqLwnvm8aUPU2QnGdJeo7gNRhvZqvj/2oH6UOwzPXoONvGcPD7p8hci78Ca2yMRQnqaYIlFWZJeiHc7wQLLoyosTyh1BBmNkHSS8DFkr5lxc6LJ0m8f6bd5diWVcL+VhYvCD+0i86pF8kJ728oI77i3zAhOE+fiKIO8ZIuzS0qP6hz+hCPU9prkIiSPuh2h/ss798nh+Bzo6yLPIpODRYpaSBn0TH2tYLM7C/h73cjQYvoNsAkjQF+ambxvrjEagysM7OD+phKeP8kHCPwQ2A+QUvyjvC2W9IbwI9LSJ7Vnp/yqlnuJPh2+TuVcKUUwTc7KPnLRkU/KBPRsniBpFpALhD7TbLow/AoM1MptyeK74/E10cvOlZJy7G2LlbvUBU9v7TXIAobgfVlvM4ys0UVOYiZPWlmxxMksG8S9Lf1B96S1LwcMTZTsFz3AUp4/xxKfHvM7G9mdhTB3+gignFA5wFvSqpTkf2nK08oNUj4relhglMct5RQrajPol3xDZK6sP+bc2U4JU7ZSQTfFKfFlE0M709OeUT7jzug+Ibww6oohnhXlCWi6PmlvQZRmAg0lXREZRzMzDaY2RtmdgPBhRTNCBJLaaYRfLbFq9ef4LWr6N8nNsZVZvaKmV1KcCq0M8HFIDWOJ5Sa5x6Cpv1dxD8F9AXBt7dBkvadFpBUD3igUiLc75eS9p1bl1QX+F348F8x9f5F8DsNldS3+E7CMQkDkhTTawRXRV0mqXhf1G0Eyfpdq/iAu5Hh/V2SmhUVFnsNovDX8P5xSYcV3xiONalQH52kgfHGprD/NFVZa8OPCO9/J2nfKPrw56IZI/5ZgfjqSOoXpzyLIOGVJ8ZqyftQahgzWyfpt8AfSti+S9LfgV8C0yS9SvA+OZ2gAzelo72LmQ18Hvb9FI3B6EwwZqTo6jHMbK2kiwmnHpH0HvA5wemsdgQdyTkEl9xWiJltkXQt8CIwRtKLBJ3vxwJnEPTJfDcJxxkv6UGCluTMYq/BeiKaXsXM3pN0B0FSmxv2GSwg+HLSgaBFNQ44qwKHeRXYImkiwVVzImj5fY3g4oF3y4jxGUmDgEsJ3j+vEbwXzidI+M/HueIwEfWAcZIKwngWEby3Tie48GOUmc2uwP7TlieUmukBgg7PvBK2DyX4hnUDMITgQ/I5gsFsxa+SSqVLCRLbFQQD4paGMdxvZgf0fYQfdEWD684k+ADaSZAA3ycYiJgUZvbv8Bvqz8NjNSZ4jR4luNQ5WUn3VoKxODcRJKm1BB+2PwdmJOkYCTOz30saT9BhfhJBkttI8PcZDjxTwUPcQfC69gHOJphIcxHBJb+PmNlBlxPHcRnBlVfXsj/Bzya4euyRCsZXdPnxQILxQucDm4F5BINbR5T81OpNxf4vnXPOuUPifSjOOeeSwhOKc865pPCE4pxzLik8oTjnnEsKTyjOOeeSwhOKc865pPCE4pxzLik8oTjnnEsKTyjOOeeS4v8BUQTpOpPcRCYAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Q7rPltFljItN",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "34928d49-d899-407a-b953-3ec7074b942f"
      },
      "source": [
        "x = range(2, len(error)+2)\n",
        "\n",
        "from kneed import KneeLocator\n",
        "kn = KneeLocator(x, error, curve='convex', direction='decreasing')\n",
        "print(kn.knee)"
      ],
      "execution_count": 20,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "6\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yHN60HACjiUm",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 296
        },
        "outputId": "607da31f-8762-479b-c316-6109bb8955ce"
      },
      "source": [
        "plt.xlabel('Number of dimensions')\n",
        "plt.ylabel('Reconstruction error')\n",
        "plt.plot(x, error, 'bx-')\n",
        "plt.vlines(kn.knee, plt.ylim()[0], plt.ylim()[1], linestyles='dashed')"
      ],
      "execution_count": 21,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.collections.LineCollection at 0x7fad89c45fd0>"
            ]
          },
          "metadata": {},
          "execution_count": 21
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debyWc/7H8denjcoSyjItKioTYz0lw1Ak1SA0yFjK+MlO9ppByIytbEMRkm1ahJFtqJxknRayFNEUOiYqjKUo1ef3x/e6nbvTOee+zzn3fa77Puf9fDyuR/f9va77uj8n6dP1XT5fc3dERETKUyfuAEREJPcpWYiISEpKFiIikpKShYiIpKRkISIiKdWLO4BsaNq0qbdu3TruMERE8sqcOXNWuHuz0s7VyGTRunVrZs+eHXcYIiJ5xcw+LeucuqFERCQlJQsREUkpa8nCzMaY2TIzez+pbYKZzY2OT8xsbtK5IWa20MwWmNlhSe09o7aFZjY4W/GKiEjZsjlmMRa4E3go0eDuxydem9kI4NvodUegH7Ar8Ctgqpm1jy69CzgUKAJmmdlkd5+fxbhFRKSErCULd59hZq1LO2dmBhwHHBw19QHGu/tqYLGZLQQ6R+cWuvui6HPjo2uVLEREqlFcYxa/A75094+j982BJUnni6K2stpFRKQaxZUsTgDGZfKGZjbQzGab2ezly5dn8tYiIrVetScLM6sHHANMSGr+HGiZ9L5F1FZW+0bcfbS7F7h7QbNmpa4pERGRSorjyaI78KG7FyW1TQb6mdkmZtYGaAfMBGYB7cysjZk1IAyCT672iEVEarlsTp0dB7wBdDCzIjM7LTrVjxJdUO4+D5hIGLj+F3COu69z97XAucALwAfAxOhaqUGmTp3K1KlT4w5DRMphNXGnvIKCAle5j/zRtWtXAKZPnx5rHCK1nZnNcfeC0s5pBbeIiKSkZCEiIikpWYiISEpKFpGbboLCwg3bCgtDu4hIbVcj97OojE6d4LjjYOJE6NYtJIrEe8mue+65J+4QRCQFJYtIt24wZgz8/vdw6qkhSSQSh2RXhw4d4g5BRFJQN1SSffeFn3+GkSPhrLOUKKrL008/zdNPPx13GCJSDiWLJPPmQf364fWdd248hiHZMWLECEaMGBF3GCJSDiWLSGKMYtw42GIL2HXX8F4JQ0REyeIXs2aFMYo+feCii+DVV+G660K7iEhtpwHuyGWXFb8eNAjuuAMmT4Znn40vJhGRXKEni1JsuSVceik89xy8/nrc0YiIxE+FBMuwciW0bQu77QbTpmUoMCnVkiVhM8SWLVumuFJEskmFBCuhcWMYMgReekmD3NnWsmVLJQqRHKdkUY4zz4TmzeHKK6EGPoDljAkTJjBhwoTUF4pIbJQsyrHppnDFFfDaa/DCC3FHU3ONGjWKUaNGxR2GiJRDySKFP/0JWrcOSUNPFyJSWylZpNCgAVx1FcyZA089FXc0IiLxULJIw8knQ7t2IWmsXx93NCIi1U/JIg316sE118B778Fjj8UdjYhI9ctasjCzMWa2zMzeL9F+npl9aGbzzOympPYhZrbQzBaY2WFJ7T2jtoVmNjhb8aZy/PGhXtTQobB2bVxR1EyTJk1i0qRJcYchIuXI5pPFWKBncoOZdQP6AHu4+67A8Ki9I9AP2DX6zEgzq2tmdYG7gF5AR+CE6NpqV6cOXHstLFgAjz4aRwQ1V9OmTWnatGncYYhIObKWLNx9BvB1ieazgBvcfXV0zbKovQ8w3t1Xu/tiYCHQOToWuvsid18DjI+ujcXRR8Nee4UuqZ9/jiuKmmfs2LGMHTs27jBEpBzVPWbRHvidmf3bzF42s05Re3NgSdJ1RVFbWe0bMbOBZjbbzGYvX748C6GDWahEu3gxPPBAVr6iVlKyEMl91Z0s6gFbA12AS4GJZmaZuLG7j3b3AncvaNasWSZuWapevaBLFxg2DH76KWtfIyKSU6o7WRQBT3gwE1gPNAU+B5KLA7WI2spqj03i6aKoCEaPjjMSEZHqU93J4p9ANwAzaw80AFYAk4F+ZraJmbUB2gEzgVlAOzNrY2YNCIPgk6s55o0cfDB07Qp/+xusWhV3NCIi2ZfNqbPjgDeADmZWZGanAWOAttF02vFA/+gpYx4wEZgP/As4x93Xufta4FzgBeADYGJ0bazMQjfUl1/CXXfFHY2ISPZpP4sq6NUrbLu6aFHYt1sqZ1X0eNaoUaOYIxGp3bSfRZZcey189RXcfnvckeS3Ro0aKVGI5Dgliyro1An69IERI+Cbb+KOJn+NHDmSkSNHxh2GiJRDyaKKrr0Wvv02JAypnIkTJzJx4sS4wxCRcihZVNHuu4e6UbfdBllaCygiEjsliwy4+mr48Ue48ca4IxERyQ4liwzYZRc46aQwjXbp0rijERHJPCWLDLnqqlC6/G9/izsSEZHMU7LIkJ12Cvt133MPfPpp3NHkl+nTpzN9+vS4wxCRcihZZNAVVxTXjhIRqUmULDKoZUs444xQvnzhwrijyR/Dhw9n+PDhcYchIuVQssiwIUOgQYOw/kLS88wzz/DMM8/EHYaIlEPJIsN22AHOPRceeQTmz487GhGRzFCyyILLLoPGjcP6CxGRmkDJIguaNoVBg+Cxx2Du3LijERGpOiWLLLn4YmjSBIYOjTuS3NewYUMaNmwYdxgiUo56cQdQUzVpApdcEqbTzpwJnTvHHVHuev755+MOQURS0JNFFp1/fuiSuvLKuCMREakaJYss2nxzuPxyePFFeOWVuKPJXcOGDWPYsGFxhyEi5Sg3WZhZXTMrrK5gaqKzz4bttw/dUTVwB9uMmDZtGtOmTYs7DBEpR7nJwt3XAevNbMuK3tjMxpjZMjN7P6ntajP73MzmRkfvpHNDzGyhmS0ws8OS2ntGbQvNbHBF44hbo0bwl7/AjBmgvw9FJF+l0w31A/Cemd1vZnckjjQ+NxboWUr7re6+Z3Q8B2BmHYF+wK7RZ0ZGTzV1gbuAXkBH4ITo2rzy3XfQrNmGTxeFhXDTTfHGJSKSrnRmQz0RHRXi7jPMrHWal/cBxrv7amCxmS0EEvOHFrr7IgAzGx9dm1dro/fbD66/Hv79b3j22bBg77jjQDuJiki+SJks3P1BM2sAtI+aFrj7z1X4znPN7BRgNnCxu38DNAfeTLqmKGoDWFKifd/SbmpmA4GBAK1atapCeJnXrRs88QT07BnKmLuHRNGtW9yR5YZtttkm7hBEJIWU3VBm1hX4mNAdNBL4yMwOrOT3jQJ2AvYElgIjKnmfjbj7aHcvcPeCZs2aZeq2GXPooeFpYvnysPeFEkWxxx9/nMcffzzuMESkHOmMWYwAerj7Qe5+IHAYcGtlvszdv3T3de6+HriX4q6mz4GWSZe2iNrKas87hYUwdSoUFITuqJtvjjsiEZH0pZMs6rv7gsQbd/8IqF+ZLzOzHZLeHg0kZkpNBvqZ2SZm1gZoB8wEZgHtzKxN1BXWL7o2rxQWFo9RzJgBrVrB4MHwz3/GHVluGDJkCEOGDIk7DBEpRzoD3HPM7D7gkej9iYTxhnKZ2TigK9DUzIqAoUBXM9sTcOAT4AwAd59nZhMJA9drgXOiabuY2bnAC0BdYIy7z0v7p8sRs2ZtOEbx1FPhCeOKK6BPn7C7Xm32xhtvxB2CiKRgnmKlmJltApwDHBA1vQKMjGYu5aSCggKfPTtlPovV8OFw6aVw331w2mlxRxOvrl27AmgfbpGYmdkcdy8o7Vy5TxbROod33H0X4JZsBFdbXXQRPP98qB/1u99B+/apPyMiEpd0VnAvMLPcmotaA9SpAw8+CJtsAiedBD9XZTKyiEiWpTPAvRUwz8ymmdnkxJHtwGqDFi3g3nvDmMY118QdTXxatGhBixYt4g5DRMqRzpjFQaW1u/vLWYkoA/JhzCLZn/4EY8fCyy+HLikRkTiUN2ZRbrKIxizmRWMWeSPfksX338Nee4WuqHfeCRsniYhUt/KShcYscsDmm8Ojj8Lnn8M558QdTfUbNGgQgwYNijsMESlHOussEmMWM4GViUZ3PzJrUdVC++4b9uu+6iro3RtOPDHuiKrP3Llz4w5BRFJIJ1loU9BqMmQI/OtfYcOk/feH1q3jjkhEJEg5GyoayP6EUPbjZUIJjreyHFetVK8ePPJIqEp78smwbl3cEYmIBOlUnT0dmATcEzU1B1TVKEvatIG77oJXX4Ubbog7GhGRIJ11FucA+wPfAbj7x8C22QyqtjvpJOjXD66+GmbOjDua7Gvfvj3ttYRdJKelM2ax2t3XWFTtzszqEQoBSpaYwahR8NprYaD77bdhs83ijip7Ro8eHXcIIpJCOk8WL5vZn4GGZnYo8BjwdHbDkiZN4OGH4T//gQsvjDsaEant0kkWg4HlwHuEkuLPAVdkMygJDjoILr88VKZ98sm4o8megQMHMnDgwLjDEJFypLMHd2JXu3uzH46UdM01MGUK/N//hbUYv/pV3BFl3kcffRR3CCKSQjpPFhKjBg3C6u6ffoL+/WH9+rgjEpHaSMkiD3ToALfeGvbwvv32uKMRkdpIySJPnH562IJ18OBQbFBEpDqlHLMws/bApcCOyde7+8FZjEtKMAt7X+y+e5hOO2sWNGwYd1SZseeee8YdgoikkM5+Fu8AdwNzgF8KULj7nOyGVnn5VqK8Il54AXr2hPPOgzvuiDsaEalJKl2iPLLW3Ue5+0x3n5M40vjSMWa2zMzeL+XcxWbmZtY0em9mdoeZLTSzd81s76Rr+5vZx9HRP414a7TDDoMLLoC//z3s4S0iUh3SSRZPm9nZZraDmW2dONL43FigZ8lGM2sJ9AA+S2ruBbSLjoHAqOjarYGhwL5AZ2ComW2VxnfXaDfcALvtBqeeCsuWxR1N1Z100kmcdNJJcYchIuVIJ1n0J4xZvE7oipoDpOzjcfcZwNelnLoVuIwNS4b0AR7y4E2giZntABwGTHH3r939G2AKpSSg2mbTTeEf/4D//S+sv0jRk5jzioqKKCoqijsMESlHOiXK25RytK3Ml5lZH+Bzdy85n6c5sCTpfVHUVlZ7afceaGazzWz28uXLKxNeXvnNb8ITxtNPwz33pL5eRKQq0ilRXt/MzjezSdFxrpnVr+gXmVkj4M/AVZUJNBV3H+3uBe5e0KxZs2x8Rc45/3zo0QMuugg+/DDuaESkJkunG2oUsA8wMjr2idoqaiegDfCOmX0CtADeMrPtgc+BlknXtojaymoXoE4d6NQJ6tcP02nXrAnthYVw003xxiYiNUs6yaKTu/d395ei41SgU0W/yN3fc/dt3b21u7cmdCnt7e5fAJOBU6JZUV2Ab919KfAC0MPMtooGtntEbRI55JDw61tvhf27CwvhuONCEskX++23H/vtt1/cYYhIOdJZZ/EWcKy7/yd63xaY5O57p/jcOKAr0BT4Ehjq7vcnnf8EKHD3FRY2y7iTMHi9CjjV3WdH1/2J0H0F8Fd3fyDVD1WT11mUprAQevcO9aO22AL++U/o1i3uqEQk35S3ziKdzY8uBQrNbBFghJXcp6b6kLufkOJ866TXTtiRr7TrxgBj0oiz1urWLay9uPFGWL26ZlamFZF4pTMbahph/cP5wHlAB3cvzHZgkr7CQrj//jDgvWYNdO8eptXmi759+9K3b9+4wxCRcpT5ZGFmB7v7S2Z2TIlTO5sZ7v5ElmOTNCTGKCZODE8YbdvCoEFw6KHw5ptQt27cEab21VdfxR2CiKRQXjfUQcBLwBGlnHNAySIHzJpVnCggdEd99hncckvYZW/48HjjE5Gaocxk4e5Do5fXuvvi5HNm1iarUUnaLrts47YRI0J31IgRYfFe/1pfUUtEqiqdqbOPl9I2KdOBSGbdcgscfDAMHBi6o0REqqK8MYtdgF2BLUuMW2wBbJrtwKRq6tcP3VOdO8PRR4fuqhYt4o6qdIckFouISM4qb8yiA3A40IQNxy2+B07PZlCSGdtsA5MnQ5cucNRR8Morublh0pVXXhl3CCKSQnljFk8BT5nZfu7+RjXGJBm0666hQm2fPnDaafDoo2HXPRGRikhnzOJMM2uSeBOV3tAiuTxyxBHw17/CuHFh4V6u6dWrF7169Yo7DBEpRzoruHd391+WeLn7N2a2VxZjkiwYPBjefRf+/OfwtHFEaROiY/Ljjz/GHYKIpJDOk0Wd5N3pot3r0kkykkPMwirvvfeGP/4R5s2LOyIRySfpJIsRwBtmNszMriPsmKcC2HmoUaNQZLBxYzjySNDCaRFJVzq1oR4CjiFUjv0COMbdH852YJIdLVrAk09CUREceyz8/HPcEYlIPkjZnWRmrYAfCHtO/NLm7p9lMzDJnv32g9GjYcCAsMve3/8ebzyHH354vAGISErpjD08S6gFBdCQsNvdAsKCPclT/fuHAe9bbgklQQYOjC+WSy65JL4vF5G0pEwW7v6b5PdmtjdwdtYikmpz441hoPucc2CXXeDAA+OOSERyVToD3Btw97eAfbMQi1SzevVg/PhQ1rxvX/j003ji6Nq1K127do3ny0UkLemMWVyU9LYOsDfw36xFJNWqSZNQEmTffcMMqddeg802izsqEck16TxZbJ50bEIYw+iTzaCkenXoABMmwPvvh0Hv9evjjkhEck25ycLM6gKbu/s10fFXd3/U3X9KdWMzG2Nmy8zs/aS2YWb2rpnNNbMXzexXUbuZ2R1mtjA6v3fSZ/qb2cfRoZ0ZsuSww+Dmm+Hxx2HYsLijEZFcU26ycPd1wP6VvPdYoGeJtpvdfXd33xN4Brgqau9F2Oe7HTAQGAW/rBYfShgj6QwMTV5NLpl14YVhltTVV4ekISKSkM7U2blmNhl4DFiZaEy1B7e7zzCz1iXavkt625jiKbl9gIfc3YE3zayJme0AdAWmuPvXAGY2hZCAxqURt1SQGdx9NyxYAKecAjvvDHvskf3vPe6447L/JSJSJekki02Br4CDk9oqvQe3mf0VOAX4Foh2jqY5sCTpsqKoraz20u47kPBUQqtWrSoTmgCbbgpPPAGdOoUB71mzYNtts/udZ5+tmdgiuS6dAe773P3U5AO4v7Jf6O5/cfeWwKPAuZW9Tyn3He3uBe5e0KxZs0zdtlbaYYdQQ2rZMvjDH8J+3tm0atUqVq1ald0vEZEqSSdZlFYMIhMFIh4F+kavPwdaJp1rEbWV1S5ZVlAAY8aE3fWOPhrci88VFsJNGSwl2bt3b3r37p25G4pIxpWZLMxsPzO7GGhmZhclHVcDdSvzZWbWLultH+DD6PVk4JRoVlQX4Ft3Xwq8APSINlzaCugRtUk1OOGEUM78uefgggtCW2EhHHdc6KYSkdqjvDGLBsBm0TWbJ7V/B/wh1Y3NbBxhgLqpmRURZjX1NrMOwHrgU+DM6PLngN7AQmAVcCqAu39tZsOAWdF11yYGu6V6PPwwLF4cig2uWAFTpsDEidCtW+rPikjNYZ7cv1DaBWY7uvun0es6wGYlZjXlnIKCAp89e3bcYdQY330H7dqFMYyuXWHqVKhbqWfL0iVKfUyfPj1zNxWRCjOzOe5eUNq5dMYsrjezLcysMfA+MN/MLs1ohJLT5syBdetgt91g+vTQBbV0adxRiUh1SidZdIyeJI4CnieUKD85q1FJzkiMUTz2WChpfuml8Pbb8OtfwwsZGj0aMGAAAwYMyMzNRCQr0kkW9c2sPiFZTHb3nyleTCc13KxZxWMUZmEW1JgxYT1Gz54weHDVd9tTshDJfekki3uATwgrrmeY2Y6EQW6pBS67bOPB7FNPDYPeAweGPTEOOqhq5c1XrFjBihUrqhaoiGRVygHuUj9kVs/d12YhnozQAHf1mTAhJI06deD+++GYYyp+Dw1wi+SG8ga409nPYhPC4rnWJa6/NiPRSV47/vgw4H388WEDpXPOgeHDQzeViNQc6XRDPUVYQLeWUEgwcYgAYae9116Diy6Cu+6CLl3go4/ijkpEMimdQoIt3L1kqXGRDTRoACNGhPGNAQNg771h1Cg4WfPmRGqEdJ4sXjez32Q9EqkRDj8c5s6FffYJZc4HDIAffog7KhGpqnSeLA4ABpjZYmA1YIC7++5ZjUzyVosWMG1a2HFv2DB4880w/Xb3Mv7EnHXWWdUboIhUWFrlPkprT5QAyUWaDZU7XnoJTjwRvvkGbr0VzjwzrNcQkdxTpXIfUVJoAhwRHU1yOVFIbjn4YHjnnVBT6uyz4dhj4X//2/CaJUuWsGTJklI/LyK5IWWyMLMLCHtPbBsdj5jZedkOTGqObbcNZc5vvBGeegr22gv+/e/i8yeffDInayRcJKelM8B9GrCvu1/l7lcBXYDTsxuW1DR16oTV4K+8EjZSOuAAuPlmWL8+7shEJB3pJAsD1iW9Xxe1iVRYly6hEOGRR4bksd9+sGbNlr+cz/QufCKSGekkiweAf5vZ1dEueW9ShT24RbbaCiZNCrvvzZwJs2Y9xDff7KVd+ERyWFq1ocxsb8IUWoBX3P3trEZVRZoNlT/uvRcGDlwH1KV+fbj++rASXDOmRKpflWZDRXtif+zud7j7HcB/zGzfTAcptdPpp0PfvouAsAr8kkugoADGj4e1OVuqUqT2SacbahSQvAb3h6hNpMoKC+Hll9tx5ZXQsCFcfHFY8X3CCWEr1zvvhFWr4o5SRNIa4Pakvip3X0961WrHmNkyM3s/qe1mM/vQzN41syfNrEnSuSFmttDMFpjZYUntPaO2hWY2OP0fTXJdYoxi+PDPOPHEBUycCA8+CCNHwpNPwvbbw3nnQatWcPXVoC0vROKTTrJYZGbnm1n96LgAWJTG58YCJQsQTgF2i0qFfAQMATCzjkA/YNfoMyPNrK6Z1QXuAnoBHYETomulBkjswvfAA6dwxhln0K1beD9nDhx1FLz+ephq+9vfwjXXhKRx7rmwKJ0/fSKSUekkizOB3wKfA0XAvsDAVB9y9xnA1yXaXkzaNOlNoEX0ug8w3t1Xu/tiYCHQOToWuvsid18DjI+ulRqgtF34unUL7RAGuQ84ACZPhnnzQtfU6NGhe6pfv5BURKR6pFPuY5m793P3bd19O3f/o7svy8B3/wl4PnrdHEiu91AUtZXVvhEzG2hms81s9vLlyzMQnuSSjh3DTnyLF4dB8OefDwPh3bvDiy+GhX4ikj3pzIZqb2bTEmMPZra7mV1RlS81s78QNlN6tCr3Sebuo929wN0LmjVrlqnbSo5p3jyUDfnss7B4b/58OOywUELkH//QDCqRbEmnG+pewtjCzwDu/i5hfKFSzGwAcDhwYtLA+edAy6TLWkRtZbVLLbfllnDppeFJY8wYWLMmVLfdeWe44w5YuTIkk8LCDT+nFeIilZNOsmjk7jNLtFXq329m1hO4DDjS3ZMnRE4G+pnZJmbWBmgHzARmAe3MrI2ZNSAkqcmV+W7JXVdccQVXXFG5h9VNNoFTT4X33w9jGy1bhpXhrVqFtmOPLU4YWiEuUnnpbH60wsx2AhzAzP4ALE31ITMbB3QFmppZETCU8ISyCTDFwhLdN939THefZ2YTgfmERHSOu6+L7nMu8AJQFxjj7vMq9iNKruvevXuV71GnDhxxRDhefz0UKXzkEahXD3r3htNOgwkTwmyrkoPqIpJaOpsftQVGE2ZEfQMsJnQh5eyeFir3kV/mzp0LwJ577pnR+y5YAMOHh26q9ethl11g7FjYV/UHREpV1c2PFrl7d6AZsAtwEMV1okSqbNCgQQwaNCjj9+3QAf74R2jSBPbfPySPLl3gwAPhmWdUHl2kIspMFma2RbSq+k4zOxRYBfQnrIE4rroCFKmsxBjFpEnw6qvw9NPQuHFIGkccAbvtFp46Vq+OO1KR3Ffek8XDQAfgPcJmR4XAscDR7q6FcZLzEivEE2MUv/99SBiDBoXxjAYNwlhGmzZwww0bb/cqIsXKHLMws/fc/TfR67qEQe1W7v5TNcZXKRqzyC9du3YFYPr06dX6ve4wdWoYDJ8yBTbbDAYODMmkZcvUnxepaSo7ZvFz4kU0M6koHxKFSLrM4NBDwwrwt9+GPn3g9tuhbVs4+WR45524IxTJHeU9WawDVibeAg0J4xYGuLtvUS0RVoKeLPLL66+/DsBvf/vbmCMJK8Nvuy3UoFq5Enr0CIv/DjlEGzJJzVfek0VaO+XlGyULqapvvoG77w6rwb/4IpQTufTSsMivXjqrk0TyUJWmzopk2+uvv/7L00Wu2GorGDIEPvkE7rsPfvwxTMPdeefQVTVsmEqJSO2iJwuJXVwD3BWxfj08+2xIBq++GgbD162Dhx+Gvn2Lp+lqhbjkMz1ZiFRRopzIK6/AG2+EsYwff4Q//CF0UR1zTCgnokQhNZWShUgFdekCjz8OH30E++wDc+eGNRrnnx8Gx7X9q9REShYilVRUBJ9+Gnb2S3RLXXhh2HPj+OPD2g2VFJGaQslCpBKSxyhuvDGUR1+xIuzmd/bZYbFfjx5hzca118KSJanvKZLLlCwkdrfddhu33XZb3GFUSMlSIt26hfcrVsCtt8J//wvjx0P79jB0KOy4YyiV/vjjYaMmkXyj2VAiWbZ4MTzwQDiKiqBZM+jfP9Sl2mWXuKMTKabZUJLTpk6dytSpU+MOI2vatAldUZ98As89B7/7XRgI//Wv4YADwh4bK1emuotIvJQsJHbXXXcd1113XdxhZF3dutCrV+iKKioKBQxXrAjbwu6wA5xxRujectf+4ZJ7lCxEYrDddnDJJfDBB2HtRt++oWx6586w555hQFz7h0suUbIQiZFZ6Ip64AFYuhTuuQc22QTuvBO++y48iRx1VFj8p0V/EiclC5EcscUWYT+NmTPDQr8zzwwrx596Cr7+GgYMgJNOCgnlgw9Cd5VIdclasjCzMWa2zMzeT2o71szmmdl6Mysocf0QM1toZgvM7LCk9p5R20IzG5yteEVyyR57wNFHh21gzzgjLPrbaSeYNi0kkY4dQ1dW376hsOHbb4dFgSLZks1iy2OBO4GHktreB44B7km+0Mw6Av2AXYFfAVPNrH10+i7gUKAImGVmk919fhbjlmp2zz33pL6olilZmPD448P7CbdvFuEAAA74SURBVBPCLn4zZoTjlVfgiSfCZ7bYInRpHXhgmHFVUBC2jhXJhKwlC3efYWatS7R9AGAb7yLTBxjv7quBxWa2EOgcnVvo7ouiz42PrlWyqEE6dOgQdwg5p6xFf7NmwcEHQ7t2YZ0GhMHwV14pTiDPPRfaGzYMdawOPDAcXbpAo0ZhRlWnThuOfxQWhntfdln1/pySP3JlG5fmwJtJ74uiNoAlJdr3Le0GZjYQGAjQqlWrLIQo2fL0008DcMQRR8QcSe4o7S/tbt1KH+Bu2TLstfHHP4b3y5aFMuqJ5HHttWF8o169kCRatYK//jXMvjriiA2fYkTKkivJosrcfTQwGsIK7pjDkQoYMWIEoGSRKdtuG0qmH3NMeP/tt/Daa8VPH088AT//DEceCdtvH85fdZWm5Ur5cmU21OdAy6T3LaK2stpFJE1bbhnqUl1/fUga//sfvPRS6Jr64gtYvTrsCrjVVmGsY+hQePnl0C6SkCvJYjLQz8w2MbM2QDtgJjALaGdmbcysAWEQfHKMcYrkvUaNwq/z58OVV8LWW4dxjIsvDgniuuuga9eQPHr0CFV1Z83SbKvaLmvdUGY2DugKNDWzImAo8DXwd6AZ8KyZzXX3w9x9nplNJAxcrwXOcfd10X3OBV4A6gJj3H1etmIWqQ1KzrTq1q34/Q03hCePl18O03RfegkGRxPWt9wyJJFDDgmD7B07hkWFUjuo6qzELh/24K5JKjob6osvwvlE8li8OLRvt11IGonk0aaNZlrlu/KqzipZSOyWRDsDtWzZMsWVkgsWLw5JI5E8vvwytLdpE0quv/pqKF/St+/GTzGS25QsRCQr3MPYRyJ5TJ8eZlcBbLMN/PADnHACHH546LbaeWeoXz/WkKUcShaS0yZMmADA8ccfH3MkUlXr1oXSI5dfHhJIkyYheST+mqlfHzp0CIlj113DoSSSO8pLFjVmnYXkr1GjRgFKFjVB3brw/ffw7rthptWoUfDss2F8Y968cMyfD7Nnw2OPVTyJaEwkPkoWIpIx5c20OvnkDa9dtSpUz50/vziRpEoi7mEsZMIEOPRQrT6vTkoWIpIx5dW0KjnA3agR7LNPOJKlk0R69AjrQ77/Hnr2DPdfsSIMsrdpE85pWm9macxCYqeps5KOlSvhww9DErn7bnj99bAd7erVYb+PZJtvHpJG27bFCSRxtG4dSr+XpC4ujVmISA3QuHF4CvnuO/joo+IxkUmTYO+94ZNPwrTeRYvCr4sXh+teeAF+/HHDe2277cZJpE4ddXGVR08WErsVK1YA0LRp05gjkVxXckwknXUc7qESbyKBlDw+/XTjUiZbbhm6w7p3h/33L35CadsWmjWruV1cmjorIjVCNrqK1q6FoqLi5PHAA2FhYfPmIYl88cWG1zduXJw4kpNI27ahiytRe6s6Ys80dUNJThs7diwAAwYMiDUOyX0V2ecjXfXqhb/kW7cO7z/8sLiLa+JE6Nx5wy6uRDfXokVhIeLKlRveb/vtS08mO+5Y9lNRPtCThcROA9ySCyrbxbV8+YYJJPn1kiWwfn3x9XXrhl933BH++98wnbh795Co2rSBpk3j7eLSk4WISAoVmfabYBYGy7fdNmxbW9KaNSFhJCeRyZPD1OCGDeHee8ORsNlmxYkjefZW4vUWW5QeR3V0cSlZiIiQnS6uBg1gp53CAeEv8PvvL+7meuyxsC3u4sXFXV2Jo7Aw1NZKtvXWpU8F3m677HdxKVmIiFSD8la39+mz8fXu8NVXxckjOZm8+254QlmzZsPPdO8eVrt/8UXmK/0qWYiIVIOKdnOZhTGMpk1L3x99/XpYunTDRPL44/DOO+HJJdMl4TXALbFbtWoVAI3KmnMoIiklnlzOOqt4JldFE4YGuCWnKUmIVE15XVyZesKok5nbiFTeyJEjGTlyZNxhiOSt8rq4MiVr3VBmNgY4HFjm7rtFbVsDE4DWwCfAce7+jZkZcDvQG1gFDHD3t6LP9AeuiG57nbs/mOq71Q2VX7TOQiQ3lNcNlc0ni7FAzxJtg4Fp7t4OmBa9B+gFtIuOgcAo+CW5DAX2BToDQ81sqyzGLCIipchasnD3GUCJwsH0ARJPBg8CRyW1P+TBm0ATM9sBOAyY4u5fu/s3wBQ2TkAiIpJl1T1msZ27L41efwFsF71uDixJuq4oaiurfSNmNtDMZpvZ7OXLl2c2ahGRWi62AW4PgyUZGzBx99HuXuDuBc2aNcvUbUVEhOqfOvulme3g7kujbqZlUfvnQMuk61pEbZ8DXUu0T6+GOKUaaWBbJPdV95PFZKB/9Lo/8FRS+ykWdAG+jbqrXgB6mNlW0cB2j6hNRESqUdaeLMxsHOGpoKmZFRFmNd0ATDSz04BPgeOiy58jTJtdSJg6eyqAu39tZsOAxGzha9295KC5iIhkmcp9iIgIEN86CxERqSGULEREJCUlCxERSUnJQkREUlKyEBGRlJQsREQkpRo5ddbMlhPWcVRWU2BFhsKpTvkaNyj2uCj2eORq7Du6e6n1kmpksqgqM5td1lzjXJavcYNij4tij0c+xq5uKBERSUnJQkREUlKyKN3ouAOopHyNGxR7XBR7PPIudo1ZiIhISnqyEBGRlJQsREQkJSWLiJm1NLNCM5tvZvPM7IK4Y6ooM6trZm+b2TNxx1IRZtbEzCaZ2Ydm9oGZ7Rd3TOkyswujPy/vm9k4M9s07pjKYmZjzGyZmb2f1La1mU0xs4+jX7eKM8aylBH7zdGfmXfN7EkzaxJnjKUpLe6kcxebmZtZ0zhiqygli2JrgYvdvSPQBTjHzDrGHFNFXQB8EHcQlXA78C933wXYgzz5GcysOXA+UODuuwF1gX7xRlWusUDPEm2DgWnu3g6YFr3PRWPZOPYpwG7uvjvwETCkuoNKw1g2jhsza0nY+fOz6g6ospQsIu6+1N3fil5/T/gLq3m8UaXPzFoAvwfuizuWijCzLYEDgfsB3H2Nu/8v3qgqpB7Q0MzqAY2A/8YcT5ncfQZQcqfJPsCD0esHgaOqNag0lRa7u7/o7mujt28CLao9sBTK+D0HuBW4DMibGUZKFqUws9bAXsC/442kQm4j/OFbH3cgFdQGWA48EHWh3WdmjeMOKh3u/jkwnPCvw6WEveNfjDeqCtsu2u8e4AtguziDqYI/Ac/HHUQ6zKwP8Lm7vxN3LBWhZFGCmW0GPA4Mcvfv4o4nHWZ2OLDM3efEHUsl1AP2Bka5+17ASnK3K2QDUf9+H0LC+xXQ2MxOijeqyvMwjz5v/qWbYGZ/IXQjPxp3LKmYWSPgz8BVccdSUUoWScysPiFRPOruT8QdTwXsDxxpZp8A44GDzeyReENKWxFQ5O6Jp7hJhOSRD7oDi919ubv/DDwB/DbmmCrqSzPbASD6dVnM8VSImQ0ADgdO9PxYNLYT4R8X70T/v7YA3jKz7WONKg1KFhEzM0K/+Qfufkvc8VSEuw9x9xbu3powwPqSu+fFv3Dd/QtgiZl1iJoOAebHGFJFfAZ0MbNG0Z+fQ8iTwfkkk4H+0ev+wFMxxlIhZtaT0PV6pLuvijuedLj7e+6+rbu3jv5/LQL2jv4/yGlKFsX2B04m/Kt8bnT0jjuoWuI84FEzexfYE/hbzPGkJXoamgS8BbxH+P8pZ8s4mNk44A2gg5kVmdlpwA3AoWb2MeFJ6YY4YyxLGbHfCWwOTIn+f7071iBLUUbceUnlPkREJCU9WYiISEpKFiIikpKShYiIpKRkISIiKSlZiIhISkoWknOiSpwjkt5fYmZXZ+jeY83sD5m4V4rvOTaqoFuYbjxRqZNqL15pZmea2SnV/b2SX+rFHYBIKVYDx5jZ9e6+Iu5gEsysXlLhulROA05391fTvb+7/1/lIqsad8+59QmSe/RkIbloLWFx24UlT5R8MjCzH6Jfu5rZy2b2lJktMrMbzOxEM5tpZu+Z2U5Jt+luZrPN7KOorlZiL5CbzWxWtD/CGUn3fcXMJlPKynIzOyG6//tmdmPUdhVwAHC/md1c4nozszvNbIGZTQW2TTo33cwKEj9XFM88M5tqZp2j84vM7Mg0Yp5uxXuEPBqtMCf6fZkfXT88arvazC6JXu9pZm9a8R4RWyXFdmP0+/mRmf0uat81apsbfaZd2v+VJa8oWUiuugs40UIJ83TtAZwJ/JqwGr+9u3cmlG0/L+m61kBnQkn3uy1sWHQaoWpsJ6ATcLqZtYmu3xu4wN3bJ3+Zmf0KuBE4mLDyvJOZHeXu1wKzCfWKLi0R49FAB6AjcApl15JqTCjbsivwPXAdcGj0+Wuja8qLeS9gUPQ9bYH9zWyb6PO7RntAXFfK9z4EXB6dfw8YmnSuXvT7OSip/UzgdnffEygglK+QGkjJQnJSVPH3IcLmQumaFe1Lshr4D5AoF/4eIUEkTHT39e7+MbAI2IWwEc0pZjaXUJp+GyDxr+SZ7r64lO/rBEyPCgkmqp4emCLGA4Fx7r7O3f8LvFTGdWuAfyXF/3JUrDD5Z0kVc5G7rwfmRp/5FviJ8MRzDLBBPaUoMTdx95ejpgdL/DyJ4ppzkmJ4A/izmV0O7OjuP6b4+SVPKVlILruN8K/n5P0t1hL9uTWzOkCDpHOrk16vT3q/ng3H50rWuHHAgPPcfc/oaJO0N8XKKv0UlfNzUhXVX36W6C//xM9SXszJvxfrCE8FawlPVJMIlVr/RcUk7rkuEYO7/wM4EvgReM7MDq7gPSVPKFlIznL3r4GJhISR8AmwT/T6SKB+JW59rJnVicYx2gILgBeAsyyUqcfM2lvqTZhmAgeZWVMzqwucALyc4jMzgOOj8YYdgG6ViD+hQjFb2KtlS3d/jjAetEfyeXf/FvgmMR5B6Mor9+cxs7bAIne/g1CxdvfK/jCS2zQbSnLdCODcpPf3Ak+Z2TuEfxlX5l/9nxH+ot8CONPdfzKz+whdK29Fg8HLSbHFqLsvNbPBQCHhX/nPunuqEt9PEsY45kdxvFGJ+BMqGvPmhN+7TaN4Lyrlmv6EcZxGhC66U1PEcBxwspn9TNhpLy8qBkvFqeqsiIikpG4oERFJSclCRERSUrIQEZGUlCxERCQlJQsREUlJyUJERFJSshARkZT+H6mxrG4WLqITAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X7V0FPcUnJYM"
      },
      "source": [
        "# Lower dimensional cell embedding using NMF"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MS6bKTnQhlE7"
      },
      "source": [
        "optDim = kn.knee"
      ],
      "execution_count": 22,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5IdCKmNLnCrd"
      },
      "source": [
        "#Mark the begining of the proposed method\n",
        "t.tic()"
      ],
      "execution_count": 23,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IGuo26QZd847"
      },
      "source": [
        "data_orig = data.copy()"
      ],
      "execution_count": 24,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eu_ZcTpqcciX"
      },
      "source": [
        "model = NMF(n_components=optDim, init='nndsvd')\n",
        "W = model.fit_transform(data)\n",
        "H = model.components_"
      ],
      "execution_count": 25,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "CcpF1opmdIuX",
        "outputId": "edf550b0-fda2-4f51-82ac-3b2db3d5d2f9"
      },
      "source": [
        "print(W.shape)\n",
        "print(H.shape)"
      ],
      "execution_count": 26,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "(1630, 6)\n",
            "(6, 1830)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nvR7jYeGgue8"
      },
      "source": [
        "data = pd.DataFrame(W)"
      ],
      "execution_count": 27,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2C6cH0nzk9Q2"
      },
      "source": [
        "# (A) Encoding labels (B) Removing p% labels (C) Labelling unknown cells as -1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oh1FGVEXMqM2"
      },
      "source": [
        "**Encoding the labels**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HABmOzUnRpCk"
      },
      "source": [
        "from sklearn import preprocessing\n",
        "le = preprocessing.LabelEncoder()\n",
        "y = pd.DataFrame(le.fit_transform(np.ravel(labels)))"
      ],
      "execution_count": 28,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Drc7NELacY-9"
      },
      "source": [
        "y1 = y.copy()"
      ],
      "execution_count": 29,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YFA1oChFrGvo"
      },
      "source": [
        "**Removing p% labels (if no label is missing)**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "J2qdQwaz31n1"
      },
      "source": [
        "p = 0.9   #input parameter p = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 0.9}\n",
        "n = int(p * data.shape[0])    \n",
        "index = random.sample(range(1, data.shape[0]), n)\n",
        "ground_truth = y.copy()"
      ],
      "execution_count": 30,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "BSCLsdJg2KE3",
        "outputId": "810dc206-f00a-41df-a166-26eb0736b0db"
      },
      "source": [
        "data.shape"
      ],
      "execution_count": 31,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(1630, 6)"
            ]
          },
          "metadata": {},
          "execution_count": 31
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "12RDczIPatOF"
      },
      "source": [
        "dataWOmissingLab = data_orig.iloc[~y.index.isin(index)]\n",
        "LabNOTmissing = y.iloc[~y.index.isin(index)]"
      ],
      "execution_count": 32,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QhBTOAgLMloY"
      },
      "source": [
        "**Labelling unknown data as -1**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WU8ItOMOWAFC"
      },
      "source": [
        "y.iloc[index] = -1           # forcibly removing p% of the labels\n",
        "indx = pd.DataFrame(np.arange(0,data.shape[0]))\n",
        "df = pd.concat([data,indx], axis=1)\n",
        "X = pd.concat([df,y], axis=1)"
      ],
      "execution_count": 33,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1l7BHf8aCxPS"
      },
      "source": [
        "#Setting k to the number of unique labels in the data excluding -1 (-1 stands for unknown labels)\n",
        "k = y[0].nunique()-1"
      ],
      "execution_count": 34,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oryceblO3i7v"
      },
      "source": [
        "# Clustering using k-Means"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pd0KBJ1X3-N6"
      },
      "source": [
        "Recursive clustering"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "I3VjKR_bCqzJ"
      },
      "source": [
        "def rec_KMeans(res):\n",
        "    if len(res) == 0:\n",
        "        return\n",
        "    \n",
        "    num_labels = set()\n",
        "    #finding no. of known labels for each cluster\n",
        "    for i in res:\n",
        "        if(i[-1] != -1):\n",
        "          num_labels.add(i[-1])\n",
        "    print(len(num_labels), num_labels)   \n",
        "\n",
        "    if(len(num_labels) == 1):     #if all labelled data belong from one class\n",
        "        final_clusters.append(res)\n",
        "        cluster_label.append(num_labels)\n",
        "        return\n",
        "    \n",
        "    if(len(num_labels) == 0):     #if no labelled data exist\n",
        "        final_clusters.append(res)\n",
        "        cluster_label.append(-1)\n",
        "        return\n",
        "    \n",
        "    #if all labelled data belong from more than one class, then cluster recursively\n",
        "    newk = len(num_labels)\n",
        "    centroid1 = []\n",
        "    rnd = random.sample(range(len(res)), newk)\n",
        "    for j in range(newk):\n",
        "        centroid1.append(res[rnd[j]])\n",
        "        \n",
        "    print('Partioning further...')\n",
        "    my_KMeans(res, centroid1, newk)"
      ],
      "execution_count": 35,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qZlzPnBo4Blc"
      },
      "source": [
        "Passing each partition to recursive k-Means function"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "khU9EUJq9uoT"
      },
      "source": [
        "def rec_KMeans_call(res):\n",
        "    for par in res:\n",
        "        print('...New cluster...')   \n",
        "        rec_KMeans(res[par])"
      ],
      "execution_count": 36,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lLu1EVVT4ZoU"
      },
      "source": [
        "Distance between a given sample and a given centroid"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lSolb_2V4tXc"
      },
      "source": [
        "def distance(a,b):\n",
        "    ans = 0\n",
        "    for i in range(len(a)):\n",
        "        ans += math.pow((a[i] - b[i]),2)\n",
        "    return math.sqrt(ans)"
      ],
      "execution_count": 37,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "R-Wmt2QchgLI"
      },
      "source": [
        "def my_KMeans(data, centroid, k):  #using user-defined KMeans\n",
        "    \n",
        "    #max_iter = 1\n",
        "    #for i in range(max_iter):\n",
        "        \n",
        "    diction = {}    #empty dictionary to store data corresponding to each cluster\n",
        "    for i in range(k):\n",
        "        diction[i]=[]\n",
        "        \n",
        "    for j in range(len(data)):\n",
        "        l = []\n",
        "        for c in centroid:\n",
        "            l.append(distance(data[j][:-2],c[:-2]))      #calculate distance of the sample from each centroid\n",
        "        diction[l.index(min(l))].append(data[j])         #assigning the sample to the nearest centroid\n",
        "    \n",
        "    for c in diction:     #update centroids\n",
        "        if(len(diction[c])!=0):\n",
        "            centroid[c] = np.average(diction[c], axis = 0)\n",
        "    rec_KMeans_call(diction)"
      ],
      "execution_count": 38,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BmRgi6qckASN"
      },
      "source": [
        "cntr = []\n",
        "final_clusters = []\n",
        "cluster_label = []"
      ],
      "execution_count": 39,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YxfW5witQ0MR"
      },
      "source": [
        "#Initializing centroids\n",
        "random.seed(0)\n",
        "rnd = random.sample(range(X.shape[0]), k)\n",
        "for i in range(k):\n",
        "        cntr.append(X.iloc[rnd[i]])"
      ],
      "execution_count": 40,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "npTH-D9YS8vw",
        "outputId": "90514e5c-3fd4-406d-9a96-318dba463812"
      },
      "source": [
        "#initial k-Means call  \n",
        "my_KMeans(X.values, cntr, k)"
      ],
      "execution_count": 41,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "...New cluster...\n",
            "1 {3.0}\n",
            "...New cluster...\n",
            "2 {3.0, 6.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {6.0}\n",
            "...New cluster...\n",
            "1 {3.0}\n",
            "...New cluster...\n",
            "1 {3.0}\n",
            "...New cluster...\n",
            "2 {10.0, 6.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {10.0}\n",
            "...New cluster...\n",
            "2 {10.0, 6.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {10.0, 6.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {10.0, 6.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {10.0, 6.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {6.0}\n",
            "...New cluster...\n",
            "1 {10.0}\n",
            "...New cluster...\n",
            "1 {10.0}\n",
            "...New cluster...\n",
            "1 {10.0}\n",
            "...New cluster...\n",
            "1 {10.0}\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {4.0}\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {5.0}\n",
            "...New cluster...\n",
            "1 {4.0}\n",
            "...New cluster...\n",
            "0 set()\n",
            "...New cluster...\n",
            "0 set()\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {4.0}\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {5.0}\n",
            "...New cluster...\n",
            "1 {4.0}\n",
            "...New cluster...\n",
            "1 {4.0}\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {5.0}\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {4.0, 5.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {5.0}\n",
            "...New cluster...\n",
            "1 {4.0}\n",
            "...New cluster...\n",
            "1 {4.0}\n",
            "...New cluster...\n",
            "2 {1.0, 7.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {7.0}\n",
            "...New cluster...\n",
            "1 {1.0}\n",
            "...New cluster...\n",
            "4 {8.0, 1.0, 2.0, 7.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {2.0, 7.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {2.0, 7.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "2 {2.0, 7.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {2.0}\n",
            "...New cluster...\n",
            "1 {7.0}\n",
            "...New cluster...\n",
            "0 set()\n",
            "...New cluster...\n",
            "1 {2.0}\n",
            "...New cluster...\n",
            "2 {8.0, 7.0}\n",
            "Partioning further...\n",
            "...New cluster...\n",
            "1 {8.0}\n",
            "...New cluster...\n",
            "1 {7.0}\n",
            "...New cluster...\n",
            "1 {2.0}\n",
            "...New cluster...\n",
            "1 {1.0}\n",
            "...New cluster...\n",
            "1 {3.0}\n",
            "...New cluster...\n",
            "1 {4.0}\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HTCGRRhEBoX7"
      },
      "source": [
        "# Results"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PtPaOvVhJzH9"
      },
      "source": [
        "Computing final centroids"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Y7jFoyv4gWK0"
      },
      "source": [
        "def get_final_centroids(final_clusters):\n",
        "    final_centroids = []\n",
        "    for i in range(len(final_clusters)):\n",
        "        fc = pd.DataFrame(final_clusters[i]).iloc[:, :-2]\n",
        "        fc_mean = np.average(fc,axis = 0)\n",
        "        final_centroids.append(fc_mean)\n",
        "    return final_centroids\n"
      ],
      "execution_count": 42,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MMAcmhYTJ45G"
      },
      "source": [
        "final_centroids = get_final_centroids(final_clusters)"
      ],
      "execution_count": 43,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "abofvkMGFEvo"
      },
      "source": [
        "def get_final_results(cluster_label):\n",
        "    no_of_points = []\n",
        "    labelled = []\n",
        "    unlabelled = []\n",
        "    cluster_lbl = []\n",
        "\n",
        "    dtf = [f'cluster{i}' for i in range(1,len(final_clusters) + 1)]\n",
        "    results = pd.DataFrame(dtf, columns = ['Clusters'])\n",
        "    \n",
        "    for i in cluster_label:\n",
        "        if( i != -1):\n",
        "            cluster_lbl.append(int(list(i)[0]))\n",
        "        else:\n",
        "            cluster_lbl.append(i)\n",
        "        \n",
        "    results['Label'] = cluster_lbl\n",
        "    \n",
        "    for i in range(len(final_clusters)):\n",
        "        no_of_points.append(len(final_clusters[i]))  #count of total samples\n",
        "        cl = 0\n",
        "        cu = 0\n",
        "        for j in final_clusters[i]:\n",
        "            if(j[-1] == -1):\n",
        "                cu += 1       #count of unlabelled samples\n",
        "            else:\n",
        "                cl += 1       #count of labelled samples\n",
        "        labelled.append(cl)\n",
        "        unlabelled.append(cu)\n",
        "\n",
        "    results['Total_Count'] = no_of_points\n",
        "    results['Count_ld'] = labelled\n",
        "    results['Count_ud'] = unlabelled\n",
        "    return results"
      ],
      "execution_count": 44,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "J5xmpi7VFIi8",
        "outputId": "e4c752f4-8c3c-43ec-d20a-13fc9691c84a"
      },
      "source": [
        "result = get_final_results(cluster_label)\n",
        "print(result)"
      ],
      "execution_count": 45,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "     Clusters  Label  Total_Count  Count_ld  Count_ud\n",
            "0    cluster1      3           99         9        90\n",
            "1    cluster2      6           85        10        75\n",
            "2    cluster3      3          308        32       276\n",
            "3    cluster4      3          236        25       211\n",
            "4    cluster5     10           83         8        75\n",
            "5    cluster6      6            5         1         4\n",
            "6    cluster7     10           17         1        16\n",
            "7    cluster8     10           20         1        19\n",
            "8    cluster9     10           27         3        24\n",
            "9   cluster10     10           32         3        29\n",
            "10  cluster11      4            4         1         3\n",
            "11  cluster12      5            4         1         3\n",
            "12  cluster13      4            7         1         6\n",
            "13  cluster14     -1            4         0         4\n",
            "14  cluster15     -1            6         0         6\n",
            "15  cluster16      4            9         1         8\n",
            "16  cluster17      5            1         1         0\n",
            "17  cluster18      4            3         1         2\n",
            "18  cluster19      4           37         3        34\n",
            "19  cluster20      5           15         1        14\n",
            "20  cluster21      5            4         1         3\n",
            "21  cluster22      4           25         1        24\n",
            "22  cluster23      4            3         1         2\n",
            "23  cluster24      7          250        25       225\n",
            "24  cluster25      1           43         3        40\n",
            "25  cluster26      2            6         1         5\n",
            "26  cluster27      7            5         1         4\n",
            "27  cluster28     -1            4         0         4\n",
            "28  cluster29      2           16         1        15\n",
            "29  cluster30      8           15         1        14\n",
            "30  cluster31      7           11         1        10\n",
            "31  cluster32      2           61         6        55\n",
            "32  cluster33      1           24         3        21\n",
            "33  cluster34      3           59         5        54\n",
            "34  cluster35      4          102        10        92\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BqFWlR1qvVCO"
      },
      "source": [
        "# Labelling clusters with both labelled and unlabelled cells"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "goZgzvQtfjJ4"
      },
      "source": [
        "pl=[]\n",
        "index=[]\n",
        "\n",
        "for i in range(len(result)):\n",
        "    cls = result['Label'][i]\n",
        "    nud = result['Count_ud'][i]\n",
        "    if cls != -1 and nud > 0:\n",
        "        for j in final_clusters[i]:\n",
        "            if(j[-1] == -1):\n",
        "                j[-1] = cls   #setting the class label\n",
        "                #marking the predicted label and corresponding index position\n",
        "                id = int(j[-2])\n",
        "                index.append(id)     \n",
        "                pl.append(cls)"
      ],
      "execution_count": 46,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U3Q8N2fxGKnj"
      },
      "source": [
        "# Setting weighted radius of clusters as threshold"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AAMpwb35CLqT"
      },
      "source": [
        "rad=[]\n",
        "for i in range(len(final_centroids)):\n",
        "    dist=[]\n",
        "    for j in final_clusters[i]:\n",
        "        dist.append(distance(final_centroids[i],j[:-2]))\n",
        "    rad.append(np.max(dist))"
      ],
      "execution_count": 47,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "a_qE4ma8Kzfd"
      },
      "source": [
        "wrad=0\n",
        "npnts=data.shape[0]\n",
        "for i in range(len(final_clusters)):\n",
        "    wrad += rad[i]*len(final_clusters[i])\n",
        "wrad/=npnts"
      ],
      "execution_count": 48,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E9VcM6nqJ_3S"
      },
      "source": [
        "# Labelling the clusters with all unlabelled cells"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_lP4YswZ6KDV"
      },
      "source": [
        "pos = [i for i in result.loc[result['Label']==-1].index]"
      ],
      "execution_count": 49,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "MGpzeKbP_g-r",
        "outputId": "2fc65093-e399-4f35-9e88-582183ee67ce"
      },
      "source": [
        "d = euclidean_distances(final_centroids, final_centroids)\n",
        "for i in pos:\n",
        "    nc = i\n",
        "    #find nearest labelled cluster\n",
        "    while result.loc[nc,'Label'] == -1:\n",
        "        d[i][nc]=999\n",
        "        nn = np.where(d[i] == np.amin(d[i]))  #finding nearest neighbor of the unlabelled cluster\n",
        "        nc = nn[0][0]\n",
        "\n",
        "    #Setting labels of unlabelled data in unlabelled clusters to the class label of the nearest cluster\n",
        "    #if distance from nearest centroid is less than or equal to threshold distance (weighted radius)\n",
        "    l = len(final_clusters[i])  #number of points in cluster i\n",
        "    print('no. of points in cluster {0}: {1}'.format(i,l))\n",
        "    c=0\n",
        "    for j in final_clusters[i]:\n",
        "        if (distance(final_centroids[nc],j[:-2]) <= wrad):\n",
        "            c+=1\n",
        "            j[-1] = result['Label'][nc]   #setting the class label\n",
        "            #marking the predicted label\n",
        "            id = int(j[-2])\n",
        "            index.append(id)\n",
        "            pl.append(int(j[-1]))\n",
        "    print('no. of points labelled now: ',c)\n",
        "    if c==l:    #implies labels of all samples in cluster i changed\n",
        "        print('Label of cluster {0} changed to the label of its nearest neighbour cluster {1}'.format(i+1,nc))\n",
        "        result.loc[i,'Label'] = result.loc[nc,'Label']"
      ],
      "execution_count": 50,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "no. of points in cluster 13: 4\n",
            "no. of points labelled now:  4\n",
            "Label of cluster 14 changed to the label of its nearest neighbour cluster 17\n",
            "no. of points in cluster 14: 6\n",
            "no. of points labelled now:  6\n",
            "Label of cluster 15 changed to the label of its nearest neighbour cluster 20\n",
            "no. of points in cluster 27: 4\n",
            "no. of points labelled now:  4\n",
            "Label of cluster 28 changed to the label of its nearest neighbour cluster 26\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RaE5-BUUo_Cw"
      },
      "source": [
        "Computing running time"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "0oJc1pdGnZDO",
        "outputId": "9d07f9e4-9c91-49b8-8dd0-9abe629735d6"
      },
      "source": [
        "#Compute time elapsed at the end of proposed algorithm\n",
        "t.toc()"
      ],
      "execution_count": 51,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Elapsed time is 3.409091 seconds.\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dissDtuEP7tN"
      },
      "source": [
        "# Storing true labels"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9KicyVZyDqI1"
      },
      "source": [
        "tl = []\n",
        "for i in index:\n",
        "    tl.append(ground_truth.iloc[i].values[0])"
      ],
      "execution_count": 52,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pHIgvGm92PMt"
      },
      "source": [
        "# Results after reassignment"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "G1alCbK337zp"
      },
      "source": [
        "def get_final_results(result):\n",
        "    labelled = []\n",
        "    unlabelled = []\n",
        "    for i in range(len(final_clusters)):\n",
        "        cl = 0\n",
        "        cu = 0\n",
        "        for j in final_clusters[i]:\n",
        "            if(j[-1] == -1):\n",
        "                cu += 1       #count of unlabelled samples\n",
        "            else:\n",
        "                cl += 1       #count of labelled samples\n",
        "        labelled.append(cl)\n",
        "        unlabelled.append(cu)\n",
        "\n",
        "    result['Count_ld_Now'] = labelled\n",
        "    result['Count_ud_Now'] = unlabelled\n",
        "    return result"
      ],
      "execution_count": 53,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Vlh_td882H32",
        "outputId": "bf0fca1e-d6c4-4506-9d8f-41e0391c244f"
      },
      "source": [
        "reassignmentResult = get_final_results(result)\n",
        "print(reassignmentResult)"
      ],
      "execution_count": 54,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "     Clusters  Label  Total_Count  ...  Count_ud  Count_ld_Now  Count_ud_Now\n",
            "0    cluster1      3           99  ...        90            99             0\n",
            "1    cluster2      6           85  ...        75            85             0\n",
            "2    cluster3      3          308  ...       276           308             0\n",
            "3    cluster4      3          236  ...       211           236             0\n",
            "4    cluster5     10           83  ...        75            83             0\n",
            "5    cluster6      6            5  ...         4             5             0\n",
            "6    cluster7     10           17  ...        16            17             0\n",
            "7    cluster8     10           20  ...        19            20             0\n",
            "8    cluster9     10           27  ...        24            27             0\n",
            "9   cluster10     10           32  ...        29            32             0\n",
            "10  cluster11      4            4  ...         3             4             0\n",
            "11  cluster12      5            4  ...         3             4             0\n",
            "12  cluster13      4            7  ...         6             7             0\n",
            "13  cluster14      4            4  ...         4             4             0\n",
            "14  cluster15      5            6  ...         6             6             0\n",
            "15  cluster16      4            9  ...         8             9             0\n",
            "16  cluster17      5            1  ...         0             1             0\n",
            "17  cluster18      4            3  ...         2             3             0\n",
            "18  cluster19      4           37  ...        34            37             0\n",
            "19  cluster20      5           15  ...        14            15             0\n",
            "20  cluster21      5            4  ...         3             4             0\n",
            "21  cluster22      4           25  ...        24            25             0\n",
            "22  cluster23      4            3  ...         2             3             0\n",
            "23  cluster24      7          250  ...       225           250             0\n",
            "24  cluster25      1           43  ...        40            43             0\n",
            "25  cluster26      2            6  ...         5             6             0\n",
            "26  cluster27      7            5  ...         4             5             0\n",
            "27  cluster28      7            4  ...         4             4             0\n",
            "28  cluster29      2           16  ...        15            16             0\n",
            "29  cluster30      8           15  ...        14            15             0\n",
            "30  cluster31      7           11  ...        10            11             0\n",
            "31  cluster32      2           61  ...        55            61             0\n",
            "32  cluster33      1           24  ...        21            24             0\n",
            "33  cluster34      3           59  ...        54            59             0\n",
            "34  cluster35      4          102  ...        92           102             0\n",
            "\n",
            "[35 rows x 7 columns]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ApHYk-t7d_Hp",
        "outputId": "9ad0d64d-3b60-42b8-f879-0590124bae14"
      },
      "source": [
        "#Finding how many samples are still unlabelled\n",
        "count = 0\n",
        "unlab = []\n",
        "for i in range(len(final_clusters)):\n",
        "  for j in final_clusters[i]:\n",
        "    if(j[-1] == -1):\n",
        "      count+=1\n",
        "      unlab.append(int(j[-2]))\n",
        "print('Count of unlabelled cells: ', count)\n",
        "print('Percentage of unlabelled cells: {0:.2f}'.format(count/data.shape[0]*100))\n",
        "if count>0:\n",
        "  print('Unknown cell indices: ',unlab)"
      ],
      "execution_count": 55,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Count of unlabelled cells:  0\n",
            "Percentage of unlabelled cells: 0.00\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NSnaLIdM9cpF"
      },
      "source": [
        "#Remove unpredicted cells from data\n",
        "index_labelled = ~data_orig.index.isin(unlab)\n",
        "data_labelled = data_orig.iloc[~data_orig.index.isin(unlab)]"
      ],
      "execution_count": 56,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sMcn54eBglRo"
      },
      "source": [
        "labWOmissing = pd.DataFrame(list(zip(LabNOTmissing.index,LabNOTmissing[0].values)), columns = ['Index', 'Predicted Label'])"
      ],
      "execution_count": 57,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kGtd9NyBapEJ"
      },
      "source": [
        "predicted = pd.DataFrame(list(zip(index, pl)), columns = ['Index', 'Predicted Label'])"
      ],
      "execution_count": 58,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SwLtK4X0gcZ2"
      },
      "source": [
        "allLabels = pd.concat([labWOmissing, predicted], axis=0).sort_values(by = 'Index')"
      ],
      "execution_count": 59,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bUm2r74aG-LG"
      },
      "source": [
        "# Validation"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_naU_y8a4cpb"
      },
      "source": [
        "actual = pd.DataFrame(list(zip(index, tl)), columns = ['Index', 'True Label'])"
      ],
      "execution_count": 60,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LNI4xiGmNMU-"
      },
      "source": [
        "merged_inner = pd.merge(left=actual, right=predicted, left_on='Index', right_on='Index')"
      ],
      "execution_count": 61,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8AB6-oumqaHq"
      },
      "source": [
        "**Classification Metrics**\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oSpxxt61OJNX",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "6fc33480-5760-45b1-c23c-7b9eef282748"
      },
      "source": [
        "from sklearn.metrics import confusion_matrix\n",
        "from sklearn.metrics import accuracy_score\n",
        "from sklearn.metrics import classification_report\n",
        "\n",
        "print('Confusion matrix: \\n', confusion_matrix(merged_inner['True Label'],merged_inner['Predicted Label']))\n",
        "\n",
        "print('Accuracy: {0:.2f}'.format(accuracy_score(merged_inner['True Label'],merged_inner['Predicted Label'])))\n",
        "\n",
        "print('Classification_report: \\n', classification_report(merged_inner['True Label'],merged_inner['Predicted Label']))"
      ],
      "execution_count": 62,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Confusion matrix: \n",
            " [[  0   1   0   0   0   0   0   0   3   0   0   0   0]\n",
            " [  0  30   0   0   0   0   0   0   8   0   0   0   0]\n",
            " [  0   1  74   1   0   0   0   0   0   0   0   0   0]\n",
            " [  0   0   0 615   0   0   4   0   0   0   3   0   0]\n",
            " [  0   0   0   0 167   3   0   0   0   0   0   0   0]\n",
            " [  0   0   0   7   5  15   0   0   0   0   0   0   0]\n",
            " [  0   0   0   0   0   8  71   0   0   0   0   0   0]\n",
            " [  0  13   1   1   0   0   0 243   2   0   0   0   0]\n",
            " [  0  12   0   0   0   0   0   0   1   0   0   0   0]\n",
            " [  0   0   0   6   0   0   0   0   0   0   0   0   0]\n",
            " [  0   0   0   0   0   0   4   0   0   0 160   0   0]\n",
            " [  0   4   0   0   0   0   0   0   0   0   0   0   0]\n",
            " [  0   0   0   1   3   0   0   0   0   0   0   0   0]]\n",
            "Accuracy: 0.94\n",
            "Classification_report: \n",
            "               precision    recall  f1-score   support\n",
            "\n",
            "           0       0.00      0.00      0.00         4\n",
            "           1       0.49      0.79      0.61        38\n",
            "           2       0.99      0.97      0.98        76\n",
            "           3       0.97      0.99      0.98       622\n",
            "           4       0.95      0.98      0.97       170\n",
            "           5       0.58      0.56      0.57        27\n",
            "           6       0.90      0.90      0.90        79\n",
            "           7       1.00      0.93      0.97       260\n",
            "           8       0.07      0.08      0.07        13\n",
            "           9       0.00      0.00      0.00         6\n",
            "          10       0.98      0.98      0.98       164\n",
            "          11       0.00      0.00      0.00         4\n",
            "          12       0.00      0.00      0.00         4\n",
            "\n",
            "    accuracy                           0.94      1467\n",
            "   macro avg       0.53      0.55      0.54      1467\n",
            "weighted avg       0.93      0.94      0.93      1467\n",
            "\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/sklearn/metrics/_classification.py:1272: UndefinedMetricWarning: Precision and F-score are ill-defined and being set to 0.0 in labels with no predicted samples. Use `zero_division` parameter to control this behavior.\n",
            "  _warn_prf(average, modifier, msg_start, len(result))\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9gKizoZWYjzb"
      },
      "source": [
        "#set(merged_inner['True Label']) - set(merged_inner['Predicted Label'])"
      ],
      "execution_count": 63,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6Nt2qaRcqqQF"
      },
      "source": [
        "**Comparing clustering results (Predicted Label verified against ground truth)**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gBWuvzz9Ylj6",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "2c89a967-d943-4e57-b6af-2b39141ad4b2"
      },
      "source": [
        "print('Adjusted Rand Score: {0:.2f}'.format(adjusted_rand_score(merged_inner['True Label'].values,merged_inner['Predicted Label'].values)))\n",
        "print('Adjusted Mutual Information Score: {0:.2f}'.format(adjusted_mutual_info_score(merged_inner['True Label'].values,merged_inner['Predicted Label'].values)))\n",
        "print('Fowlkes-Mallows score: {0:.2f}'.format(fowlkes_mallows_score(merged_inner['True Label'].values,merged_inner['Predicted Label'].values)))"
      ],
      "execution_count": 64,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Adjusted Rand Score: 0.94\n",
            "Adjusted Mutual Information Score: 0.89\n",
            "Fowlkes-Mallows score: 0.95\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iWiPdGtCqVvq"
      },
      "source": [
        "**Cohesion and separation of clusters with newly predicted predicted cells added**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tFTdARfsg54C",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "1d3498d0-4912-48b2-8438-b37323c21783"
      },
      "source": [
        "#for data with newly predicted cells added\n",
        "print('Internal indices for data including newly predicted cells (After Clustering):  \\n')\n",
        "print('Calinski Harabasz Score: {0:.2f}'.format(calinski_harabasz_score(data_labelled, allLabels['Predicted Label'].values)))"
      ],
      "execution_count": 67,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Internal indices for data including newly predicted cells (After Clustering):  \n",
            "\n",
            "Calinski Harabasz Score: 296.14\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dlpanaNnFPI-"
      },
      "source": [
        "# -----------THE END----------"
      ]
    }
  ]
}