simulateEpidemic.m 5.61 KB
Newer Older
Argyris Kalogeratos's avatar
Argyris Kalogeratos committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
%simulateEpidemic - This function performs a batch of simulations given 
%some epidemic parameters and a set of control strategies to test.
%
% function results = simulateEpidemic (graphPars, infectionSeedPars, ...
%                    epidemicPars, simulationPars, strategies, varargin)
%
% Input:
%   - graphPars: a structure with settings regarding the network to use 
%     for the simulations. 
%     > .graphType: available are 'preferentialAttachment', 'erdosRenyi', 
%       'smallWorld', 'adjacencyMatrix'. The first three are well-known 
%       random network models, while the last one allows for the graph to 
%       be an input. For each case, additional parameters should be provided
%       in this structure.
%     > .unweight  : if true, it considers a binary graph
%     > .symmetrize: if true, symmetrizes the graph
%     > .regularize: if true, regularizes the graph weights
%   - infectionSeedPars: a structure with settings regarding the initial
%     infection state. Uniform random initialization is used.
%     > .size : the percentage of initial infected nodes (seeds)
%   - epidemicPars: a structure with the SIS epidemic parameters.
%     > .beta : the independent infection rate for each edge
%     > .delta: the self-recovery rate for each infected node
%     > .rho  : the increase in recovery rate for each treated node
%   - simulationPars: a structure with more specific parameters related to
%     the performed simulations. 
%     > .T       : total simulation time (stopping criterion)
%     > .noSlopeT: characteristic time for the non stationarity test used
%       to determine when an experiment reached a saturation state. Set an 
30
%       Inf value to ignore this feature and simulate up to total time T.
Argyris Kalogeratos's avatar
Argyris Kalogeratos committed
31 32 33 34 35 36 37 38 39
%     > .timeStep: the time interval for the recording the infection network 
%       state and other results (low values wrt T will reduce simulation speed).
%     > .stateSnapshots: if true, then the treatment distribution and network 
%       state will be stored for each change.
%     > .consistentRNG: if true, it forces the preservation of the consistency
%       or the random number generator for producing reproducible results.
%   - strategies: a structure with the strategies to use in the simulations.
%     See more details in the file createStrategies.m.
%   - varargin: a list of other optional arguments, here 'verbosityLevel' 
40
%     indicates the detail of reported information to the terminal.
Argyris Kalogeratos's avatar
Argyris Kalogeratos committed
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
%
% Output:
%   - results: it is a structure with all the recorded information during
%     the simulation process.
%---
% This is part of the DRA Simulator package for Dynamic Resource Allocation
% strategies aiming to suppress SIS epidemics. See the comments header 
% of the main script DRAsimulator.m and the README.md for more details.
% This package is distributed under GNU GPL v.3 or later.
%
% Copyright (c) by Argyris Kalogeratos and Kevin Scaman, 2015.
%---

function results = simulateEpidemic (graphPars, infectionSeedPars, epidemicPars, simulationPars, strategies, varargin)

[found, verbosityLevel, ~]  = parsepar(varargin, 'verbosityLevel');
if (~found), verbosityLevel = 1; end

%% 1. NETWORK SPECIFICATION
    switch graphPars.graphType
61 62 63 64 65 66 67
    case 'erdosRenyi'
        G = erdosRenyiGraph(graphPars.N, graphPars.edge_prob);
    case 'preferentialAttachment'
        G = preferentialAttachmentGraph(graphPars.N, graphPars.edge_num, graphPars.delta);
    case 'smallWorld'
        G = smallWorldGraph(graphPars.N, graphPars.edge_num, graphPars.edge_prob);
    case 'adjacencyMatrix'
Argyris Kalogeratos's avatar
Argyris Kalogeratos committed
68
        G = graphPars.adjacencyMatrix; 
69
        otherwise
Argyris Kalogeratos's avatar
Argyris Kalogeratos committed
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
        error('The graphType field of graphPars should be either ''erdosRenyi'', ''preferentialAttachment'', ''smallWorld'', ''adjacencyMatrix''');
    end
    
    graphPars.N = size(G, 1);
    
    % Some modifications to the graph, if needed
    if (isfield(graphPars, 'unweight') && graphPars.unweight)
        G = sparse(1 * (G ~= 0));
    end
    if (isfield(graphPars, 'symmetrize') && graphPars.symmetrize)
        G = (G + G') / 2;
    end
    if (isfield(graphPars, 'regularize') && graphPars.regularize)
        G = G * (nnz(G) / sum(G(:)));
    end
        
    % Show the properties of the generated graph
    if (verbosityLevel >= 2)
        graphProperties(G);
    end
        
%% 2. INFECTION SEEDS
    infectionSeedPars.size = max(1/graphPars.N, infectionSeedPars.size);
    infectionSeeds = randomSeed(graphPars.N, infectionSeedPars.size);
    if (verbosityLevel >= 1)
        fprintf('Seeding: %d infections uniformly distributed in network.\n', nnz(infectionSeeds));
    end
    
%% 3. EVOLUTION OF THE SYSTEM (performs the epidemic and control simulation)
    % Intervals of updating the Dynamic Resource Allocation strategy
    if ~(isfield(simulationPars, 'updateInterval') && ~isempty(simulationPars.updateInterval))
        simulationPars.updateInterval = @() -1;     % will always update the resource allocation
    end

    methods = fieldnames(strategies);
    numOfMethods = numel(methods);
    res = cell(numOfMethods, 1);
    
    % use simple "for" loop or "parfor" loop to run each method in parallel
    if (simulationPars.runInParallel || simulationPars.consistentRNG) % this approach respects the random number generator
        for m=1:numOfMethods,
            res{m} = epidemicSimulation(G, infectionSeeds, epidemicPars, strategies.(methods{m}), simulationPars);
        end
    else % run in parallel using many Matlab workers 
        parfor m=1:numOfMethods,
            res{m} = epidemicSimulation(G, infectionSeeds, epidemicPars, strategies.(methods{m}), simulationPars);
        end
    end
    for m=1:numOfMethods,  % gather results
        results.(methods{m}) = res{m};
    end
end