194 lines
6.0 KiB
Python
194 lines
6.0 KiB
Python
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()
|