function [MLE,HST] = HoopStats_EstimMixtureModel(Responses,Targets,IncUniform,nVonmComps,MuVonmComps)
% "HoopStats_EstimMixtureModel"
% R01.03
% 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 = HoopStats_EstimMixtureModel(Responses,Targets,IncUniform,nVonM,tVonM)
%
% 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;
%
% OUTPUT:
%   - "MLE" : The best fitting model;
%   - "HST" : Iteration history;

%% Clear, define, and Shuffle RNG:
rng('shuffle');
global Responses_X Responses_Y ScatterPointRadius LineThickness TerminateOndNll Tolerance_MaxDeltaParams;

Responses_X = real(exp(Responses.*1i));
Responses_Y = imag(exp(Responses.*1i));

if ~exist('MuVonmComps','var')
    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

%% Estimation options:
EstimatePriors = true;
PriorBreakIn = Inf;
TerminateOndNll = true;
Tolerance_dNll = 1e-3;
Tolerance_MaxDeltaParams = 1e-3;
MaxIterations = 10^3;
NStarts = 3;
WNLLWPOST = 1; % Weight NLL with posterior probability: 0 or 1;

%%% Plot appearance:
ShowFigure = true;
ScreenSize = get(0,'ScreenSize');
% FigureHight = round(ScreenSize(1,4).*0.95);
% FigureSize = [ScreenSize(3)-(FigureHight.*16/9),(ScreenSize(4)-FigureHight)+20,FigureHight.*(16/9).*0.85,FigureHight.*0.85];
FigureSize = [1,ScreenSize(1)+40,ScreenSize(3),ScreenSize(4)-120];
ScatterPointRadius = 30;
LineThickness  = 3;

