#!/usr/bin/env python

import sys
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#~ if len(sys.argv) < 2:
	#~ print 'Usage:\npython <script> <monitorFile> \n'
	#~ print 'Output:\n a more human readbale file'
	#~ exit()

class Coords:
	def __init__(self, x, y, z):
		self.x = x
		self.y = y
		self.z = y

	def __repr__(self):
		return "Coords: {}, {}, {}".format(self.x, self.y, self.z)

FeBio={}
Sofa={}
for file in sys.argv:
	if 'BoxPositions_clean' in file:
		f = open(file, 'r')
		#read data
		line=f.readline()				# *Title:
		line=f.readline()				# *Step = 1
		while line:
			# these line repeat at each time
			t=float(f.readline().split()[-1])	# *Time = 0.2
			line=f.readline()			# *Data = x;y;z
			FeBio[t]={}
			line=f.readline()
			while not line.startswith('*'):	# *Step = N   that is the next timestep
				if line:
					line=line.split()
					#~ print line
					#~ print tuple(float(i) for i in line[1:4])
					FeBio[t][int(line[0][:-3])-319] = Coords(*tuple(float(i) for i in line[1:4]))
					line=f.readline()
				else:	break
		#
		f.close()
	elif 'BoxMonitorPosition' in file:
		f = open(file, 'r')
		line=f.readline().split()	# # Gnuplot File : positions of N particle(s) Monitored
		N=int(line[6])
		line=f.readline()		# # 1st Column : time, others : particle(s) number 0 1 2...(N-1) 
		line=f.readline()		# First line of data
		while line:
			line=line.split()
			t=float(line[0])
			Sofa[t]={}
			for i in xrange(len(line)/3):
				Sofa[t][i+1]=Coords( *tuple([float(j) for j in line[i*3+1:i*3+4]]) )
			line=f.readline()
		f.close()
	elif '.py' in file:
		pass
	else:
		print 'Error'
		exit()

#~ # Check difference between first and last values
#~ for n in FeBio[0.02].iterkeys():
	#~ diff = ((FeBio[0.02][n].x-FeBio[1.0][n].x)**2+(FeBio[0.02][n].y-FeBio[1.0][n].y)**2+(FeBio[0.02][n].z-FeBio[1.0][n].z)**2)**0.5
	#~ sDiff = ((Sofa[0.02][n].x-Sofa[1.0][n].x)**2+(Sofa[0.02][n].y-Sofa[1.0][n].y)**2+(Sofa[0.02][n].z-Sofa[1.0][n].z)**2)**0.5
	#~ print '{}	{}	{}'.format(n, diff, sDiff)




highTime=min( max(FeBio.keys()), max(Sofa.keys()) )	# highest time in both lists...hopefully

histData={}
highest=0

for t,d in sorted(Sofa.iteritems()):
	if t<=highTime:
		histData[t]=[]
		for n,S in sorted(d.iteritems()):
				F=FeBio[t][n]
				diff = ((F.x-S.x)**2+(F.y-S.y)**2+(F.z-S.z)**2)**0.5
				if diff > highest:	highest=diff
				histData[t].append(diff)

histNum = len(histData)
print 'histNum', histNum

order=1
if highest-round(highest,order) < 5*10**-(order+1):
		highest = highest+5*10**-(order+1)
highest=round(highest,order)

mult=10

threeD = False
if threeD:
	bins = np.linspace(0, highest, mult*(highest*10)+1)
	print 'bins'
	print bins
	#
	# mgrid in this case works backwards...the first dimension will be y; the second dimension wil be x
	grid = np.mgrid[0:1:histNum*1j, 0:(highest-10**-order/mult):mult*(highest*10)*1j]
	x_data= grid[1].flatten()
	y_data= grid[0].flatten()
	print x_data
	print y_data
	#
	z_data=np.empty([0])
	for i in sorted(histData.keys()):
		hist, edges = np.histogram(histData[i], bins)
		z_data=np.append(z_data, hist)
	#
	print 'x_data', len(x_data)
	print 'y_data', len(y_data)
	print 'z_data', len(z_data)
	#
	colorList = ['red', 'green', 'blue', 'aqua', 'burlywood']
		  #~ , 'cadetblue', 'chocolate', 'cornflowerblue',
		  #~ 'crimson', 'darkcyan', 'darkgoldenrod', 'darkgreen',
		  #~ 'purple', 'darkred', 'darkslateblue', 'darkviolet']
	barColors=[]
	for i in range(histNum):
		for n in (np.linspace(0, highest-10**-order/mult, mult*(highest*10))):
			barColors.append(colorList[i%5])
	#
	fig = plt.figure()
	ax = fig.add_subplot(111, projection='3d')
	#
	ax.bar3d(	x_data, y_data, np.zeros(len(x_data)), 0.1/mult, 0.02, z_data, color=barColors, alpha = 0.3)
	ax.set_title('3D Histogram of Differences') 
	ax.set_xlabel('Difference')
	ax.set_ylabel('Time')
	ax.set_zlabel('Frequency of Difference')
	
	#~ plt.xlabel('Smarts')
	#~ plt.ylabel('Probability')
	#~ plt.title(r'$\mathrm{Histogram\ of\ IQ:}\ \mu=100,\ \sigma=15$')
	#~ plt.axis([0, highest, 0, 1])
	#~ plt.legend( sorted(histData.keys()) )
	plt.grid(True)
else:
	bins = np.linspace(0, highest, mult*(highest*10)+1)
	print mult*(highest*10)+1
	#~ print bins
	print 1.0*(histNum-1)/histNum
	#~ grid = np.mgrid[0:1:histNum*1j, 0:(highest):mult*(highest*10)*1j]
	#~ grid = np.mgrid[0:(1.0*(histNum-1)/histNum):histNum*1j, 0:(highest-10**-order/mult):mult*(highest*10)*1j]
	#~ print grid
	#~ x_data = grid[0].flatten()
	#~ other_data= grid[1].flatten()
	#~ print 'xData'
	#~ print 'len', len(x_data)
	#~ print x_data
	#~ print 'other_data'
	#~ print other_data
	x_data = np.linspace(0, 1, histNum*(mult*highest*10))
	#
	#~ y_data=np.empty([0])
	y_data=np.ones( (histNum,mult*highest*10) )
	#~ print y_data
	
	#~ print sorted(histData.keys())
	#~ print sorted(histData.keys(), reverse=True)
	counter=0
	for i in sorted(histData.keys(), reverse=True):
		hist, edges = np.histogram(histData[i], bins)
		#~ print int(i*histNum)
		y_data[counter]=hist
		counter+=1
	#
	print y_data
	#~ print len(x_data), len(y_data)
	#~ plt.hist2d(x_data, y_data, bins=(len(histData)-1 , mult*(highest*10)+1) )#, norm=LogNorm())
	#~ plt.hist2d(other_data, y_data, bins=(len(histData)-1 , mult*(highest*10)+1) )#, norm=LogNorm())
	#~ plt.colorbar()
	#~ plt.xlabel('Time')
	#~ H, xedges, yedges = np.histogram2d(x_data, y_data, bins=(histNum,mult*(highest*10)+1) )
	#~ print H
	#~ print 'len H', len(H)
	#~ print xedges
	#~ print yedges
	#~ plt.imshow(y_data)
	plt.imshow(y_data, extent=[0, highest, 0, 1])
	#~ plt.imshow(y_data, extent=[0, 50, 0, 78])
	plt.colorbar()
	plt.xlabel('Difference')
	plt.ylabel('Time')
	#~ plt.axis([0, highest, 0, 1])

plt.show()