function fig12

% The computation engine is in the functions opDyn and boundConf

fig11data
fig11animation

function fig11animated

load fig11data 

n = 101;
E = 1:30;

T = [];
for i=E
    T = [T size(data{i},2)-2];
end
maxT=max(T);

for t = 1:maxT+1
    clf
    hold on
    axis([ 0.5 30.5 -0.5 n-0.5 ])
    set(gca,'Layer','Top','XTick',2:2:30,'YTick',0:20:100,'Box','on','Position',[0.1 0.16 0.8 202/300]);
    for e=E
        plot(e*ones(size(data{e}(:,min([t size(data{e},2)])))),data{e}(:,min([t size(data{e},2)])),'r.');
    end
    xlabel('\epsilon')
    ylabel('opinions')
    title(['t = ' num2str(t-1)])
    set(gcf,'Position',[232   100   600   300],'PaperPositionMode','auto')
    print('-dtiff','-r96',['fig11t' num2str(t)]);
end

for e=3:18
    X=data{e}(:,size(data{e},2));X(:,2)=1;
    X=reduceOp(X);
    for i=1:size(X,1)
        text(e+0.2,X(i,1)+1.5,num2str(X(i,2)),'FontSize',7)    
    end
end
print('-dtiff','-r96',['fig11t' num2str(maxT+1)]);

function fig11data

% creates data for a bifurcation diagram for the HK agent-based with n =
% 101

n = 101;
E = 1:30;
X = [0:n-1]'; X(:,2)=1

data = cell(length(E),1);

for e = E
    data{e} = opDyn(X,e+10^-10);% e+10^-10 must be chosen due to numerical vulnerablity 
                                % of abs(X-Xnew)<=e when abs(X-Xnew)==e
    save fig11data data
end


function Y = opDyn(X,e)

Y = X(:,1);
Xnew = X;
X=X-1; % only to make X ~= Xnew

while not(all(X(:,1) == Xnew(:,1)))
    X = reduceOp(Xnew);
    Xnew = boundConf(X,e);
	Y = [Y allOp(Xnew)];
end


function X = boundConf(X,e)

% X = boundConf(X,e)
% X    : Matrix n x 2, 
%        first column : n opinions, 
%        second column : number of agents with this opinion
% e    : bound of confidence
%
% boundConf calculates one step in the HK process of OD
% 
% 02/2003 Jan Lorenz

n = size(X,1); % number of opinions

A=[]; % will get the confidence matrix

% produce rows of the confidence matrix
for i = 1:n
   d = abs(X(:,1)-X(i,1));
   a = [d <= e ]'; % a is confidence set of opinion i
   a = a.*X(:,2)'; % a weighted with number of agents of opinions
   a = a/sum(a); 
   A = [A; a]; 
end
X(:,1)=A*X(:,1);


function X = reduceOp(X)

% X = reduceOp(X)
% X    : Matrix n x 2, 
%        first column : n opinions, 
%        second column : number of agents with this opinion
%
% reduceOp collects all opinions with the same value in one row
% and adds the numbers of agents with this opinion
%
% 09/2002 Jan Lorenz
% Version 2.0

k=1; % row counter

while size(X,1) > k
   I = ( find(X(:,1) == X(k,1)) ); % row indices equal to X(k,1)
   X(k,2) = sum(X(I,2));
   X = X( sort([k setdiff([1:size(X,1)],I)]) ,:); % delete rows
   k=k+1; 
end


function X = allOp(X)

% X = allOp(X)
%
% X (in): Matrix n x 2, 
%        first column : n opinions, 
%        second column : number of agents with this opinion
% X (out): column Vector
%
% allOp blows up an opinion profile with collected equal opinions
% to a vector with single opinions, for plotting reasons

M = [];
for i = 1:size(X,1)
    m = zeros(X(i,2),size(X,1));
    m(:,i) = ones(X(i,2),1);
    M = [M;m];
end
X = M*X;
X = X(:,1);