%SW_STEADY.M
% 
% SIMULATIONS ON RING/REWIRED RING SMALL WORLD MODEL 
% FIXED I.C. IS SIGMOIDAL
% COMPUTES TIME TO (APPROXIMATE) STEADY STATE AS REWIRING PARAM VARIES
% 0 TO 1
%
%
% DJH Sep 2019

rand('state',100)
r = 0.5;
theta = 0.5;

N = 100;

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;

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

pvals = [0.1:0.025:0.9];             %rewiring probabilities

pruns = 1e4;    % number of Monte Carlo simulations for each choice of p

for s = 1:length(pvals)
        s
    
        p = pvals(s);
        
        for count = 1:pruns    % these are independent simulations for each choice of p
        
           %%%%% Rewiring phase 
           A = Aring;
           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);  % replace 2s with 1s
  
           % Need to rescale so that all rows sum to one
           deg = A*ones(N,1);   % out degree vector
           Ahat = diag(1./deg)*A;
            
           maxsteps = 10000;
           Udiff = 1; % record change in solution over each time step
                           % as part of numerical definition of steady state
           Uchangesum = 0; % record changes over time
                           % to test for promiscuity 
           steps = 0;
           Uk = Uzero; 
           while Udiff > 1e-8 & steps < maxsteps
              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
               Udiff = norm(Ukp1-Uk)/norm(Uk);
               Uchangesum = Uchangesum + sum(abs(Ukp1-Uk));
               Uk = Ukp1;
               steps = steps + 1;
           end
           psteps(count) = steps;
           Uchange(count) = Uchangesum/steps;
        end
        stepmean(s) = mean(psteps);
        stepstd(s) = std(psteps); 
        stepvar(s) = var(psteps);
        changemean(s) = mean(Uchange);
        changestd(s) = std(Uchange);
        changevar(s) = var(Uchange);
end


figure(1)
clf
subplot(1,2,1)
plot(pvals,stepmean,'b-','LineWidth',4)
hold on
plot(pvals,stepmean+stepstd/sqrt(pruns),'b--','LineWidth',2)
plot(pvals,stepmean-stepstd/sqrt(pruns),'b--','LineWidth',2)
grid on
xlabel('p')
ylabel('E[time]')
set(gca,'FontSize',14,'FontWeight','Bold')

subplot(1,2,2)
plot(pvals,changemean,'b-','LineWidth',4) 
hold on
plot(pvals,changemean+changestd/sqrt(pruns),'b--','LineWidth',2)
plot(pvals,changemean-changestd/sqrt(pruns),'b--','LineWidth',2)
grid on
xlabel('p')
ylabel('E[vacillation]')
set(gca,'FontSize',14,'FontWeight','Bold')

print -dpdf pic_steps.pdf 
print -dpng pic_steps.png

