
close all
clear all

if 1
    rng('shuffle');
    seed=rng;
    save seed;
else
    load seed;
    rng(seed)
end
ones(10)*ones(10);

% Parameters
N=1000;             % Population size  
n=20;               % group size
groups=N/n;
G=1;                % authority optimum
FF=(0:0.25:0.75);   % authority effort
s=0.15;             % probability of success in hunting
m=0.1;              % migration rate

T=1000;
Runs=10;


% Initial conditions: lognormal distribution of y, specific value of p 
Mu=0.2;
Sigma=0.1;
mu=log(Mu^2/sqrt(Sigma^2+Mu^2));        % for lognormal distribution of y
sigma=sqrt(log(Sigma^2/Mu^2+1));        % for lognormal distribution of y

% Psychological parameters

uy=0.5;         % beliefs update rate
lambda=10;      % precision in x
e=0.1;          % error in y
b=1.;

% edges=0:0.02:1.;
P=zeros(Runs,T);
Y=zeros(Runs,T);
X=zeros(N,T);
str_color=["darkgreen","orange","blue","red","magenta","cyan"];

figure
tic
for i_f=1:length(FF)
    F=FF(i_f);
    
    for run=1:Runs
        
        % Initial conditions
        p=zeros(N,1);               % no sharing initially
        y=lognrnd(mu,sigma,N,1);
        y0=y;                       % save initial distribution of y
        
        % individual parameters: broken stick distibution
        U=rand(N,4);
        v=U(:,1);
        k=(1-v).*U(:,2);
        d=1-v-k;
        alpha=U(:,3);
        beta=(1-alpha).*U(:,4);
        gama=1-alpha-beta;
        
        % Dynamics
        for t=1:T
            delta_U=y-b+v.*(2*y-1)+k.*(2*p-1)+d*F*(2*G-1);      % difference in utilities
            x=double(1./(1+exp(-lambda*delta_U))>rand(N,1));    % decide on action QRE
            X(:,t)=x;
            
            % same p for all groupmates
            h=double(rand(N,1)<s);                         % hunter: yes (1) or no (0)
            H=sum(reshape(h,n,groups),1);                  % number of hunters per group
            h1=reshape(x.*h,n,groups);                     % hunters who share
            S=sum(h1,1);                                   % their number per group
            p=S./(H+realmin);                              % actual frequency of sharing per group
            p=reshape(repmat(p,n,1),N,1);                  % vector form      

            % update y
            y_new=y+h.*alpha.*(x-y)+beta.*(p-y)+gama.*(G-y)*F+e*randn(N,1); % if update y            
            update_y=rand(N,1)<uy;              % random updating of y
            y(update_y)=y_new(update_y);
            y(y<0)=0;                  
            y(y>1)=1; 


            % migration
            emmigrants=find(rand(N,1)<m);
            immigrants=emmigrants(randperm(numel(emmigrants)));
            
            p(immigrants)=p(emmigrants);        % personal knowledge about p
            y(immigrants)=y(emmigrants);
            
            v(immigrants)=v(emmigrants);
            k(immigrants)=k(emmigrants);
            d(immigrants)=d(emmigrants);       
            alpha(immigrants)=alpha(emmigrants);
            beta(immigrants)=beta(emmigrants); 
            gama(immigrants)=gama(emmigrants);
            
            P(run,t)=mean(p);
            Y(run,t)=mean(y);           
        end
    end
    
    % Plots
    subplot(length(FF),2,2*(i_f-1)+1)
    plot(P','Color',rgb(str_color(5)))
    hold on
    plot(Y','b')
    hold off
    xlim([0 T])
    xlabel('time')
    ylabel({['F=', num2str(FF(i_f))];'y,p'});
    subplot(length(FF),2,2*(i_f-1)+2)
    edges=0:0.025:max(y)+0.025;
    histogram(y,edges)
    hold on
    histogram(y0,edges,'FaceColor','w')
    xlim([0 Inf])
    xlabel('y')
    ylabel('frequency')
    M=ylim;
    plot([P(end,end) P(end,end)],[0 M(2)],'color',rgb(str_color(5)),'LineWidth',2)
    plot([Y(end,end) Y(end,end)],[0 M(2)],'b','LineWidth',2)
    hold off
end
toc

