function [f,b,a] = genalg(x, y)
%[f,b,a] = genalg(x, y)
%genetic algorithm for determining the filter coefficients in the filtering network between x and y
% Outputs:
% f: estimate for y
% b, a: filter coefficients (see filter)
%
%Copyright (C) 2010 Andras Hartmann <hdbandi@gmail.com>
%
%This file is part of GapReconst, a software for
%reconstructing short loss of physiological signals.
%
%GapReconst is free software: you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation, either version 3 of the License.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program. If not, see <http://www.gnu.org/licenses/>.
bdim = 4;
adim = 4;
zdim = size(x,1);
population = 100;
%Choosing the first population from 100*population to have an appropriate prior
bs = randn(100*population,bdim,zdim);
as = randn(100*population,adim,zdim);
%The actual iterations of the genetic algorithm
for iter = 1:2000;
c = [];
st1 = size(as,1);
for i = 1:st1
f = multifilter(squeeze(bs(i,:,:)),squeeze(as(i,:,:)),x);
ccf = corrcoef(y,f);
c(i) = mse(y-f)*(1-ccf(2,1));
end
% Setting up the next generation
newas = [];
newbs = [];
%population/5: survival of the fittest
pop5 = population/5;
for i = 1:pop5
[value,index] = min(c);
% uncomment this if you would like to have some info printed about the generations
% if i == 1
% value
% index
% end
c(index) = [];
if ~isnan(value)
newas = [newas; as(index,:,:)];
newbs = [newbs; bs(index,:,:)];
end
end
%The rest of the population is set up as recombination of two random specimen
st1 = size(newas,1);
for i = st1+1:population
mo = unidrnd(st1);%mother
fa = unidrnd(st1);%father
newas(i,:,:) = mean([newas(mo,:,:); newas(fa,:,:)]);
newbs(i,:,:) = mean([newbs(mo,:,:); newbs(fa,:,:)]);
end
%Adding random (Gauss) mutation in three different way
%some mutation to all all coeffitients
newas(pop5+1:2*pop5,:,:) = newas(pop5+1:2*pop5,:,:)+randn(pop5,adim,zdim)*0.001;
newbs(pop5+1:2*pop5,:,:) = newbs(pop5+1:2*pop5,:,:)+randn(pop5,bdim,zdim)*0.001;
%some more mutation to all coeffitients, this time with larger amplitude
newas(2*pop5+1:3*pop5,:,:) = newas(2*pop5+1:3*pop5,:,:)+randn(pop5,adim,zdim)*0.01;
newbs(2*pop5+1:3*pop5,:,:) = newbs(2*pop5+1:3*pop5,:,:)+randn(pop5,bdim,zdim)*0.01;
%some more mutation, this time with even larger amplitude, but adding only with 0.1 probability to the coeffitients
R1 = rand(pop5,adim,zdim);
R = zeros(pop5,adim,zdim);
R (find(R1>0.9)) = 1;
newas(3*pop5+1:4*pop5,:,:) = newas(3*pop5+1:4*pop5,:,:)+R.*randn(pop5,adim,zdim)*0.1;
R1 = rand(pop5,bdim,zdim);
R = zeros(pop5,bdim,zdim);
R (find(R1>0.9)) = 1;
newbs(3*pop5+1:4*pop5,:,:) = newbs(3*pop5+1:4*pop5,:,:)+R.*randn(pop5,bdim,zdim)*0.1;
%Replacing the actual generation with the new one
as = newas;
bs = newbs;
%This would be the condition to stop, but not used
%setting limit
%if (c(1) < 0.5)
% break;
%end
end;
% the best filters to return
a = squeeze(as(1,:,:));
b = squeeze(bs(1,:,:));