cuckoo search, latin hypercube sampling

diff --git a/cs.py b/cs.py
new file mode 100644 (file)
index 0000000..7c91119
--- /dev/null
+++ b/cs.py
@@ -0,0 +1,71 @@
+"""Cuckoo search optimization"""
+
+import random
+import lhs
+from math import sin, gamma, pi
+
+class Space(object):
+    def __init__(self, dimension, bound, d_function, n_nest):
+        self.pa = 0.25 #parameter
+        self.dimension = dimension
+        self.bound = bound
+        self.d_function = d_function
+        self.nests = [(d_function(*p), p) for p in lhs.latin_hypercube(dimension, bound, n_nest)]
+        self.best_value, self.best = max(self.nests)
+
+def new_nest(space):
+        position = [2 * space.bound * random.random() 
+                    - space.bound for _ in xrange(space.dimension)]
+        value = space.d_function(*position)
+        return (value, position)
+
+def get_cuckoos(space):
+    beta = 1.5
+    sigma = (gamma(1. + beta) * sin(pi * beta / 2.) / (gamma((1. + beta) / 2.) *
+             beta * 2. ** ((beta - 1.) / 2))) ** (1. / beta)
+    u_a = [[random.gauss(0, 1) * sigma for _ in xrange(space.dimension)] for _ in
+         xrange(len(space.nests))]
+    v_a = [[random.gauss(0, 1) for _ in xrange(space.dimension)] for _ in
+         xrange(len(space.nests))]
+    r_a = [[random.gauss(0, 1) for _ in xrange(space.dimension)] for _ in
+         xrange(len(space.nests))]
+    step = [[u / abs(v) ** (1. / beta) for (u, v) in zip(u_l, v_l)]
+            for (u_l, v_l) in zip(u_a, v_a)]
+    stepsize = [[0.01 * st * (n_e - be) for (st, n_e, be) 
+                in zip(step_l, n_l, space.best)]
+                for (step_l, (_, n_l)) in zip(step, space.nests)]
+    s = [[n + st * r for (n, st, r) in zip(n_l, st_l, r_l)]
+         for ((_, n_l), st_l, r_l) in zip(space.nests, stepsize, r_a)]
+    cuckoos = [[min(max(st, - space.bound), space.bound) for st in st_l]
+         for st_l in s]
+    return [(space.d_function(*c), c) for c in cuckoos]
+
+def get_empty(space):
+    r = random.random()
+    r_arr = [[random.random() for _ in xrange(space.dimension)] for _ in
+             xrange(len(space.nests))]
+    perm1 = [n for (_, n) in space.nests]
+    random.shuffle(perm1)
+    perm2 = [n for (_, n) in space.nests]
+    random.shuffle(perm2)
+    stepsize = [[p1 - p2 for (p1, p2) in zip (p1l, p2l)] for (p1l, p2l) in
+                zip(perm1, perm2)]
+    step = [[(r * p * (1 if random.random() > space.pa else 0)) for p in n] for n in stepsize]
+    empty = [[(p + s) for (p, s) in zip(sl, n)] 
+            for (sl, (_, n)) in zip(step, space.nests)]
+    empty = [[min(max(st, - space.bound), space.bound) for st in st_l]
+         for st_l in empty]
+    return [(space.d_function(*e), e) for e in empty]
+
+def next_turn(space):
+    cuckoos = get_cuckoos(space)
+    space.nests = [max(n, m) for (n, m) in zip(space.nests, cuckoos)]
+    nests = get_empty(space)
+    space.nests = [max(n, m) for (n, m) in zip(space.nests, nests)]
+    space.best_value, space.best = max(space.nests)
+
+def optimize(dimension, boundary, function_d, n_nest, n_turns):
+    space = Space(dimension, boundary, function_d, n_nest)
+    for _ in xrange(n_turns):
+        next_turn(space)
+    return space.best

diff --git a/gridf.py b/gridf.py
index e957880..03e3014 100644 (file)
--- a/gridf.py
+++ b/gridf.py
@@ -10,7 +10,7 @@ from manual import lines as g_grid
 from intrsc import intersections_from_angl_dist
 from linef import line_from_angl_dist
 import pcf
 from intrsc import intersections_from_angl_dist
 from linef import line_from_angl_dist
 import pcf

-import pso
+import cs as Optimizer

 
 class GridFittingFailedError(Exception):
     pass
 
 class GridFittingFailedError(Exception):
     pass


