From: Tomas Musil Date: Tue, 25 Dec 2012 23:14:39 +0000 (+0100) Subject: cuckoo search, latin hypercube sampling X-Git-Url: http://git.tomasm.cz/imago.git/commitdiff_plain/fb4d31c62744561af4bc3941d1cb9da430c5c64a?ds=sidebyside cuckoo search, latin hypercube sampling --- diff --git a/cs.py b/cs.py new file mode 100644 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 --- 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 -import pso +import cs as Optimizer 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) - 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() @@ -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) - 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) diff --git a/lhs.py b/lhs.py new file mode 100644 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 --- a/pso.py +++ b/pso.py @@ -4,8 +4,13 @@ import random 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 @@ -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 - 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,