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.ones([ndata,mcol+1],dtype=np.float)
yd=np.zeros([ndata],dtype=np.float)
for i in range(0,ndata):
    text=fin.readline()
    text=text.strip()
    text=text.split(',')
    yd[i]=float(text[0])
    for j in range(0,mcol):
        xd[i,j+1]=float(text[j+1])
fin.close()

# Calculation
a=np.dot(xd.T,xd)
b=np.dot(yd,xd)
c=np.linalg.solve(a,b)
ye=np.dot(xd,c)
rr=np.corrcoef(yd,ye)[0][1]

# Data output
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,mcol+1):
    fout.write('c({0:d}),{1:g}\n'.format(k,c[k]))
fout.write('i,Obs.,Est.\n')
for i in range(0,ndata):
    fout.write('{0:d},{1:g},{2:g}\n'.format(i+1,yd[i],ye[i]))
fout.close()

# Draw graph(1)
xmin=0
ymin=0
xmax=120000
ymax=120000
fig=pl.figure()
ax = pl.gca()
pl.subplot(1,1,1)
pl.xlabel('Observed data',fontsize=12)
pl.ylabel('Estimated data',fontsize=12)
pl.xlim(xmin,xmax)
pl.ylim(ymin,ymax)
aspect = 1.0*(ax.get_xlim()[1] - ax.get_xlim()[0]) / (ax.get_ylim()[1] - ax.get_ylim()[0])
ax.set_aspect(aspect)
pl.plot(yd,ye,'ro')
pl.plot([xmin,xmax],[ymin,ymax],'k-',lw=0.5)
txt='r={0:.3f}'.format(rr)
pl.text(0.05*(xmax-xmin),0.95*(ymax-ymin),txt,fontsize=16,ha='left',va='top') 
pl.grid(True,which='both',ls='-', color='0.65')
pl.savefig("fig_sta_MRA1.png", dpi=100)
pl.show()

# Draw graph(2)
xd[:,0]=yd[:]
mc=mcol+1
ag = np.arange(1,mc*mc+1,1).reshape(mc,mc)
fig=pl.figure(figsize=(8,8))
for i in range(0,mc):
    for j in range(0,mc):
        ax=fig.add_subplot(mc,mc,ag[i,j])
        pl.xticks([])
        pl.yticks([])
        if i==j:
            if i==0:
                txt='y'
            else:
                txt='x'+str(i)
            ax.text(0.5,0.5,txt,fontsize=20,fontweight='bold',color='black',ha='center',va='center',transform=ax.transAxes)
        elif j<i:
            xmin=min(xd[:,j])
            xmax=max(xd[:,j])
            ymin=min(xd[:,i])
            ymax=max(xd[:,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(xd[:,j],xd[:,i],'ro')
        else:
            rr=np.corrcoef(xd[:,i],xd[:,j])[0][1]
            txt='{0:.3f}'.format(rr)
            ax.text(0.5,0.5,txt,fontsize=16,fontweight='bold',color='black',ha='center',va='center',transform=ax.transAxes)
pl.subplots_adjust(hspace=0)
pl.subplots_adjust(wspace=0)
pl.savefig("fig_sta_MRA2.png", dpi=100)
pl.show()