{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "11197577",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 33600 characters]\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np \n",
    "import numpy.matlib \n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colors\n",
    "\n",
    "# np.random.rand = uniform random generation, [0, 1]\n",
    "# np.random.randn = normal random generation, N(0, 1)\n",
    "r1 = 10*(np.random.rand(2,2) - 0.5) # [-1, 1]\n",
    "r2 = 10*(np.random.rand(2,2) - 0.5) # [-1, 1]\n",
    "r3 = 10*(np.random.rand(2,2) - 0.5) # [-1, 1]\n",
    "r4 = 10*(np.random.rand(2,2) - 0.5) # [-1, 1]\n",
    "\n",
    "t1 = 50*(np.random.rand(2) - 0.5)\n",
    "t2 = 50*(np.random.rand(2) - 0.5)\n",
    "t3 = 50*(np.random.rand(2) - 0.5)\n",
    "t4 = 50*(np.random.rand(2) - 0.5)\n",
    "\n",
    "X1 = np.random.randn(1000,2)\n",
    "X1 = X1.dot(r1) + t1\n",
    "\n",
    "X2 = np.random.randn(500,2)\n",
    "X2 = X2.dot(r2) + t2\n",
    "\n",
    "X3 = np.random.randn(500,2)\n",
    "X3 = X3.dot(r3) + t3\n",
    "\n",
    "X4 = np.random.randn(500,2)\n",
    "X4 = X4.dot(r4) + t4\n",
    "\n",
    "# unlabel data\n",
    "X = np.vstack((X1, X2, X3, X4))\n",
    "\n",
    "# normalization, [0, 1]\n",
    "Xmin = np.min(X, axis = 0)\n",
    "Xmax = np.max(X, axis = 0)\n",
    "X = (X - Xmin) / (Xmax - Xmin)\n",
    "\n",
    "plt.figure(figsize = (5,5))\n",
    "plt.plot(X[:,0],X[:,1],\".\",markersize = 1, color = [0.75, 0.75, 0.75])\n",
    "\n",
    "N,D = X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "29697d3c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jacob\\AppData\\Local\\Temp\\ipykernel_73348\\2930623055.py:24: RuntimeWarning: Mean of empty slice\n",
      "  mu[:,k] = np.nanmean( X[label == k,:], axis = 0)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x20ea7413130>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 22752 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def computeDistance2(samples, point):\n",
    "    N,D = samples.shape\n",
    "    D, = point.shape\n",
    "    dist2 = np.sum( (samples - np.matlib.repmat(point.reshape(1,-1),N,1))**2, axis = 1)\n",
    "    \n",
    "    return dist2.reshape(-1)\n",
    "\n",
    "inertia = np.zeros(20)\n",
    "inertia[:] = np.nan\n",
    "\n",
    "for K in range(2, 8):\n",
    "    dist2 = np.zeros((N,K))\n",
    "    dist2[:] = np.nan\n",
    "    \n",
    "    while (np.sum(np.isnan(dist2)) > 0):\n",
    "        mu = np.random.rand(D,K)\n",
    "        iterMax = 10\n",
    "        for iter in range(iterMax):\n",
    "            for k in range(K):\n",
    "                dist2[:,k] = computeDistance2(X, mu[:,k])\n",
    "            label = np.argmin(dist2, axis = 1)\n",
    "\n",
    "            for k in range(K):\n",
    "                mu[:,k] = np.nanmean( X[label == k,:], axis = 0)            \n",
    "\n",
    "    inertia[K] = np.mean( np.min(dist2, axis = 1) )\n",
    "    \n",
    "plt.figure()\n",
    "plt.plot(inertia)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7b1ff52a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jacob\\AppData\\Local\\Temp\\ipykernel_73348\\1345307982.py:24: RuntimeWarning: Mean of empty slice\n",
      "  mu[:,k] = np.nanmean( X[label == k,:], axis = 0)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x20ea88f8fa0>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 28256 characters]\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 32668 characters]\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 36020 characters]\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 43020 characters]\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 41908 characters]\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 49976 characters]\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 25456 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def computeDistance2(samples, point):\n",
    "    N,D = samples.shape\n",
    "    D, = point.shape\n",
    "    dist2 = np.sum( (samples - np.matlib.repmat(point.reshape(1,-1),N,1))**2, axis = 1)\n",
    "    \n",
    "    return dist2.reshape(-1)\n",
    "\n",
    "s_score = np.zeros(20)\n",
    "s_score[:] = np.nan\n",
    "\n",
    "for K in range(2, 8):\n",
    "    dist2 = np.zeros((N,K))\n",
    "    dist2[:] = np.nan\n",
    "    \n",
    "    while (np.sum(np.isnan(dist2)) > 0):\n",
    "        mu = np.random.rand(D,K)\n",
    "        iterMax = 10\n",
    "        for iter in range(iterMax):\n",
    "            for k in range(K):\n",
    "                dist2[:,k] = computeDistance2(X, mu[:,k])\n",
    "            label = np.argmin(dist2, axis = 1)\n",
    "\n",
    "            for k in range(K):\n",
    "                mu[:,k] = np.nanmean( X[label == k,:], axis = 0)            \n",
    "\n",
    "    s_coeff = np.zeros(N)\n",
    "    a = np.zeros(N) # intra-cluster distance\n",
    "    b = np.zeros(N) # nearest-cluster distance\n",
    "    \n",
    "    for n in range(N):\n",
    "        # 우리가 속한 cluster의 index (intra)\n",
    "        # 그 index에 해당되는 sample들을 추림\n",
    "        # 그 샘플들이랑, 나 = x[n,:]랑 거리를 구해요\n",
    "        # 그 거리들의 평균을 냅니다 (a[n])\n",
    "        label_intra = label[n]\n",
    "        mask_intra = (label == label_intra)\n",
    "        a[n] = np.mean(np.sqrt(computeDistance2(X[mask_intra,:], X[n,:])))\n",
    "        \n",
    "        # 나랑 제일 가까운 cluster의 index (nearest)\n",
    "        # 그 index에 해당되는 sample들을 추림\n",
    "        # 그 샘플들이랑, 나 = x[n,:]랑 거리를 구해요\n",
    "        # 그 거리들의 평균을 냅니다 (b[n])\n",
    "        dist2_sort = np.sort(dist2[n,:])\n",
    "        label_nearest = np.where(dist2[n,:] == dist2_sort[1])[0]\n",
    "        mask_nearest = (label == label_nearest)\n",
    "        b[n] = np.mean(np.sqrt(computeDistance2(X[mask_nearest,:], X[n,:])))\n",
    "        \n",
    "        # s_coeff[n] = (b[n] - a[n]) / max([a[n], b[n]])\n",
    "        s_coeff[n] = (b[n] - a[n]) / np.max([a[n], b[n]])\n",
    "        \n",
    "    s_score[K] = np.nanmean(s_coeff)\n",
    "    \n",
    "    plt.figure(figsize = (5,5))\n",
    "    for k in range(0,K):\n",
    "        plt.plot(np.sort(s_coeff[label == k]))\n",
    "    plt.xlabel('samples')\n",
    "    plt.ylabel('silhouette coefficient')\n",
    "    \n",
    "        \n",
    "plt.figure()\n",
    "plt.plot(s_score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4c1f9dd1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x20eaba959d0>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 42828 characters]\n",
      "text/plain": [
       "<Figure size 800x800 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "from sklearn.metrics import silhouette_score \n",
    "\n",
    "max_cluster = 10\n",
    "\n",
    "score_inertia = np.full(max_cluster + 1, np.nan)\n",
    "score_silhouette = np.full(max_cluster + 1, np.nan)\n",
    "\n",
    "for n_cluster in range(2, max_cluster + 1):\n",
    "    kmeans = KMeans(n_clusters = n_cluster, random_state = 0).fit(X)\n",
    "    score_inertia[n_cluster] = kmeans.inertia_\n",
    "    score_silhouette[n_cluster] = silhouette_score(X, kmeans.labels_)\n",
    "        \n",
    "plt.figure(figsize = (8,8))\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(score_inertia,'o-')\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(score_silhouette,'o-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d039261f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
