%SW_SIM.M
% 
% SIMULATIONS ON RING/REWIRED RING SMALL WORLD MODEL 
% FIXED I.C. IS SIGMOIDAL 
%
%
% DJH Sep 2019

%rand('state',100)
rand('state',500)   % used this one for the other random number seed
r = 0.5;
theta = 0.5;

N = 100;

% Set up A for nearest neighbour ring network
A = zeros(N,N);
A = diag(ones(N-1,1),1) + diag(ones(N-1,1),-1);
A(1,N) = 1;
A(N,1) = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Rewiring phase 
p = 0.1;

for i = 1:N
    for j = 1:N
        if A(i,j) == 1
            %% rewire the directed edge with probability p
            if rand < p
                v = randperm(N);
                newj = v(1);
                A(i,j) = 0;
                A(i,newj) = 2;  % use 2 instead of 1 to avoid double rewire later
            end
        end
    end
end
A = sign(A);   % convert any 2s to 1s

deg = A*ones(N,1);        %out-degree vector
Ahat = diag(1./deg)*A;   %scale so that row sums = 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 %%%%
x = linspace(-10,10,N)'; 
sig1 = exp(x)./(1 + exp(x));
Uzero = 0.3 + 0.225*sig1; 
Uk = Uzero;
steps = 500-1;
steps = 10000-1;
%steps = 1500-1;
steps = 2000-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(1)
clf
waterfall(Usol')
ylabel('time')
xlabel('U component')
zlabel('opinion')
xlim([1 N])
set(gca,'FontSize',14,'FontWeight','Bold')
print -dpdf pic_500swsima.pdf 
print -dpng pic_500swsima.png

figure(2)
clf
semilogy(Udiff,'LineWidth',4)
xlim([1 steps])
xlabel('Time')
ylabel('Relative change in U') 
set(gca,'FontSize',14,'FontWeight','Bold')
print -dpdf pic_500swsimb.pdf 
print -dpdf pic_500swsimb.png 

figure(3)
clf
plot(Usol(:,end),'k-v','LineWidth',4,'MarkerSize',10)
xlim([1 N])
xlabel('U component')
ylabel('opinion')
legend('Final','Location','Southwest') 
set(gca,'FontSize',14,'FontWeight','Bold')
print -dpdf pic_500swsimc.pdf 
print -dpdf pic_500swsimc.png 