%PREF_SIM.M
% 
% SIMULATIONS ON Barabasi/Albert MODEL
% USES PREF.M FROM THE CONTSET PACKAGE
%
% I.C. is piecewise constant
% Idea is to compare two different I.C.'s
%
% Plot evolution over 50 steps
%        profile after 50 steps
%        increments in U over 50 steps
%        increments in U over 5,000 steps   
%
% DJH Oct 2019

rand('state',100)

r = 0.5;
theta = 0.5;

N = 100;

% Set up preferential attachment network %
% Default parameters in pref.m % 
A = pref(N);

A = full(A);   % don't need to exploit sparsity

%Look at the matrix
figure(1)
subplot(2,1,1)
spy(A)
xlabel('')
ylabel('')
set(gca,'xtick',[])
set(gca,'ytick',[])
title('A')
subplot(2,1,2)
plot(sum(A),'r*-')
title('degree')
set(gca,'FontSize',14,'FontWeight','Bold')

print -dpdf pic_prefsima.pdf 
print -dpng pic_prefsima.png

% Need to rescale so that all rows sum to one
deg = sum(A);   % A is symmetric
Ahat = diag(1./deg)*A;

% Check that we have rescaled properly
Achk = norm(ones(N,1) - Ahat*ones(N,1))

epsi = 2.5;
u1 = 0.5 - 0.5*sqrt(1 - 4*r/epsi)
u2 = 0.5 + 0.5*sqrt(1 - 4*r/epsi)

%%%%% Initial values %%%%
Uzero = 0.8*ones(N,1);
Uzero(1:21) = 0.1;
Uk = Uzero;

steps = 5000-1;
Usol = zeros(N,steps+1);
Usol(:,1) = Uzero;
Udiff = zeros(steps,1); % record change in solution over each time step
for k = 1:steps   
   Ukp1 = Uk + r*(theta - Uk) + (epsi)*(Uk.*(1-Uk)).*(Ahat*Uk - theta);      
        %%%%%% Reflection---could make this better? %%%%
   [a1,b1] = find(Ukp1<0);
   [a2,b2] = find(Ukp1>1);
   Ukp1(a1) = 0;        % just perturb back to 0 
   Ukp1(a2) = 1;        % just perturb back to 1
   Usol(:,k+1) = Ukp1;
   Udiff(k) = norm(Ukp1-Uk)/norm(Uk);
   Uk = Ukp1;
end

figure(2)
clf
subplot(2,1,1)
waterfall(Usol(:,1:50)')
ylabel('time')
xlabel('U component')
zlabel('opinion')
xlim([1 N])
set(gca,'FontSize',14,'FontWeight','Bold')

%%%% Now repeat with a slightly different initial condition %%%%%
%%%%% Initial values %%%%
Uzero = 0.8*ones(N,1);
Uzero(1:22) = 0.1;
Uk = Uzero;

steps = 5000-1;
Usol = zeros(N,steps+1);
Usol(:,1) = Uzero;
Udiff = zeros(steps,1); % record change in solution over each time step
for k = 1:steps   
   Ukp1 = Uk + r*(theta - Uk) + (epsi)*(Uk.*(1-Uk)).*(Ahat*Uk - theta);      
        %%%%%% Reflection---could make this better? %%%%
   [a1,b1] = find(Ukp1<0);
   [a2,b2] = find(Ukp1>1);
   Ukp1(a1) = 0;        % just perturb back to 0 
   Ukp1(a2) = 1;        % just perturb back to 1
   Usol(:,k+1) = Ukp1;
   Udiff(k) = norm(Ukp1-Uk)/norm(Uk);
   Uk = Ukp1;
end

subplot(2,1,2)
waterfall(Usol(:,1:50)')
ylabel('time')
xlabel('U component')
zlabel('opinion')
xlim([1 N])
set(gca,'FontSize',14,'FontWeight','Bold')


print -dpdf pic_prefsimb.pdf 
print -dpng pic_prefsimb.png




