spread(net, netLabels, alpha)

% This function will create the following activation spreading
% matrix('activ'):
% Each ROW is a successive time step
% Each COLUMN is a brain area, corresponding to the layout in 'net' matrix
% ('netLabels' aswell)
% The function uses the heaviside function to calculate the activation of
% each brain area at successive time steps
% Passed to this function:
% net               -   connectivity matrix
% netLabels         -   node labels for net
% listOfStartNodes  -   cell array from gui list box of user-selected starting nodes
% listOfRemoveNodes -   cell array from gui of user-selected nodes for removal (will NOT be included in spreading activation)
% timeSteps         -   no. of time steps user wants to model (from gui slider control)

function [activ, activLabels] = spread(net, netLabels, alpha)

listOfStartNodes = get(findobj('tag', 'list_sel_startNodes'), 'string');
listOfRemoveNodes = get(findobj('tag', 'list_sel_removal'), 'string');
timeSteps = str2num(get(findobj('tag', 'time_step_editbox'), 'String'));
alpha = str2num(get(findobj('tag', 'alpha_editbox'), 'String'));
log_thresh = [];

%-----------------------------------------------------
% Checks
if isempty(net)
    errordlg('Cant run simulation until network data has been loaded!')
    activ = [];
    activLabels = [];
    return;
end

if isempty(listOfStartNodes)
    errordlg('Cant run simulation until 1 or more starting nodes have been selected')
    activ = [];
    activLabels = [];
    return;
end
%-----------------------------------------------------

set(findobj('tag', 'label_simulation_working'), 'visible', 'on');
drawnow;

noOfStartNodes = length(listOfStartNodes);

%-----------------------------------------------------
% Perform a check to make sure list of nodes for removal does not
% contain any nodes selected for starting

for checkCounter = 1: noOfStartNodes
    if isempty( strmatch( listOfStartNodes{checkCounter}, listOfRemoveNodes, 'exact')) == 0
        errordlg(['Warning - you have selected a starting node, ' listOfStartNodes{checkCounter} ', which is also a node selected for removal! Aborting simulation.'])
        activ = [];
        activLabels = [];
        set(findobj('tag', 'label_simulation_working'), 'visible', 'off');
        return;
    end
end
%------------------------------------------------------


% -----------------------------------------------------
% REMOVE nodes selected for removal from 'net' and 'netLabels'

noOfRemovals = length(listOfRemoveNodes);

for removeCounter = 1: noOfRemovals
    removeIndex = strmatch(listOfRemoveNodes{removeCounter}  , netLabels, 'exact');

    if isempty(removeIndex)
        errordlg(['Error in spread.m - couldnt remove node ' listOfRemoveNodes{removeCounter}])
        activ = [];
        set(findobj('tag', 'label_simulation_working'), 'visible', 'off');
        return;
    end

    net(:, removeIndex) = [];
    net(removeIndex, :) = [];
    netLabels(removeIndex) = [];

end
%----------------------------------------------------------

%-----------------------------------------------
% From netLabels and listOfStartNodes, work out which indices of activ are
% the starting nodes.

startIndices = [];

for startCounter = 1: noOfStartNodes
    startIndices(end+1) = strmatch(listOfStartNodes{startCounter}  , netLabels, 'exact');
end
%----------------------------------------------

noOfNodes = length(net);

% --------------------------------------------
% Set all start nodes in activ to 1 at time step 1 (1st row of activ), others to 0
% 1 = active, 0 = inactive

activ(1,1:noOfNodes) = 0; % Set all in 1st row to 0

for startActivCounter = 1: noOfStartNodes
    activ(1,startIndices(startActivCounter)) = 1; % Set all start nodes in 1st row to 1
end
%-----------------------------------------------

%start at t=2
for timeStepCnt = 2 : timeSteps % For every time step t

    %For every node nodeCnt store 0 (not active) or 1 (active)
    % at every timestep in net 'activ'
    % Determined by heaviside function - area is active if
    % the sum of w(ij) x(r)i - thet > 0    --  if the sum of (coupling
    % strength between this node and each other node times the 'activ'
    % activation value of the other node) at the previous time step
    % (timeStepCnt - 1), minus theta (threshold value), is greater than 0

    % Calculate threshold
    %Go round in a while loop until the correct value for thresh is found
    %- initially set at noOfNodes, each time decrement by 1, test to see if
    %sum of activities is greater than at last time step, and less than
    %alpha. Each subsequent time round the loop the last row of activ will
    %be over-written - the last time it is overwritten will be the correct
    %values (using correct thresh value)

    % After some tests, threshold values seemed to be generally in the region of 10 -
    % 30. So we shall start at 50 and work our way down. A warning will be
    % given if it looks like threshold should be above 50.
    % Additional condition that threshold is greater than 0

    thresh = 50;
    while testThresh(thresh, noOfNodes, activ, timeStepCnt, alpha, net) == 0 & thresh >= 1
        thresh = thresh - 1;
    end

    thresh = thresh + 1; %Because we have gone 1 step too far
    log_thresh = [log_thresh, ' ', num2str(thresh)]; %To print in the log at the end

    if thresh >= 49
        errordlg('Warning - threshold was calculated as >= 49. May need to readjust threshold calculating routine because the search starts at 50!');
    end

    for nodeCnt = 1:noOfNodes

        %calculate sum of inputs ("other nodes")
        sumOfInputs = 0;
        for otherNodeCnt = 1:noOfNodes
            if (otherNodeCnt ~= nodeCnt)
                couplStr = net(nodeCnt, otherNodeCnt) ^ 2; % Coupling strength = fibre density squared between startnode nodeCnt and endNode otherNodeCnt
                sumOfInputs = sumOfInputs + couplStr * activ(timeStepCnt-1, otherNodeCnt);
            end
        end

        sumInputsMinusThresh = sumOfInputs - thresh;

        if sumInputsMinusThresh >= 0
            activ(timeStepCnt, nodeCnt) = 1;
        else
            activ(timeStepCnt, nodeCnt) = 0;
        end
    end
end

% Return the labels cell array 'activLabels' which corresponds to the activ number array.
% Note that if there were nodes selected for removal, activLabels will no
% longer be the same as netLabels - it has some removed (see beginning of function).
activLabels = netLabels;

set(findobj('tag', 'label_simulation_working'), 'visible', 'off');

% Write an entry in the GUI_log
fc_saveLog('------------ Simulation run --------------');
log_startNodes = 'Start node(s): ';
for counter = 1: length(listOfStartNodes)
    log_startNodes = [log_startNodes, ' ' listOfStartNodes{counter}];
end

fc_saveLog(log_startNodes);

if isempty(listOfRemoveNodes)
    log_removeNodes = 'No nodes were selected for removal.';
else
    log_removeNodes = 'Nodes selected for removal: ';
    for counter = 1: length(listOfRemoveNodes)
        log_removeNodes = [log_removeNodes, ' ' listOfRemoveNodes{counter}];
    end
end

fc_saveLog(log_removeNodes);
fc_saveLog(['No. of time steps: ' num2str(timeSteps)]);
fc_saveLog(['Alpha value: ' num2str(alpha)]);
fc_saveLog(['Calculated threshold values : ' log_thresh]);

return;