function [MLE,HST] = HoopStats_EstimMixtureModel(Responses,Targets,IncUniform,nVonmComps,MuVonmComps,ShowFigure)
% "HoopStats_EstimMixtureModel"
% R02.00
% By Sam Berens (sam.berens@york.ac.uk).
%
% Estimates a mixture model for circularly distributed data. Components
% within the model can include a uniform distribution and arbitrarily many
% target/fixed position von Mises distributions.
%
% [MLE,HST] = HoopStats_EstimMixtureModel(...
%   Responses,Targets,IncUniform,nVonmComps,[MuVonmComps],[ShowFigure]);
%
% INPUT:
%   - "Responses" : an nx1 numeric denoting n angles (in radians) of
%                   circular responses;
%   - "Targets" : an nxt numeric denoting n angles (in radians) of t
%                 targets distributions. If sum(isnan(Targets),2) ==
%                 size(Targets,1), then response to target mappings are
%                 assumed to be mutually exclusive;
%   - "IncUniform" : a 1x1 logical specifying whether model should include
%                    a uniform distribution;
%   - "nVonmComps" : a 1x1 numeric denoting the number of non-target von
%                    Mises distributions to be modelled;
%   - "MuVonmComps" : [Optional] a cx1 numeric denoting the preferred
%                     direction (in radians) for each of c non-target von
%                     Mises distributions;
%   - "ShowFigure" : [Optional] a 1x1 numeric denoting the figure
%                    preference. If ShowFigure == 0, no figures will be
%                    drawn. If ShowFigure == 1, a figure of the final EM
%                    step will be drawn. If ShowFigure == 2, figures will
%                    be drawn to show each EM step;
%
% OUTPUT:
%   - "MLE" : The best fitting model;
%   - "HST" : Iteration history;

%% Initial set-up

% Define global variables:
global NoVonmC;

% Error check inputs and set MuVonmComps if not provided:
if ~exist('MuVonmComps','var') || isempty('MuVonmComps')
    MuVonmComps = nan(nVonmComps,1);
elseif numel(MuVonmComps) ~= nVonmComps
    error('MuVonmComps should be an nx1 double where n = nVonmComps;');
elseif size(MuVonmComps,1) ~= nVonmComps
    MuVonmComps = MuVonmComps(:);
end

% If ShowFigure input is not provided, draw nothing:
if ~exist('ShowFigure','var')
    ShowFigure = 0;
end

% NoVonmC is a bool that determines whether the scatter plots depicts
% target errors (when NoVonmC==1) or absolute response angles (when
% NoVonmC==1).
NoVonmC = nVonmComps == 0;

%% Try estimating the model using EM:
[MLE,HST] = HoopStats_RunEM(...
    Responses,Targets,IncUniform,nVonmComps,MuVonmComps,ShowFigure);


%% Alternative fitting procedures:
% Try other fitting procedures if the EM does not provide a good fit,
% or did not converge. Here, model fit is assessed by whether the
% deltaBIC statistic is larger or smaller than an arbitrary cutoff.
% Larger (less-negative) deltaBIC values indicate a comparatively poor
% fit to the data given that deltaBIC = BIC_H1 - BIC_H0.
if isempty(MLE) || (MLE.deltaBIC > -10)
    
    % In the case where there is one target, no von Mises components
    % and a uniform:
    if (size(Targets,2) == 1) && IncUniform && NoVonmC
        
        % Try hard clustering:
        MLE = HoopStats_HardCluster(Responses,Targets);
        
        if isempty(MLE) || MLE.deltaBIC > -10
            % If the model is still no good, discard it:
            MLE = [];
        end
        HST = [];
        
    end
    
end

return