function [bound, bestpop, elapsedtime, POPULATION] = jsr_genetic(M, params)

% JSR_GENETIC Computes a lower bound on the jsr using genetic algorithm.
%
%    [bound, bestpop, elapsedtime, population] = jsr_genetic(M, params)
%
%    Inputs:
%      M      - (required) a non-empty cell array containing square matrices of
%               identical dimensions.
%      params - (optional) a structure containing custom parameters (see below).
%
%    Outputs:
%      bound       - a lower bound on the jsr of M.
%      bestpop     - a product of matrices corresponding to this lower bound:
%                    the product A1*A2*A3 is expressed as [1 2 3] in this order.
%      elapsedtime - cpu time used to run the algorithm.
%      population  - the complete population at the end of the algorithm.
%
%    Parameters:
%      verbose     - (default 1) defines the printing level:
%                    <  1 - no output is printed on the screen;
%                    >= 1 - prints the final lower bound and the associated
%                           product at the end of the algorithm;
%                    >= 2 - prints the current lower bound and the associated
%                           product at each generation of the algorithm;
%                    >= 3 - prints the current maximum product length and the
%                           number of stalling iterations at this current
%                           maximum product length at each generation.
%      popsize     - (default 150, minimum 10) population size.
%      maxgen      - (default 1000) maximum number of generations.
%      maxstall    - (default 15) maximum number of stalling iterations before
%                    increasing the maximum product length.
%      maxtotstall - (default 100) maximum number of stalling iterations before
%                    terminating the algorithm.
%      mutantprop  - (default 0.3) mutation probability of a given product.
%      muteprop    - (default 0.2) mutation proportion of a given product.

%% Default parameters
error(nargchk(1, 2, nargin, 'struct'));

verbose = 1;
POPSIZE = 150;
MAXGEN = 1000;
MAXSTALLING = 15;
MAXTOTALSTALLING = 100;
MUTANTPROP = 0.3;
MUTEPROP = 0.2;

if (nargin == 2) && isstruct(params),
    if isfield(params, 'verbose'), verbose = params.verbose; end;
    if isfield(params, 'popsize'), POPSIZE = params.popsize; end;
    if isfield(params, 'maxgen'), MAXGEN = params.maxgen; end;
    if isfield(params, 'maxstall'), MAXSTALLING = params.maxstall; end;
    if isfield(params, 'maxtotstall'), MAXTOTALSTALLING = params.maxtotstall; end;
    if isfield(params, 'mutantprop'), MUTANTPROP = params.mutantprop; end;
    if isfield(params, 'muteprop'), MUTEPROP = params.muteprop; end;
end

POPSIZE = max(10, POPSIZE);

%% Initialization
starttime = cputime;
m = numel(M);
n = size(M{1}, 1);
k = ceil(log(POPSIZE)/log(m+1));
scaler = ((m+1).^(k-1:-1:0))';

% Cache
cache = genCache(M, k, m, n);
ncache = size(cache, 2);
stalling = 0;
totalstalling = 0;
bound = 0;
bestpop = [];
for i = 2:ncache,
    value = max(abs(eig(mat(cache(:, i)))));
    if (value > bound),
        bound = value;
        bestpop = zeros(1, k);
        key = i-1;
        for j = 1:k,
            bestpop(j) = floor(key/scaler(j));
            key = key - bestpop(j)*scaler(j);
        end
    end
end
bestpop = bestpop(bestpop ~= 0);
bound = bound^(1/k);
X = eye(n);
testbestpop = bestpop;
for i = 1:length(testbestpop)-1,
    X = X * M{testbestpop(i)};
    testval = max(abs(eig(X)))^(1/i);
    if (testval >= bound),
        bestpop = testbestpop(1:i);
        bound = testval;
    end
end
if (verbose >= 2),
    fprintf(' \n');
    fprintf('Starting population: init lower bound on the JSR = %.15g with product: %s\n', bound, num2str(deperiod(bestpop)));
end

% Initial population
CURLENGTH = 2*k;
POPULATION = floor((m+1)*rand(POPSIZE, CURLENGTH));

