import csc
import numpy as np
import os

thetas = np.arange(0,91,1)

files=[]
true_thetas = []
for theta in thetas:
	filename = 'Ba2NaOsO6-001-cant-010-'+str(int(theta))+'-Mag-histo.tsv'
	if os.path.exists(filename):
		meta,properties,f = csc.begin_read(filename)
		files.append(f)
		true_thetas.append(theta)
		del f

i=0
is_eof = False
meanfile = open('Ba2NaOsO6-001-cant-010-3d-Mag-means.tsv','w')
medianfile = open('Ba2NaOsO6-001-cant-010-3d-Mag-medians.tsv','w')
modefile = open('Ba2NaOsO6-001-cant-010-3d-Mag-modes.tsv','w')

for i in range(len(true_thetas)):
	mu = []
	pdf = []
	while True:
		values,error = csc.readline(files[i],properties)
		if error == 'EOF':
			break
		elif error == None:
			mu.append(values['mu'])
			pdf.append(values['pdf'])
	mu = np.array(mu,np.float)
	pdf = np.array(pdf,np.float)
	binwidth=mu[1]-mu[0] #could extract from metadata
	print thetas[i],np.sum(mu),np.sum(pdf)
	mean = np.sum(mu*pdf*binwidth)
	stdev = np.sum(mu**2*pdf*binwidth-mean**2)
	total_prob = 0.
	for j in range(len(mu)):
		prev_total_prob = total_prob
		total_prob += pdf[j]*binwidth
		if total_prob > 0.5:
			median = mu[j-1]+(mu[j]-mu[j-1])*(0.5-prev_total_prob)/(total_prob-prev_total_prob)
			break
	total_prob = 0.
	for j in range(len(mu)):
		prev_total_prob = total_prob
		total_prob += pdf[j]*binwidth
		if total_prob > 0.25:
			lq = mu[j-1]+(mu[j]-mu[j-1])*(0.25-prev_total_prob)/(total_prob-prev_total_prob)
			break
	total_prob = 0.
	for j in range(len(mu)):
		prev_total_prob = total_prob
		total_prob += pdf[j]*binwidth
		if total_prob > 0.75:
			uq = mu[j-1]+(mu[j]-mu[j-1])*(0.75-prev_total_prob)/(total_prob-prev_total_prob)
			break
	max_val=0.
	mode=0.
	for j in range(len(mu)):
		if pdf[j] > max_val:
			max_val = pdf[j]
			mode = mu[j]
	meanfile.write(str(true_thetas[i])+'\t'+str(mean)+'\t'+str(stdev)+'\n')
	medianfile.write(str(true_thetas[i])+'\t'+str(median)+'\t'+str(lq)+'\t'+str(uq)+'\n')
	modefile.write(str(true_thetas[i])+'\t'+str(mode)+'\n')