2023/seminars/seminar10/sem10_clustering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7nDp0NrBscGu"
   },
   "source": [
    "# Кластеризация\n",
    "\n",
    "Кластеризация &mdash; это метод машинного обучения, который включает группировку данных в пространстве признаков. Теоретически, точки, находящиеся в одной группе (кластере), должны иметь схожие свойства, в то время как точки в разных группах должны иметь (сильно) отличающиеся свойства. \n",
    "\n",
    "Кластеризация является методом обучения без учителя (вопрос: чем отличаются методы обучения с учителем и без учителя?) и распространенным методом статистического анализа данных, используемым во многих областях. В частности используется при составлении портретов пользователей, поиске аномалий, анализе геоданных, представлении рыночной информации.\n",
    "\n",
    "В анализе данных часто прибегают к кластеризации, чтобы получить ценную информацию из данных, наблюдая, в какие группы попадают точки при применении алгоритма кластеризации.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KL5ExmINEJUF"
   },
   "source": [
    "# K-means\n",
    "\n",
    "Напомним, что сам алгоритм можно схематически представить в виде следующих шагов:\n",
    "\n",
    "1. Инициализируем центры кластеров случайным образом (должно быть задано количество кластеров).\n",
    "2. Относим точки к соответствующим кластерам (с минимальным расстоянием до их центра).\n",
    "3. Производится пересчет центров кластеров по формуле центра масс для всех точек, принадлежащих кластеру.\n",
    "4. Пункты 2-3 повторяются до тех пор пока центры кластеров перестанут меняться (сильно)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fj3oLr6CscGx"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "J4Y7Cw6AscGy"
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [11, 8]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "k7FPBb8VscGz"
   },
   "source": [
    "Посмотрим на то, как работает метод."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "xqUP7QlFscGz"
   },
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "\n",
    "np.random.seed(123)\n",
    "X1 = np.random.randn(100,2)\n",
    "X2 = np.random.randn(100,2) - np.array([10,1])\n",
    "X3 = np.random.randn(100,2) - np.array([1,10])\n",
    "X = np.vstack((X1,X2,X3))\n",
    "y = np.array([1]*100 + [2]*100 + [3]*100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 699
    },
    "id": "iZ8HlxCwscG0",
    "outputId": "090ecd7f-ae7e-4a46-a525-6598b19a0718"
   },
   "outputs": [],
   "source": [
    "k_means = KMeans(n_clusters = 3)\n",
    "k_means = k_means.fit(X)\n",
    "clusters = k_means.predict(X)\n",
    "plt.scatter(X[:,0], X[:,1], c = clusters)\n",
    "plt.title('K-means clustering')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dQ_KZmy7scG1"
   },
   "source": [
    "Посмотрим, что будет происходить, если мы не угадали с числом кластеров."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 699
    },
    "id": "dEwj1QjYscG1",
    "outputId": "c3f34e71-ea29-4baf-befd-072eef20b585"
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize= (15,8))\n",
    "for n_c in range(2,8):\n",
    "    k_means = KMeans(n_clusters = n_c)\n",
    "    k_means = k_means.fit(X)\n",
    "    clusters = k_means.predict(X)\n",
    "    plt.subplot(2,3,n_c - 1)\n",
    "    plt.scatter(X[:,0], X[:,1], c = clusters)\n",
    "    plt.title('n_clusters = {}'.format(n_c))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2YXIZZn5scG2"
   },
   "source": [
    "Как мы видим, k-means обязательно пытается отдать каждому кластеру какие-то объекты и, как большинство алгоритмов кластеризации, зависит от заданного числа кластеров. Есть огромное количество вариаций, как выбирать количество кластеров автоматически — например, ввести вероятностный подход к выбору числа кластеров, но данные методы мы не будем рассматривать в этом курсе. \n",
    "\n",
    "Одним из главных недостатков k-means является случайная инициализация центров кластеров, что может привести к различным результатам кластеризации.\n",
    "\n",
    "Главным же достоинством является скорость алгоритма. На каждой итерации требуется пересчет только расстояний до центров кластеров."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hw5KSPFSscG2"
   },
   "source": [
    "Также есть вариация k-medians, которая использует медиану вместо среднего. Это позволяет алгоритму стать более устойчивым к выбросам."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sYqP-iN8scG3"
   },
   "source": [
    "[Визуализация работы K-Means](https://www.naftaliharris.com/blog/visualizing-k-means-clustering/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vHZuQCXGD9LK"
   },
   "source": [
    "# Оценка качества кластеризации\n",
    "\n",
    "В случае, когда неизвестны истинные лейблы (внешняя информация), могут быть использованы внутренние меры качества кластеризации, которые оценивают качество структуры кластеров: \n",
    "\n",
    "* компактность (cohesion) - чем ближе (более похожи) объекты в кластере, тем лучше\n",
    "* отделимость (separation) - чем дальше друг от друга (менее похожи)  находятся объекты из разных кластеров, тем лучше"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6m0hUX8sGD1X"
   },
   "source": [
    "Например, **индекс Силуэта** (Silhouette Coefficient):\n",
    "\n",
    "$$Sil(X, C) = \\frac{1}{|C|}\\sum_{c_k \\in C}\\frac{1}{|c_k|}\\sum_{x_i \\in c_k}\\frac{separation(x_i, c_k) - cohension(x_i, c_k)}{max\\{separation(x_i, c_k) , cohension(x_i, c_k)\\}}$$\n",
    "\n",
    "где\n",
    "\n",
    "$separation(x_i, c_k) = min_{c_l \\in C ∖ \\{c_k\\}} \\{\\frac{1}{|c_l|}\\sum_{x_j \\in c_l}||x_i - x_j||\\}$ - среднее расстояние до объектов ближайшего кластера\n",
    "\n",
    "$cohension(x_i, c_k) = \\frac{1}{|c_k|-1}\\sum_{x_j \\in c_k}||x_i - x_j||$ - среднее расстояние до объектов того же кластера\n",
    "\n",
    "\n",
    "$-1 \\le Sil(X, C) \\le 1$\n",
    "\n",
    "Чем ближе к 1, тем лучше. Значения около 0 - пересекающиеся кластера. Отрицательные - обычно объектам присовены не те кластеры, так как объекты других кластеров более похожи, чем объекты из одного."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "4LPXx0sQCS9P",
    "outputId": "2f398a6f-b95a-4cf9-db23-f3aa669e7ba9"
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import silhouette_score\n",
    "\n",
    "best_k, best_score = None, -1\n",
    "for k in range(2,15):\n",
    "    k_means = KMeans(n_clusters = k)\n",
    "    k_means = k_means.fit(X)\n",
    "    clusters = k_means.predict(X)\n",
    "    score = np.round(silhouette_score(X=X,\n",
    "                             labels=clusters), 2)\n",
    "    if score > best_score:\n",
    "      best_score = score\n",
    "      best_k = k\n",
    "print('Best score {}, k = {}'.format(best_score, best_k))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zSa9dnj_scG3"
   },
   "source": [
    "# DBSCAN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zNv6x4rXscG3"
   },
   "source": [
    "(Density-based spatial clustering of applications with noise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IYJVg1n2scG4"
   },
   "source": [
    "Это алгоритм, основанный на плотности — если дан набор объектов в некотором пространстве, алгоритм группирует вместе объекты, которые расположены близко и помечает как выбросы (шум) объекты, которые находятся в областях с малой плотностью (ближайшие соседи которых лежат далеко)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PL4q2InwscG4"
   },
   "source": [
    "Алгоритм имеет два основных гиперпараметра:\n",
    "1. `eps` &mdash; радиус рассматриваемой окрестности\n",
    "2. `min_samples` &mdash; число соседей в окрестности"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gL_GQyUVscG4"
   },
   "source": [
    "Для выполнения кластеризации DBSCAN точки делятся на основные точки, достижимые по плотности точки и выпадающие следующим образом:\n",
    "\n",
    "- Точка $p$ является основной точкой, если по меньшей мере `min_samples` точек находятся на расстоянии, не превосходящем \n",
    "`eps` от неё. Говорят, что эти точки достижимы прямо из $p$.\n",
    "\n",
    "-  Точка $q$ прямо достижима из $p$, если точка $q$ находится на расстоянии, не большем `eps`, от точки $p$, и $p$ — основная точка.\n",
    "Точка $q$ достижима из $p$, если имеется путь \n",
    "$p_1,…,p_n$ где $p_1=p$ и $p_n=q$ , а каждая точка $p_{i+1}$ достижима прямо из $p_i$ (все точки на пути должны быть основными, за исключением $q$).\n",
    "\n",
    "Все точки, не достижимые из основных точек, считаются шумом.\n",
    "\n",
    "Теперь, если $p$ является основной точкой, то она формирует кластер вместе со всеми точками (основными или неосновными), достижимыми из этой точки. Каждый кластер содержит по меньшей мере одну основную точку. Неосновные точки могут быть частью кластера, но они формируют его «край», поскольку не могут быть использованы для достижения других точек.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rYu77WsRscG5"
   },
   "source": [
    "Рассмотрим диаграмму, параметр `min_samples`=4.\n",
    "\n",
    "Точка $A$ и другие красные точки являются основными точками, поскольку область с радиусом \n",
    "`eps` , окружающая эти точки, содержит по меньшей мере 4 точки (включая саму точку). Поскольку все они достижимы друг из друга, точки образуют один кластер. Точки $B$ и $C$ основными не являются, но достижимы из $A$ (через другие основные точки), и также принадлежат кластеру. Точка $N$ является точкой шума, она не является ни основной точкой, ни доступной прямо."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "OEO2I-VkscG5"
   },
   "source": [
    "<p><a href=\"https://commons.wikimedia.org/wiki/File:DBSCAN-Illustration.svg#/media/Файл:DBSCAN-Illustration.svg\"><img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/a/af/DBSCAN-Illustration.svg/1200px-DBSCAN-Illustration.svg.png\" alt=\"DBSCAN-Illustration.svg\" width=\"450\" height=\"450\"> </a><br>Автор: <a href=\"//commons.wikimedia.org/wiki/User:Chire\" title=\"User:Chire\">Chire</a> &mdash; <span class=\"int-own-work\" lang=\"ru\">собственная работа</span>, <a href=\"https://creativecommons.org/licenses/by-sa/3.0\" title=\"Creative Commons Attribution-Share Alike 3.0\">CC BY-SA 3.0</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-6JPcn_cscG6"
   },
   "source": [
    "Посмотрим на результаты кластеризации при разном выборе параметров `eps` и `min_samples`.\n",
    "\n",
    "Шумовые точки помечаются -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "6IPA8WYzsnZm",
    "outputId": "b8c59570-c924-438e-cbad-4fe726323b93"
   },
   "outputs": [],
   "source": [
    "clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "rZVbiA1jscG6",
    "outputId": "9297e13e-fa75-471c-9325-0623ad61fd02"
   },
   "outputs": [],
   "source": [
    "from sklearn.cluster import DBSCAN\n",
    "import matplotlib as mpl\n",
    "\n",
    "def get_colors(clusters):\n",
    "  n_clusters = len(np.unique(clusters[clusters != -1]))\n",
    "  \n",
    "  palette = mpl.colormaps['viridis'].resampled(n_clusters)(np.linspace(0, 1, n_clusters))\n",
    "  colors = np.zeros((clusters.shape[0], 4))\n",
    "  colors[clusters != -1] = palette[clusters[clusters != -1]]\n",
    "  \n",
    "  #Шум будет серо-прозрачным\n",
    "  colors[clusters == -1] = [0.5, 0.5, 0.5, 0.3]\n",
    "  \n",
    "  return colors\n",
    "\n",
    "plt.figure(figsize= (15,20))\n",
    "i = 1\n",
    "\n",
    "for e in [0.2, 1, 3, 5, 10]:\n",
    "    for samples in [2, 5, 40]:\n",
    "        dbscan = DBSCAN(eps=e, min_samples=samples)\n",
    "        clusters = dbscan.fit_predict(X)\n",
    "\n",
    "        n_clusters = len(np.unique(clusters[clusters != -1]))\n",
    "        colors = get_colors(clusters)\n",
    "        \n",
    "        plt.subplot(5, 3, i)\n",
    "\n",
    "        plt.scatter(X[:,0], X[:,1], c=colors)\n",
    "        plt.title('eps = {}, min_samples = {} -> {} clusters'.format(e, samples, n_clusters))\n",
    "        i += 1\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZkNiAtjTscG6"
   },
   "source": [
    "[Визуализация работы DBSCAN](https://www.naftaliharris.com/blog/visualizing-dbscan-clustering/)\n",
    "(Стоит посмотреть)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5kyjmOnwscG7"
   },
   "source": [
    "Результаты работы методов на других датасетах:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "pmnyRfkwscG7"
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_moons, make_circles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "gPU-9FBtscG7"
   },
   "outputs": [],
   "source": [
    "X, y = make_moons(n_samples=500, noise=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 676
    },
    "id": "7dXHIjIkscG7",
    "outputId": "7dc81244-2fbc-4ead-8c4f-08bfc3207a6c"
   },
   "outputs": [],
   "source": [
    "plt.scatter(X[:,0], X[:,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 676
    },
    "id": "LI_g06zqscG8",
    "outputId": "4074bbb5-e28f-4315-d2f4-1a2a778477f5"
   },
   "outputs": [],
   "source": [
    "k_means = KMeans(n_clusters = 2)\n",
    "k_means = k_means.fit(X)\n",
    "clusters = k_means.predict(X)\n",
    "\n",
    "colors = get_colors(clusters)\n",
    "plt.scatter(X[:,0], X[:,1], c=colors)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 699
    },
    "id": "8I175QxQscG8",
    "outputId": "f97cd2b9-8447-4338-f624-0b96c9204bbc"
   },
   "outputs": [],
   "source": [
    "dbscan = DBSCAN(eps=0.2, min_samples=10)\n",
    "clusters = dbscan.fit_predict(X)\n",
    "\n",
    "colors = get_colors(clusters)\n",
    "plt.scatter(X[:,0], X[:,1], c=colors)\n",
    "plt.title('eps = {}, min_samples = {} -> {} clusters'.format(dbscan.eps, dbscan.min_samples, len(np.unique(clusters[clusters != -1]))))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ablfUF7Z5RRS"
   },
   "source": [
    "Как вы думаете, почему так выходит? (Пояснения ниже)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bht8Hg3KscG8"
   },
   "outputs": [],
   "source": [
    "X, y = make_circles(n_samples=500, noise=0.05, factor=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 676
    },
    "id": "EytHhWvCscG9",
    "outputId": "f9ecc245-cc1f-47e2-9dba-85483ff5625b"
   },
   "outputs": [],
   "source": [
    "plt.scatter(X[:,0], X[:,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 676
    },
    "id": "L33MtpupscG9",
    "outputId": "9d92093d-3740-4f10-857d-ce6c4ecc924a"
   },
   "outputs": [],
   "source": [
    "k_means = KMeans(n_clusters = 2)\n",
    "k_means = k_means.fit(X)\n",
    "clusters = k_means.predict(X)\n",
    "plt.scatter(X[:,0], X[:,1], c=clusters)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 699
    },
    "id": "-0F1AHHAscG9",
    "outputId": "328c49fd-7bb8-4461-e253-1648f1ff0db8"
   },
   "outputs": [],
   "source": [
    "dbscan = DBSCAN(eps=0.2, min_samples=10)\n",
    "clusters = dbscan.fit_predict(X)\n",
    "colors = get_colors(clusters)\n",
    "plt.scatter(X[:,0], X[:,1], c=colors)\n",
    "plt.title('eps = {}, min_samples = {} -> {} clusters'.format(dbscan.eps, dbscan.min_samples, len(np.unique(clusters[clusters != -1]))))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0ekVSSmU3Hqf"
   },
   "source": [
    "k-Means выделяет только выпуклые кластеры, сходного размера. Не очень хорошо работает с шумными данными. Хуже видит аномалии. С другой стороны, k-means лучше справляется с большими данными и размерностью, а также разреженными данными с переменной плотностью (DBSCAN скорее относит разные плотности к разным кластерам)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ftKxRB6xscG-"
   },
   "source": [
    "# Иерархическая кластеризация"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eyjbSsXxscG-"
   },
   "source": [
    "Другим вариантом к построению кластеров является иерархический подход, в котором алгоритм жадным образом строит кластеры. Существует два варианта иерархической кластеризации:\n",
    "\n",
    "1. аггломеративная, в которой алгоритм на каждой итерации объединяет два меньших кластера в один\n",
    "2. дивизивная, в которой алгоритм на каждой итерации разбивает один кластер на два более мелких\n",
    "\n",
    "Мы рассмотрим аггломеративный подход к кластеризации (дивизивный можно рассмотреть по аналогии).\n",
    "\n",
    "Опишем схематически алгоритм аггломеративной иерархической кластеризации:\n",
    "\n",
    "- Инициализируем наше множество кластеров, каждая точка считается своим кластером. То есть для выборки размера $N$ у нас на первой итерации будет $N$ кластеров. Также в качестве входного параметра алгоритму подается метрика расстояния между двумя кластерами. Одной из популярных метрик является расстояние Уорда.\n",
    "\n",
    "- На каждой итерации  мы объединяем два кластера в один. Объединяющиеся кластера выбираются в соответствии с наименьшим расстоянием Уорда. То есть в соответствии с выбранным нами расстоянием эти два кластера будут наиболее похожи и поэтому объединяются.\n",
    "\n",
    "- Предыдущий шаг повторяется вплоть до объединения всех точек один кластер.\n",
    "\n",
    "Расстояние Уорда — между кластерами берётся прирост суммы квадратов расстояний объектов до центра кластера, получаемого в результате их объединения.\n",
    "\n",
    "$$ \n",
    "\\Delta = \\sum_{x_i \\in A \\cup B}{(x_i-\\bar{x})^2} - \\sum_{x_i \\in A}(x_i - \\bar{a})^2 - \\sum_{x_i \\in B}(x_i - \\bar{b})^2\n",
    "$$\n",
    "\n",
    "В результате в данном подходе мы можем выбрать любое количество кластеров после завершения процедуры, просто остановив алгоритм на нужном нам шаге. К тому же данный алгоритм гораздо менее чувствителен к выбору метрики между точками, тогда как другие алгоритмы сильно зависят от этого.\n",
    "\n",
    "Для визуализации иерархической кластеризации удобно строить дендрограммы, в которых разница между уровнями равна выбранному расстоянию объединяющихся на данном этапе кластеров.\n",
    "\n",
    "Посмотрим на иерархическую кластеризацию на примере:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "au_BcuLp1oGB",
    "outputId": "0f75f43f-72cf-4996-93f1-cb37e50ed951"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.cluster.hierarchy import dendrogram\n",
    "from sklearn.datasets import load_iris\n",
    "from sklearn.cluster import AgglomerativeClustering\n",
    "\n",
    "\n",
    "def plot_dendrogram(model, **kwargs):\n",
    "    # Create linkage matrix and then plot the dendrogram\n",
    "\n",
    "    # create the counts of samples under each node\n",
    "    counts = np.zeros(model.children_.shape[0])\n",
    "    n_samples = len(model.labels_)\n",
    "    for i, merge in enumerate(model.children_):\n",
    "        current_count = 0\n",
    "        for child_idx in merge:\n",
    "            if child_idx < n_samples:\n",
    "                current_count += 1  # leaf node\n",
    "            else:\n",
    "                current_count += counts[child_idx - n_samples]\n",
    "        counts[i] = current_count\n",
    "\n",
    "    linkage_matrix = np.column_stack(\n",
    "        [model.children_, model.distances_, counts]\n",
    "    ).astype(float)\n",
    "\n",
    "    # Plot the corresponding dendrogram\n",
    "    dendrogram(linkage_matrix, **kwargs)\n",
    "\n",
    "\n",
    "random_state = 170\n",
    "X, _ = make_blobs(n_samples=n_samples, random_state=random_state)\n",
    "\n",
    "# setting distance_threshold=0 ensures we compute the full tree.\n",
    "model = AgglomerativeClustering(distance_threshold=0, n_clusters=None)\n",
    "\n",
    "model = model.fit(X)\n",
    "plt.title(\"Hierarchical Clustering Dendrogram\")\n",
    "# plot the top three levels of the dendrogram\n",
    "plot_dendrogram(model, truncate_mode=\"level\", p=5)\n",
    "plt.xlabel(\"Number of points in node (or index of point if no parenthesis).\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "H6EDC5voscHC"
   },
   "source": [
    "Если нет каких-то специализированных условий (например, известно, что кластеров должно быть не более $K$), то число кластеров можно выбирать по резкому скачку дендрограммы. Кроме того, в некоторых задачах важно понимать, для чего делается кластеризация и доменную область задачи — исходя из этого можно сильно сократить искомое количество кластеров."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4QO-cn2kscHE"
   },
   "source": [
    "# Сравнение работы алгоритмов"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-a0LlHmoscHE"
   },
   "source": [
    "Сгенерируем кластеры разной формы и посмотрим на результаты работы алгоритмов"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "sribQsXiscHF",
    "outputId": "f2467871-5df7-43c5-fdf8-6cbf01c722a6"
   },
   "outputs": [],
   "source": [
    "print(__doc__)\n",
    "\n",
    "import time\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn import cluster, datasets, mixture\n",
    "from sklearn.neighbors import kneighbors_graph\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from itertools import cycle, islice\n",
    "warnings.simplefilter('ignore')\n",
    "\n",
    "np.random.seed(0)\n",
    "\n",
    "# ============\n",
    "# Generate datasets. We choose the size big enough to see the scalability\n",
    "# of the algorithms, but not too big to avoid too long running times\n",
    "# ============\n",
    "n_samples = 1500\n",
    "noisy_circles = datasets.make_circles(n_samples=n_samples, factor=.5,\n",
    "                                      noise=.05)\n",
    "noisy_moons = datasets.make_moons(n_samples=n_samples, noise=.05)\n",
    "blobs = datasets.make_blobs(n_samples=n_samples, random_state=8)\n",
    "no_structure = np.random.rand(n_samples, 2), None\n",
    "\n",
    "# Anisotropicly distributed data\n",
    "random_state = 170\n",
    "X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)\n",
    "transformation = [[0.6, -0.6], [-0.4, 0.8]]\n",
    "X_aniso = np.dot(X, transformation)\n",
    "aniso = (X_aniso, y)\n",
    "\n",
    "# blobs with varied variances\n",
    "varied = datasets.make_blobs(n_samples=n_samples,\n",
    "                             cluster_std=[1.0, 2.5, 0.5],\n",
    "                             random_state=random_state)\n",
    "\n",
    "# ============\n",
    "# Set up cluster parameters\n",
    "# ============\n",
    "plt.figure(figsize=(3 * 2 + 3, 10.5))\n",
    "plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05,\n",
    "                    hspace=.01)\n",
    "\n",
    "plot_num = 1\n",
    "\n",
    "default_base = {'quantile': .3,\n",
    "                'eps': .3,\n",
    "                'damping': .9,\n",
    "                'preference': -200,\n",
    "                'n_neighbors': 10,\n",
    "                'n_clusters': 3,\n",
    "                'min_samples': 20,\n",
    "                'xi': 0.05,\n",
    "                'min_cluster_size': 0.1}\n",
    "\n",
    "datasets = [\n",
    "    (noisy_circles, {'damping': .77, 'preference': -240,\n",
    "                     'quantile': .2, 'n_clusters': 2,\n",
    "                     'min_samples': 20, 'xi': 0.25}),\n",
    "    (noisy_moons, {'damping': .75, 'preference': -220, 'n_clusters': 2}),\n",
    "    (varied, {'eps': .18, 'n_neighbors': 2,\n",
    "              'min_samples': 5, 'xi': 0.035, 'min_cluster_size': .2}),\n",
    "    (aniso, {'eps': .15, 'n_neighbors': 2,\n",
    "             'min_samples': 20, 'xi': 0.1, 'min_cluster_size': .2}),\n",
    "    (blobs, {}),\n",
    "    (no_structure, {})]\n",
    "\n",
    "for i_dataset, (dataset, algo_params) in enumerate(datasets):\n",
    "    # update parameters with dataset-specific values\n",
    "    params = default_base.copy()\n",
    "    params.update(algo_params)\n",
    "\n",
    "    X, y = dataset\n",
    "\n",
    "    # normalize dataset for easier parameter selection\n",
    "    X = StandardScaler().fit_transform(X)\n",
    "    connectivity = kneighbors_graph(\n",
    "        X, n_neighbors=params['n_neighbors'], include_self=False)\n",
    "    \n",
    "    # make connectivity symmetric\n",
    "    connectivity = 0.5 * (connectivity + connectivity.T)\n",
    "\n",
    "    # ============\n",
    "    # Create cluster objects\n",
    "    # ============\n",
    "    \n",
    "    k_means = cluster.KMeans(n_clusters=params['n_clusters'])\n",
    "\n",
    "\n",
    "    dbscan = cluster.DBSCAN(eps=params['eps'])\n",
    "    \n",
    "    average_linkage =  cluster.AgglomerativeClustering(\n",
    "        linkage=\"average\", affinity=\"cityblock\",\n",
    "        n_clusters=params['n_clusters'], connectivity=connectivity)\n",
    "    \n",
    "\n",
    "    clustering_algorithms = (\n",
    "        ('KMeans', k_means),\n",
    "        ('AgglomerativeClustering', average_linkage),\n",
    "        ('DBSCAN', dbscan),\n",
    "    )\n",
    "\n",
    "    for name, algorithm in clustering_algorithms:\n",
    "        t0 = time.time()\n",
    "\n",
    "        algorithm.fit(X)\n",
    "\n",
    "        t1 = time.time()\n",
    "        if hasattr(algorithm, 'labels_'):\n",
    "            y_pred = algorithm.labels_.astype(np.int)\n",
    "        else:\n",
    "            y_pred = algorithm.predict(X)\n",
    "\n",
    "        plt.subplot(len(datasets), len(clustering_algorithms), plot_num)\n",
    "        if i_dataset == 0:\n",
    "            plt.title(name, size=18)\n",
    "\n",
    "        colors = np.array(list(islice(cycle(['#377eb8', '#ff7f00', '#4daf4a',\n",
    "                                             '#f781bf', '#a65628', '#984ea3',\n",
    "                                             '#999999', '#e41a1c', '#dede00']),\n",
    "                                      int(max(y_pred) + 1))))\n",
    "        # add black color for outliers (if any)\n",
    "        colors = np.append(colors, [\"#000000\"])\n",
    "        plt.scatter(X[:, 0], X[:, 1], s=10, color=colors[y_pred])\n",
    "\n",
    "        plt.xlim(-2.5, 2.5)\n",
    "        plt.ylim(-2.5, 2.5)\n",
    "        plt.xticks(())\n",
    "        plt.yticks(())\n",
    "        plt.text(.99, .01, ('%.2fs' % (t1 - t0)).lstrip('0'),\n",
    "                 transform=plt.gca().transAxes, size=15,\n",
    "                 horizontalalignment='right')\n",
    "        plot_num += 1\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MTOCn-Fh5vA7"
   },
   "source": [
    "# Бонус: кластеризуем координаты российских населенных пунктов и яблоки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fi8rLDFE6VDG"
   },
   "outputs": [],
   "source": [
    "import pandas\n",
    "data = pandas.read_csv('https://raw.githubusercontent.com/nadiinchi/HSE_minor_DataAnalysis_seminars_iad16/master/materials/sem_b3_coord.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 206
    },
    "id": "uFv51hD96XPN",
    "outputId": "5170fc88-4aa8-49ef-9966-3766cbda635b"
   },
   "outputs": [],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 519
    },
    "id": "Lg0Nf7YM6bQi",
    "outputId": "0b8672c7-dfcb-4d96-82e5-8ea7785dc5ff"
   },
   "outputs": [],
   "source": [
    "# Просто показываем точки\n",
    "\n",
    "plt.figure(figsize= (30,10))\n",
    "plt.scatter(data['LONG'], data['LAT'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cQYb_fru6yRW"
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import pairwise_distances\n",
    "\n",
    "class k_means:\n",
    "    def __init__(self, n_clusters):\n",
    "        # гиперпараметр - количество кластеров\n",
    "        self.n_clusters = n_clusters\n",
    "        \n",
    "    def fit(self, X, max_iter=10, visualize=False):\n",
    "        \"\"\"\n",
    "        X - выборка размера количество объектов x количество признаков\n",
    "        max_iter - максимальное число итераций\n",
    "        visualize - визуализировать точки (\"карту\") на каждой итерации \n",
    "                    (реализуйте поддержку этого флага, только если уже реализовали алгоритм)\n",
    "        \"\"\"\n",
    "        # Цель: обучить центры кластеров self.centers, форма: число кластеров x число признаков\n",
    "        \n",
    "        # Шаг 0: сколько объектов и признаков?\n",
    "  \n",
    "        ### ╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ\n",
    "\n",
    "        # Шаг 1.1: инициализируйте y - вектор длины число объектов, \n",
    "        # хранящий, к какому кластеру относится каждый объект\n",
    "        # каждый элемент - случайно выбранный номер кластера (используйте np.random.randint)\n",
    "\n",
    "        ### ╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ\n",
    "        \n",
    "        # Шаг 1.2: создайте матрицу для хранения центров кластеров (используйте np.zeros)\n",
    "        \n",
    "        ### ╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ\n",
    "        \n",
    "        for _ in range(max_iter):\n",
    "            # Шаг 2: обновите центры кластеров (используйте цикл по номерам кластеров, индексацию, np.mean)\n",
    "\n",
    "            ### ╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ\n",
    "            \n",
    "            # Шаг 3: обновите y. Для этого вычислите расстояния между всеми объектами \n",
    "            # и центрами кластеров (используйте pairwise_distances), \n",
    "            # затем найдите ближайший к каждому объекту кластер (используйте np.argmin с указанием axis)\n",
    "            \n",
    "            ### ╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ\n",
    "            \n",
    "            # (опциональный пункт)\n",
    "            # Здесь можно реализовать поддержку флага visualize и отрисовать карту, если visualize=True\n",
    "            # Отрисовывать точки нужно с помощью plt.scatter, в конце отрисовки графика используйте plt.show()\n",
    "            if visualize:\n",
    "                \n",
    "                ### ╰( ͡° ͜ʖ ͡° )つ──☆*:・ﾟ\n",
    "    \n",
    "    def transform(self, X):\n",
    "        \"\"\"\n",
    "        X - выборка размера количество объектов x количество признаков\n",
    "        \"\"\"\n",
    "        # Цель: вернуть y - номера кластеров, к которым относится каждый объект\n",
    "        # вектор длины число объектов в X\n",
    "        dists = pairwise_distances(X, self.centers)\n",
    "        return np.argmin(dists, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "_taBG5Ii62u1",
    "outputId": "0c5ad171-101e-4ecc-b741-c67e06f84a94"
   },
   "outputs": [],
   "source": [
    "cls = k_means(10)\n",
    "cls.fit(data.values, visualize=True) # visualize=True выводит график на каждой итерации\n",
    "y = cls.transform(data.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 519
    },
    "id": "HH9ghqzy67UN",
    "outputId": "0ddbd174-fe1b-4ea0-a432-f54bff3f9153"
   },
   "outputs": [],
   "source": [
    "# Итог\n",
    "\n",
    "plt.figure(figsize= (30,10))\n",
    "plt.scatter(data['LONG'], data['LAT'], c=y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cwCmcrgG7OkK"
   },
   "source": [
    "Теперь очередь DBSCAN. Какие ожидания?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "adbcwu4T7Sq8",
    "outputId": "07d2b23e-a2eb-49c2-e7b9-a95f658b893b"
   },
   "outputs": [],
   "source": [
    "# Пробуем менять min_samples\n",
    "from sklearn.cluster import DBSCAN\n",
    "for k in [1, 2, 4, 8, 16]:\n",
    "    cls = DBSCAN(min_samples=k)\n",
    "    y = cls.fit_predict(data.values)\n",
    "    plt.figure(figsize= (30,10))\n",
    "    colors = get_colors(y)\n",
    "   \n",
    "    plt.scatter(data['LONG'], data['LAT'], c=colors)\n",
    "    ratio_noise = np.round(len(y[y == -1]) / len(y), 2)\n",
    "    plt.title(\"Min samples: {} -> clusters: {}, noise ratio: {}\".format(k, len(np.unique(y[y != -1])), ratio_noise), fontsize=20)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "Obcprq4I7Zhf",
    "outputId": "4cca8b3d-0c29-4904-d99c-32b516d79357"
   },
   "outputs": [],
   "source": [
    "# Выбрали min_samples = 1\n",
    "# Пробуем менять eps\n",
    "from sklearn.cluster import DBSCAN\n",
    "for eps in 0.8**np.arange(7):\n",
    "    cls = DBSCAN(min_samples=1, eps=eps)\n",
    "    y = cls.fit_predict(data.values)\n",
    "    plt.figure(figsize= (30,10))\n",
    " \n",
    "    colors = get_colors(y)\n",
    "    plt.scatter(data['LONG'], data['LAT'], c=colors)\n",
    "    ratio_noise = np.round(len(y[y == -1]) / len(y), 2)\n",
    "    plt.title(\"Eps: {} -> clusters: {}, noise ratio: {}\".format(np.round(eps, 3), len(np.unique(y[y != -1])), ratio_noise), fontsize=20)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 538
    },
    "id": "qWCwrbxE7gBs",
    "outputId": "acd9971b-99a3-4068-e7f7-8ea96c8fba61"
   },
   "outputs": [],
   "source": [
    "# Итог\n",
    "\n",
    "cls = DBSCAN(min_samples=1, eps=0.512)\n",
    "y = cls.fit_predict(data.values)\n",
    "plt.figure(figsize= (30,10))\n",
    "colors = get_colors(y)\n",
    "plt.scatter(data['LONG'], data['LAT'], c=colors)\n",
    "plt.title(\"Eps: {} Min samples: {} -> clusters: {}\".format(np.round(cls.eps, 3), cls.min_samples, len(np.unique(y[y != -1]))), fontsize=20)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cSX7huMr7zJ8"
   },
   "source": [
    "Наконец, можно кластеризовать картинки (пиксели)!\n",
    "\n",
    "Применим наш метод кластеризации к сегментации изображений. Для этого загрузим картинку и будем каждый ее пиксель рассматривать как отдельный объект с тремя признаками R G B, отвечающими за цвет пикселя. Нам понадобится изменить форму тензора с изображением, вытянув все пиксели вдоль одной оси, чтобы получить матрицу объекты-признаки."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "YrK9n38b8Y2m",
    "outputId": "8c612375-eeef-4a7d-b5f1-3e27cd0e9cc7"
   },
   "outputs": [],
   "source": [
    "import PIL\n",
    "import urllib\n",
    "img = np.array(PIL.Image.open(urllib.request.urlopen('https://raw.githubusercontent.com/iad34/seminars/master/materials/apple.png')))\n",
    "img.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 697
    },
    "id": "6VKBjUAY8di9",
    "outputId": "915f4db4-0a39-4a5d-beb3-3232f5a51dd9"
   },
   "outputs": [],
   "source": [
    "plt.imshow(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 715
    },
    "id": "KFX0OvZ28fzT",
    "outputId": "74db19ca-e41d-49f5-8883-cbb62c6ea863"
   },
   "outputs": [],
   "source": [
    "X = img.reshape(-1, 3)\n",
    "cls = k_means(3)\n",
    "cls.fit(X)\n",
    "y = cls.transform(X)\n",
    "new_img = y.reshape(img.shape[:2])\n",
    "print(new_img.shape)\n",
    "plt.imshow(new_img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pBSocCGJ8mH7"
   },
   "source": [
    "Обратите внимание, что при кластеризации никак не учитывалась информация о положении пикселя на изображении, а учитывался только цвет."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lcInzugh9J-N"
   },
   "source": [
    "# Credits:\n",
    "\n",
    "Ковалев Евгений, Чиркова Надежда, Коган Александра"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}