from math import *
	

#--------------------------------------------------------------------------------------------
def minichi(method,Meas,Upper,Eavg):	
	nm=len(Meas)
	nu=len(Upper)
	
	#no useful results of any kind:
	if nm == 0 and nu == 0: return None
		
	#Upper limits only:  choose the lowest one.	
	if nm == 0 and nu >0: 
		ul=99999.9
		#print "In minichi UL section"
		for E in Upper:
			#print "type(E.data)=",type(E.data[0])
			#print "E.data = ",E.data[0]
			#print "Considering E=%s in minichi; current UL threshold is %6.2f, while E.data = %6.2f"%(E.name,ul,E.data[0])
			if E.data[0] < ul: 
				ul=E.data[0]
				Eavg.tex = E.tex
				Eavg.avg = E.data[0]
				Eavg.newtex = Eavg.tex  # May 13th, 2004
				#print "minichi: Eavg.newtex=",Eavg.newtex
		return None
		
	#A single Measurement available: choose it.	
	if nm == 1: 
		E = Meas[0]
		Eavg.avg = E.data
		return 1
	
	#All remaining cases have more than one central value measurement and need to be averaged.

	#initialize the interpolated error (at the center); (and positivize the neg err)
	Measurements=[E.data for E in Meas]
	for measurement in Measurements:
		measurement[2]=fabs(measurement[2]) # make the low-side errors positive numbers
		evaluateItHere=measurement[0]
		error,weight=InterpolateError(evaluateItHere,measurement)
		measurement.append(error)
		measurement.append(weight)
		#change any 0 error into 1/10 its companion error
		if measurement[2] < 0.0001: measurement[2] = 0.1*measurement[1]
		if measurement[1] < 0.0001: measurement[1] = 0.1*measurement[2]
		for m in measurement:
			if type(m) != type(0.0): print "type(measurement)=",type(measurement), measurement
	print "mean = %7.3f      sigR = %7.3f       sigL = %7.3f"%tuple(measurement[:3])


	#Below each averaging algorithm is implemented in turn:

	print "Averages:",50*"-"
	print "Original tex=     ", Eavg.tex

	#0th algorithm
	avg,sigL,sigR = method0(Measurements)
	err = sigR  # symmetric errors
	if method == 0: Eavg.avg=[avg,sigR,sigL]
	print "Method 0: %10.2f +/- %6.2f"%(avg,err)#, "     Chisq=%6.2f"%chisq
	
	#1st algorithm
	#newavg,sigL,sigR = method1(avg, Measurements)
	#if method == 1: Eavg.avg=[newavg,sigR,sigL]
	#format = "Method 1: %10.2f (%+6.2f %+6.2f)     "
	#print    format%(newavg,sigR,sigL)#,"     Chisq=%6.2f"%chisq
		
	#2nd algorithm
	newavg,sigL,sigR=method2(avg,Measurements)	
	if method == 2: Eavg.avg=[newavg,sigR,sigL]
	format = "Method 2: %10.2f (%+6.2f %+6.2f)     "
	print    format%(newavg,sigR,sigL)#,"     Chisq=%6.2f"%chisq
	
	#Look for significant differences w/r 0th algrithm.
	try:  #identify cases where this algo is 5% or more different from previous
		da = abs(newavg-avg)/avg
		dl = abs(abs(sigL)-err)/err
		dr = abs(sigR-err)/err
		if da>0.05 or dl > 0.05 or dr > 0.05: print "<<<<<<<<< noticeable diff w/r method 0 >>>>>>>>>"
	except: pass

	return 1
			
#--------------------------------------------------------------------------------------------
def method0(Measurements):
	for i in range(10):
		avg,err,chisq = getavg(Measurements)
		for measurement in Measurements:
			error,weight = InterpolateError(avg,measurement)
			measurement[3]=error
			measurement[4]=weight		
	return avg,-err,err
			
#--------------------------------------------------------------------------------------------
def method1(avg, Measurements):
	Left,Right=defineSearchWindow(avg,Measurements)
	#print id(myLikelihood)
	Normalize(Left,Right,myLikelihood)
	newavg,maxval = findMean(Left, Right, myLikelihood)
	sigL, sigR = findSigmas(Left,Right, Mean, myLikelihood)
	return newavg, sigL, sigR

#--------------------------------------------------------------------------------------------
def method2(avg, Measurements):
	Left,Right=defineSearchWindow(avg,Measurements)
	Normalize(Left,Right,myLikelihood)
	Mean,maxval = findMean(Left,Right,myLikelihood)
	sigL=SearchValueL(Left,Mean,Mean,myLikelihood)
	sigR=SearchValueR(Mean,Right,Mean,myLikelihood)
	return Mean,sigL-Mean,sigR-Mean
	
