% helperBLEMeshMonteCarloSimulations Helper script for running Monte Carlo
% simulations for a Bluetooth mesh network to observe the following
% results:
%
%   1. Probability of a message to reach from source to destination in case
%   of relay nodes dropping randomly
%
%   2. Average hop count required between source and destination
%
%   3. Critical relays in having a path between source and destination

% Copyright 2019 The MathWorks, Inc.

% Number of nodes in the mesh network
totalNodes = 50;

% Initialize meshNodes vector with objects of type helperBLEMeshNetworkNode
meshNodes(1, totalNodes) = helperBLEMeshNetworkNode;
% Configure mesh nodes
for idx = 1:totalNodes
    meshNode = helperBLEMeshNetworkNode;
    meshNode.Identifier = idx;
    meshNode.NetworkLayer.ElementAddresses = dec2hex(idx, 4);
    meshNodes(idx) = meshNode;
end

% Load node positions
load('bleMeshNetworkNodePositions.mat');

% Create network
meshNetworkPlot = helperBLEMeshVisualizeNetwork();
meshNetworkPlot.NumberOfNodes = totalNodes;
meshNetworkPlot.NodePositionType = 'UserInput';
meshNetworkPlot.Positions = bleMeshNetworkNodePositions; % In meters
meshNetworkPlot.VicinityRange = 25; % In meters

% Simulation time, in milliseconds
simulationTime = 600;

% Enable/disable visualization
enableVisualization = false;

% Source destination
srcDstPair = [16 12];

% TTL value in message publication
ttl = 25;

% Relay nodes
relayNodeIDs = [21 15 25 11 38 19 46 8 39 20 37 32 30 5 45 49 43 3 28 36 47 ...
    34 23 48 41 44 42 10 4];

% Pre-allocate
numSeeds = 10000;
relayDropIdx = 1;
criticalRelaysIdx = 1;
simResultsIdx = 1;

% Simulation results
simResults = table('Size', [numel(relayNodeIDs)*numSeeds, 4], ...
    'VariableTypes', {'double','cell','cell', 'double'}, 'VariableNames', ...
    {'SimulationID', 'DisabledRelays', 'Path', 'NumberOfHops'});
criticalRelays = zeros(1, numSeeds);

% Run Monte-Carlo simulations
for idx = 1:numSeeds
    rng('shuffle');
    selectedRelays = relayNodeIDs;
    % Run simulation by disabling the relay feature for random nodes
    for j = 1:numel(selectedRelays)
        relayDropIdx = randi([1 numel(selectedRelays)]);
        droppedRelay = selectedRelays(relayDropIdx);
        selectedRelays = setdiff(selectedRelays, droppedRelay);
        pathSrcToDstScenario = helperBLEMeshFloodingSimulation(meshNodes, ...
            meshNetworkPlot, simulationTime, srcDstPair, ttl, ...
            selectedRelays, [], enableVisualization);
        
        % Update simResults
        simResults{simResultsIdx, 1} = simResultsIdx;
        simResults{simResultsIdx, 2} = {setdiff(relayNodeIDs, selectedRelays)};
        simResults{simResultsIdx, 3} = pathSrcToDstScenario;
        simResults{simResultsIdx, 4} = (numel(pathSrcToDstScenario{1}) - 1);
        simResultsIdx = simResultsIdx+1;
        
        % Run these simulations until there is no path between source and
        % destination. Then mark the dropped relay node as a critical node
        % for this route
        if pathSrcToDstScenario{1} == -1
            criticalRelays(criticalRelaysIdx) = droppedRelay;
            criticalRelaysIdx = criticalRelaysIdx+1;
            break;
        end
    end
end

% Simulation results
simResults = simResults(1:simResultsIdx-1, :);

% Total simulations
totalSims = size(simResults, 1);

% Failures (number of simulations having no path between source and
% destination)
failureSims = numel(find(cell2mat(simResults.Path') == -1));

% Success (number of simulations having path between source and
% destination)
successfulSims = totalSims - failureSims;

% Probability
probabilityOfSuccess = successfulSims/totalSims*100;

% Hop count
averageHopCount = ceil(mean(simResults.NumberOfHops));

% Critical nodes
uniqueCriticalRelays = unique(criticalRelays);
criticalRelaysInfo = table('Size', [numel(uniqueCriticalRelays), 2], ...
    'VariableTypes', {'double','double'}, 'VariableNames', ...
    {'Relay', 'Probability of failure'});
for idx = 1:numel(uniqueCriticalRelays)
    criticalnessPercentage = ...
        nnz(criticalRelays == uniqueCriticalRelays(idx))/numel(criticalRelays);
    criticalRelaysInfo{idx, 1} = ceil(uniqueCriticalRelays(idx));
    criticalRelaysInfo{idx, 2} = criticalnessPercentage*100;
end
criticalRelaysInfo = sortrows(criticalRelaysInfo, 'Probability of failure', 'descend');
vicinityNodesForSrc = helperBLEMeshVicinityNodes(srcDstPair(1), bleMeshNetworkNodePositions, 25);
vicinityNodesForDst = helperBLEMeshVicinityNodes(srcDstPair(2), bleMeshNetworkNodePositions, 25);
delNodes = [vicinityNodesForSrc' vicinityNodesForDst'];
delIdx = any(criticalRelaysInfo.Relay == delNodes, 2);
criticalRelaysInfo(delIdx, :) = [];

% Probability of having a route between source to destination for the above
% network, while randomly disabling the relay feature is stored in variable
% 'probabilityOfSuccess'
disp(['Probability of having a route between nodes ' ...
    num2str(srcDstPair(1)) ' and ' num2str(srcDstPair(2)) ' is: ' ...
    num2str(probabilityOfSuccess) '%']);

% Average hop count required between source to destination is stored in
% variable 'averageHopCount'
disp(['Average hop count between nodes ' num2str(srcDstPair(1)) ' and ' ...
    num2str(srcDstPair(2)) ' is: ' ...
    num2str(averageHopCount)]);

% Critical relays and their importance percentage in having a path between
% source to destination is stored in variable 'criticalRelaysInfo'
disp(['Apart from the vicinity nodes of node ' num2str(srcDstPair(1)) ...
    ' and node ' num2str(srcDstPair(2)) ', nodes [' num2str(criticalRelaysInfo{1:5, 1}') ...
    '] are the top five crucial relays in having a path between them.']);
