gatk-3.8/python/theoPost.py

import math
import sys

#if ( sys.version_info < (3,0) ):
#	raise "Must use python version 3 or later. See /broad/software/free/Linux/redhat_5_x86_64/pkgs/python_3.1.2/bin/python3.1"

class controls:
	def __init__(self,l_x,l_y,c_x,c_y,r_x,r_y):
		self.x_a = l_x - 2*c_x + r_x
		self.y_a = l_y - 2*c_y + r_y
		self.x_b = -2*l_x + 2*c_x
		self.y_b = -2*l_y + 2*l_x
		self.x_c = l_x
		self.y_c = l_y
		
class IntegrationCollection:
	def __init__(self,start,stop,err,num_ints):
		self.start = start
		self.stop = stop
		self.err = err
		self.num_ints = num_ints

lfmap = dict()
def logfact(a):
	global lfmap
	if ( a < 2 ):
		return 0.0
	if ( a not in lfmap ):
		lfmap[a] = math.log10(a) + logfact(a-1)
	return lfmap[a]

def logchoose(a,b):
	return logfact(a)-logfact(b)-logfact(a-b)

def logbinomial(success,trials,prob):
	return logchoose(trials,success) + success*math.log10(prob) + (trials-success)*math.log10(1-prob)

def quad(a,b,c,x):
	return a*x*x + b*x + c

def qformula(a,b,c,equivVal):
	return (-b + math.sqrt(b*b-4*a*c))/(2*a)

def cbezierf(cts,pt):
	t = qformula(cts.x_a,cts.x_b,cts.x_c,pt)
	y = quad(cts.y_a,cts.y_b,cts.y_c,t)
	return y

bez_cts = controls(-7,10.99919,-2.849154,0.1444735,-0.0043648054,-1.559080) # based on previous gradient descent


def simpson(f,ic):
	class DeprecationError(Exception):
		def __init__(self,val):
			self.value = val
		def __str__(self):
			return repr(self.value)
	raise DeprecationError("Simpson is deprecated. Do not use it.")

def simpAux(f,a,b,eps,s,fa,fb,fc,cap):
	if ( s == 0 ):
		return []
	c = ( a + b )/2
	h = b-a
	d = (a + c)/2
	e = (c + b)/2
	fd = f(d)
	fe = f(e)
	s_l = (h/12)*(fa + 4*fd + fc)
	s_r = (h/12)*(fc + 4*fe + fb)
	s_2 = s_l + s_r
	if ( cap <= 0 or abs(s_2 - s) <= 15*eps ):
		try:
			return [math.log10(s_2 + (s_2 - s)/15.0)]
		except OverflowError:
			print(s_2)
			print(s_2-s)
			return [-350]
	return simpAux(f,a,c,eps/2,s_l,fa,fc,fd,cap-1) + simpAux(f,c,b,eps/2,s_r,fc,fb,fe,cap-1)

def adaptiveSimpson(f,start,stop,error,cap):
	mid = (start + stop)/2
	size = stop - start
	fa = f(start)
	fb = f(mid)
	fc = f(stop)
	s = (size/6)*(fa + 4*fc + fb)
	h = simpAux(f,start,stop,error,s,fa,fb,fc,int(cap))
	h.sort()
	#print("first: "+str(h[0]))
	#print("last: "+str(h[len(h)-1]))
	return sum(map(lambda x: 10**x,h))

def neutral(x):
	return -1.0*math.log10(x)

def twoState(x):
	if ( x < 0.04 ):
		return -1.5*math.log10(x) + 0.5*math.log10(0.04)
	else:
		return -1.0*math.log10(x)

def bezier(x):
	return cbezierf(bez_cts,math.log10(x))

norm_cache = (None,None)
def resampleProbability(logshape,ic,ac,ns,ac_new,ns_new):
	global norm_cache
	logpost = lambda x: logshape(x) + logbinomial(ac,2*ns,x)
	if ( norm_cache[1] == None or norm_cache[0] != (ac,ns,logshape) ):
		print("Caching posterior norm")
		norm_cache = ((ac,ns,logshape),math.log10(adaptiveSimpson( lambda v: math.pow(10,logpost(v)), ic.start,ic.stop,ic.err,ic.num_ints)))
	logpost_normed = lambda v: logpost(v) - norm_cache[1]
	newshape = lambda y: math.pow(10,logpost_normed(y) + logbinomial(ac_new, 2*ns_new, y))
	return adaptiveSimpson(newshape,ic.start,ic.stop,ic.err,ic.num_ints)