%% Start Loop:
try
    UniqueStarts = struct;
    for iStart = 1:1:NStarts
        
        %%% Specify starting parametres in ParamStruct:
        ParamStruct = struct;
        ParamStruct(1,1).Iteration = 0;
        for iTarget = 1:1:size(Targets,2)
            TargetId = sprintf('T%i',iTarget-1);
            ParamStruct(1,1).(TargetId) = struct;
            ParamStruct(1,1).(TargetId).DistributionType = 'vonMises';
            ParamStruct(1,1).(TargetId).Mu = 0;
            ParamStruct(1,1).(TargetId).K = (rand(1,1)*1)+0.5;
            ParamStruct(1,1).(TargetId).Prior = 1/(nVonmComps+double(IncUniform)+size(Targets,2));
        end
        if IncUniform
            ParamStruct(1,1).C0 = struct;
            ParamStruct(1,1).C0.DistributionType = 'Uniform';
            ParamStruct(1,1).C0.Prior = 1/(nVonmComps+1+size(Targets,2));
            %%% Here the p(Responses_i|C0) is always = 1/(2*pi);
            %%% Only the prior is estimated;
            %%% The posterior is proportunal to p(Responses_i|C0) * Prior;
        end
        StartMus = sort(rand(nVonmComps,1)*2*pi);
        for iVonmComp = 1:1:nVonmComps
            ComponentId = sprintf('C%i',iVonmComp);
            ParamStruct(1,1).(ComponentId) = struct;
            ParamStruct(1,1).(ComponentId).DistributionType = 'vonMises';
            if isnan(MuVonmComps(iVonmComp,1))
                ParamStruct(1,1).(ComponentId).Mu = StartMus(iVonmComp,1);
            else
                ParamStruct(1,1).(ComponentId).Mu = MuVonmComps(iVonmComp,1);
            end
            ParamStruct(1,1).(ComponentId).K = (rand(1,1)*1)+0.5;
            ParamStruct(1,1).(ComponentId).Prior = 1/(nVonmComps+double(IncUniform)+size(Targets,2));
        end
        for iTarget = 1:1:size(Targets,2)
            TargetId = sprintf('T%i',iTarget-1);
            ParamStruct(1,1).(['W_',TargetId]) = nan(size(Responses,1),1);
        end
        if IncUniform
            ParamStruct(1,1).W_C0 = nan(size(Responses,1),1);
        end
        for iVonmComp = 1:1:nVonmComps
            ComponentId = sprintf('C%i',iVonmComp);
            ParamStruct(1,1).(['W_',ComponentId]) = nan(size(Responses,1),1);
            % W is for Posterior;
        end
        ParamStruct(1,1).MaxDeltaParams = NaN;
        ParamStruct(1,1).NLL = NaN;
        
        %% Get the numbers and types of components being fitted:
        AFun_CellFun00 = @(s) strcmp(s(1,1),'C');
        AFun_CellFun01 = @(s) s(1,2);
        Fields = fieldnames(ParamStruct);
        DistNames = Fields(cellfun(AFun_CellFun00,Fields));
        DistNumsPre = cellfun(AFun_CellFun01,DistNames);
        DistNums = nan(1,size(DistNumsPre,1));
        for i_DistNumsPre = 1:1:size(DistNumsPre,1)
            DistNums(1,i_DistNumsPre) = str2double(DistNumsPre(i_DistNumsPre,1));
        end
        
        %% Draw initial figure:
        if ShowFigure
            FigureH = figure;
            set(FigureH,'Position',FigureSize);
            set(FigureH,'color','white');
            PlotLastEmStep(ParamStruct,FigureH);
            pause(1);
        end
        
        %% Start EM:
        Finished = false;
        iIteration = 0;
        CurrentState = struct;
        while ~Finished
            
            %%% Incrament iteration count and set variables:
            iIteration = iIteration + 1;
            %%% Targets:
            for iTarget = 1:1:size(Targets,2)
                TargetId = sprintf('T%i',iTarget-1);
                CurrentState.(TargetId).K = ParamStruct(iIteration,1).(TargetId).K;
                CurrentState.(TargetId).Prior = ParamStruct(iIteration,1).(TargetId).Prior;
            end
            if IncUniform
                CurrentState.C0.Prior = ParamStruct(iIteration,1).C0.Prior;
                CurrentState.C0.Pdents = ones(size(Responses,1),1) .* (1/(2*pi));
            end
            for iVonmComp = 1:1:nVonmComps
                % Remember that the first entry of 'ParamStruct' (index (i) = 1)
                % contains model starting parameters. As such, the data for
                % iteration 'a' is contained in index i = a + 1. Below, we read
                % in model values from the previous iteration (i = a);
                ComponentId = sprintf('C%i',iVonmComp);
                CurrentState.(ComponentId).Mu = ParamStruct(iIteration,1).(ComponentId).Mu;
                CurrentState.(ComponentId).K = ParamStruct(iIteration,1).(ComponentId).K;
                CurrentState.(ComponentId).Prior = ParamStruct(iIteration,1).(ComponentId).Prior;
            end
            
            %% E-Step:
            %%% Fixed Targets:
            for iTarget = 1:1:size(Targets,2)
                TargetId = sprintf('T%i',iTarget-1);
                Selector = ~isnan(Targets(:,iTarget));
                Acc = wrapToPi(Responses(Selector) - Targets(Selector,iTarget));
                [CurrentState.(TargetId).Pdents,~] = HoopStats_VonmFit(0,CurrentState.(TargetId).K,Acc);
            end
            
            %%% Componants:
            for iVonmComp = 1:1:nVonmComps
                ComponentId = sprintf('C%i',iVonmComp);
                [CurrentState.(ComponentId).Pdents,~] = HoopStats_VonmFit(CurrentState.(ComponentId).Mu,CurrentState.(ComponentId).K,Responses);
            end
            
            %%% Make Denominator:
            Denominator = zeros(size(Responses));
            % Targets:
            for iTarget = 1:1:size(Targets,2)
                TargetId = sprintf('T%i',iTarget-1);
                Selector = ~isnan(Targets(:,iTarget));
                TermToAdd = zeros(size(Responses));
                TermToAdd(Selector) = CurrentState.(TargetId).Prior .* CurrentState.(TargetId).Pdents;
                Denominator = Denominator + TermToAdd;
            end
            % Components:
            for nComp = DistNums % This includes the Uniform if present;
                ComponentId = sprintf('C%i',nComp);
                TermToAdd = CurrentState.(ComponentId).Prior .* CurrentState.(ComponentId).Pdents;
                Denominator = Denominator + TermToAdd;
            end
            
            %%% Compute weights:
            % Targets:
            for iTarget = 1:1:size(Targets,2)
                TargetId = sprintf('T%i',iTarget-1);
                Selector = ~isnan(Targets(:,iTarget));
                CurrentState.(TargetId).W = zeros(size(Responses));
                CurrentState.(TargetId).W(Selector) = (CurrentState.(TargetId).Prior .* CurrentState.(TargetId).Pdents) ./ Denominator(Selector);
            end
            % Components:
            for nComp = DistNums
                ComponentId = sprintf('C%i',nComp);
                CurrentState.(ComponentId).W = (CurrentState.(ComponentId).Prior .* CurrentState.(ComponentId).Pdents) ./ Denominator;
            end
            
            %% M-Step:
            % Targets:
            for iTarget = 1:1:size(Targets,2)
                TargetId = sprintf('T%i',iTarget-1);
                Selector = ~isnan(Targets(:,iTarget));
                Acc = wrapToPi(Responses(Selector) - Targets(Selector,iTarget));
                [~,CurrentState.(TargetId).K] = HoopStats_GetParams(Acc,CurrentState.(TargetId).W(Selector));
                CurrentState.(TargetId).Prior = mean(CurrentState.(TargetId).W(Selector),1);
            end
            % Components:
            for nComp = DistNums
                ComponentId = sprintf('C%i',nComp);
                if nComp ~= 0
                    if isnan(MuVonmComps(nComp,1))
                        [CurrentState.(ComponentId).Mu,CurrentState.(ComponentId).K] = HoopStats_GetParams(Responses,CurrentState.(ComponentId).W);
                    else
                        [~,CurrentState.(ComponentId).K] = HoopStats_GetParams(Responses,CurrentState.(ComponentId).W);
                    end
                end
                CurrentState.(ComponentId).Prior = mean(CurrentState.(ComponentId).W,1);
            end
            
            %% Update DeltaParams:
            iEntry = 0;
            DeltaParams = nan(size(Targets,2)+numel(DistNums),3);
            for iTarget = 1:1:size(Targets,2)
                iEntry = iEntry + 1;
                TargetId = sprintf('T%i',iTarget-1);
                DeltaParams(iEntry,2) = abs((1/CurrentState.(TargetId).K) - (1/ParamStruct(iIteration,1).(TargetId).K)); % Here, DeltaParams for K are difined as [1/K_i - 1/K_(i+1)];
                DeltaParams(iEntry,3) = abs(CurrentState.(TargetId).Prior - ParamStruct(iIteration,1).(TargetId).Prior);
            end
            for nComp = DistNums
                iEntry = iEntry + 1;
                ComponentId = sprintf('C%i',nComp);
                if nComp ~= 0
                    DeltaParams(iEntry,1) = abs(CurrentState.(ComponentId).Mu - ParamStruct(iIteration,1).(ComponentId).Mu);
                    DeltaParams(iEntry,2) = abs(CurrentState.(ComponentId).K - ParamStruct(iIteration,1).(ComponentId).K);
                end
                DeltaParams(iEntry,3) = abs(CurrentState.(ComponentId).Prior - ParamStruct(iIteration,1).(ComponentId).Prior);
            end
            MaxDeltaParamsExcludingPriors = DeltaParams(:,1:2);
            MaxDeltaParamsExcludingPriors = nanmax(MaxDeltaParamsExcludingPriors(:));
            MDPEP = MaxDeltaParamsExcludingPriors;
            clear MaxDeltaParamsExcludingPriors;
            MaxDeltaParams = nanmax(DeltaParams(:));
            
            %% Calculate NLL:
            % Here, likelihood is calculated to be the sum of the probability
            % densities for each component (T&C) and is weighted by the
            % posterior terms P(C*|Responses) [W*].
            iDist = 0;
            Ps = nan(size(Responses,1),size(Targets,2)+numel(DistNums));
            for iTarget = 1:1:size(Targets,2)
                iDist = iDist + 1;
                TargetId = sprintf('T%i',iTarget-1);
                Selector = ~isnan(Targets(:,iTarget));
                Acc = wrapToPi(Responses(Selector) - Targets(Selector,iTarget));
                [CurrentState.(TargetId).Pdents,~] = HoopStats_VonmFit(0,CurrentState.(TargetId).K,Acc);
                Ps(Selector,iDist) = CurrentState.(TargetId).Pdents .* (CurrentState.(TargetId).W(Selector).^WNLLWPOST);
            end
            for iVonmComp = 1:1:nVonmComps
                ComponentId = sprintf('C%i',iVonmComp);
                [CurrentState.(ComponentId).Pdents,~] = HoopStats_VonmFit(CurrentState.(ComponentId).Mu,CurrentState.(ComponentId).K,Responses);
            end
            for nComp = DistNums
                iDist = iDist + 1;
                ComponentId = sprintf('C%i',nComp);
                Ps(:,iDist) = CurrentState.(ComponentId).Pdents .* (CurrentState.(ComponentId).W.^WNLLWPOST);
            end
            NLL = -sum(log(nansum(Ps,2)),1);
            
            %% Update ParamStruct:
            ParamStruct(iIteration+1,1).Iteration = iIteration;
            %%% Targets:
            for iTarget = 1:1:size(Targets,2)
                TargetId = sprintf('T%i',iTarget-1);
                ParamStruct(iIteration+1,1).(TargetId).Mu = 0;
                ParamStruct(iIteration+1,1).(TargetId).K = CurrentState.(TargetId).K;
                ParamStruct(iIteration+1,1).(['W_',TargetId]) = CurrentState.(TargetId).W;
                if EstimatePriors && (MDPEP < PriorBreakIn)
                    ParamStruct(iIteration+1,1).(TargetId).Prior = CurrentState.(TargetId).Prior;
                else
                    ParamStruct(iIteration+1,1).(TargetId).Prior = ParamStruct(iIteration,1).(TargetId).Prior;
                end
            end
            %%% Components:
            for nComp = DistNums
                ComponentId = sprintf('C%i',nComp);
                if nComp ~= 0
                    ParamStruct(iIteration+1,1).(ComponentId).Mu = CurrentState.(ComponentId).Mu;
                    ParamStruct(iIteration+1,1).(ComponentId).K = CurrentState.(ComponentId).K;
                end
                ParamStruct(iIteration+1,1).(['W_',ComponentId]) = CurrentState.(ComponentId).W;
                if EstimatePriors && (MDPEP < PriorBreakIn)
                    ParamStruct(iIteration+1,1).(ComponentId).Prior = CurrentState.(ComponentId).Prior;
                else
                    ParamStruct(iIteration+1,1).(ComponentId).Prior = ParamStruct(iIteration,1).(ComponentId).Prior;
                end
            end
            ParamStruct(iIteration+1,1).MaxDeltaParams = MaxDeltaParams;
            ParamStruct(iIteration+1,1).NLL = NLL;
            
            %% Plot:
            if ShowFigure% && (nanmin([ParamStruct.NLL],[],2) == NLL)
                PlotLastEmStep(ParamStruct,FigureH)
            end
            
            %% Terminate if MaxDetlaParams is less than tollerence:
            if TerminateOndNll
                dNll = ParamStruct(iIteration,1).NLL - ParamStruct(iIteration+1,1).NLL;
                if abs(dNll) < Tolerance_dNll
                    Finished = true;
                elseif (abs(dNll) < 0.01) && (MaxDeltaParams < 0.001)
                    Finished = true;
                end
            else
                if EstimatePriors %#ok<UNRCH>
                    if MaxDeltaParams < Tolerance_MaxDeltaParams
                        Finished = true;
                    end
                else
                    if MDPEP < Tolerance_MaxDeltaParams
                        Finished = true;
                    end
                end
            end
            
            %% Terminate if reached Iteration limit.
            if iIteration == MaxIterations
                Finished = true;
            end
        end
        
        %%% Record Result:
        UniqueStarts(iStart,1).ParamStruct = ParamStruct;
        UniqueStarts(iStart,1).NLL_Last = ParamStruct(end,1).NLL;
        [UniqueStarts(iStart,1).NLL_Min,UniqueStarts(iStart,1).NLL_MinIndex] = nanmin([ParamStruct.NLL],[],2);
        close all;
    end
    
    UsedBestNll = true;
    if UsedBestNll
        [~,BestStart_Index] = min([UniqueStarts.NLL_Min],[],2);
        BestIte_Index = UniqueStarts(BestStart_Index,1).NLL_MinIndex;
    else
        [~,BestStart_Index] = min([UniqueStarts.NLL_Last],[],2); %#ok<UNRCH>
        BestIte_Index = size(UniqueStarts(BestStart_Index,1).ParamStruct,1);
    end
    HST = UniqueStarts(BestStart_Index,1).ParamStruct(1:BestIte_Index,1);
    MLE = HST(end,1);
    
    close all;
    FigureH = figure;
    set(FigureH,'Position',FigureSize);
    set(FigureH,'color','white');
    PlotLastEmStep(HST,FigureH);
    Notification_A01;
    
