diff options
Diffstat (limited to 'hawkes')
| -rw-r--r-- | hawkes/main.py | 96 | ||||
| -rw-r--r-- | hawkes/refine.py | 47 |
2 files changed, 87 insertions, 56 deletions
diff --git a/hawkes/main.py b/hawkes/main.py index a922812..cddc85b 100644 --- a/hawkes/main.py +++ b/hawkes/main.py @@ -23,7 +23,9 @@ def gss(f, a, b, tol=1e-20): a = c c = d d = a + GR * (b - a) - return (b + a) / 2 + sys.stderr.write("gss:" + " ".join(map(str, [a, b, fc, fd])) + "\n") + sys.stderr.flush() + return (b + a) / 2, fc def iter_events(events): @@ -55,42 +57,63 @@ def ll(lamb, alpha, mu): return -(r1 - r2 - r3) -def sa(x, y, z, sigma=0.1, niter=100): - T = 1 - e = 1.01 - for _ in xrange(niter): +def sa(x, y, z, sigma=0.5, niter=70, fc=None): + T = 0.1 + e = 1.1 + if fc: + fo = fc + else: fo = ll(x, y, z) - #print T, sigma, x, y, z, fo + for _ in xrange(niter): + sys.stderr.write("sa: " + " ".join(map(str, [T, sigma, x, y, z, fo])) + + "\n") + sys.stderr.flush() yn = max(y + gauss(0, sigma * y + 1e-10), 0) zn = max(z + gauss(0, sigma * z + 1e-10), 0) fn = ll(x, yn, zn) if fn < fo or exp((fo - fn) / T) > random(): y = yn z = zn - sigma *= e + sigma *= 2 + fo = fn else: sigma /= e if sigma < 1e-5: break T *= 0.99 - return y, z + return y, z, fo -def optimize_with_sa(x, y, z): +def optimize_with_sa(x, y, z, niter=200): + + def f(x): + return ll(x, y, z) + # y, z = sa(x, y, z) - while True: - x = gss(f, 0, 1, tol=1e-10) - y, z = sa(x, y, z) - print x, y, z, ll(x, y, z) + y, z, fc = sa(x, y, z) + for _ in xrange(niter): + x, fc = gss(f, 0, 1e-3, tol=1e-10) + y, z, fc = sa(x, y, z, fc=fc) + print x, y, z, fc sys.stdout.flush() -def optimize_with_gss(x, y, z): - while True: - x = gss(f, 0, 1, tol=1e-20) - y = gss(g, 0, 1, tol=1e-20) - z = gss(h, 0, 1, tol=1e-20) - print x, y, z, ll(x, y, z) +def optimize_with_gss(x, y, z, niter=100): + + def f(x): + return ll(x, y, z) + + def g(y): + return ll(x, y, z) + + def h(z): + return ll(x, y, z) + + for _ in xrange(niter): + y, fc = gss(g, 0, 1, tol=1e-10) + z, fc = gss(h, 0, 1, tol=1e-10) + x, fc = gss(f, 0, 1e-3, tol=1e-10) + print x, y, z, fc sys.stdout.flush() @@ -100,7 +123,7 @@ def coarse_search(): v = map(float, line.strip().split()) d[tuple(v[:3])] = v[3] p = Pool(5) - lamb = [20. ** i for i in range(-15, 15)] + lamb = [20. ** i for i in range(-10, 0)] alpha = [20. ** i for i in range(-15, 15)] mu = [20. ** i for i in range(-15, 15)] @@ -112,32 +135,21 @@ def coarse_search(): sys.stdout.flush() p.map(aux, product(lamb, alpha, mu)) - # for l, a, m in product(lamb, alpha, mu): - # print l, a, m, ll(*(l, a, m)) if __name__ == "__main__": nodes, edges, events, event_edges = load(open("data2.pickle", "rb")) - def f(x): - global y, z - return ll(x, y, z) - - def g(y): - global x, z - return ll(x, y, z) + x = 1.2e-5 + y = 0.002 + z = 0.004 + # sa(x, y, z) - def h(z): - global x, y - return ll(x, y, z) + with open(sys.argv[1]) as fh: + l = [map(float, line.strip().split()[:3]) for line in fh] + for e in l: + optimize_with_sa(*e) - # x, y, z = map(float, sys.argv[1:]) - x = 1.8e-5 - y = 1. - z = 1. - #x, y, z = 1.8387870804e-05, 0.0383954874502, 0.00108402386463 - x, y, z = 1.83e-5, 0.00125, 0.00125 - #optimize_with_sa(x, y, z) - #optimize_with_gss(x, y, z) - print ll(x, y, z) - #coarse_search() + # optimize_with_gss(x, y, z) + # print ll(x, y, z) + # coarse_search() diff --git a/hawkes/refine.py b/hawkes/refine.py index 675bee3..3dfcfcc 100644 --- a/hawkes/refine.py +++ b/hawkes/refine.py @@ -1,6 +1,8 @@ from cPickle import load from itertools import product -from math import exp, log +from math import exp, log, sin +import sys +from multiprocessing import Pool def iter_events(events): @@ -13,31 +15,48 @@ def inprod(a, b): return tuple(float(x * y) for x, y in zip(a, b)) -def ll(a): - lamb, alpha, mu = a - - r1 = sum(log(lamb + sum(alpha * w * mu * exp(-mu * (t1 - t2)) - for (n2, t2, w) in s)) +def ll(lamb, alpha, mu): + r1 = sum(log(lamb * (1 + 0.43 * sin(0.0172 * t1 + 4.36)) + + sum(alpha / d * mu * exp(-mu * (t1 - t2)) + for (n2, t2, d) in s)) for ((n1, t1), s) in event_edges.iteritems()) - r2 = sum(sum(alpha * w * (1 - exp(-mu * (nodes[n2] - t1))) - for n2, w in edges[n1].iteritems() if nodes[n2] > t1) + r2 = sum(sum(alpha / d * (1 - exp(-mu * (nodes[n2] - t1))) + for n2, d in edges[n1].iteritems() if nodes[n2] > t1) for (n1, t1) in iter_events(events)) r3 = lamb * sum(nodes.itervalues()) return -(r1 - r2 - r3) -if __name__ == "__main__": - nodes, edges, events, event_edges = load(open("data.pickle", "rb")) +def get_values(): d = {} - for line in open("test.txt"): + for line in open(sys.argv[1]): v = map(float, line.strip().split()) d[tuple(v[:3])] = v[3] l = d.items() l.sort(key=lambda x: x[1]) - for a, _ in l[:10]: - t = [1. / i for i in range(2, 6)] + [float(i) for i in range(1, 6)] + for line in open("refine.txt"): + v = map(float, line.strip().split()) + d[tuple(v[:3])] = v[3] + for a, _ in l[:100]: + t = [1. / i for i in range(2, 4)] + [float(i) for i in range(1, 4)] for b in product(t, repeat=3): l, al, m = inprod(a, b) if (l, al, m) in d: continue - print l, al, m, ll((l, al, m)) + yield (l, al, m) + + +def refine(): + p = Pool(5) + + def aux(x): + l, a, m = x + print l, a, m, ll(l, a, m) + sys.stdout.flush() + + p.map(aux, get_values()) + + +if __name__ == "__main__": + nodes, edges, events, event_edges = load(open("data2.pickle", "rb")) + refine() |