repmut.m

%function repmut
% 12/08/14

function repmut

clear
format compact


f = 1;  % this is a dummy init to make f global
N = 64;     % 64
p = 1/sqrt(N);
time_expand = N;

mutype = 1;     % 0 = Hamming   1 = rand
pay = 0;        % 0 = Hamming   1 = 1/sqrt(N) 
ep = .06;      % average mutation rate: 0.1 to 0.01 typical  (0.4835)

%%%%% Set original population
x0temp = rand(1,N);    % Initial population
sx = sum(x0temp);
y0 = x0temp/sx;
Pop0 = sum(y0);


%%%%% Set Adjacency

%node = makeglobal(N);
node = makeER(N,0.1);
%node = makeSF(N,4);
%node = makeSW(N,4,0.5);
[Adj,degree,Lap] = adjacency(node);



%%%%%% Set Hamming distance
for yloop = 1:N
    for xloop = 1:N
        H(yloop,xloop) = hamming(yloop-1,xloop-1);
    end
end



%%%%%%% Set Mutation matrix
if mutype == 0
    Qtemp = 1./(1+H/ep);    %Mutation matrix on Hamming
    %Qtemp = exp(-H/(ep*50));
    Qsum = sum(Qtemp,2);
    
    % Normalize mutation among species
    for yloop = 1:N
        for xloop = 1:N
            Q(yloop,xloop) = Qtemp(yloop,xloop)/Qsum(xloop);
        end
    end
    
elseif mutype == 1  
    S = stochasticmatrix(N);
    Stemp = S - diag(diag(S));
    Qtemp = ep*Stemp;
    sm = sum(Qtemp,2)';
    Q = Qtemp + diag(ones(1,N) - sm);
end


figure(1)
imagesc(Q)
title('Mutation Matrix')
colormap(jet)

%%%%%%% Set payoff matrix
if pay == 1
    payoff = zeros(N,N);
    for yloop = 1:N
        payoff(yloop,yloop) = 1;
        for xloop = yloop + 1:N
            payoff(yloop,xloop) = p;
            payoff(xloop,yloop) = p;
            %payoff(yloop,xloop) = p*2*(0.5 - randbin(1,0.5));
            %payoff(xloop,yloop) = payoff(yloop,xloop);
            %payoff(xloop,yloop) = -payoff(yloop,xloop);
        end
    end
elseif pay == 0
    payoff = exp(-1*H);
end


figure(2)
imagesc(payoff)
title('Payoff Matrix')
colormap(jet)


% Run time evolution
tspan = [0 1000];
[t,x] = ode45(@quasispec,tspan,y0);

Pop0
[sz,dum] = size(t);
Popend = sum(x(sz,:))

phistar = sum(f.*x(sz,:))       % final average fitness



figure(3)
clf
h = colormap(lines);
for loop = 1:N
    plot(t,x(:,loop),'Color',h(round(loop*64/N),:))
    hold on
end
hold off

figure(4)
clf
for loop = 1:N
    semilogx(t,x(:,loop),'Color',h(round(loop*64/N),:))
    hold on
end
hold off

figure(5)
clf
for loop = 1:N
    semilogy(t,x(:,loop),'Color',h(round(loop*64/N),:))
    hold on
end
hold off

%keyboard

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function yd = quasispec(~,y)
        
        for loop = 1:N
            ftemp(loop) = sum(payoff(:,loop).*y);
        end
        f = time_expand*ftemp/sum(ftemp);
        
        
        % Transition matrix
        for yloop = 1:N
            for xloop = 1:N
                W(yloop,xloop) = f(yloop)*(Adj(yloop,xloop)*Q(yloop,xloop));
            end
        end
        
        phi = sum(f'.*y);   % Average fitness of population
        
        yd = W*y - phi*y;
        
        
    end     % end quasispec

end