function ModelParams = samplemodelparametersfromfullconditional( ModelParams, PriorParams, ModelSpec, DataSpec, PositionData, updateStepTime, stateTrajectory, maxNTransitions, stateOccupancyDurations, transitionCountsMat, firstArrivals, spikeCountsMat, stateEachObsTime, currentPosCountEachState, previousPosCountEachState, nObsEachState, currentTimeParticleInds, stateParticleIndsEachObsTime )
    % This function is used to sample parameter values from the full
    % conditional distribution. Sampling statistics have already been
    % computed from the sampled MJP paths and observations, although some
    % work may still need to be done in functions called from here.
    %
    % This function operates on multiple particles of the same
    % subpopulation (and same block, in fact), i.e. the max number of
    % states is the same for all particles and parameter structures.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ModelParams - struct, not a struct array, i.e. it cannot be indexed.
    % Contains all model parameters for particles of a single block/subpop.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 20.10.2014 - created, coded functionality for spike train model
    % paramaters only.
    % 07.11.2014 - coded for position model parameters.
    % 11.12.2014 - replaced coordsTrajectory with posTimeVec.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 10,
        error('Insufficient arguments!');
    end
    if ~isfield(ModelSpec, 'modelIndicator'),
        error('modelIndicator must be a field of ModelSpec!')
    end
    if ~isfield(PriorParams, 'transitionRatesDirichletHyperparams'),
        error('transitionRatesDirichletHyperparams must be a field of PriorParams!')
    end
    if ModelSpec.modelIndicator == 2,
        if ~isfield(PriorParams, 'jumpRateGammaHyperparams'),
            error('jumpRateGammaHyperparams must be a field of PriorParams!')
        end
        if ~isfield(ModelSpec, 'dominatingRateScalingFactor'),
            error('dominatingRateScalingFactor must be a field of ModelSpec!')
        end
    end
    if ~isfield(ModelSpec, 'useAugStateModel'),
        error('useAugStateModel must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'spikeTrainModelIndicator'),
        error('spikeTrainModelIndicator must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'positionModelIndicator'),
        error('positionModelIndicator must be a field of ModelSpec!')
    end
    if ~isequal(size(ModelParams), [1, 1]),
        error('ModelParams must be a singleton struct (not a struct array)!')
    end
    if ~isfield(ModelParams, 'nTrueStatesThisBlock'),
        error('nTrueStatesThisBlock must be a field of ModelParams!')
    end
    if ~isfield(ModelParams, 'nStatesThisBlock'),
        error('nStatesThisBlock must be a field of ModelParams!')
    end
    if ~isfield(ModelParams, 'nParticlesThisBlock'),
        error('nParticlesThisBlock must be a field of ModelParams!')
    end
    if updateStepTime < 0,
        error('updateStepTime must be positive!')
    end
    if ModelSpec.useAugStateModel,
        if nargin < 11,
            error('Insufficient arguments!');
        end
        if ~isfield(ModelSpec, 'augStateComponents'),
            error('augStateComponents must be a field of ModelSpec!')
        end
    else
        if ~isfield(ModelSpec, 'initialStateModelIndicator'),
            error('initialStateModelIndicator must be a field of ModelSpec!')
        end
        if ModelSpec.initialStateModelIndicator == 3,
            if ~isfield(PriorParams, 'initialStateDistDirichletHyperparams'),
                error('initialStateDistDirichletHyperparams must be a field of PriorParams!')
            end
        end
    end
    if ModelSpec.spikeTrainModelIndicator > 0,
        if nargin < 12,
            error('Insufficient arguments!');
        end
        if ~isfield(PriorParams, 'spikeRateGammaHyperparams'),
            error('spikeRateGammaHyperparams must be a field of PriorParams!')
        end
        if ~isfield(DataSpec, 'nSpikeTrains'),
            error('nSpikeTrains must be a field of DataSpec!')
        end
    end
    if ModelSpec.positionModelIndicator > 0,
        if nargin < 18,
            error('Insufficient arguments!');
        end
        if ~isfield(ModelSpec, 'positionJitterModelIndicator'),
            error('positionJitterModelIndicator must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'positionObsErrorPriorIndicator'),
            error('positionObsErrorPriorIndicator must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'positionObsErrorTolerance'),
            error('positionObsErrorTolerance must be a field of ModelSpec!')
        end
        if ModelSpec.positionModelIndicator == 1,
            if ~isfield(PriorParams, 'driftRateUniformHyperparams'),
                error('driftRateUniformHyperparams must be a field of PriorParams!')
            end
%             if ~isfield(PriorParams, 'obsErrorPrecisionGammaHyperparams'),
%                 error('obsErrorPrecisionGammaHyperparams must be a field of PriorParams!')
%             end
%             if ~isfield(ModelParams, 'posObsVar'),
%                 error('posObsVar must be a field of ModelParams!')
%             end
        end
        if ~isfield(PriorParams, 'asympMeanPriorDist'),
            error('asympMeanPriorDist must be a field of PriorParams!')
        end
        if (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
            if ~isfield(PriorParams, 'jitterPrecisionGammaHyperparams'),
                error('jitterPrecisionGammaHyperparams must be a field of PriorParams!')
            end
        elseif (ModelSpec.positionJitterModelIndicator == 3) && (ModelSpec.positionModelIndicator == 2)
            if ~isfield(PriorParams, 'jitterCovPriorDof'),
                error('jitterCovPriorDof must be a field of PriorParams!')
            end
            if ~isfield(PriorParams, 'jitterCovPriorScaleMat'),
                error('jitterCovPriorScaleMat must be a field of PriorParams!')
            end
        end
        if ~isfield(DataSpec, 'posObsInterval'),
            error('posObsInterval must be a field of DataSpec!')
        end
        if ~isfield(DataSpec, 'spatialBinWidth'),
            error('spatialBinWidth must be a field of DataSpec!')
        end
        if ~isfield(PositionData, 'posTimeVec'),
            error('posTimeVec must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'linearisedPositionTrajectory'),
            error('linearisedPositionTrajectory must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'gridCentreCoordsTrajectory'),
            error('gridCentreCoordsTrajectory must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'spaceTransformationFactor'),
            error('spaceTransformationFactor must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'positionIndToGridCentreCoordsMap'),
            error('positionIndToGridCentreCoordsMap must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'nValidDiscretePositions'),
            error('nValidDiscretePositions must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'validDiscretePositionsVector'),
            error('validDiscretePositionsVector must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'sandMats'),
            error('sandMats must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'resShifterMats'),
            error('resShifterMats must be a field of PositionData!')
        end
        if ~isfield(ModelParams, 'asympMeanPosInd'),
            error('asympMeanPosInd must be a field of ModelParams!')
        end
        if ~isfield(ModelParams, 'posCov'),
            error('posCov must be a field of ModelParams!')
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % State process parameters.

    if ModelSpec.useAugStateModel,
        % Generator and jump rate.
        if ModelSpec.modelIndicator == 2,
            [ModelParams.generatorMat, ModelParams.jumpRate] = samplegeneratormatrix( ModelSpec, PriorParams.transitionRatesDirichletHyperparams(1:(ModelParams.nTrueStatesThisBlock - 1)), PriorParams.jumpRateGammaHyperparams, maxNTransitions, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, stateOccupancyDurations, transitionCountsMat, firstArrivals, ModelParams.nStatesThisBlock );
        end

        % Initial state.
        ModelParams.logInitialStateDist = -inf([ModelParams.nStatesThisBlock, ModelParams.nParticlesThisBlock]);
        ModelParams.logInitialStateDist(1, :) = 0;
    else
        % Generator and jump rate.
        if ModelSpec.modelIndicator == 2,
            [ModelParams.generatorMat, ModelParams.jumpRate] = samplegeneratormatrix( ModelSpec, PriorParams.transitionRatesDirichletHyperparams(1:(ModelParams.nTrueStatesThisBlock - 1)), PriorParams.jumpRateGammaHyperparams, maxNTransitions, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, stateOccupancyDurations, transitionCountsMat );
        end
        
        % Initial state.
        if ModelSpec.initialStateModelIndicator == 1,
            ModelParams.logInitialStateDist = -log(repmat(ModelParams.nTrueStatesThisBlock, [ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock]));
        elseif ModelSpec.initialStateModelIndicator == 2,
            ModelParams.logInitialStateDist = computemjpstationarydistribution( ModelParams.generatorMat, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock );
        elseif ModelSpec.initialStateModelIndicator == 3,
            ModelParams.logInitialStateDist = sampleinitialstatedistribution( PriorParams.initialStateDistDirichletHyperparams(1:ModelParams.nTrueStatesThisBlock), maxNTransitions, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, stateTrajectory );
        else
            error('Invalid initialStateModelIndicator!')
        end
    end

    if ModelSpec.modelIndicator == 2,
        % Dominating rate.
        ModelParams.dominatingRate = max(abs(ModelParams.jumpRate), [], 1)' * ModelSpec.dominatingRateScalingFactor;

        % Dominating process transition matrix.
        ModelParams.logTransitionMat = formdominatingprocesstransitionmat( ModelParams.generatorMat, ModelParams.dominatingRate );
    elseif ModelSpec.modelIndicator == 1,
        ModelParams.logTransitionMat = samplehmmtransitionmatrix( ModelSpec, transitionCountsMat, firstArrivals, PriorParams.transitionRatesDirichletHyperparams(1:ModelParams.nTrueStatesThisBlock), ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, ModelParams.nStatesThisBlock );
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Spike train model parameters.
    
    if ModelSpec.spikeTrainModelIndicator > 0,
        if ModelSpec.modelIndicator == 1,
            ModelParams.logSpikeRate = samplediscretetimespikerate( ModelSpec, PriorParams.spikeRateGammaHyperparams, updateStepTime, DataSpec.nSpikeTrains, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, stateOccupancyDurations, spikeCountsMat );
        elseif ModelSpec.modelIndicator == 2,
            ModelParams.logSpikeRate = samplespikerate( PriorParams.spikeRateGammaHyperparams, updateStepTime, DataSpec.nSpikeTrains, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, stateOccupancyDurations, spikeCountsMat );
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Position model parameters.
    if ModelSpec.positionModelIndicator > 0,

        validObsTimeBinInds = PositionData.posTimeVec < updateStepTime;

        % Asmyptotic mean position.
        if ModelSpec.modelIndicator == 2,
            ModelParams.asympMeanPosInd = sampleasymptoticmeanpos( PriorParams, ModelSpec, PositionData, updateStepTime, ModelParams.driftRate, ModelParams.posCov, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, ModelSpec.randomWalkState, stateEachObsTime, nObsEachState, currentPosCountEachState, previousPosCountEachState, validObsTimeBinInds );
        else
            [ModelParams.asympMeanPosInd, nObsEachSateAll] = sampleasymptoticmeanposdiscretetime( PriorParams, ModelSpec, PositionData, updateStepTime, ModelParams.driftRate, ModelParams.posCov, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, stateEachObsTime, validObsTimeBinInds );
        end

        % Another sampling statistic common to more than one full
        % conditional.
        distancesFromMean = computedistancesfrommean(PositionData, ModelParams.asympMeanPosInd, stateParticleIndsEachObsTime, validObsTimeBinInds, ModelParams.nParticlesThisBlock);

        % Drift rate.
        % Only difference between continuous time and discrete time models:
        % sampling drift rate. In discrete time model this is kept as 0.
        if ModelSpec.modelIndicator == 1,
            ModelParams.driftRate = zeros([ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock]);
        else
%             ModelParams.driftRate = sampledriftrate( PositionData, PriorParams.driftRateUniformHyperparams, updateStepTime, ModelParams.posCov, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, ModelSpec.randomWalkState, stateEachObsTime, currentTimeParticleInds, distancesFromMean, validObsTimeBinInds );
        end
        
        % 'Jitter' model.
        if (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
            % Independent spatial dimensions. Covariances are stored as
            % matrices; covariance parameters in columns with 2 rows (Y
            % dim, X dim) rather than a 2x2 matrix (off-diagonals not
            % needed).

            % We sample the jitter covariance of the AR(1) model with fixed
            % time interval. It has dependence on drift rate, and is
            % proportional to a uniform times a gamma.
%             if (ModelSpec.modelIndicator == 2) && (ModelSpec.positionObsErrorPriorIndicator > 0),
%                 % Also observation error variance, specific to continuous
%                 % observations model.
%                 [ModelParams.jitterCov, ModelParams.posObsVar] = samplecovariances( PositionData, ModelSpec, PriorParams, updateStepTime, ModelParams.driftRate, ModelParams.posObsVar, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, stateEachObsTime, nObsEachState, currentTimeParticleInds, stateParticleIndsEachObsTime, distancesFromMean, validObsTimeBinInds );
%             else
                ModelParams.jitterCov = samplecovariances( PositionData, ModelSpec, PriorParams, updateStepTime, ModelParams.driftRate, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, stateEachObsTime, nObsEachState, currentTimeParticleInds, stateParticleIndsEachObsTime, distancesFromMean, validObsTimeBinInds );
%             end

            % Total covariance.
            ModelParams.posCov = ModelParams.jitterCov;
        elseif (ModelSpec.modelIndicator == 1) && (ModelSpec.positionJitterModelIndicator == 3),
            ModelParams.jitterCov = samplecovariancematrix( PriorParams, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, nObsEachState, distancesFromMean, stateEachObsTime );

            % Total covariance.
            ModelParams.posCov = ModelParams.jitterCov;
        else
            error('Invalid jitter model!')
        end

        if ModelSpec.positionModelIndicator == 1,
%             if ModelSpec.positionObsErrorPriorIndicator == 0,
%                 expectedPrecision = (norminv(1 - ModelSpec.positionObsErrorTolerance) ./ DataSpec.spatialBinWidth) ./ 2;
%                 ModelParams.posObsVar = repmat(1 ./ expectedPrecision, [ModelParams.nParticlesThisBlock, 1]);
%             elseif ModelSpec.positionObsErrorPriorIndicator == -1,
%                 ModelParams.posObsVar = ones([ModelParams.nParticlesThisBlock, 1]);
%             end

            % Total covariance - incorporate observations variance.
%             ModelParams.posCov = bsxfun(@times, ModelParams.jitterCov, permute(ModelParams.posObsVar, [3, 2, 1]));
            ModelParams.posCov = ModelParams.jitterCov;

            %             % The instantaneous jitter covariance, i.e. that which features
            %             % in the SDE. In the continuously observed model the AR(1)
            %             % jitter covariance is a function of the observation error.
            %             driftRateTerm = 2 .* log(ModelParams.driftRate) ./ (ModelParams.driftRate .^ 2 - 1);
            %             ModelParams.instJitterCov = bsxfun(@times, (bsxfun(@minus, ModelParams.jitterCov, permute(ModelParams.driftRate .^ 2, [3, 1, 2])) - 1), permute(ModelParams.posObsVar, [3, 2, 1]));
            %             ModelParams.instJitterCov = bsxfun(@times, ModelParams.instJitterCov, permute(driftRateTerm, [3, 1, 2])) ./ DataSpec.posObsInterval;
            %         elseif ModelSpec.positionModelIndicator == 2,
            %             driftRateTerm = log(ModelParams.driftRate) ./ (ModelParams.driftRate .^ 2 - 1);
            %             ModelParams.instJitterCov = bsxfun(@times, ModelParams.jitterCov, permute(driftRateTerm, [3, 1, 2])) ./ DataSpec.posObsInterval;
        end

%         % Drift rate.
%         ModelParams.driftRate = sampledriftrate( PositionData, PriorParams.driftRateUniformHyperparams, updateStepTime, ModelParams.asympMeanPosInd, ModelParams.posCov, ModelParams.nTrueStatesThisBlock, ModelParams.nParticlesThisBlock, stateEachObsTime, currentTimeParticleInds, distancesFromMean, validObsTimeBinInds );

%         % Instantaneous drift rate.
%         ModelParams.instDriftRate = -log(ModelParams.driftRate) ./ DataSpec.posObsInterval;

%         if ModelSpec.positionModelIndicator == 2,
%             ModelParams.stationaryPosCov = ModelParams.posCov;
%         elseif ModelSpec.positionModelIndicator == 1,
%             [ ~, ~, ~, ModelParams.stationaryPosCov ] = computeinstantaneousparameters( ModelParams.driftRate, ModelParams.jitterCov, ModelParams.posObsVar, ModelSpec.positionModelIndicator, DataSpec.posObsInterval );
%         end
%         if ModelSpec.modelIndicator == 2,
%             [ ~, ~, ~, ModelParams.stationaryPosCov ] = computeinstantaneousparameters( ModelParams.driftRate, ModelParams.jitterCov, ModelParams.posObsVar, ModelSpec.positionModelIndicator, DataSpec.posObsInterval, ModelSpec.randomWalkState );
%         end
                
    end


    
    
    function distancesFromMean = computedistancesfrommean(PositionData, asympMeanPosInd, stateParticleIndsEachObsTime, validObsTimeBinInds, nParticles)

        % Strip data.
        PositionData.linearisedPositionTrajectory = PositionData.linearisedPositionTrajectory(validObsTimeBinInds);
        PositionData.gridCentreCoordsTrajectory = PositionData.gridCentreCoordsTrajectory(validObsTimeBinInds, :);
        nObs = size(PositionData.linearisedPositionTrajectory, 1);
%         if nObs == 1,
%             distancesFromMean = [];
%             return;
%         end
        
        % Get position transformation factors for distance from positions
        % to asymptotic mean...
        asympMeanPosIndEachTime = asympMeanPosInd(stateParticleIndsEachObsTime);
        posMeanFactorMatInds = bsxfun(@plus, PositionData.linearisedPositionTrajectory, (asympMeanPosIndEachTime - 1) .* PositionData.nValidDiscretePositions);
        posMeanTransformationFactors = reshape(PositionData.spaceTransformationFactor(posMeanFactorMatInds(:)), [nObs, 1, nParticles]);
        clear posMeanFactorMatInds;

        % ...and real space coords of asymptotic means.
        asympMeanGridCentreCoords = permute(reshape(PositionData.positionIndToGridCentreCoordsMap(asympMeanPosIndEachTime(:), :), [nObs, nParticles, 2]), [1, 3, 2]);

        % Now we can compute the squared distances.
        distancesFromMean = bsxfun(@times, bsxfun(@minus, PositionData.gridCentreCoordsTrajectory, asympMeanGridCentreCoords), posMeanTransformationFactors);
%         distancesFromMean = bsxfun(@times, bsxfun(@minus, asympMeanGridCentreCoords, PositionData.gridCentreCoordsTrajectory(2:end, :)), posMeanTransformationFactors);
    end

%     function distancesFromMean = computedistancesfrommean(PositionData, asympMeanPosInd, stateParticleIndsEachObsTime, validObsTimeBinInds, nParticles)
% 
%         % Strip data.
%         PositionData.linearisedPositionTrajectory = PositionData.linearisedPositionTrajectory(validObsTimeBinInds);
%         PositionData.gridCentreCoordsTrajectory = PositionData.gridCentreCoordsTrajectory(validObsTimeBinInds, :);
%         nObs = size(PositionData.linearisedPositionTrajectory, 1);
%         if nObs == 1,
%             distancesFromMean = [];
%             return;
%         end
%         
%         % Get position transformation factors for distance from positions
%         % to asymptotic mean...
%         asympMeanPosIndEachTime = asympMeanPosInd(stateParticleIndsEachObsTime(2:end, :));
%         posMeanFactorMatInds = bsxfun(@plus, PositionData.linearisedPositionTrajectory(2:end), (asympMeanPosIndEachTime - 1) .* PositionData.nValidDiscretePositions);
%         posMeanTransformationFactors = reshape(PositionData.spaceTransformationFactor(posMeanFactorMatInds(:)), [(nObs - 1), 1, nParticles]);
%         clear posMeanFactorMatInds;
% 
%         % ...and real space coords of asymptotic means.
%         asympMeanGridCentreCoords = permute(reshape(PositionData.positionIndToGridCentreCoordsMap(asympMeanPosIndEachTime(:), :), [(nObs - 1), nParticles, 2]), [1, 3, 2]);
% 
%         % Now we can compute the squared distances from 'current' position
%         % to asymptotic mean.
%         distancesFromMean = bsxfun(@times, bsxfun(@minus, PositionData.gridCentreCoordsTrajectory(2:end, :), asympMeanGridCentreCoords), posMeanTransformationFactors);
% %         distancesFromMean = bsxfun(@times, bsxfun(@minus, asympMeanGridCentreCoords, PositionData.gridCentreCoordsTrajectory(2:end, :)), posMeanTransformationFactors);
%     end

end
