import numpy as np
import sys
from math import *
import matplotlib.pyplot as pl

def INITC(kk,xd):
    nn=np.zeros([kk],dtype=np.int)
    cli=np.zeros([ndata],dtype=np.int)
    d1=np.zeros([ndata],dtype=np.float)
    d2=np.zeros([ndata],dtype=np.float)
    dr=np.zeros([ndata],dtype=np.float)
    di=np.zeros([kk,ndata],dtype=np.float)
    for i in range(0,ndata):
        d1[i]=1.0e10
    nn[0]=np.random.randint(0,ndata,1) # ndata nos random number from 0 to kk-1
    for k in range(1,kk):
        ck=xd[nn[k-1],:]
        for i in range(0,ndata):
            d2[i]=sqrt(np.sum((xd[i,:]-ck[:])**2))
        for i in range(0,ndata):
            dr[i]=np.amin([d1[i],d2[i]])
        nn[k]=np.argmax(dr)
        d1=dr
    for k in range(0,kk):
        for i in range(0,ndata):
            di[k,i]=sqrt(np.sum((xd[i,:]-xd[nn[k],:])**2))
    for i in range(0,ndata):
        cli[i]=np.argmin(di[:,i])
    return cli

def MAIN(mcol,ndata,kk,itmax,xd):
    cl=np.zeros([ndata],dtype=np.int)
    clr=np.zeros([ndata],dtype=np.int)
    ck=np.zeros([mcol],dtype=np.float)
    xk=np.zeros([ndata,mcol],dtype=np.float)
    dis=np.zeros([kk,ndata],dtype=np.float)
    kkk=0
    cl=INITC(kk,xd)
#   cl=np.random.randint(0,kk,ndata) # ndata nos random number from 0 to kk-1
    for iii in range(0,itmax):
        kkk=kkk+1
        iflag=0
        for k in range(0,kk):
            n=-1
            for i in range(0,ndata):
                if cl[i]==k:
                    n=n+1
                    xk[n,:]=xd[i,:]
            if 0<n:
                wxi=np.linalg.inv(np.cov(xk[0:n+1,0:mcol].T))
                for j in range(0,mcol):
                    ck[j]=np.sum(xk[0:n+1,j])/float(n+1)
                for i in range(0,ndata):
                    w0=xd[i,:]-ck[:]
                    w1=np.dot(wxi,w0)
                    w2=np.dot(w0,w1)
                    dis[k,i]=sqrt(w2) # Mahalanobis distance
#                    dis[k,i]=sqrt(np.sum((xd[i,:]-ck[:])**2)) # Euclidean distance
            else:
                iflag=1
                break
        if iflag==0:
            for i in range(0,ndata):
                cl[i]=np.argmin(dis[:,i])
            if 0<kkk and np.sum(cl-clr)==0:
                break
        else:
            kkk=0
            cl=INITC(kk,xd)
#            cl=np.random.randint(0,kk,ndata) # ndata nos random number from 0 to kk-1
    return cl

#Main routine
param=sys.argv
kk    =int(param[1])
mds   =int(param[2])
itmax =int(param[3])
fnameR=param[4]
fnameW=param[5]

# Data input
fin=open(fnameR,'r')
text=fin.readline()
text=text.strip()
strcom=text

text=fin.readline()
text=text.strip()
text=text.split(',')
mcol=int(text[0])
ndata=int(text[1])

xd=np.zeros([ndata,mcol],dtype=np.float)
for i in range(0,ndata):
    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    for j in range(0,mcol):
        xd[i,j]=float(text[j])
fin.close()

num=np.zeros([kk],dtype=np.int)
for jjj in range(0,int(itmax/10)):
    cl=MAIN(mcol,ndata,kk,itmax,xd)
    for k in range(0,kk):
        num[k]=0
        for i in range(0,ndata):
            if cl[i]==k:
                num[k]=num[k]+1
    print(num)
    if mds<=min(num):
        break

# Data output
iw=mcol
if 10<iw:
    iw=10
fout=open(fnameW,'w')
fout.write('{0:s}\n'.format(strcom))
fout.write('mcol,{0:d}\n'.format(mcol))
fout.write('ndata,{0:d}\n'.format(ndata))
for k in range(0,kk):
    fout.write('Cl({0:d}),{1:d}\n'.format(k,num[k]))
text='No.'
for j in range(0,iw):
    text=text+',x{0:d}'.format(j+1)
text=text+',cl\n'
fout.write(text)
for i in range(0,ndata):
    text='{0:d}'.format(i+1)
    for j in range(0,iw):
        text=text+',{0:g}'.format(xd[i,j])
    text=text+',{0:d}\n'.format(cl[i])
    fout.write(text)
fout.close()

# Calculation
a=np.corrcoef(xd.T)     # correlation matrix
ev,vec=np.linalg.eig(a) # Eigenvalue analysis
xx=np.dot(xd,vec)       # score

# Draw graph
icol=['red','yellow','blue','cyan','lime']
ag = np.arange(1,mcol*mcol+1,1).reshape(mcol,mcol)
fig=pl.figure(figsize=(8,8))
for i in range(0,mcol):
    for j in range(0,mcol):
        ax=fig.add_subplot(mcol,mcol,ag[i,j])
        pl.xticks([])
        pl.yticks([])
        if i==j:
            s='S'+str(i+1)
            pl.hist(xx[:,i],color='r',alpha=0.3)
            ax.text(0.95,0.95,s,fontsize=20,fontweight='bold',color='black',ha='right',va='top',transform=ax.transAxes)
        elif j<i:
            xmin=min(xx[:,j])
            xmax=max(xx[:,j])
            ymin=min(xx[:,i])
            ymax=max(xx[:,i])
            sx=xmax-xmin
            sy=ymax-ymin
            pl.xlim(xmin-0.1*sx,xmax+0.1*sx)
            pl.ylim(ymin-0.1*sy,ymax+0.1*sy)
            for ii in range(0,ndata):
                pl.plot(xx[ii,j],xx[ii,i],'o',color=icol[cl[ii]])
        else:
            xmin=min(xx[:,i])
            xmax=max(xx[:,i])
            ymin=min(xx[:,j])
            ymax=max(xx[:,j])
            sx=xmax-xmin
            sy=ymax-ymin
            pl.xlim(xmin-0.1*sx,xmax+0.1*sx)
            pl.ylim(ymin-0.1*sy,ymax+0.1*sy)
            for ii in range(0,ndata):
                pl.plot(xx[ii,i],xx[ii,j],'o',color=icol[cl[ii]])
pl.subplots_adjust(hspace=0)
pl.subplots_adjust(wspace=0)
pl.savefig("fig_sta_KMP.png", dpi=100)
pl.show()