function fig5

% These are graphic display commands
% The computation engine is in the functions markovHKsymodd

clf;
n=101;
k=16;

%hk
p=ones(1,n);
t=0;

plotting=[0 1 2 4 6 10];

subplot(length(plotting),2,1);
colormap([0 0 0])
bar(0:n-1,p,1);
axis([-1 n 0 40 ]);
text(4,35,['t = 0']);
title('k = 16');
a=get(gca,'Position');
set(gca,'XTick',[0:20:100],'XTickLabel',[],'YTick',[0 20 40],'Layer','top')

for i=2:length(plotting)
    for j=t+1:plotting(i)
        p=markovHKsymodd(p,k);
    end
    t=plotting(i);
    subplot(length(plotting),2,2*i-1);
    colormap([0 0 0])
    bar(p,1);
    axis([-1 n 0 40]);
    a=get(gca,'Position');
    set(gca,'XTick',[0:20:100],'XTickLabel',[],'YTick',[0 20 40],'Layer','top')
    if i==6
        axis([0 n+1 0 1 ]);
        set(gca,'YColor',[1 0 0],'XTick',[0:20:100],'XTickLabel',[0:20:100],'YTick',[0 0.5 1],'Layer','top')
        text(4,0.875,['t = ' num2str(t)]);
        text(2,0.15,'zoom!','color',[1 0 0])
        xlabel('101 opinion classes')
    else
        text(5,35,['t = ' num2str(t)]);
    end    
end

plotting = [374 375 376 377 378 379];

for i=1:length(plotting)
    for j=t+1:plotting(i)
        p=markovHKsymodd(p,k); % due to numerical vulnerability the symmetry of P can break
                               % thus we have a special function for
                               % symmetric initial profiles
    end
    t=plotting(i);
    subplot(length(plotting),2,i*2);
    colormap([0 0 0])
    bar(0:n-1,p,1);
    axis([-1 n 0 40 ]);
    a=get(gca,'Position');
    set(gca,'XTick',[0:20:100],'XTickLabel',[],'YTick',[0 20 40],'Layer','top')
    if i==1
        axis([-1 n 0 1 ]);
        set(gca,'YColor',[1 0 0],'XTick',[0:20:100],'XTickLabel',[],'YTick',[0 0.5 1])
        title('k = 16')
        text(4,0.875,['t = ' num2str(t)]);
        text(2,0.15,'zoom!','color',[1 0 0])
    else
        text(4,35,['t = ' num2str(t)]);
    end
    if i == 6 
        set(gca,'XTickLabel',[0:20:100])
        xlabel('101 opinion classes')
    end
end

set(gcf,'Position',[232    49   600   500],'PaperPositionMode','auto')
print -dtiff -r96 fig4


function P = markovHKsymodd(P,k)

% P = markovHKsymodd(P,k)
% P : 1 x n vector of the opinion distribution
% k : discrete opinion dynamics
%
% computes one step in an HK interactive Markov chain for symmetric P

k = floor(k);
n = length(P);

I = [];
J = [];
S = []; % I, J, and S are to construct the transition matrix in sparse form
for i = 1:n
    minindP = max([1 i-k]);
    maxindP = min([n i+k]);
    inds = minindP:maxindP;
    mass = sum(P(inds));
    if mass > 10*eps
        mean = sum(P(inds).*inds)/mass;
    else
        mean = i;
    end
    j1 = floor(mean+10*eps);
    j2 = ceil(mean-10*eps);
    if j1 ~= j2
        I=[I;i;i];
        J=[J;j1;j2];
        S=[S;j2-mean;mean-j1];
    else
        I=[I;i];
        J=[J;j1];
        S=[S;1];
    end
end
A = sparse(I,J,S,n,n);
P = P*A;
P = (P+fliplr(P))/2;