import numpy as np
import sys
from math import *
import matplotlib.pyplot as pl
 
#Main routine
param=sys.argv
fnameR=param[1]
fnameW=param[2]

# 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()

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

# 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))
text='Number'
for j in range(0,iw-1):
    text=text+',{0:d}'.format(j+1)
text=text+',{0:d}\n'.format(iw)
fout.write(text)
text='Eigenvalue'
for j in range(0,iw-1):
    text=text+',{0:g}'.format(ev[j])
text=text+',{0:g}\n'.format(ev[iw-1])
fout.write(text)
text='Contribution'
for j in range(0,iw-1):
    text=text+',{0:g}'.format(ev[j]/np.sum(ev))
text=text+',{0:g}\n'.format(ev[iw-1]/np.sum(ev))
fout.write(text)
for i in range(0,mcol):
    text='Eigenvector({0:d})'.format(i+1)
    for j in range(0,iw-1):
        text=text+',{0:g}'.format(vec[i,j])
    text=text+',{0:g}\n'.format(vec[i,iw-1])
    fout.write(text)
fout.write('Score\n')
for i in range(0,ndata):
    text='{0:d}'.format(i+1)
    for j in range(0,iw-1):
        text=text+',{0:g}'.format(xx[i,j])
    text=text+',{0:g}\n'.format(xx[i,iw-1])
    fout.write(text)
fout.close()

# Draw graph
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)
            pl.plot(xx[:,j],xx[:,i],'ro')
        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)
            pl.plot(xx[:,i],xx[:,j],'ro')
pl.subplots_adjust(hspace=0)
pl.subplots_adjust(wspace=0)
pl.savefig("fig_sta_PCA.png", dpi=100)
pl.show()