catch
    Notification_A02;
    warning('Mixture model failed to converge! %04d%02d%02d-%02d%02d%02d;',fix(clock));
end

return

function [] = PlotLastEmStep(ParamStruct,FigureH)

global Responses_X Responses_Y ScatterPointRadius LineThickness TerminateOndNll Tolerance_MaxDeltaParams;
%VectorTextOffset = 0.6 + (-0.1*1i);

%% Open figure:
set(0,'CurrentFigure',FigureH);
clf;
subplot(2,2,[1,3]);
DrawCircle;
hold on;

%% Get Colours and scatter:
Fields = fieldnames(ParamStruct);
M = [];
for iFields = 1:1:size(Fields,1)
    if strcmp(Fields{iFields,1}(1,1),'W')
        M = [M,ParamStruct(end,1).(Fields{iFields,1})]; %#ok<AGROW>
    end
end
[PointColours_RGB,Stereotypes_RGB] = MakeColours(M);
scatter(Responses_X,Responses_Y,ScatterPointRadius,PointColours_RGB,'filled',...
    'MarkerEdgeColor',[0,0,0],'LineWidth',1);

%%
axis equal;
axis off;
text(0.5,1,sprintf('Iteration: %i;',ParamStruct(end,1).Iteration),'FontSize',10);