%% Genetic evolution
for gen = 1:MAXGEN,
    % Evaluation
    [nbeyes, idx] = sort(sum(POPULATION == 0, 2), 'descend');
    POPULATION = POPULATION(idx, :);
    POPULATION((nbeyes == CURLENGTH), 1) = ceil(m*rand(1));
    nbeyes = min(nbeyes, CURLENGTH-1);
    fitness = zeros(POPSIZE, 1);
    cached = floor(CURLENGTH/k);
    cachekey = zeros(POPSIZE, cached);
    for j = 1:cached,
        cachekey(:, j) = POPULATION(:, k*j-k+1:k*j)*scaler;
    end
    cachekey = cachekey + 1;
    for i = 1:POPSIZE,
        X = eye(n);
        for j = 1:cached,
            X = X * mat(cache(:, cachekey(i, j)));
        end
        for j = k*cached+1:CURLENGTH,
            if (POPULATION(i, j) ~= 0),
                X = X * M{POPULATION(i, j)};
            end
        end
        fitness(i) = max(abs(eig(X)))^(1/(CURLENGTH-nbeyes(i)));
    end
    [fitness, idx] = sort(fitness, 'descend');
    POPULATION = POPULATION(idx, :);
    
    % Local optimization
    if (fitness(1) > bound),
        stalling = 0;
        totalstalling = 0;
        bound = fitness(1);
        bestpop = POPULATION(1, :);
        bestpop = bestpop(bestpop ~= 0);
        testbestpop = bestpop;
        lenbestpop = length(testbestpop);
        XX = cell(lenbestpop + 1, 1);
        X = eye(n);
        XX{1} = X;
        for i = 1:lenbestpop,
            X = X * M{testbestpop(i)};
            XX{i+1} = X;
            testval = max(abs(eig(X)))^(1/i);
            if (testval > bound) || ((testval == bound) && (i < length(bestpop))),
                bestpop = testbestpop(1:i);
                bound = testval;
            end
        end
        lenbestpop = length(testbestpop);
        X = eye(n);
        localimprove = 0;
        for i = lenbestpop:-1:1,
            testval = max(abs(eig(XX{i}*X)))^(1/(lenbestpop-1));
            if (testval > bound),
                bestpop = testbestpop([1:i-1 i+1:end]);
                bound = testval;
                localimprove = localimprove + 1;
            end
            for j = 1:m,
                testval = max(abs(eig(XX{i}*M{j}*X)))^(1/lenbestpop);
                if (testval > bound),
                    bestpop = testbestpop;
                    bestpop(i) = j;
                    bound = testval;
                    localimprove = localimprove + 1;
                end
            end
            if (lenbestpop < CURLENGTH),
                for j = 1:m,
                    testval = max(abs(eig(XX{i+1}*M{j}*X)))^(1/(lenbestpop+1));
                    if (testval > bound),
                        bestpop = [testbestpop(1:i)  j  testbestpop(i+1:end)];
                        bound = testval;
                        localimprove = localimprove + 1;
                    end
                end
            end
            X = M{testbestpop(i)}*X;
        end
        if localimprove,
            POPULATION = [bestpop, zeros(1, CURLENGTH-length(bestpop)); POPULATION(1:end-1, :)];
        end
    else
        stalling = stalling + 1;
        totalstalling = totalstalling + 1;
    end
    
    % Stopping criterion
    if (verbose >= 3),
        fprintf('Generation #%3d: STA = %2d, LEN = %2d, lower bound = %.15g with product: %s\n', gen, stalling, CURLENGTH, bound, num2str(deperiod(bestpop)));
    elseif (verbose >= 2),
        fprintf('Generation #%3d:  current lower bound on the JSR = %.15g with product: %s\n', gen, bound, num2str(deperiod(bestpop)));
    end
    if (totalstalling >= MAXTOTALSTALLING),
        break;
    end
    if (stalling >=  MAXSTALLING),
        stalling = 0;
        CURLENGTH = CURLENGTH + 1;
        POPULATION = [POPULATION, zeros(POPSIZE, 1)]; %#ok<AGROW>
    end
    
    % Selection and crossover
    NB_ELITE = max(2, min(3, floor(POPSIZE/50)));
    NB_SPAWN = max(4, floor(POPSIZE/50));
    NB_SWAP = min(POPSIZE - NB_ELITE - NB_SPAWN, floor(POPSIZE/2));
    NB_MIX = POPSIZE - NB_ELITE - NB_SPAWN - NB_SWAP;
    ID_SWAP = ceil((POPSIZE/2)*rand(NB_SWAP, 2));
    ID_MIX = ceil(POPSIZE*rand(NB_MIX, 2));
    
    POP_ELITE = POPULATION(1:NB_ELITE, :);
    POP_SPAWN = ceil(m*rand(NB_SPAWN, CURLENGTH));
    POP_SWAP = zeros(NB_SWAP, CURLENGTH);
    for i = 1:NB_SWAP,
        cut = ceil(CURLENGTH*rand(1));
        POP_SWAP(i, :) = [POPULATION(ID_SWAP(i, 1), 1:cut) POPULATION(ID_SWAP(i, 2), cut+1:end)];
    end
    POP_MIX = zeros(NB_MIX, CURLENGTH);
    for i = 1:NB_MIX,
        mixer = floor(2*rand(CURLENGTH, 1))';
        POP_MIX(i, :) = POPULATION(ID_MIX(i, 1), :) .* mixer + POPULATION(ID_MIX(i, 2), :) .* (1-mixer);
    end
    
    POPULATION = [POP_ELITE ; POP_SPAWN ; POP_SWAP ; POP_MIX];
    
    % Mutation
    MUTESTR = ceil(CURLENGTH * MUTEPROP);
    for i = 2:POPSIZE,
        if (rand(1) < MUTANTPROP),
            POPULATION(i, ceil(CURLENGTH * rand(1, MUTESTR))) = floor((m+1)*rand(1, MUTESTR));
        end
    end
end

%% Termination
elapsedtime = cputime - starttime;
bestpop = deperiod(bestpop);
if (verbose >= 1),
    fprintf('Algorithm terminated with lower bound on the JSR = %.15g with product: %s\n', bound, num2str(bestpop));
end


%% %%%%%%%%%%%%%%%% %%
%  CACHE GENERATION  %
%%%%%%%%%%%%%%%%%%%%%%

function cache = genCache(M, k, m, n)

if (k <= 1),
    cache = zeros(n*n, m+1);
    cache(:, 1) = vec(eye(n));
    for i = 1:m,
        cache(:, i+1) = vec(M{i});
    end
    return;
end

cache = repmat(genCache(M, k-1, m, n), 1, m+1);
gsize = (m+1)^(k-1);
for i = 1:m,
    cache(:, gsize*i+1:gsize*(i+1)) = kron(eye(n), M{i}) * cache(:, gsize*i+1:gsize*(i+1));
end


%% %%%%%%%%%%%%%%% %%
%  DEPERIODIZATION  %
%%%%%%%%%%%%%%%%%%%%%

function [x, k] = deperiod(X)

n = length(X);

for k = 1:n/2,
    if (mod(n, k) == 0),
        ok = 1;
        for t = 1:n-k,
            if X(t) ~= X(t+k),
                ok = 0;
                break;
            end
        end
        if ok,
            x = X(1:k);
            return;
        end
    end
end

x = X;
k = n;