innoculate.m

% innoculate.m
% 10-6-14
% Successively removes highest-degree nodes (or randomly)
% SF, ER and SW graphs
% finite difference, node-average (but not homogeneous model)

clear

close all

N = 100;
beta = 0.1;     % 0.1
mu = 0.5;       % 0.5

removetype = 2;     % 1 == high-degree    2 == random

netype = 3; % 1 == SF   2 == ER    3 == SW

if netype == 1   % SF
    disp('SF Graph')
    m = 3;
    node = makeSF(N,m);
    
elseif netype == 2  % ER
    disp('ER Graph')
    p = 0.06;
    node = makeER(N,p);
    
elseif netype == 3   %SW
    disp('SW Graph')
    m = 3;
    p = 0.1;
    node = makeSW(N,m,p);
end

[A,degree,Lap] = adjacency(node);

[V,D] = eig(Lap);

for loop = 1:N
    eigval(loop) = D(loop,loop);
end

figure(1)
plot(eigval)
title('Eigenvalues')

[N,e,avgdegree,maxdegree,mindegree,numclus,meanclus,Lmax,L2,LmaxL2] = clusterstats(node);
displine('Average degree = ',avgdegree)
thr = mu/(avgdegree*beta);
displine('Thresh = ',thr)
pause(1)

% initial values
do = 1; ind = 0;
while do == 1
    ind = ind + 1;
    if degree(ind) == ceil(avgdegree);
        do = 0;
        initind = ind;
    end
end


newnode = node;
for loop = 1:20             % 20 for removetype == 1     40 for removetype == 2
    
    [N,e,avgdegree,maxdegree,mindegree,numclus,meanclus,Lmax,L2,LmaxL2] = clusterstats(newnode);
    
    disp(' ')
    displine('Number of nodes = ',N)
    disp(strcat('Number of edges = ',num2str(e)))
    disp(strcat('Mean degree = ',num2str(avgdegree)))
    displine('Maximum degree = ',maxdegree)
    disp(strcat('Number of clusters = ',num2str(numclus)))
    disp(strcat('mean cluster coefficient = ',num2str(meanclus)))
    disp(' ')
    disp(strcat('Lmax = ',num2str(Lmax)))
    disp(strcat('L2 = ',num2str(L2)))
    disp(strcat('Lmax/L2 = ',num2str(LmaxL2)))
    disp(' ')
    displine('Average degree = ',avgdegree)
    thr = mu/(avgdegree*beta);
    displine('Thresh = ',thr)
    disp(' ')
    
    Th(loop) = thr;
    
    
    [A,degree,Lap] = adjacency(newnode);
    
    [V,D] = eig(Lap);
    
    
    % initial values
    a = zeros(N,1);
    a(initind) = 1;

    
    dt = .5;
    c = a;
    for timeloop = 1:200
                
        c = eye(N,N)*(1-mu*dt)*c + beta*dt*(A*c).*(1-c);
        
        concentration(timeloop,:) = c';
        
    end
    
    
    x = 0:199;
    h = colormap(jet);
    figure(3)
    for nodeloop = 1:N
        rn = round(nodeloop*63/N + 1);
        y =  concentration(:,nodeloop)+0.001;
        semilogy(x,y,'Color',h(rn,:))
        hold on
    end
    hold off
    pause(0.75)
    
    if removetype == 1
        [Y,I] = max(degree);
    else
        newI = 0;
        while newI == 0
            I = randint(1,N);
            if (newnode(I).numlink ~=0)&(initind ~=I)
                newI = 1;
            end
        end
    end
    
    newnode = subnode(I,newnode);
    
    
    con = concentration(100,:);
    
    [y mask] = maskbilevel(con,1/(N-loop),1/(N-loop),0,1);
    
    Num = sum(y);
    Nm(loop) = Num;
    
    
end


figure(4)
plot(Th)
title('Threshold')
xlabel('Innoculation Iteration')

figure(5)
plot(Nm)
title('# Conc. Greater than 1/N')
xlabel('Innoculation Iteration')

% figure(6)
% DynamicDrawNet(newnode)