iDist = 0;
Priors = nan(0,1);
DistNames = cell(0,1);
ParamsToDisplay = nan(0,2);
InfoContent = nan(0,1);
for iFields = 1:1:size(Fields,1)
    if strcmp(Fields{iFields,1}(1,1),'T')
        iDist = iDist + 1;
        DistNames{iDist,1} = Fields{iFields,1};
        Prior = ParamStruct(end,1).(DistNames{iDist,1}).Prior;
        K = ParamStruct(end,1).(DistNames{iDist,1}).K;
        Priors(iDist,1) = Prior;
        ParamsToDisplay(iDist,2) = K;
        InfoContent(iDist,1) = HoopStats_InfoContent(Prior,K);
    elseif strcmp(Fields{iFields,1}(1,1),'C')
        iDist = iDist + 1;
        DistNum = str2double(Fields{iFields,1}(1,2));
        DistNames{iDist,1} = Fields{iFields,1};
        if DistNum == 0
            %%% Uniform distrobution:
            Prior = ParamStruct(end,1).C0.Prior;
            ParamsToDisplay(iDist,1) = NaN;
            ParamsToDisplay(iDist,2) = NaN;
            Thetas = linspace(0,2*pi);
            X = cos(Thetas).*Prior;
            Y = sin(Thetas).*Prior;
            plot(X,Y,'Color',Stereotypes_RGB(iDist,:),'LineWidth',LineThickness./2);
            InfoContent(iDist,1) = Prior * log(2*pi);
        else
            %% vonM:
            Mu = ParamStruct(end,1).(DistNames{iDist,1}).Mu;
            K = ParamStruct(end,1).(DistNames{iDist,1}).K;
            Mag = HoopStats_K2R(K);
            Prior = ParamStruct(end,1).(DistNames{iDist,1}).Prior;
            ParamsToDisplay(iDist,1) = Mu;
            ParamsToDisplay(iDist,2) = K;
            InfoContent(iDist,1) = HoopStats_InfoContent(Prior,K);
            %%% Maths:
            Vector = exp(Mu*1i) * Mag;
            Vector_X = real(Vector);
            Vector_Y = imag(Vector);
            %             TextPos = Vector * VectorTextOffset;
            %             Text_X = real(TextPos);
            %             Text_Y = imag(TextPos);
            %%% Draw:
            plot([0;Vector_X],[0;Vector_Y],...
                'Color',Stereotypes_RGB(iDist,:),...
                'LineStyle','-',...
                'LineWidth',LineThickness);
            %             text(Text_X,Text_Y,...
            %                 ['\mu = ',num2str(Mu),';',10,...
            %                 '\it{k} = ',num2str(K),';']);
        end
        Priors(iDist,1) = Prior;
    end
