function gamma = func_empiricalLoss(M,A,kaux)
% ------ Empirical Loss Function -------
% A contains the emission densities and the transition matrix.
% M:  m*m*m empirical joint distribution
% kaux: order of current model
% Used only in "func_CMAES.m" and "func_CMAESchrono.m"

maux=size(M,1);
Aaux=reshape(A,[maux+kaux-1,kaux]);
Oaux = Aaux(1:maux,:);

Qaux = zeros(kaux);
Qaux(:, 1:kaux-1) = (Aaux(maux+1:maux+kaux-1,:))';
Qaux(:, kaux) = ones(kaux,1) - sum(Qaux, 2);

piaux = func_stationaryLaw(Qaux);

Mhataux = zeros(maux,maux,maux);
for b=1:maux
    Mhataux(:,b,:) = Oaux * diag(piaux) * Qaux * diag(Oaux(b,:)) * Qaux * Oaux';
end
gamma = sum(sum(sum( Mhataux .* (Mhataux - 2 .* M) )));
