function [ stateTrajectory, trajectoryMatBinInds, jumpTimes, interjumpIntervals, maxNTransitions ] = sampleinitialstate( ModelParams, AlgorithmSpec, nParticles )
    % This function initialises the Markov jump process state trajectory
    % structures by sampling an initial state (for each particle).
    %
    % This is for a single block, i.e. ModelParams is a scalar struct and
    % all particles represent the same subpopulation (max number of
    % states).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Note: if using the augmented state model, every process has to start
    % in state 1, and the initial state distribution should reflect this.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 21.10.2014 - created, coded functionality.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 2,
        error('Insufficient arguments!')
    end
    if ~isfield(AlgorithmSpec, 'updateStepTimes'),
        error('updateStepTimes must be a field of ModelParams!')
    end
    if ~isfield(ModelParams, 'logInitialStateDist'),
        error('logInitialStateDist must be a field of ModelParams!')
    end
    if length(AlgorithmSpec.updateStepTimes) < 2,
        error('Must have at least two update step times (including start time)!')
    end
    if nargin < 3,
        nParticles = size(ModelParams.logInitialStateDist, 2);
    else
        if size(ModelParams.logInitialStateDist, 2) ~= nParticles,
            error('Size of logInitialStateDist does not match nParticles!')
        end
    end

    % Straightforward to initialise most of these structures.
    maxNTransitions = 1;
    trajectoryMatBinInds = true([1, nParticles]);
    [jumpTimes{1:nParticles}] = deal([AlgorithmSpec.updateStepTimes(1)]);
    jumpTimes = jumpTimes';
    [interjumpIntervals{1:nParticles}] = deal([AlgorithmSpec.updateStepTimes(2) - AlgorithmSpec.updateStepTimes(1)]);
    interjumpIntervals = interjumpIntervals';
    
    % We may need to sample initial states.
    stateTrajectory = samplefromcategoricaldistributions(exp(ModelParams.logInitialStateDist), 1);
    
end