end
ParamsToDisplay = [ParamsToDisplay,Priors];

%%
subplot(2,2,2);
for iInfoContent = 1:1:size(InfoContent,1)
    H = bar(iInfoContent,InfoContent(iInfoContent,1));
    if iInfoContent == 1
        hold on;
    end
    set(H,'FaceColor',Stereotypes_RGB(iInfoContent,:));
end

ParamsToDisplay(ParamsToDisplay==0) = NaN;
for iParamsToDisplay = 1:1:size(ParamsToDisplay,1)
    Text_X = iParamsToDisplay;
    Text_Y = -0.5;
    Mu = ParamsToDisplay(iParamsToDisplay,1);
    K = ParamsToDisplay(iParamsToDisplay,2);
    Prior = ParamsToDisplay(iParamsToDisplay,3);
    if isnan(Mu) && isnan(K)
        text(Text_X,Text_Y,...
            ['\it{p} = ',sprintf('%0.2f',Prior),';'],...
            'HorizontalAlignment','center');
    elseif isnan(Mu)
        text(Text_X,Text_Y,...
            ['\it{k} = ',sprintf('%0.2f',K),';',10,...
            '\it{p} = ',sprintf('%0.2f',Prior),';'],...
            'HorizontalAlignment','center');
    else
        text(Text_X,Text_Y,...
            ['\mu = ',sprintf('%0.2f',Mu),';',10,...
            '\it{k} = ',sprintf('%0.2f',K),';',10,...
            '\it{p} = ',sprintf('%0.2f',Prior),';'],...
            'HorizontalAlignment','center');
    end
