#! python

"""This computes the "Fraction significant" column in the
Speech Prosody 2010 paper.    It returns the fraction of
each variety of linear regression that is statistically
significant, and also prints some interesting details that
are not used in the paper.
"""

import sys
import math
import numpy
import random

from gmisclib import die
from gmisclib import dictops
from gmisclib import weighted_percentile as WP

sys.path.insert(0, '.')
import cross_eval
import model5



def print_N_statistics(d, prefix=''):
	if len(d) < 50:
		return
	ys = 0
	yss = 0
	ns = 0
	nss = 0
	Ns = []
	Nns = []
	for datum in d:
		n = int(datum['N'])
		if int(datum['__marked__']):
			ys += n
			yss += n*n
			Ns.append(n)
		else:
			ns += n
			nss += n*n
			Nns.append(n)
	yn = float(len(Ns))
	if yn < 10:
		return
	ybar = float(ys)/yn
	ysigma = math.sqrt((yss - ys*(ys/yn))/(yn-1))
	ybarsigma = ysigma / math.sqrt(yn)
	nn = float(len(Nns))
	if nn < 10:
		return
	nbar = float(ns)/float(nn)
	nsigma = math.sqrt((nss - ns*(ns/nn))/(nn-1))
	nbarsigma = nsigma / math.sqrt(nn)
	if abs(ybar-nbar)/math.hypot(ybarsigma, nbarsigma) > 3.0:
		print "N: %sSig Nbar=%.1f+-%.1f Nsig Nbar=%.1f+-%.1f" % (
			prefix, ybar, ybarsigma, nbar, nbarsigma
			)


def lang_statistics(d, overall):
	sig = 0
	n = 0
	for datum in d:
		sig += int(datum['__marked__'])
		n += 1
	s = float(sig)/float(n)
	sigma_s = math.sqrt(s*(1.0-s)/n)
	if s >= overall**0.2 and s-overall>sigma_s*3:
		return (1, s)
	elif s <= overall**5 and overall-s>sigma_s*3:
		return (-1, s)
	return (0, s)


def filterlist(x, lbl, target):
	return [q for q in x if q['lbl']==lbl and q['target']==target ]


def filtersplit(dlist, lbl, target):
	d = filterlist(dlist, lbl, target)
	per_lang = dictops.dict_of_lists()
	for info in d:
		per_lang.add(info['language'], info)
	return per_lang



def process(per_lang):
	alld = []
	for (lang, dl) in per_lang.items():
		alld.extend(dl)
	print_N_statistics(alld, '')
	overall = 0
	n = 0
	for (lang, dl) in per_lang.items():
		for datum in dl:
			overall += int(datum['__marked__'])
			n += 1
	print 'overall=%.3f' % (float(overall)/float(n))
	overall = float(overall)/float(n)
	for (lang, dl) in per_lang.items():
		side, s = lang_statistics(dl, overall)
		if side:
			print 'Lang', lang, side, s, overall
		print_N_statistics(dl, '%s ' % lang)
	



def replica_langrank(esLN):
	ll = []
	for (lang, esN) in esLN.items():
		esA = [random.choice(elist) for (n,elist) in esN.items()]
		ll.append( (max(esA), lang) )
	ll.sort()
	# return ''.join([l[0] for (e,l) in ll])
	return "%s%s" % (ll[0][1][0], ll[-1][1][0])


def lang_rank(esLN):
	ll = []
	for (lang, esN) in esLN.items():
		esA = [numpy.mean(elist) for (n,elist) in esN.items()]
		ll.append( (max(esA), lang) )
	ll.sort()
	# return ''.join([l[0] for (e,l) in ll])
	return "%s%s" % (ll[0][1][0], ll[-1][1][0])


def mark_effect_size(per_lang):
	ll = []
	esLN = {}
	for (lang, dl) in per_lang.items():
		esN = dictops.dict_of_lists()
		for d in dl:
			esN.add(int(d['N']), float(d['effect_size']))
		esLN[lang] = esN
		esA = [numpy.mean(elist) for (n,elist) in sorted(esN.items()) ]
		ll.append( (max(esA), lang)  )
		es = list([float(d['effect_size']) for d in dl])
		print '\t%s max_effect_size=%.3f' % (lang, max(esA))
		pctl = WP.wp(es, None, [0.95])[0]
		eslim = min(0.95*max(es), pctl)
		for d in dl:
			d['__marked__'] = float(d['effect_size']) > eslim
		ll.sort()
	
	NREP = 1000
	langrank = lang_rank(esLN)
	# langrank = ''.join([l[0] for (e,l) in ll])
	matchrank = 0
	for i in range(NREP):
		if replica_langrank(esLN) == langrank:
			matchrank += 1
	mark = "%.3f" % (float(matchrank)/float(NREP))
	# mark = '*' * ((matchrank>0.25*NREP)+(matchrank>0.75*NREP)+(matchrank>0.9*NREP))
	mean_effect_size = numpy.mean([e for (e,l) in ll])
	stdev_effect_size = numpy.std([e for (e,l) in ll])*math.sqrt(float(len(ll))/float(len(ll)-1))
	print 'EffectSize: %.3f LangRank: %s%s' % (mean_effect_size, langrank, mark)
	print 'fractional_stdev_of_effect_size_over_langs:', stdev_effect_size/mean_effect_size
	print 'effect_size: max/min:', ll[-1][0]/ll[0][0]


def mark(per_lang, marktype):
	if marktype == 'significant':
		for (lang, dl) in per_lang.items():
			for d in dl:
				d['__marked__'] = d['significant']
	elif marktype == 'effect_size':
		mark_effect_size(per_lang)
	else:
		die.die('Whoops')
	return per_lang



def loop(root, what_to_do):
	dlist = cross_eval.read_models(root)
	for (lbl, target) in model5.LTlist:
		print 'lbl=', lbl, 'target=', target
		process(mark(filtersplit(dlist, lbl, target), what_to_do))
		print


if __name__ == '__main__':
	loop('.', sys.argv[1])