def getPost(logshape,ic,ac,ns):
	global norm_cache
	logpost = lambda x: logshape(x) + logbinomial(ac,2*ns,x)
	if ( norm_cache[1] == None or norm_cache[0] != (ac,ns,logshape) ):
		print("Caching posterior norm")
		norm_cache = ((ac,ns,logshape),math.log10(adaptiveSimpson(lambda v: math.pow(10,logpost(v)),ic.start,ic.stop,ic.err,ic.num_ints)))
	return lambda v: logpost(v) - norm_cache[1]

sim_ic = IntegrationCollection(5e-8,0.999,1e-2000,22)
sys.setrecursionlimit(int(2e6))
#neutral_post = map( lambda v: resampleProbability(neutral,sim_ic,1,900,v,900), range(0,21) ) 
#twostate_post = list(map( lambda v: resampleProbability(twoState,sim_ic,1,900,v,900), range(0,21) ))
#g = open("n_ts.txt",'w')
#idx = 0
#for e in neutral_post:
#	g.write(str(idx))
#	g.write("\t"+str(e)+"\t"+str(twostate_post[idx])+"\n")
#	idx += 1

DO_1 = False
if ( DO_1 ):
	eomiautism_ac_1 = 317763
	eomiautism_ac_2 = 78844
	eomiautism_ac_3p = 239526 # all of these go on chip by default

	new_set = 10000-917-998

	num_unseen_sites = 125*new_set

	unseen_unseen = resampleProbability(twoState,sim_ic,0,917+998,0,new_set)
	unseen_1 = resampleProbability(twoState,sim_ic,0,917+998,1,new_set)/(1-unseen_unseen)
	unseen_2 = resampleProbability(twoState,sim_ic,0,917+998,2,new_set)/(1-unseen_unseen)
	ac1_unseen = resampleProbability(twoState,sim_ic,1,917+998,0,new_set)
	ac1_ac1 = resampleProbability(twoState,sim_ic,1,917+998,1,new_set)
	ac2_unseen = resampleProbability(twoState,sim_ic,2,917+998,0,new_set)

	total = 636133 + num_unseen_sites
	ac1 = unseen_1*num_unseen_sites + ac1_unseen*eomiautism_ac_1
	ac2 = unseen_2*num_unseen_sites + ac1_unseen*eomiautism_ac_1 + ac2_unseen*eomiautism_ac_2

	print("\t".join(map(lambda u: str(u), [unseen_unseen,unseen_1,unseen_2])))
	print("\t".join(map(lambda u: str(u), [total,ac1,ac2])))

	ea_ns = 343877
	ea_ns_ac1 = 204223
	ea_ns_ac2 = 42280

	ns_new_ac1 = ea_ns_ac1*ac1_unseen + unseen_1*num_unseen_sites*(1.7/(1+1.7))
	ns_new_ac2 = ea_ns_ac2*ac2_unseen + unseen_2*num_unseen_sites*(1.4/(1+1.4))
	ns_new_total = ea_ns + ns_new_ac1 + ns_new_ac2 + (num_unseen_sites*(1-unseen_1-unseen_2))*(0.6/(1+0.6))

	print("\t".join(map(lambda u: str(u), [ns_new_total,ns_new_ac1,ns_new_ac2])))
	print(1-resampleProbability(twoState,sim_ic,2,1000,0,10000)-resampleProbability(twoState,sim_ic,2,1000,1,10000))
	print(1-resampleProbability(twoState,sim_ic,1,100,0,2000)-resampleProbability(twoState,sim_ic,1,100,1,2000))
	print(1-resampleProbability(twoState,sim_ic,2,100,0,2000)-resampleProbability(twoState,sim_ic,2,100,1,2000))
	print(1-resampleProbability(twoState,sim_ic,20,1000,0,2000)-resampleProbability(twoState,sim_ic,20,1000,1,2000))

def emitPosterior(ac):
	post_ac = getPost(twoState,sim_ic,ac,10000)
	o = open("post_%d.txt" % ac,'w')
	pt = sim_ic.start
	while ( pt < 0.2 ):
		o.write("%e\t%e\n" % (pt,post_ac(pt)))
		pt = 1.015*pt
	o.close()

emitPosterior(2)
emitPosterior(3)
emitPosterior(10)
emitPosterior(25)

#o = open("test2s.txt",'w')
#pt = sim_ic.start
#while ( pt<0.4 ):
#	o.write("%e\t%e\n" % (pt,twoState(pt)))
#	pt = 1.015*pt
#o.close()