end
title('Parameter estimates');
ylabel('Information content (nats)');
ylim([-1,2.2]);
set(gca, 'XTick',1:size(DistNames,1));
set(gca,'XTickLabel',DistNames');
plot(xlim,[log(2*pi),log(2*pi)],'LineStyle','--','LineWidth',LineThickness./4,'Color','k');
    
%%
subplot(2,2,4);
XaxisData = (1:size(ParamStruct,1))';
Y2Add = nan(0,1);
if size(ParamStruct,1) < 100
    N2Add = 100 - size(ParamStruct,1);
    XaxisData = (1:1:100)';
    Y2Add = nan(N2Add,1);
end
[HFit_Ax,HFit_L1,HFit_L2] = plotyy(...
    XaxisData,[[ParamStruct.NLL]';Y2Add],...
    XaxisData,[log10([ParamStruct.MaxDeltaParams]');Y2Add]); %#ok<ASGLU>
title('Iteration history');
ylabel(HFit_Ax(1),'Negative loig-likelihood');
ylabel(HFit_Ax(2),'Log_1_0 max \Delta params');
xlabel('Iteration number');
if TerminateOndNll
    % Nothing
else
    HFit_Ax(1,2).YLim = [log10(Tolerance_MaxDeltaParams),2];
end
HFit_Ax(1,2).YTick = flip([0,-2,-4,-6,-8,-10]);
if size(ParamStruct,1) < 100
    HFit_Ax(1,1).XLim = [0,100];
    HFit_Ax(1,2).XLim = [0,100];
end

pause(0.05);
hold off;

return

function [] = DrawCircle()
global LineThickness;
Thetas = linspace(0,2*pi);
X = cos(Thetas);
Y = sin(Thetas);
plot(X,Y,'k-','LineWidth',LineThickness./3);
return

function [PointColours_RGB,Stereotypes_RGB,FSDisp] = MakeColours(M)
if sum(sum(isnan(M))) == numel(M)
    M = ones(size(M)) ./ size(M,2);
else
    M(isnan(M)) = 0;
end
nColours = size(M,2);
Range = 1;
Offset = 0;
%%
Stereotypes_Angles = mod((((Range/nColours).*((1:1:nColours)-1))+Offset),1)' .* (2*pi);
PointColours_Complex = M * exp(Stereotypes_Angles*1i);
PointColours_Angles = angle(PointColours_Complex);
PointColours_Angles(PointColours_Angles<0) = (2*pi) + PointColours_Angles(PointColours_Angles<0);
PointColours_Angles  = PointColours_Angles ./ (2*pi);
PointColours_Mags = abs(PointColours_Complex);
PointColours_HSV = [PointColours_Angles,PointColours_Mags,ones(size(PointColours_Complex,1),1)];
PointColours_RGB = hsv2rgb(PointColours_HSV);
%%
Stereotypes_HSV = [Stereotypes_Angles./(2*pi),ones(size(Stereotypes_Angles)),ones(size(Stereotypes_Angles))];
Stereotypes_RGB = hsv2rgb(Stereotypes_HSV);
FSDisp = @() image(reshape(Stereotypes_RGB,size(Stereotypes_RGB,1),1,3));
return

function [] = Notification_A01()
Fs = 48000;
Vol = 0.4;
L1 = 0.1; % Seconds;
P1 = 1000;
L2 = 0.1; % Seconds;
P2 = 1200;

%
[T1] = MakeSound(Fs,L1,P1,Vol);
[T2] = MakeSound(Fs,L2,P2,Vol);
sound(T1,Fs);
pause(L1);
sound(T2,Fs);
pause(L2);
sound(T2,Fs);
pause(L1);
sound(T1,Fs);
return

function [] = Notification_A02()
Fs = 48000;
Vol = 0.4;
L1 = 0.1; % Seconds;
P1 = 1200;
L2 = 0.1; % Seconds;
P2 = 900;
%
[T1] = MakeSound(Fs,L1,P1,Vol);
[T2] = MakeSound(Fs,L2,P2,Vol);
sound(T1,Fs);
pause(L1);
sound(T2,Fs);
pause(L2);
sound(T2,Fs);
pause(L1);
sound(T2,Fs);
return

function [fT] = MakeSound(Fs,Length,Pitch,Vol)
FadeTime = 0.01; % Seconds;
fT = 0:(Pitch/Fs):(Pitch*Length);
fT = sin(fT.*2.*pi);
InEnv = 0:(1/(Fs*FadeTime)):1;
OutEnv = 1:-(1/(Fs*FadeTime)):0;
fT(1:numel(InEnv)) = fT(1:numel(InEnv)) .* InEnv;
fT(end-numel(InEnv)+1:end) = fT(end-numel(InEnv)+1:end) .* OutEnv;
fT = fT .* Vol;
return