import numpy as np
from sys import argv
import matplotlib.pyplot as plt
import seaborn as sns
from math import ceil

if len(argv) < 2:
	print 'expects fnames'
	quit()
SKIP = 0
fnames = []

print argv
for i in range(len(argv)-1):
	if '=' not in argv[i+1]:
		fnames.append(argv[i+1])


col1 = 1
col2 = 14 # 13-17 for 100,200,300,400,500
col3 = 12 # sigma
col4 = 11 # balanced state frequency
xpr=3
ypr=3


output = 'test.eps'
remove_cont = 0
MAXROWS =50

ysc = 1
xsc = 0.01
offset = 0
for things in argv:
	if len(things.split('='))>1:
		temp = things.split('=')
		try:
			exec(things)
		except:
			temp[1] = str(temp[1])
			exec('%s=\'%s\'' % (temp[0],temp[1]))

i=0
sa = []
xa=[]
ya=[]
za=[]
qa=[]

fig, ax = plt.subplots()
from scipy.signal import savgol_filter
i=0
color =['red','green','red','black']
for fname in fnames:
	xd, yd, zd, qd = np.loadtxt(fname, usecols=(col1,col2,col3,col4),unpack=True)
	dln = len(xd)
	dhalf = dln#int(dln/2)+1
	xa.append(xd[0:dhalf]*xsc)
	#ya.append(np.multiply(yd[0:dhalf],qd[0:dhalf])) # actual area per 50 ms
	ya.append(yd[0:dhalf]) # area/mean
	za.append(zd[0:dhalf]) # std of balanced state
	qa.append(qd[0:dhalf]) # mean of balanced state
	
	#sd = zd[0:dhalf]
	za[i] = np.divide(za[i],qa[i]) # normalized std
	#za[i]=np.multiply(za[i],qa[i])
	#np.multiply(sd,qa[i])
	print fname
	#yhat = savgol_filter(yd,21,3)
	#shat = savgol_filter(sd,21,3)
	dl = len(xd)
	#plt.fill_between(xd,np.subtract(yhat,shat),np.add(yhat,shat),alpha=0.50)
	plt.plot(xa[i],ya[i],color=color[i])
	plt.fill_between(xa[i],np.subtract(ya[i],za[i]),np.add(ya[i],za[i]),alpha=0.5,color=color[i])
	#plt.fill_between(xa[i],np.subtract(ya[i],za[i]),np.add(ya[i],za[i]),alpha=0.5,color=color[i])

	#plt.errorbar(xa[i],ya[i],yerr=sd,elinewidth=0.5,fmt='.',markersize=2)
	i+=1
	#plt.plot(xd,yd,'.')
#plt.subplots_adjust(left=0.2,bottom=0.3,right=0.8,top=0.9)
ax.set(xlim=(1,50),ylim=(0,1.2))
plt.savefig(output, format='eps')
plt.show()