{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "*c* argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with *x* & *y*.  Please use the *color* keyword-argument or provide a 2D array with a single row if you intend to specify the same RGB or RGBA value for all points.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x24045096d90>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 29144 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",
    "\n",
    "# np.random.rand = uniform random generation, [0, 1]\n",
    "# np.random.randn = normal random generation, N(0, 1)\n",
    "\n",
    "# r = 2 x 2 matrix, 찌그러트린다\n",
    "r1 = 10*(np.random.rand(2,2) - 0.5)\n",
    "r2 = 10*(np.random.rand(2,2) - 0.5)\n",
    "r3 = 10*(np.random.rand(2,2) - 0.5)\n",
    "\n",
    "# t = 1 x 2 vector, 평균을 이동\n",
    "t1 = 50*(np.random.rand(1,2) - 0.5)\n",
    "t2 = 50*(np.random.rand(1,2) - 0.5)\n",
    "t3 = 50*(np.random.rand(1,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",
    "# unlabel data\n",
    "X = np.vstack([X1, X2, X3])\n",
    "np.random.shuffle(X)\n",
    "\n",
    "# normalization [0, 1] 사이로 분포\n",
    "Xmin = np.min(X, axis = 0)\n",
    "X = (X - Xmin)\n",
    "\n",
    "Xmax = np.max(X, axis = 0)\n",
    "X = X / Xmax\n",
    "\n",
    "plt.figure(figsize = (5,5))\n",
    "plt.scatter(X[:,0], X[:,1], s = 1, c = [0.75, 0.75, 0.75])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 34248 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34500 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34800 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34768 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34812 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34788 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34788 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34788 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34788 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34788 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# k-means\n",
    "\n",
    "# 전체 클러스터의 개수\n",
    "K = 3\n",
    "maxIter = 10\n",
    "\n",
    "# 초기화\n",
    "N,D = X.shape # N개의 sample, D차원\n",
    "mu = np.random.rand(K, D)\n",
    "# k번째 centroid = mu[k,:]\n",
    "\n",
    "dist2 = np.zeros((N,K))\n",
    "#dist2[n,k] = || x[n,:] - mu[k,:] ||**2\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# iteration\n",
    "for iteration in range(maxIter):\n",
    "\n",
    "    # 주어진 sample, 주어진 centroid => 거리 => label\n",
    "    for k in range(K):\n",
    "        dist2[:,k] = np.sum( ( X - np.matlib.repmat( mu[k,:], N, 1) )**2, axis = 1)\n",
    "\n",
    "    # 정수형 label\n",
    "    c = np.argmin(dist2, axis = 1)\n",
    "\n",
    "    # one-hot coding\n",
    "    r = np.zeros((N,K), dtype = bool)\n",
    "    r[np.arange(N), c] = 1\n",
    "\n",
    "    # 주어진 sample, 주어진 label => centroid 갱신\n",
    "    for k in range(K):\n",
    "        mu[k,:] = np.mean(X[r[:,k],:], axis = 0)\n",
    "        \n",
    "    plt.figure()\n",
    "    for k in range(0,K):\n",
    "        plt.scatter(X[r[:,k],0],X[r[:,k],1], s = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 36200 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34952 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34784 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34652 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import multivariate_normal\n",
    "\n",
    "# \n",
    "K = 3\n",
    "maxIter = 20\n",
    "\n",
    "N,D = X.shape\n",
    "\n",
    "# 파라미터 초기화\n",
    "pi = np.random.rand(K)\n",
    "pi = pi / np.sum(pi)\n",
    "mu = np.random.rand(K,D) # mu[k,:] = k번째 cluster의 mean\n",
    "sigma = np.zeros((D,D,K)) # sigma[:,:,k] = k번째 cluster의 covariance\n",
    "for k in range(K):\n",
    "    sigma[:,:,k] = 1e-2*np.identity(D)\n",
    "    \n",
    "# iteration\n",
    "for iteration in range(maxIter):\n",
    "    # 주어진 sample, 주어진 param => responsibility\n",
    "    gamma = np.zeros((N,K))\n",
    "    for k in range(K):\n",
    "        gamma[:,k] = pi[k]*multivariate_normal.pdf(X, mean = mu[k,:], cov = sigma[:,:,k])\n",
    "    for n in range(N):\n",
    "        gamma[n,:] = gamma[n,:] / np.sum(gamma[n,:])\n",
    "\n",
    "    c = np.argmax(gamma, axis = 1)\n",
    "\n",
    "    # 주어진 sample, 주어진 responsibility => param 업데이트\n",
    "    Nk = np.zeros(K)\n",
    "    mu = np.zeros((K,D))\n",
    "    sigma = np.zeros((D,D,K))\n",
    "    for k in range(K):\n",
    "        Nk[k] = np.sum(gamma[:,k])\n",
    "        pi[k] = Nk[k] / N\n",
    "\n",
    "        for n in range(N):\n",
    "            mu[k,:] += gamma[n,k] * X[n,:]\n",
    "        mu[k,:] = mu[k,:] / Nk[k]\n",
    "\n",
    "        for n in range(N):\n",
    "            error = X[n,:] - mu[k,:]\n",
    "            sigma[:,:,k] += gamma[n,k]*(error.reshape(D,1)).dot(error.reshape(1,D))        \n",
    "        sigma[:,:,k] = sigma[:,:,k] / Nk[k]\n",
    "\n",
    "    plt.figure()\n",
    "    for k in range(K):\n",
    "        plt.scatter(X[c == k,0],X[c == k,1],s = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 351228 characters]\n",
      "text/plain": [
       "<Figure size 1500x500 with 10 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.datasets import make_blobs\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "# Generate sample data\n",
    "n_samples = 1000\n",
    "random_state = 0\n",
    "\n",
    "X, _ = make_blobs(n_samples=n_samples, centers=2, random_state=random_state)\n",
    "\n",
    "n_clusters_list = [2, 3, 4, 5, 6]\n",
    "\n",
    "# Algorithms to compare\n",
    "clustering_algorithms = {\n",
    "    \"K-Means\": KMeans,\n",
    "    \"MoG model\" : GaussianMixture\n",
    "}\n",
    "\n",
    "# Make subplots for each variant\n",
    "fig, axs = plt.subplots(\n",
    "    len(clustering_algorithms), len(n_clusters_list), figsize=(15, 5)\n",
    ")\n",
    "\n",
    "axs = axs.T\n",
    "\n",
    "for i, (algorithm_name, Algorithm) in enumerate(clustering_algorithms.items()):\n",
    "    for j, n_clusters in enumerate(n_clusters_list):\n",
    "        \n",
    "        if (algorithm_name == \"K-Means\"):\n",
    "            algo = Algorithm(n_clusters=n_clusters, random_state=random_state)\n",
    "        elif (algorithm_name == \"MoG model\"):\n",
    "            algo = Algorithm(n_components=n_clusters, random_state=random_state)\n",
    "            \n",
    "        algo.fit(X)\n",
    "        \n",
    "        if (algorithm_name == \"K-Means\"):\n",
    "            centers = algo.cluster_centers_    \n",
    "            labels = algo.labels_\n",
    "        elif (algorithm_name == \"MoG model\"):            \n",
    "            centers = algo.means_   \n",
    "            labels = algo.fit_predict(X)\n",
    "\n",
    "        axs[j, i].scatter(X[:, 0], X[:, 1], s=10, c=labels)\n",
    "        axs[j, i].scatter(centers[:, 0], centers[:, 1], c=\"r\", s=20)\n",
    "\n",
    "        axs[j, i].set_title(f\"{algorithm_name} : {n_clusters} clusters\")\n",
    "\n",
    "\n",
    "# Hide x labels and tick labels for top plots and y ticks for right plots.\n",
    "for ax in axs.flat:\n",
    "    ax.label_outer()\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