@@ -52,7 +52,7 @@ def find(lines, size, l1, l2, bounds, hough, do_something, im_h):
     for line in sum(lines, []):
         dr_l.line(line_from_angl_dist(line, size), width=1, fill=255)
 
     for line in sum(lines, []):
         dr_l.line(line_from_angl_dist(line, size), width=1, fill=255)
 

-    im_l = im_l.filter(MyGaussianBlur(radius=5))
+    im_l = im_l.filter(MyGaussianBlur(radius=3))

     #GaussianBlur is undocumented class, may not work in future versions of PIL
     im_l_s = im_l.tostring()
 
     #GaussianBlur is undocumented class, may not work in future versions of PIL
     im_l_s = im_l.tostring()
 


@@ -62,7 +62,7 @@ def find(lines, size, l1, l2, bounds, hough, do_something, im_h):
     f_dist = partial(job_4, im_l=im_l_s, v1=v1, v2=v2, h1=h1, h2=h2,
                      dv=delta_v, dh=delta_h, size=size)
 
     f_dist = partial(job_4, im_l=im_l_s, v1=v1, v2=v2, h1=h1, h2=h2,
                      dv=delta_v, dh=delta_h, size=size)
 

-    x_v, y_v, x_h, y_h = pso.optimize(4, 30, f_dist, 32, 1028)
+    x_v, y_v, x_h, y_h = Optimizer.optimize(4, 30, f_dist, 128, 256)

 
     v1 = (v1[0] + x_v * delta_v, v1[1] + x_v)
     v2 = (v2[0] + y_v * delta_v, v2[1] + y_v)
 
     v1 = (v1[0] + x_v * delta_v, v1[1] + x_v)
     v2 = (v2[0] + y_v * delta_v, v2[1] + y_v)

diff --git a/lhs.py b/lhs.py
new file mode 100644 (file)
index 0000000..c890dfe
--- /dev/null
+++ b/lhs.py
@@ -0,0 +1,25 @@
+import random
+
+def test():
+    bound = 10.
+    m = 25
+    l = latin_hypercube(2, bound, m)
+    import matplotlib.pyplot as pyplot
+    fig = pyplot.figure()
+    fig.add_subplot(121)
+    pyplot.plot([v[0] for v in l], [v[1] for v in l], 'o')
+    fig.add_subplot(122)
+    pyplot.plot([random.random() * 2 * bound - bound for _  in xrange(m)], 
+                [random.random() * 2 * bound - bound for _  in xrange(m)], 'o')
+    pyplot.show()
+    import sys
+    sys.exit()
+
+def latin_hypercube(dim, bound, m):
+    dv = (2 * bound) / float(m)
+    dim_p = [range(m) for _ in xrange(dim)] 
+    for p in dim_p:
+        random.shuffle(p)
+    points = [list(l) for l in zip(*dim_p)]
+    points = [[(float(l) + random.random()) * dv - bound  for l in p] for p in points]
+    return points

diff --git a/pso.py b/pso.py
index c30bcfd..6d6fada 100644 (file)
--- a/pso.py
+++ b/pso.py
@@ -4,8 +4,13 @@ import random
 import multiprocessing
 from functools import partial
 
 import multiprocessing
 from functools import partial
 

-def particle(dimension, bound, v_max, func_d):
-    position = [2 * bound * random.random() - bound for _ in xrange(dimension)]
+import lhs
+
+def particle(dimension, bound, v_max, func_d, pos=None):
+    if not pos:
+        position = [2 * bound * random.random() - bound for _ in xrange(dimension)]
+    else:
+        position = pos

     velocity = [2 * v_max * random.random() - v_max for _ in xrange(dimension)]
     value = func_d(*position)
     return value, position, velocity, value, position
     velocity = [2 * v_max * random.random() - v_max for _ in xrange(dimension)]
     value = func_d(*position)
     return value, position, velocity, value, position


@@ -27,8 +32,8 @@ def move(particle, omega, phi_p, phi_g, v_max, global_best, func_d):
 def optimize(dimension, boundary, function_d, n_parts, n_turns):
     pool = multiprocessing.Pool(None)
     v_max = boundary
 def optimize(dimension, boundary, function_d, n_parts, n_turns):
     pool = multiprocessing.Pool(None)
     v_max = boundary

-    particles = [particle(dimension, boundary, v_max, function_d)
-                 for _ in xrange(n_parts)]
+    particles = [particle(dimension, boundary, v_max, function_d, pos)
+                 for pos in lhs.latin_hypercube(dimension, bound, n_parts)]

     gl_best = max(particles)
     for _ in xrange(n_turns):
         move_p = partial(move, 
     gl_best = max(particles)
     for _ in xrange(n_turns):
         move_p = partial(move,