function []=Simulation_Mixed_Games(W,r,N,MaxLoops)
%
% Mixed_Games(W, range, N, MaxLoops)
%    where W is a 2x2 matrix of the payoff to the two pure strategies.
%    range is the range of strategy values to consider
%    N is the population size
%    MaxLoops is the number of loops
%

if nargin<4
    MaxLoops=1000;
    if nargin<3
        N=1000;
        if nargin<2
            r=[0 1];
            if nargin<1
                W=[1 0; 2 -1];  % hawk-dove model
            end
        end
    end
end
r=[min(r) max(r)];
N=2*floor(N/2); % make sure N is even.

P=r(1) + (r(2)-r(1))*rand(N,1); % pick random strategies.

for L=1:MaxLoops
    
    %calculate fitnesses against the average population
    %
    for i=N:-1:1
        Q=[P(i); 1-P(i)]; Pave=[mean(P); 1-mean(P)];
        F(i)=(Q'*W)*Pave;
    end
    
%     %calculate fitnesses against a random pick
%     
%     for i=N:-1:1
%         Q=[P(i); 1-P(i)]; R=ceil(N*rand(1,1)); Prnd=[P(R); 1-P(R)];
%         F(i)=(Q'*W)*Prnd;
%     end
    
    
    %largest half breed
    [y i]=sort(F,'descend'); i=i(1:(N/2));
    P=[P(i); P(i)];
    % include mutation
    P=P+0.01*randn(N,1);
    % reflect from boundaries
    P(P<0)=-P(P<0);
    P(P>1)=2-P(P>1);
    
    hist(P,r(1)+(r(2)-r(1))*[0.005:0.01:1]); drawnow;
end
end