function [discreteSpikeTrainsMat, discreteTimeBinWidth, nTimeBins] = discretisespiketrains( spikeTrainsArray, modelIndicator, nTimeBins, discreteTimeBinWidth, startTime, endTime )
    % Take the cell array of spike trains (as column vectors of spike times
    % - in seconds) and return a matrix of binned spike train data. Can use
    % Bernoulli or discrete Poisson model.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % spikeTrainsArray - cell array with one cell for each neuron in the
    % sample. Each cell is a column vector of times, in seconds, of spikes.
    %
    % nTimeBins - the number of time bins of width discreteTimeBinWidth
    % that cover the duration of the entire data interval for analysis. Can
    % be [] if discreteTimeBinWidth is not [], which implies we want to
    % derive nTimeBins.
    %
    % discreteTimeBinWidth - the duration in seconds of a time bin for
    % discretisation. Assumes all bins are of equal duration. Can be [] if
    % nTimeBins is not [], which implies we want to derive
    % discreteTimeBinWidth.
    %
    % startTime is the lower time limit for the analysis, in seconds. A
    % suggested value is the time of the first position observation;
    % default is first spike observation time. Likely to be 0.
    %
    % endTime is the upper time limit for the analysis, in seconds. A
    % default value is the time of the last spike observation.
    %
    % modelIndicator is a positive integer: 1 for Bernoulli model, 2 for
    % discrete Poisson model.
    %
    % Note on discretisation time bins: if the ith bin starts at t_i, it
    % includes all spike observations in the interval [t_i, t_(i + 1)),
    % i.e. observations at the interval start is included but not at its
    % end.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 09.04.2013 - created. Binning is similar to discretisepositiondata,
    % using histc, but we 'chop off' the last bin which histc creates for
    % data beyond endTime.
    % 10.04.2013 - removed the 'chopping off' the last bin created by histc
    % binning vector. Now we always count the last bin even if it only
    % contains one time point (happens when discreteTimeBinWidth divides
    % endTime - startTime exactly). Reason: we were potentially losing
    % spikes at the end unless we manually set the endTime beyond the last
    % spike time.
    % - CORRECTION: reversed this last change. If we were to do that, the
    % data in the last bin could be over a shorter interval than other time
    % bins, which will introduce bias.
    % 18.04.2013 - added bin width and nTimeBins to output, made the
    % equivalent inputs optional (but not both; one can be derived from the
    % other), made startTime and endTime optional and provided defaults.
    % - changed so that if we supply nTimeBins, it overrides any value of
    % discreteTimeBinWidth also supplied.
    % 19.04.2013 - changed so that empty spike trains are allowed.
    % 22.04.2013 - added some code that converts empty spike trains to
    % [NaN], which ensures they lead to a column of all zeros/false in the
    % discretised spikes matrix.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 4,
        error('Insufficient arguments!')
    end
    if nargin < 5,
        startTime = min(cellfun(@min, spikeTrainsArray));
    end
    if nargin < 6,
        endTime = max(cellfun(@max, spikeTrainsArray));
    end
%     if any(cellfun('isempty', spikeTrainsArray)),
%         error('There are empty spike trains!')
%     end
    if startTime >= endTime,
        error('startTime must be before endTime!')
    end
    if (rem(modelIndicator, 1) ~= 0) || modelIndicator < 1,
        error('modelIndicator must be a positive integer!')
    end
%     if modelIndicator > 2,
%         error('modelIndicator must be 1 (Bernoulli) or 2 (discrete Poisson)!')
%     end
    if isempty(nTimeBins) && isempty(discreteTimeBinWidth),
        error('At least one of nTimeBins or discreteTimeBinWidth must be supplied (non-empty)!')
    end
    if isempty(nTimeBins),
        nTimeBins = floor((endTime - startTime) / discreteTimeBinWidth);
    else
        if isempty(discreteTimeBinWidth),
            discreteTimeBinWidth = (endTime - startTime) / nTimeBins;
        end
    end
    
    nSpikeTrains = length(spikeTrainsArray);
    
    % Replace empty spike trains with [NaN] to get histc to count all
    % zeros.
    emptySpikeTrainsLog = cellfun('isempty', spikeTrainsArray);
    if any(emptySpikeTrainsLog),
        spikeTrainsArray(emptySpikeTrainsLog) = cellfun(@(iSpikeTrain) NaN, spikeTrainsArray(emptySpikeTrainsLog), 'UniformOutput', false);
    end

    % These boundary times cover the entire session, so that endTime falls
    % in the last bin.
    discreteTimeBoundaries = startTime:discreteTimeBinWidth:(startTime + nTimeBins * discreteTimeBinWidth);

    discreteSpikeTrainsMat = cellfun(@(spikeTrain) histc(spikeTrain, discreteTimeBoundaries, 1), spikeTrainsArray, 'UniformOutput', false);
    discreteSpikeTrainsMat = reshape(cell2mat(discreteSpikeTrainsMat), [nTimeBins + 1, nSpikeTrains]);

    % Chop off the last time bin. These are counts beyond endTime.
    discreteSpikeTrainsMat = discreteSpikeTrainsMat(1:nTimeBins, :);
    
    if modelIndicator == 1,
        discreteSpikeTrainsMat = discreteSpikeTrainsMat > 0;
    end
end