#!/usr/local/bin/pythonw
# coding: iso-8859-1
# Copyright 2009 by Amos Newcombe
# 
# This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License along with this program.  If not, see <http://www.gnu.org/licenses/>.

'''Kuiper\'s variant of the Kolmogorov-Smirnov test

_Numerical Recipes in C_ pp 626-7
'''

from Statistics import Statistics
from mathlib.Probability.Uniform import Uniform
from math import exp, sqrt
from operator import sub

class KSKp1(Statistics):
	'Compare a sample with a theoretical distribution.'

	def __init__(self, lxData, prob=None):
		if prob == None: prob = Uniform()
		Statistics.__init__(self, lxData, prob)
	
	def Calculate(self):
		prob = self.tpArgs[0]
		cData = len(self.getData())
		self['n'] = cData
		self.getData().sort()
		lxP = map(prob.Distribution, self.getData())
		lxY = [float(i)/float(cData) for i in range(cData)]
		lxD = map(sub, lxP, lxY)
		xV = 1.0/cData - min(lxD) + max(lxD)
		self['maxdev'] = xV
		xRootN = sqrt(cData)
		self['sig'] = Q(xV * (xRootN + 0.155 + 0.24/xRootN))
		return self

def Q(xL):
	if xL < 0.4: return 1.0
	if 9.0 < xL: return 0.0
	epsTerm, epsSum = 1.0e-3, 1.0e-8
	xLT = 2.0 * xL * xL
	i = 1
	xSum, xTermAbs = 0.0, 0.0
	while i < 100:
		xLTi = xLT * i * i
		xTerm = 2.0 * (2.0*xLTi - 1) * exp(-xLTi)
		xSum += xTerm
		if xTerm != 0.0:
			# i*xL == 0.5 => xTerm == 0, but don't cut it off for that
			if abs(xTerm) <= epsTerm * xTermAbs: break
			if abs(xTerm) <= epsSum  * xSum    : break
		xTermAbs = abs(xTerm)
		i = i+1
	else: xSum = 1.0
	return xSum

if __name__ == '__main__':
	
	from mathlib.Probability.Gaussian import StdGaussian
	from StatTest import Test
	
	print 'Q:', map(Q, [0.0, 0.5, 1.0, 2.0, 4.0, 8.0, 16.0])
	
	def getStat(prob, cSample):
		return KSKp1([prob.Sample() for n in range(cSample)], prob)
	
	cSample = 128
	prob = StdGaussian()
	cReps = 32
	Test(getStat(prob, cSample), globals(), locals())
	print 'Summary significance:', KSKp1([
		getStat(prob, cSample).Calculate()['sig']
		for i in range(cReps)
	]).Calculate()['sig']