def SearchValueL(Left, Right, Mean, myLikelihood):
	target = myLikelihood(Mean) * exp(-0.5)  # equiv to finding -2lnL up +1
	x = (Left+Right)/2.0
	thisValue = myLikelihood(x)
	#print "Left = ",x, " thisValue = ", thisValue
	if thisValue > target*1.005: x=SearchValueL(Left, x,Mean,myLikelihood)
	if thisValue < target*0.995: x=SearchValueL(x,Right,Mean,myLikelihood)
	return x
def SearchValueR(Left, Right, Mean, myLikelihood):
	target = myLikelihood(Mean) * exp(-0.5)  # equiv to finding -2lnL up +1
	x = (Left+Right)/2.0
	thisValue = myLikelihood(x)
	#print "Right = ",x, " thisValue = ", thisValue
	if thisValue > target*1.005: x=SearchValueR(x,Right,Mean,myLikelihood)
	if thisValue < target*0.995: x=SearchValueR(Left, x,Mean,myLikelihood)
	return x
	

#--------------------------------------------------------------------------------------------
def defineSearchWindow(avg,Measurements):
	global gMeasurements,Anorm
	gMeasurements = Measurements
	largestsigma=-999.
	for meas in Measurements:
		if meas[1] > largestsigma: largestsigma=meas[1]
		if meas[2] > largestsigma: largestsigma=meas[2]
	Left  = avg - 5.0 * largestsigma
	Right = avg + 5.0 * largestsigma
	return Left,Right

#--------------------------------------------------------------------------------------------
def findMean(Left, Right, myLikelihood):
	global Mean
	maximum = -9999.
	npoints = 50
	dx = (Right-Left)/float(npoints)
	x = Left
	while x < Right:
		funcvalue = myLikelihood(x) 
		if funcvalue > maximum: 
			maximum = funcvalue
			xmax = x
		x += dx
	Mean = xmax
	#recursive loop to zero in on peak:
	Right= Mean + 5.0*dx
	Left = Mean - 5.0*dx
	if dx > 0.001*fabs(Mean): findMean(Left,Right,myLikelihood)
	return xmax,maximum
	
#--------------------------------------------------------------------------------------------
def getavg(Measurements):
	#Returns (centralvalue, error, chisq)
	#Measurements = [ m1, m2, m3, ...]
	#m1 = [x, dpx, dnx, dix, weight]
	#note "dix" = the interpolated error; weight = 1/dix^2
	#dix and weight depend on where you evaluate for them 
	#x, dpx, dnx are fixed results from an expt.
	#dnx is a positive number, |lo-side error|
	
	#Take weighted average of the measurements
	sum1 =0.0
	sum2 =0.0
	for measurement in Measurements:
		sum1 += measurement[4]*measurement[0]
		sum2 += measurement[4]
	avg = sum1/sum2
	err = 1.0/sqrt(sum2)
	
	#Get Chisq corresponding to this average
	chisq=0.0
	for measurement in Measurements:
		error,weight = InterpolateError(avg,measurement) #Use this error 
		chisq += weight * pow(measurement[0]-avg, 2) 
	chisq  = chisq/float(len(Measurements)-1)
	
	#assemble and return the result
	return avg,err,chisq
	

#--------------------------------------------------------------------------------------------
def InterpolateError(evaluateItHere,measurement):
	#returns error that linearly grades between lo and hi side errs
	x  =measurement[0]
	dpx=measurement[1]
	dnx=measurement[2] # note I have made this a positive number
	averageError = (dpx+dnx)/2.0
	AsymmetricError = fabs(dpx-dnx)/averageError
	if evaluateItHere < x-dnx: return dnx, pow(dnx, -2)
	if evaluateItHere > x+dpx: return dpx, pow(dpx, -2)
	#otherwise, compute the linearly interpolated error
	slope=(dpx-dnx)/(dpx+dnx)
	#print "InterpolateError:",x,dpx,dnx,slope,evaluateItHere
	error = dnx + slope*(evaluateItHere - (x-dnx))
	if error < 1.0E-4:
		print "Replacing zero error by max(dpx,dnx)",error,max(dpx,dnx) 
		error = max(dpx,dnx)
	print error
	weight = pow(error,-2)
	#print "Result of Interpolation of Error: ", error
	return error,weight


#--------------------------------------------------------------------------------------------
def SimpsonsRule(a,b,function, m=100):
	h = (b-a)/(2.0*float(m))
	integral = 0.0
	for i in range(m):
		fi=float(i)
		term1 = function(a + (2.0*fi+0.0)*h)
		term2 = function(a + (2.0*fi+1.0)*h)
		term3 = function(a + (2.0*fi+2.0)*h)
		integral += (h/3.0)* (term1 + 4.0*term2 + term3)
	return integral

