from Gnuplot import Gnuplot, Data
from math import sqrt
from random import gauss

#HERE THE DATA.DAT file FOR THE EXERCISE IS CREATED
#Create Data Points
data_file=open("data3.dat","w")
data_file.write("### First Data Set\n")
for i in range(6450):
	num=gauss(14.5,2.5)
	data_file.write("%4.2f\n"%(num))
data_file.write("### Second Data Set\n")
for i in range(20000):
	num=gauss(12.5,1.5)
	data_file.write("%4.2f\n"%(num))
data_file.write("### Thrid Data Set\n")
for i in range(300):
	num=gauss(8.5,0.5)
	data_file.write("%4.2f\n"%(num))
data_file.close()

#EXERCISE 1
### Read Data Sets from File
data_points=[]
data_file=open("data.dat","r")
for line in data_file:
	if "###" in line:
		data_points.append([])
	else:
		num=float(line)
		data_points[-1].append(num)

# Plot Data Points
gp=Gnuplot()
gpdata1=Data(range(len(data_points[0])),data_points[0],with_="points",title="Set 1")
gpdata2=Data(range(len(data_points[1])),data_points[1],with_="points",title="Set 2")
gpdata3=Data(range(len(data_points[2])),data_points[2],with_="points",title="Set 3")
gp("set xlabel 'x'; set ylabel 'y'")
gp.plot(gpdata1,gpdata2,gpdata3)


#EXERCISE 2
#Create Bins for Data and empty histogram
Bins=[]
bin_data=[[], [], []]
for N in range(0,300):
	Bins.append(N/10.)
	for k in range(3): #This creates an empty histogram filled with 0's
		bin_data[k].append(0)
Bin_width=Bins[1]-Bins[0]

#Fill histogram
for k in range(3):
	for point in data_points[k]:
		#Calculate the Bin Number for the data_point corresponding 
		#to Bins
		bin_number=long(point/Bin_width+0.5)
		bin_data[k][bin_number]+=1

#Plot Histogram
gp2=Gnuplot()
gpdata1=Data(Bins,bin_data[0],with_="histeps",title="Binned Set 1")
gpdata2=Data(Bins,bin_data[1],with_="histeps",title="Binned Set 2")
gpdata3=Data(Bins,bin_data[2],with_="histeps",title="Binned Set 3")
gp("set xlabel 'x'; set ylabel 'y'")
gp2.plot(gpdata1,gpdata2,gpdata3)


#EXERCISE 3
#Derive mu and sigma for histogram
for k in range(3):
	total_A=sum(bin_data[k])
	mu=sum(data_points[k])/len(data_points[k])
	sig=0
	for point in data_points[k]:
		sig+=(point-mu)**2
	sig=sqrt(1./(len(data_points[k])-1)*sig)
	gp2("f%s(x)=%s/(%s*sqrt(2*pi)) * exp(-0.5*((x-%s)/%s)**2)*%s"%(k,total_A,sig,mu,sig,Bin_width))
	gp2.replot("f%s(x) title 'Gauss Set %i'"%(k,k+1))
raw_input()