Final tweaks to the optimization

1 files changed, 54 insertions, 42 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()