127 lines
3.4 KiB
Python
127 lines
3.4 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)
|
|
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)
|
|
|
|
sim_ic = IntegrationCollection(4e-7,0.999,1e-200,16)
|
|
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
|