遺伝的アルゴリズムの簡単なテスト

以前から遺伝的アルゴリズムというものに興味を持っていたのですが、本を読むだけじゃなく実際に手を動かさなければ、ということでプログラムを書いてみました。勉強中のLISPで書きたいところでしたが、いちばん慣れているMatlabを使ってみました。

y=e^k\sin\left( 2\pi x\right)

の最大値を出すx∈[0,1)の範囲で求めるものです。下のコードは未完成なので、優良な遺伝子を選び出す部分がぜんぜんダメですが、とりあえず何回かループさせるとx=0.2695, y=1.2995という値が得られます。こんな適当なプログラムなのに、ある程度の結果が出ているのが感動的ですらあります。

% test of genetic algorithm
function test_ga
numGene = 10;
numBits = 8;

% rows: number of individuals
% columns: length of a gene
genes = randn(numGene, numBits) >= 0;

clc;

disp([decode(genes) genes]);
fprintf(1,'\n');

for n=1:5
  genes = selection(genes, 0.2);
  genes = crossover(genes, 0.8);
  genes = mutate(genes, 0.8);  
end

fprintf(1,'\n');
disp([decode(genes) genes]);

end



%% CROSSOVER within 'gene' farm with probability 'rate'
function newgenes = crossover(genes, rate)

len = size(genes,1);
A = (rand(len,1) <= rate) .* randperm(len)';
A = A(A~=0);

if length(A)>1
  for n=1:floor(length(A)/2)
    cutpoint = floor(rand() * (length(genes(A(n),:))-1))+1;
    newgene1 = [genes(A(n),   1:cutpoint) genes(A(n+1), cutpoint+1:end)];
    newgene2 = [genes(A(n+1), 1:cutpoint) genes(A(n),   cutpoint+1:end)];
    genes(A(n),  :) = newgene1;
    genes(A(n+1),:) = newgene2;
  end
end
newgenes = genes;
end



%% MUTATE
function newgenes = mutate(genes, rate)
for n=1:size(genes,1)
  if rand <= rate
    m = floor(rand() * length(genes(n,:))) + 1;
    if genes(n,m) == 1
      genes(n,m) = 0;
    else
      genes(n,m) = 1;
    end
  end
end
newgenes = genes;
end



%% DECODE binary genes into decimal number
function dec = decode(genes)
dec = zeros(size(genes,1),1);

for n=1:size(genes, 1)
  A = '';
  for m=1:size(genes, 2)
    A = [A num2str(genes(n,m))];
  end
  dec(n) = bin2dec(A);
end
end



%% compute score for arbitrary function (this is the objective function!)
function score = objective_function(genes)
k = decode(genes) / 256;
score = sin(k*2*pi) .* exp(k);
disp([score k])
end



%% SELECTION
function newgenes = selection(genes, rate)
score = objective_function(genes);
[tmp, index] = sort(score, 1, 'descend');
disp(tmp(1));
for n=1:floor(size(genes,1) * rate)
  newgenes(n,:) = genes(index(n),:);
end
for n=floor(size(genes,1) * rate)+1:size(genes,1)
  newgenes(n,:) = genes(index(n-floor(size(genes,1)*rate)),:);
end
end