From ab0b1f3cefedb35327a19ec1b6afd560bfdf802d Mon Sep 17 00:00:00 2001 From: Thibaut Horel Date: Mon, 14 Sep 2015 23:08:02 -0400 Subject: Import supplements and repo reorganization --- hawkes_experiments/main.py | 151 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 151 insertions(+) create mode 100644 hawkes_experiments/main.py (limited to 'hawkes_experiments/main.py') diff --git a/hawkes_experiments/main.py b/hawkes_experiments/main.py new file mode 100644 index 0000000..8f30dcf --- /dev/null +++ b/hawkes_experiments/main.py @@ -0,0 +1,151 @@ +from cPickle import load +from math import log, exp, sqrt, sin +from random import random, gauss +import sys +from itertools import product +from multiprocessing import Pool + +GR = (sqrt(5) - 1) / 2 + + +def gss(f, a, b, tol=1e-20): + """golden section search to minimize along one coordinate""" + c = b - GR * (b - a) + d = a + GR * (b - a) + while abs(c - d) > tol: + fc = f(c) + fd = f(d) + if fc < fd: + b = d + d = c + c = b - GR * (b - a) + else: + a = c + c = d + d = a + GR * (b - a) + sys.stderr.write("gss:" + " ".join(map(str, [a, b, fc, fd])) + "\n") + sys.stderr.flush() + return (b + a) / 2, fc + + +def iter_events(events): + for n, s in events.iteritems(): + for t in s: + yield (n, t) + + +def approx(x): + if x > 1e-10: + return 1 - exp(-x) + else: + return x + + +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 / d * approx(mu * (nodes[n2][0] - t1)) + for n2, d in edges[n1].iteritems() + if nodes[n2][0] > t1) + for (n1, t1) in iter_events(events)) + r3 = lamb * sum(node[1] for node in nodes.itervalues()) + print r1, r2, r3 + return -(r1 - r2 - r3) + + +def sa(x, y, z, sigma=0.5, niter=1000, fc=None): + T = 0.1 + e = 1.1 + if fc: + fo = fc + else: + fo = ll(x, y, z) + 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), 0) + zn = max(z + gauss(0, sigma * z), 0) + fn = ll(x, yn, zn) + if fn < fo or exp((fo - fn) / T) > random(): + y = yn + z = zn + sigma *= 2 + fo = fn + else: + sigma /= e + if sigma < 1e-5: + break + T *= 0.99 + return y, z, fo + + +def optimize_with_sa(x, y, z, niter=10): + + def f(x): + return ll(x, y, z) + + # y, z = sa(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, niter=5): + + 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, 1e50, tol=1e-10) + z, fc = gss(h, 0, 1e50, tol=1e-10) + x, fc = gss(f, 0, 1e-3, tol=1e-10) + print x, y, z, fc + sys.stdout.flush() + + +def coarse_search(): + d = {} + for line in open("values.txt"): + v = map(float, line.strip().split()) + d[tuple(v[:3])] = v[3] + p = Pool(5) + 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)] + + def aux(x): + if x in d: + return + l, a, m = x + print l, a, m, ll(l, a, m) + sys.stdout.flush() + + p.map(aux, product(lamb, alpha, mu)) + + +if __name__ == "__main__": + nodes, edges, events, event_edges = load(open("data-dist1.pickle", "rb")) + x = 1.25e-5 + y = 1.2e16 + z = 1.5e-20 + sa(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_gss(*e) + + # print ll(x, y, z) + # coarse_search() -- cgit v1.2.3-70-g09d2