function fig10

% The computation engine is in the function markovHKsymodd

fig9data
fig9animation



function fig9animation

% Produces 778 bifurcation diagrams

load fig9data

n = size(data{1},2);
K = 1:length(data);
maxt=0;
for i=K
    maxt = max([maxt size(data{i},1)]);
end

for t=1:778
    clf
    hold on
	for k = K
        for i = 1:n
            tt = min([size(data{k},1) t]);
            if  data{k}(tt,i) > 10^(-14)
                fill([k-.5 k+.5 k+.5 k-.5],[(i-1)-.5 (i-1)-.5 (i-1)+.5 (i-1)+.5],(1 - ((data{k}(tt,i))/n)^(1/7))*[1 1 1],'EdgeColor','none')
            end
        end
	end
	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]);
	xlabel('k')
	ylabel('101 opinion classes')
    title(['t = ' num2str(t-1)])
	set(gcf,'Position',[232   100   600   300],'PaperPositionMode','auto')
	print('-dtiff','-r96',['fig9_' num2str(t-1)]);
end




function fig9data

% saves data for bifurcation diagrams

n = 101;
K = 1:30;
data = cell(length(K),1);
for i = 1:length(K);
    data{i} = ones(1,n); % data-dims: t,n
end
T = zeros(length(K),1);
t = 1;

while any(T==0)
    for k = K(find(T==0))
        p = markovHKsymodd(data{k}(t,:),k);
        if abs(p - data{k}(t,:)) <= 10^-12
            T(k) = t-1;
        else
            data{k}(t+1,:) = p;
        end
    end
    t = t+1;
    save fig9data data T
end


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;