#--------------------------------------------------------------------------------------------
def myLikelihood(x):
	global gMeasurements, Anorm	
	func = 1.0
	#                                               mean     sigL     sigR
	for meas in gMeasurements: func *= bifgauss(x, meas[0], meas[2], meas[1] )
	return Anorm * func

#--------------------------------------------------------------------------------------------
def bifgauss(x,mean,sigleft, sigright):
	'''normalized asymmetric gaussian'''
	factor = 2.0 / (sqrt(2.0*pi) * (sigleft + sigright))
	sigma = sigleft
	if x > mean: sigma=sigright
	return factor * exp(-(x-mean)*(x-mean)/ (2.0 * sigma*sigma))

#--------------------------------------------------------------------------------------------
def findSigmas(Left, Right, Mean, myLikelihood):
	l= findSigmaL(Left, Right, Mean, myLikelihood)
	r= findSigmaR(Left, Right, Mean, myLikelihood)
	#print "sigL: l, l-Mean = %6.2f, %6.2f"%(l, l-Mean)
	#print "sigR: r, r-Mean = %6.2f, %6.2f"%(r, r-Mean)
	return (l-Mean),(r-Mean)

def findSigmaL(Left, Right, Mean, myLikelihood):
	Normalize(Left, Mean, myLikelihood)
	return SearchL(Left, Mean, Mean, myLikelihood)

def findSigmaR(Left, Right, Mean, myLikelihood):
	Normalize(Mean, Right, myLikelihood)
	return SearchR(Mean, Right, Mean, myLikelihood)

def SearchL(Left, Right, Mean, myLikelihood):
	x = (Left+Right)/2.0
	thisArea = SimpsonsRule(x,Mean,myLikelihood)
	#print "Left = ",x, " thisArea = ", thisArea
	if thisArea > 0.684: x=SearchL(x,Right,Mean,myLikelihood)
	if thisArea < 0.682: x=SearchL(Left, x,Mean,myLikelihood)
	return x

def SearchR(Left, Right, Mean, myLikelihood):
	x = (Left+Right)/2.0
	thisArea = SimpsonsRule(Mean,x,myLikelihood)
	#print "Right = ",x, " thisArea = ", thisArea
	if thisArea > 0.684: x=SearchR(Left, x,Mean,myLikelihood)
	if thisArea < 0.682: x=SearchR(x,Right,Mean,myLikelihood)
	return x


#--------------------------------------------------------------------------------------------
def Normalize(Left, Right, myLikelihood):
	global Anorm
	Anorm = 1.0
	Anorm = 1.0/SimpsonsRule(Left,Right,myLikelihood)
	return Anorm



#--------------------------------------------------------------------------------------------
def plotfunc(f, xlo, xhi, nsteps, fp):
	dx = (xhi-xlo)/float(nsteps)
	x=xlo
	format="%12.4g       %12.4g\n"
	while x <= xhi:  
		fp.write(format%(x,f(x)))
		x += dx
def newplot(fp):
	fp.write("  \n   ")


#--------------------------------------------------------------------------------------------
def testplot():
	global gMeasurements, Anorm, Mean
	filename = '/Users/jima/LEPP/ANALYSIS/HFAG/JULY31_2003/test.txt'
	format = "%8.3f"
	print
	print

	fp = open(filename, "w")
	gMeasurements = [ [0.0, 1.0, 1.0] ]
	m = gMeasurements[0]
	Normalize(-5.0,10.0,myLikelihood)
	plotfunc(myLikelihood, -10.0, 10.0, 800, fp)
	Mean,Max = findMean(-2.0, 9.0, myLikelihood)
	sigLeft, sigRight = findSigmas(-5.0, 10.0, Mean, myLikelihood)
	print "Mean = ",format%Mean, " +/- ", format%sigLeft, format%sigRight
	
	newplot(fp)
	del gMeasurements
	gMeasurements= [ [1.0, 2.0, 2.0] ]
	Normalize(-5.0,10.0,myLikelihood)
	plotfunc(myLikelihood, -10.0, 10.0, 800, fp)
	Mean,Max = findMean(-4.0, 9.0, myLikelihood)
	sigLeft, sigRight = findSigmas(-5.0, 10.0, Mean, myLikelihood)
	print "Mean = ",format%Mean, " +/- ", format%sigLeft, format%sigRight
	
	newplot(fp)
	gMeasurements.append( m )
	Normalize(-5.0,10.0,myLikelihood)
	plotfunc(myLikelihood, -10.0, 10.0, 800, fp)
	Mean,Max = findMean(-4.0, 9.0, myLikelihood)
	sigLeft, sigRight = findSigmas(-5.0, 10.0, Mean, myLikelihood)
	print "Mean = ",format%Mean, " +/- ", format%sigLeft, format%sigRight
	
	fp.close()