{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# clear all variables and import various modules\n",
    "%reset -f\n",
    "import os\n",
    "from datetime import date\n",
    "from math import exp, pi, sin, sqrt, floor, ceil\n",
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "def eda_kmeans(S, mu):\n",
    "    # k-mean cluster analsis of NxM sample natrix S\n",
    "    # starting from initial guess KxM mu of cluster means\n",
    "    # returns:\n",
    "    # newmu: new cluster means\n",
    "    # k: Nx1 vector of cluster assignments\n",
    "    # dist: Nx1 vector of distance to cluster center\n",
    "    # count: Kx1 vector of samples in cluster\n",
    "    # Kdist: KxK matrix of inter-cluster distances\n",
    "    N, M = np.shape(S);   # N samples with M elements\n",
    "    K, M2 = np.shape(mu); # K cluster meabs\n",
    "    k = np.zeros((N,1),dtype=np.int); # assoc sample w cluster\n",
    "    dist = np.zeros((N,1)); # distance from cluster mean\n",
    "    Kdist = np.zeros((K,K)); # cluste to nearest cluster distance\n",
    "    oldmu = np.copy(mu);  # old cluster means\n",
    "    \n",
    "    MAXP = 2*K*M;  # max iterations\n",
    "    for p in range(MAXP):  # iterate\n",
    "    \n",
    "        # cluster assignment step\n",
    "        changes = 0; # number of changes\n",
    "        for i in range(N):     # loop over samples\n",
    "            for j in range(K): # loop over cluster means\n",
    "                # deviation of sample from mean\n",
    "                delta = (S[i:i+1,0:M] - oldmu[j:j+1,0:M]).T;\n",
    "                # squared distance\n",
    "                r2 = np.matmul(delta.T,delta);\n",
    "                if( j==0 ): # always accept first distance\n",
    "                    r2min = r2;\n",
    "                    newk = j;\n",
    "                elif r2<r2min: # accept if smaller\n",
    "                    r2min = r2;\n",
    "                    newk = j;\n",
    "            if k[i,0] != newk: # count reassignements\n",
    "                changes=changes+1;\n",
    "            k[i,0] = newk;     # reassigned cluster\n",
    "            dist[i,0] = r2min; # sq distance to that cluster\n",
    "\n",
    "        # end loop over iterations on no further changes\n",
    "        if( (p!=0) and (changes==0) ):\n",
    "            break;\n",
    "\n",
    "        # means step\n",
    "        newmu = np.zeros((K,M));  # new cluster means\n",
    "        count = np.zeros((K,1),dtype=np.int); # samples in cluster\n",
    "        for i in range(N): # loop over samples\n",
    "            for j in range(M): # loop over elements\n",
    "                # cluster sum\n",
    "                newmu[k[i,0],j] = newmu[k[i,0],j] + S[i,j];\n",
    "            # count samples in cluster\n",
    "            count[k[i,0],0] = count[k[i,0],0] + 1;\n",
    "        for j in range(K): # loop over clusters\n",
    "            if( count[j,0] == 0 ): # cluser with no samples\n",
    "                # no change in cluster mean\n",
    "                newmu[j,0:M] = oldmu[j,0:M];\n",
    "            else:\n",
    "                # new cluster mean\n",
    "                newmu[j,0:M] = newmu[j,0:M]/count[j,0]\n",
    "        oldmu=newmu; # update cluster mean\n",
    "        \n",
    "    # inter-cluster distances\n",
    "    for i in range(K-1):     # loop over cluster means\n",
    "        for j in range(1,K): # loop over cluster means\n",
    "            delta = (newmu[i:i+1,0:M] - newmu[j:j+1,0:M]).T;\n",
    "            r2 = np.matmul(delta.T,delta);\n",
    "            Kdist[i,j] = r2; # sq distance between clusters\n",
    "            Kdist[j,i] = r2; # sq distance between clusters\n",
    "    return newmu, k, dist, count, Kdist;\n",
    "\n",
    "def eda_kmeans_mod(S, mu):\n",
    "    # k-mean cluster analsis of NxM sample natrix S\n",
    "    # modified per Bradley & Fayyad (1998) never to\n",
    "    # return cluster with zero samples.\n",
    "    # Starts with initial guess KxM mu of cluster means\n",
    "    # returns:\n",
    "    # newmu: new cluster means\n",
    "    # k: Nx1 vector of cluster assignments\n",
    "    # dist: Nx1 vector of distance to cluster center\n",
    "    # count: Kx1 vector of samples in cluster\n",
    "    initmu = np.copy(mu);\n",
    "    mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "    while np.min(mycount)==0:  # while cluster with no samples\n",
    "        ic=np.argmin(mycount); # index of cluster\n",
    "        ir=np.argmax(mydist);  # index of most distant sample\n",
    "        ik=myk[ir,0];          # cluster of this sample\n",
    "        mymu[ic:ic+1,0:M] = S[ir:ir+1,0:M]; # new mean\n",
    "        mycount[ic,0]=1;       # reset cluster count\n",
    "        myk[ir,0]=ic;          # reset sample membership\n",
    "        mydist[ir,0]=0.0;      # reset sample distance\n",
    "        # note: in rare cases, next line could create\n",
    "        # infinite loop, so don't do it.\n",
    "        # mycount[ik,0]=mycount[ik,0]-1;\n",
    "    # redo k-means\n",
    "    mymu2, myk2, mydist2, mycount2, myKdist2 = eda_kmeans(S, mymu);  \n",
    "    return mymu2, myk2, mydist2, mycount2, myKdist2;\n",
    "\n",
    "def eda_forgy_mean(S,K):\n",
    "    # create initial k-mean cluster means, from samples\n",
    "    # randomly drawn NxM sample matrix S\n",
    "    # returns KxM matrix mu of means\n",
    "    N,M=np.shape(S);\n",
    "    k = np.random.randint(low=0,high=N,size=K,dtype=np.int);\n",
    "    mu = S[k,0:M];\n",
    "    return mu;\n",
    "\n",
    "def eda_random_partition_mean(S,K):\n",
    "    # create initial k-mean cluster means, by K random\n",
    "    # partitions of NxM sample matrix S\n",
    "    # returns KxM matrix mu of means\n",
    "    N,M=np.shape(S);\n",
    "    k = np.random.randint(low=0,high=K,size=N,dtype=np.int);\n",
    "    mu = np.zeros((K,M));\n",
    "    for i in range(K):\n",
    "        j = np.where(k==i);\n",
    "        mu[i:i+1,0:M] = np.mean(S[j,0:M],axis=1);\n",
    "    return mu;\n",
    "\n",
    "def eda_bradley_fayyad_mean(S,K,Nsets,fraction):\n",
    "    # create initial k-mean cluster means, by\n",
    "    # the Bradley & Fayyad (1998) method, which\n",
    "    # is based on clustering clusters. The NxM\n",
    "    # sample matrix S is randomly subsampled Nsets\n",
    "    # times, with each set containing (fraction*N+2*K)\n",
    "    # samples.\n",
    "    # returns KxM matrix mu of means\n",
    "    \n",
    "    initmu = eda_random_partition_mean(S,K);\n",
    "\n",
    "    KiMAX = K*Nsets;  # maximum number of means mui\n",
    "    Ni = floor(fraction*N)+2*K; # size of snall sets\n",
    "\n",
    "    # Step 1, create KixM table CM of means of subsampled data\n",
    "    CM = np.zeros((KiMAX,M));\n",
    "    Ki=0; # counter for rows of CM\n",
    "    for j in range(Nsets):\n",
    "        k = np.random.randint(low=0,high=N,size=Ni,dtype=np.int);\n",
    "        Si = S[k,0:M];\n",
    "        mymu, myk, mydist, mycount, myKdist = eda_kmeans_mod(Si, initmu);\n",
    "        CM[Ki:Ki+K,0:M] = mymu[0:K,0:M];\n",
    "        Ki=Ki+K;\n",
    "    CM = np.copy(CM[0:Ki,0:M]);\n",
    "\n",
    "    # Step 2, create KixM table FMS of means of cluster means\n",
    "    FMS = np.zeros((KiMAX,M));\n",
    "    Ki=0; # counter for rows of CM\n",
    "    for j in range(Nsets):\n",
    "        mymu, myk, mydist, mycount, myKdist = eda_kmeans(CM, CM[3*j:3*j+K,0:M]);\n",
    "        FMS[Ki:Ki+K,0:M] = mymu[0:K,0:M];\n",
    "        Ki=Ki+K;\n",
    "    FMS = np.copy(FMS[0:Ki,0:M]);\n",
    "\n",
    "    # Step 3A: associate each KxM member FMi of Fm\n",
    "    #          with all individual 1xM elements of CM\n",
    "    #          and compute the overall error\n",
    "    error = np.zeros((Nsets,1));\n",
    "    Ki=0;\n",
    "    for i in range(Nsets):\n",
    "        FMi = FMS[Ki:Ki+K,0:M] # one group of means\n",
    "        Ki=Ki+K;\n",
    "        for j in range(Ki):\n",
    "            CMj = CM[j:j+1,0:M]; # one member of CM\n",
    "            for m in range(K):\n",
    "                delta = (CMj[0:1,0:M] - FMi[m:m+1,0:M]).T;\n",
    "                r2 = np.matmul(delta.T,delta);\n",
    "                if( m==0 ): # always acceot first distance\n",
    "                    r2min = r2;\n",
    "                elif r2<r2min: # accept if smaller\n",
    "                    r2min = r2;\n",
    "        error[i,0] = error[i,0] + r2min;\n",
    "\n",
    "    # Step 3B: select the FMi with lowest error\n",
    "    i = np.argmin(error);\n",
    "    mu = np.copy( FMS[K*i:K*i+K,0:M] );\n",
    "    return(mu);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13508 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# set up for cluster analysys\n",
    "\n",
    "K=3; # three clusters\n",
    "# 100 samples per cluster\n",
    "N1=100;\n",
    "N2=100;\n",
    "N3=100;\n",
    "N=N1+N2+N3;\n",
    "\n",
    "M=2; # two elements\n",
    "\n",
    "# means set by hand\n",
    "mu = np.zeros((K,M));\n",
    "mu[0,0] = 1.0;\n",
    "mu[0,1] = 2.0;\n",
    "mu[1,0] = 3.0;\n",
    "mu[1,1] = 4.0;\n",
    "mu[2,0] = 4.0;\n",
    "mu[2,1] = 1.0;\n",
    "\n",
    "# sample matrix randomly assigned\n",
    "S = np.zeros((N,M));\n",
    "sc = 0.5;\n",
    "S[0:N1,0:1] = np.random.normal( loc=mu[0,0], scale=sc, size=(N1,1) );\n",
    "S[0:N1,1:2] = np.random.normal( loc=mu[0,1], scale=sc, size=(N1,1) );\n",
    "S[N1:N1+N2,0:1] = np.random.normal( loc=mu[1,0], scale=sc, size=(N2,1) );\n",
    "S[N1:N1+N2,1:2] = np.random.normal( loc=mu[1,1], scale=sc, size=(N2,1) );\n",
    "S[N1+N2:N,0:1] = np.random.normal( loc=mu[2,0], scale=sc, size=(N3,1) );\n",
    "S[N1+N2:N,1:2] = np.random.normal( loc=mu[2,1], scale=sc, size=(N3,1) );\n",
    "\n",
    "# cluster membership\n",
    "k1 = np.zeros((N1,1),dtype=np.int);\n",
    "k2 = np.ones((N2,1),dtype=np.int);\n",
    "k3 = 2*np.ones((N3,1),dtype=np.int);\n",
    "myk = np.concatenate((k1,k2,k3),axis=0);\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ko', lw=1 );\n",
    "plt.plot( mu[0:K,0:1], mu[0:K,1:2], 'b^', lw=2 );\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();\n",
    "\n",
    "Dm = np.concatenate( (S,myk),axis=1);\n",
    "np.savetxt(\"../Data/threeclusters.txt\",Dm, delimiter=\"\\t\");\n",
    "np.savetxt(\"../Data/threeclustermeans.txt\",mu, delimiter=\"\\t\");\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster means\n",
      "[[1.03045877 2.09680122]\n",
      " [3.03078748 4.06452326]\n",
      " [3.97327896 1.02852854]]\n",
      "mean sample distance\n",
      "0.6229797879025173\n",
      "cluster counts\n",
      "[[100]\n",
      " [100]\n",
      " [100]]\n",
      "mean inter-cluster distance\n",
      "2.025682677938847\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13776 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# cluster analysis using prior mean initialization\n",
    "initmu = np.zeros((K,M));\n",
    "initmu[0,0] = 1.2;\n",
    "initmu[0,1] = 2.2;\n",
    "initmu[1,0] = 3.4;\n",
    "initmu[1,1] = 4.4;\n",
    "initmu[2,0] = 4.1;\n",
    "initmu[2,1] = 1.1;\n",
    "mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "\n",
    "print('cluster means');\n",
    "print(mymu);\n",
    "print('mean sample distance');\n",
    "print(np.mean(np.sqrt(mydist)));\n",
    "print('cluster counts');\n",
    "print(mycount);\n",
    "print('mean inter-cluster distance');\n",
    "print(np.mean(np.sqrt(myKdist)));\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "plt.plot( initmu[0:K,0:1], initmu[0:K,1:2], 'r^', lw=2 )\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ko', lw=1 );\n",
    "plt.plot( mymu[0:K,0:1], mymu[0:K,1:2], 'b^', lw=2 )\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster means\n",
      "[[2.99849561 4.06511189]\n",
      " [4.00682973 0.95783841]\n",
      " [1.03566731 2.00330593]]\n",
      "mean sample distance\n",
      "0.6462561767205608\n",
      "cluster counts\n",
      "[[ 98]\n",
      " [ 99]\n",
      " [103]]\n",
      "mean inter-cluster distance\n",
      "2.058494097555173\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13752 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# cluster analysis using random-partition initialization\n",
    "initmu = eda_random_partition_mean(S,K);\n",
    "mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "\n",
    "print('cluster means');\n",
    "print(mymu);\n",
    "print('mean sample distance');\n",
    "print(np.mean(np.sqrt(mydist)));\n",
    "print('cluster counts');\n",
    "print(mycount);\n",
    "print('mean inter-cluster distance');\n",
    "print(np.mean(np.sqrt(myKdist)));\n",
    "\n",
    "fig3 = plt.figure(3,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "plt.plot( initmu[0:K,0:1], initmu[0:K,1:2], 'r^', lw=2 )\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ko', lw=1 );\n",
    "plt.plot( mymu[0:K,0:1], mymu[0:K,1:2], 'b^', lw=2 );\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster means\n",
      "[[4.00682973 0.95783841]\n",
      " [2.99849561 4.06511189]\n",
      " [1.03566731 2.00330593]]\n",
      "mean sample distance\n",
      "0.6462561767205608\n",
      "cluster counts\n",
      "[[ 98]\n",
      " [ 99]\n",
      " [103]]\n",
      "mean inter-cluster distance\n",
      "2.058494097555173\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13684 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# cluster analysis using forgy initialization\n",
    "initmu = eda_forgy_mean(S,K);\n",
    "mymu, myk, mydist, mycountm, myKdist = eda_kmeans(S, initmu);\n",
    "\n",
    "print('cluster means');\n",
    "print(mymu);\n",
    "print('mean sample distance');\n",
    "print(np.mean(np.sqrt(mydist)));\n",
    "print('cluster counts');\n",
    "print(mycount);\n",
    "print('mean inter-cluster distance');\n",
    "print(np.mean(np.sqrt(myKdist)));\n",
    "\n",
    "fig4 = plt.figure(4,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "plt.plot( initmu[0:K,0:1], initmu[0:K,1:2], 'r^', lw=2 )\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ko', lw=1 );\n",
    "plt.plot( mymu[0:K,0:1], mymu[0:K,1:2], 'b^', lw=2 )\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster means\n",
      "[[1.03566731 2.00330593]\n",
      " [4.00682973 0.95783841]\n",
      " [2.99849561 4.06511189]]\n",
      "mean sample distance\n",
      "0.6462561767205608\n",
      "cluster counts\n",
      "[[103]\n",
      " [ 99]\n",
      " [ 98]]\n",
      "mean inter-cluster distance\n",
      "2.058494097555173\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 13516 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# cluster analysis using bradley-fayyad initialization\n",
    "Nsets = 100;     # number of small sets\n",
    "fraction = 0.05; # size of small sets = fraction of N\n",
    "initmu = eda_bradley_fayyad_mean(S,K,Nsets,fraction);\n",
    "mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "\n",
    "print('cluster means');\n",
    "print(mymu);\n",
    "print('mean sample distance');\n",
    "print(np.mean(np.sqrt(mydist)));\n",
    "print('cluster counts');\n",
    "print(mycount);\n",
    "print('mean inter-cluster distance');\n",
    "print(np.mean(np.sqrt(myKdist)));\n",
    "\n",
    "fig5 = plt.figure(5,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "plt.plot( initmu[0:K,0:1], initmu[0:K,1:2], 'r^', lw=2 )\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ko', lw=1 );\n",
    "plt.plot( mymu[0:K,0:1], mymu[0:K,1:2], 'b^', lw=2 )\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster means\n",
      "[[3.07599548 3.91479333]\n",
      " [3.8587548  4.00093578]\n",
      " [4.00682973 0.95783841]\n",
      " [2.83015838 4.62503362]\n",
      " [0.99341924 1.94257799]\n",
      " [2.20755041 3.443305  ]]\n",
      "mean sample distance\n",
      "0.553983781014169\n",
      "cluster counts\n",
      "[[38]\n",
      " [15]\n",
      " [99]\n",
      " [27]\n",
      " [97]\n",
      " [24]]\n",
      "mean inter-cluster distance\n",
      "1.92351186935618\n",
      "inter-cluster distances\n",
      "[[0.         0.78748502 3.10000562 0.75158309 2.86823241 0.98817919]\n",
      " [0.78748502 0.         3.04669785 1.20312456 3.52802843 1.74282186]\n",
      " [3.10000562 3.04669785 0.         3.85134732 3.17022946 3.06837909]\n",
      " [0.75158309 1.20312456 3.85134732 0.         3.25102736 1.33571076]\n",
      " [2.86823241 3.52802843 3.17022946 3.25102736 0.         1.93036163]\n",
      " [0.98817919 1.74282186 3.06837909 1.33571076 1.93036163 0.        ]]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 15656 characters]\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# cluster analysis using forgy initialization\n",
    "# six clusers\n",
    "K6=6;\n",
    "initmu = eda_forgy_mean(S,K6);\n",
    "mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "\n",
    "print('cluster means');\n",
    "print(mymu);\n",
    "print('mean sample distance');\n",
    "print(np.mean(np.sqrt(mydist)));\n",
    "print('cluster counts');\n",
    "print(mycount);\n",
    "print('mean inter-cluster distance');\n",
    "print(np.mean(np.sqrt(myKdist)));\n",
    "print('inter-cluster distances');\n",
    "print(np.sqrt(myKdist));\n",
    "\n",
    "fig4 = plt.figure(4,figsize=(5,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, 5, 0, 5] );\n",
    "for i in range(N):\n",
    "    if myk[i,0]==0:\n",
    "        plt.plot( S[i,0], S[i,1], 'k.', lw=1 );\n",
    "    elif myk[i,0]==1:\n",
    "        plt.plot( S[i,0], S[i,1], 'k+', lw=1 );\n",
    "    elif myk[i,0]==2:\n",
    "        plt.plot( S[i,0], S[i,1], 'k^', lw=1 );\n",
    "    elif myk[i,0]==3:\n",
    "        plt.plot( S[i,0], S[i,1], 'kx', lw=1 );\n",
    "    elif myk[i,0]==4:\n",
    "        plt.plot( S[i,0], S[i,1], 'kv', lw=1 );\n",
    "    else:\n",
    "        plt.plot( S[i,0], S[i,1], 'ks', lw=1 );\n",
    "plt.plot( mymu[0:K6,0:1], mymu[0:K6,1:2], 'r^', lw=2 )\n",
    "plt.xlabel('x');\n",
    "plt.ylabel('y');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K mean dist\n",
      "1 2.450872469766111\n",
      "2 2.159369027284355\n",
      "3 1.87730362348076\n",
      "4 1.6844183908052115\n",
      "5 1.5978670360338794\n",
      "6 1.5205777930497673\n",
      "7 1.4512789214760604\n",
      "8 1.406422437353292\n",
      "9 1.3793476892712468\n",
      "10 1.3477254467136452\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 12952 characters]\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "S = np.genfromtxt('rocks.txt', delimiter='\\t')\n",
    "N, M = np.shape(S);\n",
    "\n",
    "KMAX = 11;\n",
    "Klist = np.zeros((KMAX,1));\n",
    "Dlist = np.zeros((KMAX,1));\n",
    "print('K','mean dist');\n",
    "for K in range(1,KMAX):\n",
    "    initmu = eda_forgy_mean(S,K);\n",
    "    mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "    Klist[K,0]=K;\n",
    "    Dlist[K,0]=np.mean(np.sqrt(mydist));\n",
    "    print(K,Dlist[K,0]);\n",
    "    \n",
    "fig1 = plt.figure(1,figsize=(10,5));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, KMAX, 0, 3] );\n",
    "plt.plot( Klist[1:KMAX,0:1], Dlist[1:KMAX,0:1], 'k-', lw=2 )\n",
    "plt.plot( Klist[1:KMAX,0:1], Dlist[1:KMAX,0:1], 'ko', lw=2 )\n",
    "plt.xlabel('K');\n",
    "plt.ylabel('mean Rij');\n",
    "plt.show();\n",
    "\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cluster\tSiO2\tTiO2\tAl203\tFeO-t\tMgO\tCaO\tNa2O\tK2O\tcount\n",
      "0\t50.66\t1.55\t14.11\t12.14\t6.89\t11.54\t2.26\t0.15\t901\n",
      "1\t51.47\t1.83\t15.01\t10.36\t6.56\t10.60\t3.00\t0.32\t729\n",
      "2\t50.98\t1.19\t15.00\t9.71\t7.80\t12.24\t2.24\t0.16\t1570\n",
      "3\t44.90\t0.06\t2.51\t8.20\t40.22\t1.52\t0.13\t0.03\t54\n",
      "4\t46.53\t0.61\t13.41\t8.72\t15.06\t10.87\t1.65\t0.16\t35\n",
      "5\t49.28\t1.04\t16.27\t8.83\t8.47\t12.33\t2.28\t0.19\t931\n",
      "6\t50.46\t1.50\t15.63\t9.54\t7.75\t11.29\t2.78\t0.19\t1879\n",
      "7\t47.85\t1.00\t17.98\t7.90\t6.54\t13.63\t2.18\t0.19\t257\n"
     ]
    }
   ],
   "source": [
    "S = np.genfromtxt('rocks.txt', delimiter='\\t')\n",
    "N, M = np.shape(S);\n",
    "\n",
    "K=8;\n",
    "initmu = eda_forgy_mean(S,K);\n",
    "mymu, myk, mydist, mycount, myKdist = eda_kmeans(S, initmu);\n",
    "names = [\"SiO2\", \"TiO2\", \"Al203\", \"FeO-t\", \"MgO\", \"CaO\", \"Na2O\", \"K2O\" ];\n",
    "print('cluster\\t', end='');\n",
    "for j in range(M):\n",
    "    print(\"%s\\t\" % (names[j]), end='');\n",
    "print('count');\n",
    "for i in range(K):\n",
    "    print(\"%d\\t\" % (i), end='');\n",
    "    for j in range(M):\n",
    "        print(\"%.2f\\t\" % (mymu[i,j]), end='');\n",
    "    print(\"%d\" % (mycount[i,0]));\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 45208 characters]\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#\n",
    "m1 = 2; # Al203\n",
    "xmin = 0.0\n",
    "xmax = 50.0;\n",
    "\n",
    "m2 = 4; # MgO\n",
    "ymin = 0.0\n",
    "ymax = 50.0;\n",
    "\n",
    "m3 = 5; # Ca0\n",
    "zmin = 0.0\n",
    "zmax = 20.0;\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(10,10));\n",
    "ax1 = fig2.add_subplot(121, projection='3d')\n",
    "ax1.axes.set_xlim3d(left=xmin, right=xmax);\n",
    "ax1.axes.set_ylim3d(bottom=ymin, top=ymax);\n",
    "ax1.axes.set_zlim3d(bottom=zmin, top=zmax);\n",
    "ax1.view_init(1, 30);\n",
    "\n",
    "ax1.scatter( xs=S[:,m1], ys=S[:,m2], zs=S[:,m3], color=(0,0,0) );\n",
    "#ax1.scatter( xs=mymu[:,m1], ys=mymu[:,m2], zs=mymu[:,m3], color=(0.8,0,0) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymin,ymin], zs=[zmin,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymin,ymax], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmax], ys=[ymin,ymin], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymax], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymax,ymax], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymin,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymin,ymin], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymax,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymax,ymax], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymin], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymax,ymax], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.set_xlabel(names[m1]);\n",
    "ax1.set_ylabel(names[m2]);\n",
    "ax1.set_zlabel(names[m3]);\n",
    "\n",
    "ax1 = fig2.add_subplot(122, projection='3d')\n",
    "ax1.axes.set_xlim3d(left=xmin, right=xmax);\n",
    "ax1.axes.set_ylim3d(bottom=ymin, top=ymax);\n",
    "ax1.axes.set_zlim3d(bottom=zmin, top=zmax);\n",
    "ax1.view_init(1, 30);\n",
    "\n",
    "#ax1.scatter( xs=S[:,m1], ys=S[:,m2], zs=S[:,m3], color=(0,0,0,0.1) );\n",
    "ax1.scatter( xs=mymu[:,m1], ys=mymu[:,m2], zs=mymu[:,m3], color=(0,0,0) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymin,ymin], zs=[zmin,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymin,ymax], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmax], ys=[ymin,ymin], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymax], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymax,ymax], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymin,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymin,ymin], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymax,ymin], zs=[zmax,zmax], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmin,xmin], ys=[ymax,ymax], zs=[zmax,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.plot( xs=[xmax,xmax], ys=[ymax,ymin], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "ax1.plot( xs=[xmax,xmin], ys=[ymax,ymax], zs=[zmin,zmin], color=(0.8,0.8,0.8) );\n",
    "\n",
    "ax1.set_xlabel(names[m1]);\n",
    "ax1.set_ylabel(names[m2]);\n",
    "ax1.set_zlabel(names[m3]);\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
