function [AlgorithmParams, AlgorithmSpec] = setupalgorithm( ModelParams, ModelSpec, DataSpec, parametersFile )
    % Initialise algorithm parameters and relevant structures.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % AlgorithmParams is a struct array indexed by block. Initially there
    % is one block per state dimension.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 01.10.2014 - copied from old version.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 3,
        error('Insufficient arguments! Have you run setupdata yet?')
    end
    if nargin < 4,
        error('No parameters file supplied!')
    end
    if exist(parametersFile, 'file') ~= 2,
        error('You must supply the filename of a parameters file (.m file) on the search path!')
    end
    
    % Load script containing parameter values.
    % Chop off any '.m' ending.
    if strcmp(parametersFile((end - 1):end) , '.m'),
        parametersFile = parametersFile(1:(end - 2));
    end
    eval(sprintf('%s;', parametersFile));
    
    if ~isfield(ModelSpec, 'useAugStateModel'),
        error('useAugStateModel must be a field of ModelSpec!')
    end
    if ModelSpec.useAugStateModel,
        if ~isfield(ModelSpec, 'maxNAugStates'),
            error('maxNAugStates must be a field of ModelSpec!')
        end
    end
    if ~isfield(ModelSpec, 'maxNStates'),
        error('maxNStates must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'spikeTrainModelIndicator'),
        error('spikeTrainModelIndicator must be a field of ModelSpec!')
    end
    if ~isfield(ModelParams, 'nParticlesThisBlock'),
        error('nParticlesThisBlock must be a field of ModelParams!')
    end
    if ~isfield(ModelParams, 'nTrueStatesThisBlock'),
        error('nTrueStatesThisBlock must be a field of ModelParams!')
    end

    % Set default values where we can.
    if ~exist('updateStepInterval', 'var'),
        updateStepInterval = 1;
    end
    if ~exist('nUpdateSteps', 'var'),
        nUpdateSteps = 1000;
    end
    if ~exist('memorySafetyFactor', 'var'),
        memorySafetyFactor = 1 / 2;
    end
    if ~exist('unixMemory', 'var'),
        unixMemory = 4e9;
    end
    if ~exist('debugMode', 'var'),
        debugMode = true;
    end
    if ~exist('inferenceAlgorithmIndicator', 'var'),
        inferenceAlgorithmIndicator = 1;
    end
    if ~exist('nSweepsForMjpPost', 'var'),
        nSweepsForMjpPost = 10;
    end
%     if inferenceAlgorithmIndicator == 1,
        if ~exist('nParticles', 'var'),
            nParticles = 1000;
        end
        if sum([ModelParams.nParticlesThisBlock]) ~= nParticles,
            error('Supplied ModelParams.nParticlesThisBlock does not match nParticles!')
        end
        if ~exist('useFastMoveKernel', 'var'),
            useFastMoveKernel = false;
        end
        if ~exist('useResampling', 'var'),
            useResampling = true;
        end
        if ~exist('alwaysResample', 'var'),
            alwaysResample = false;
        end
        if ~exist('alwaysSampleParameters', 'var'),
            alwaysSampleParameters = false;
        end
        if ~exist('essThreshold', 'var'),
            essThreshold = 10;
        end
        if ~exist('nGibbsSweepsPerResampleMove', 'var'),
            nGibbsSweepsPerResampleMove = 1;
        else
%             if nGibbsSweepsPerResampleMove > 1,
%                 error('Multiple Gibbs sweeps per resample-move is not implemented!')
%             end
        end
        if ~exist('usePositiveDiscrimination', 'var'),
            usePositiveDiscrimination = true;
        end
        if ~exist('positiveDiscriminationThreshold', 'var'),
            positiveDiscriminationThreshold = 0.1;
        end
        if ~exist('allowSmallSubpopulationsToGoExtinct', 'var'),
            allowSmallSubpopulationsToGoExtinct = true;
        end
        if ~exist('smallSubpopulationBuffer', 'var'),
            smallSubpopulationBuffer = 2;
        end
%     elseif inferenceAlgorithmIndicator == 2,
%         if ~exist('nStreams', 'var'),
%             nStreams = 1;
%         end
%         if ~exist('nGibbsSweeps', 'var'),
%             nGibbsSweeps = 1000;
%         end
%         if ~exist('nBurnInSweeps', 'var'),
%             nBurnInSweeps = 0;
%         end
%     end
    if ~exist('pauseTime', 'var'),
        pauseTime = [];
    end
%     if ModelSpec.modelIndicator == 1,
%         if ~exist('nResampleLeadTimeBins', 'var'),
%             nResampleLeadTimeBins = 10;
%         end
%     end
       
    if ModelSpec.useAugStateModel,
        nStates = ModelSpec.maxNAugStates;
    else
        nStates = ModelSpec.maxNStates;
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % If we have an updateStepInterval, use it to determine number of
    % steps. If we cannot fit integer number of intervals of the same
    % duration, we make the first step briefer than the others.
    if ModelSpec.modelIndicator == 1,
        AlgorithmSpec.updateStepInterval = DataSpec.updateStepInterval;
        AlgorithmSpec.nUpdateSteps = DataSpec.nUpdateSteps;
        AlgorithmSpec.updateStepTimes = (DataSpec.startTime:AlgorithmSpec.updateStepInterval:(DataSpec.startTime + AlgorithmSpec.updateStepInterval * AlgorithmSpec.nUpdateSteps))';
        
    elseif (~isnan(updateStepInterval)) && (updateStepInterval > 0),
        AlgorithmSpec.updateStepInterval = updateStepInterval;
        AlgorithmSpec.nUpdateSteps = ceil((DataSpec.endTime - DataSpec.startTime) / AlgorithmSpec.updateStepInterval);
        firstIntervalDuration = DataSpec.endTime - DataSpec.startTime - (AlgorithmSpec.nUpdateSteps - 1) * AlgorithmSpec.updateStepInterval;
        AlgorithmSpec.updateStepTimes = [DataSpec.startTime, (DataSpec.startTime + firstIntervalDuration):AlgorithmSpec.updateStepInterval:(DataSpec.startTime + firstIntervalDuration + AlgorithmSpec.updateStepInterval * (AlgorithmSpec.nUpdateSteps - 1))]';
        
    % Otherwise we determine updateStepInterval from the number of steps.
    elseif (~isnan(nUpdateSteps)) && (nUpdateSteps > 0) && (mod(nUpdateSteps, 1) == 0),
        AlgorithmSpec.nUpdateSteps = nUpdateSteps;
        AlgorithmSpec.updateStepInterval = (DataSpec.endTime - DataSpec.startTime) / AlgorithmSpec.nUpdateSteps;
        AlgorithmSpec.updateStepTimes = (DataSpec.startTime:AlgorithmSpec.updateStepInterval:(DataSpec.startTime + AlgorithmSpec.updateStepInterval * AlgorithmSpec.nUpdateSteps))';
        
    else
        error('Either updateStepInterval must be a positive real or nUpdateSteps must be a positive integer!')
    end
    AlgorithmSpec.debugMode = debugMode;
    AlgorithmSpec.inferenceAlgorithmIndicator = inferenceAlgorithmIndicator;
    AlgorithmSpec.nSweepsForMjpPost = nSweepsForMjpPost;
%     if AlgorithmSpec.inferenceAlgorithmIndicator == 1,
        AlgorithmSpec.nParticles = nParticles;
        AlgorithmSpec.useFastMoveKernel = useFastMoveKernel;
        AlgorithmSpec.useResampling = useResampling;
        AlgorithmSpec.alwaysResample = alwaysResample;
        AlgorithmSpec.alwaysSampleParameters = alwaysSampleParameters;
        AlgorithmSpec.essThreshold = essThreshold;
        AlgorithmSpec.nGibbsSweepsPerResampleMove = nGibbsSweepsPerResampleMove;
        AlgorithmSpec.usePositiveDiscrimination = usePositiveDiscrimination;
        AlgorithmSpec.positiveDiscriminationThreshold = positiveDiscriminationThreshold;
        AlgorithmSpec.allowSmallSubpopulationsToGoExtinct = allowSmallSubpopulationsToGoExtinct;
        AlgorithmSpec.smallSubpopulationBuffer = smallSubpopulationBuffer;
%     elseif AlgorithmSpec.inferenceAlgorithmIndicator == 2,
%         AlgorithmSpec.nStreams = nStreams;
%         AlgorithmSpec.nGibbsSweeps = nGibbsSweeps;
%         AlgorithmSpec.nBurnInSweeps = nBurnInSweeps;
%     end
    AlgorithmSpec.memorySafetyFactor = memorySafetyFactor;

    AlgorithmSpec.availableMemory = getmemoryavailable( memorySafetyFactor, unixMemory );
    if ModelSpec.spikeTrainModelIndicator > 0,
        largestMatrixSizePerJumpPerParticle = max([ModelSpec.maxNStates * DataSpec.nSpikeTrains, ...
                                                    nStates ^ 2]);
    else
        largestMatrixSizePerJumpPerParticle = nStates ^ 2;
    end
    maxMaxNJumpsLimit = floor(AlgorithmSpec.availableMemory / largestMatrixSizePerJumpPerParticle / 8);
    if (~isnan(maxNJumpsLimit)) && (maxNJumpsLimit > 0) && (mod(maxNJumpsLimit, 1) == 0),
        AlgorithmSpec.maxNJumpsLimit = min([maxNJumpsLimit, maxMaxNJumpsLimit]);
        if maxNJumpsLimit > maxMaxNJumpsLimit,
            sprintf('Warning: maxNJumpsLimit is greater than memory resources can permit. Using %d instead.', maxMaxNJumpsLimit)
        end
    else
        AlgorithmSpec.maxNJumpsLimit = maxMaxNJumpsLimit;
    end
    AlgorithmSpec.pauseTime = pauseTime;
%     if ModelSpec.modelIndicator == 1,
%         AlgorithmSpec.nResampleLeadTimeBins = nResampleLeadTimeBins;
%     end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % These are indexed by block.
    AlgorithmParams.nBlocks = size(ModelParams, 1);
    AlgorithmParams.nBlocksEachSubpop = ones([AlgorithmParams.nBlocks, 1]);
    AlgorithmParams.nParticlesEachBlock = [ModelParams.nParticlesThisBlock]';
    AlgorithmParams.nTrueStatesEachBlock = [ModelParams.nTrueStatesThisBlock]';
    [AlgorithmParams.logParticleWeights{1:AlgorithmParams.nBlocks}] = deal([]);
    [AlgorithmParams.logNormalisedParticleWeights{1:AlgorithmParams.nBlocks}] = deal([]);
    % Initialise weights.
    for iBlock = 1:AlgorithmParams.nBlocks,
        AlgorithmParams.logNormalisedParticleWeights{iBlock} = -log(AlgorithmSpec.nParticles) * ones([AlgorithmParams.nParticlesEachBlock(iBlock), 1]);
        AlgorithmParams.logParticleWeights{iBlock} = zeros([AlgorithmParams.nParticlesEachBlock(iBlock), 1]);
    end
    AlgorithmParams.logParticleWeights = AlgorithmParams.logParticleWeights';
    AlgorithmParams.logNormalisedParticleWeights = AlgorithmParams.logNormalisedParticleWeights';
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end
