function fig4

% These are graphic display commands
% The computation engine is in the functions onesidedOpDyn and boundConf

box on
X = [0;1];X(:,2)=1; % this initial profile will be extended in function onesidedOpDyn
e = 100;
t = 40;
Y = onesidedOpDyn(X,e+10^-10,t); % e+10^-10 must be chosen due to numerical vulnerablity 
                                 % of abs(X-Xnew)<=e when abs(X-Xnew)==e
n = size(Y,1);
c = colormap(hsv(n));
set(gca,'Layer','top','nextplot','replacechildren','ColorOrder',c)
axis([0 t 0 26*e])
plot(0:t,Y')
title(['\epsilon = ' num2str(e)])
xlabel('t')
ylabel('agents')


% Compute Statistics for the first 10 Clusters
t = size(Y,2);
for i=1:10
    clLoc(i) = Y(1,t);
    clSize(i) = sum(Y(:,t) == clLoc(i));
    inds = find(Y(:,t) == clLoc(i));
    k = t;
    while all(Y(inds,k) == clLoc(i))
        k = k-1;
    end
    T(i) = k;
    Y = Y(max(inds)+1:size(Y,1),:);
end
dclLoc = clLoc(2:length(clLoc))-clLoc(1:length(clLoc)-1);

save(['fig3statisticse' num2str(e)],  'T', 'clSize', 'clLoc', 'dclLoc')

% add data to figure
for i=1:length(T)
    text(T(i),clLoc(i),num2str(T(i)));
end
for i=1:length(dclLoc)
    text(34,(clLoc(i)+clLoc(i+1))/2,['\Delta = ' num2str(round(100*dclLoc(i))/100)]);
end
for i=1:length(clSize)
    text(40.5,clLoc(i),['# = ' num2str(clSize(i))]);
end
    
set(gcf,'Position',[232   258   560   317],'PaperPositionMode','auto')
print -dtiff -r96 fig3




function Y = onesidedOpDyn(X,e,t)

% Y = onesidedOpDyn(X,e)
% X    : Matrix n x 2, 
%        first column : n opinions, 
%        second column : number of agents with this opinion
% e    : bound of confidence
% t    : time step until we compute the dynamic
%
% Y(:,k) is the opinion profile of time step k
%
% The opinion profile will be extended in the process

Y = X(:,1);

for i=1:t
	% pou = point (index) of uniformity
    n = size(X,1);
	pou = max(find(X(2:n,1)-X(1:n-1,1) ~= 1))+1;
	if isempty(pou)
        pou=0;
	end
	
	X(pou+1:pou+2+2*e,1) = [[pou+1:pou+2+2*e]-1]'; 	% add enough length
	X(pou+1:pou+2+2*e,2) = 1;
	X = X(1:pou+2+2*e,:);
    X = reduceOp(X);
    
	X = boundConf(X,e);
    X = allOp(X);
    X(:,2) = 1;
	X(pou+1+e:pou+2+2*e,1) = [[pou+1+e:pou+2+2*e]-1]'; % flat tail to uniformit
    
    Y = [[Y; [size(Y,1):size(X,1)-1]'*ones(1,size(Y,2))] X(:,1)];
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);
    